Critical endpoint of QCD in a finite volume

We investigate the impact of finite volume and the corresponding restrictions on long-range correlations on the location of the critical endpoint in the QCD phase diagram. To this end, we employ a sophisticated combination of lattice Yang--Mills theory and a (truncated) version of Dyson--Schwinger equations in Landau gauge for $2 + 1$ quark flavors that has been studied extensively in the past. In the infinite-volume limit, this system predicts a critical endpoint at moderate temperature and large chemical potential. We study this system at small and intermediate volumes and determine the dependence of the location of the critical endpoint on the boundary conditions and the volume of a three-dimensional cube with edge length $L$. We demonstrate that noticeable volume effects of more than five percent only occur for $L \lesssim 5 \, \text{fm}$ and that volumes as large as $L^3 \gtrsim (8 \, \text{fm})^3$ are very close to the infinite-volume limit.


I. INTRODUCTION
There are a number of reasons why finite-volume studies of the location and properties of a putative chiral critical endpoint (CEP) in the phase diagram of QCD are interesting. First, the search of such a CEP in heavy-ion collisions produces an initial fireball of finite extent with typical scales of several femtometers in each direction. Volume effects on important observables such as fluctuations of conserved charges may be considerable and need to be taken into account [1,2]. Second, the location of all phase boundaries including the crossover at zero chemical potential, the CEP, and the chiral first-order transition will certainly become volume dependent for small enough volumes. At small chemical potential, these are accessible by various extrapolation methods applied to lattice QCD. The corresponding simulations are carried out at finite boxes with (anti)periodic boundary conditions, and a thorough understanding of the volume dependence of the results is mandatory. Third, effects due to changes in the volume are not only annoying artifacts but interesting in themselves because they serve to probe the reaction of a physical system on one of its external parameters. This is particularly interesting when it comes to small quark masses and important questions on the properties of the upper and lower left (chiral) corners in the Columbia plot.
Theoretical studies of finite-volume effects on the QCD phase diagram have been carried out in a number of approaches besides lattice QCD. Model studies in the Nambu-Jona-Lasinio or the quark-meson model and their Polyakov-loop enhanced versions serve to highlight important general aspects of small-volume physics; see, e.g., Refs. [3][4][5][6] and Ref. [7] for a review. In particular, the renormalization-group treatment of these models allows * julian.bernhardt@physik.uni-giessen.de † christian.fischer@theo.physik.uni-giessen.de ‡ philipp.isserstedt@physik.uni-giessen.de § bernd-jochen.schaefer@theo.physik.uni-giessen.de to study the reaction of fluctuations on changes of volume [3,8,9]. Functional approaches to QCD via Dyson-Schwinger equations (DSEs) or the functional renormalization group (FRG) have the additional benefit of treating the Yang-Mills sector explicitly, which allows for quantitative studies on the location of the CEP [10][11][12][13][14].
In this work, we employ the DSE framework of Ref. [11] and treat the system in a finite three-dimensional box of equal edge lengths L. We contrast periodic and antiperiodic boundary conditions and vary the box size between L = 3 fm and L = 8 fm. We trace the location of the CEP in the QCD phase diagram for various volumes and determine the box size necessary to approach the infinitevolume results. Furthermore, we discuss volume effects on the curvature of the crossover line at small chemical potential and compare with lattice QCD.
The paper is organized as follows. In Sec. II, we discuss all technical details of our framework. We briefly summarize the truncation scheme for the DSEs and explain our implementation of finite volume. In particular, we explain a method to get rid of cubic artifacts at large momenta in the ultraviolet (UV) that drastically improves the approach to the infinite-volume limit. In Sec. III, we discuss our results for the volume dependence of the chiral crossover, its curvature, and the CEP and compare UV improved with unimproved calculations. Finally, we summarize in Sec. IV.

A. In-medium propagators and finite volume
The quantities of interest in this work are the dressed quark and gluon propagators S f and D νσ , respectively, where f labels the quark flavor. In Landau gauge and at nonzero temperature T and quark chemical potential µ f , they are given by where p = (ω n , p) is the four-momentum with Matsubara frequencies ω n , n ∈ Z. 1 An additional quark tensor structure γ 4 γ · p is in principle possible but negligible [15]. The dressing functions C f , A f , and B f of the quark and Z T and Z L of the gluon have a nontrivial momentum dependence and contain all nonperturbative information. Furthermore, they depend implicitly on temperature and quark chemical potential. The structure of Eqs.
(1) and (2) reflects the breaking of the O(4) symmetry in vacuum down to O(3) at nonzero temperature due to the coupling to the heat bath. 2 In particular, the gluon splits into a part transversal (T) and longitudinal (L) with respect to the heat bath. The corresponding projectors read The fermionic Matsubara frequencies include the chemical potential, ω n = (2n + 1)πT + iµ f , and the bosonic frequencies read ω n = 2nπT .

Finite volume without UV improvement
The introduction of a finite, uniform, three-dimensional spatial volume with (anti)periodic boundary conditions has similarities with the introduction of nonzero temperature. Within the Matsubara formalism, nonzero temperature is implemented in the QCD action by compactifying the imaginary-time integration to the finite interval [0, 1/T ], i.e., where L is the QCD Lagrangian. Restricting this action further to a finite spatial cube with edge length L amounts to the additional replacement As a consequence, in momentum space, the spatial threemomentum becomes discrete in addition to the already discrete Matsubara frequencies ω n in (imaginary) time direction. Hence, each generic momentum integral with kernel K is replaced by a corresponding sum over spatial modes: with q z = 3 i=1 ω L zi e i and the Cartesian basis {e 1 , e 2 , e 3 } that spans the three-dimensional momentum space. While the Matsubara frequencies in temporal direction are fixed by the spin-statistics theorem-odd (even) multiples of πT for quarks (gluons)-there is no similar constraint in spatial direction. Although it turns out that kinematics necessitates periodic spatial boundary conditions for the gluon, the quarks are completely free. In this work, we will choose periodic spatial boundary conditions (PBC) or antiperiodic spatial boundary conditions (ABC) for the quarks according to where n ∈ Z. Note that PBC contain a zero mode q zero = 0 corresponding to the lowest possible momentum. Because its inclusion in both the quark self-energy and the quark-loop contribution to the gluon self-energy requires special attention, we will discuss it in more detail below. Moreover, in order to compare with other works, we also perform calculations where the PBC zero mode is omitted; we refer to this setup as PBC * . From a numerical point of view, it is beneficial to rearrange the three summations in Eq. (7) such that they resemble a spherical coordinate system [16], The index j labels the spheres with constant radius |q z |, and m = m(j) denotes the multiplicity of the individual momentum vectors on a given sphere j. The corresponding vectors are denoted by q jm ; see Fig. 1 for a two-dimensional sketch.
Within the spherical momentum summation, it is straightforward to implement an O(3) invariant cutoff Λ. To this end, we distinguish between "complete" and "incomplete" spheres: while the number of points on complete spheres remains the same when the number of grid points in each Cartesian direction is increased, the incomplete ones obtain additional points. Thus, on each (sufficiently large) grid we have a large number of complete spheres in the innermost region up to some distance from the origin, whereas quite a large number of the outermost spheres are necessarily incomplete. We sum only over complete spheres (red, solid circles in Fig. 1) and discard the incomplete ones (blue, dashed circle in Fig. 1). This procedure helps to reduce so-called cubic artifacts [16,17]. However, as we will discuss in the next section, additional efforts are necessary to get rid of these completely.

Finite volume with UV improvement
The spherical finite-volume summation introduced above has a severe drawback in terms of numerical cost. The larger the grid the more dense are the complete outer spheres and the more points are on every one of these. On the other hand, since the dressing functions run logarithmically as functions of large squared momenta in the UV, they do not change much from outer sphere to outer sphere. Therefore, a lot of numerical effort is spent to integrate an almost constant function. This effort can be drastically reduced by the following procedure. We consider discrete spheres up to some matching cutoff Λ vol and replace the spheres with radii between Λ 2 vol < q 2 < Λ 2 with a continuous momentum integral. As a consequence, the original replacement Eq. (7) is modified such that an additional integral over the continuous momenta is added to the sum over the spatial modes, This modification, called "UV improvement" in the following, allows for much larger values for the ultraviolet cutoff Λ, which should mitigate artifacts that are caused by a potentially too small cutoff. In addition, this improvement also allows to renormalize the system at a large subtraction point where medium and volume effects are negligible. Thus, the renormalization procedure can be carried out identically as in infinite-volume calculations. Recently, a similar treatment has been used in Ref. [18] within a simpler truncation of the corresponding DSEs.
In practice, we observe that our results do not depend on Λ vol as long as Λ vol is larger than any other typical scale of the system like temperature, chemical potential, and quark masses. This gives us the freedom to choose different matching points for different box sizes L and essentially work with a fixed number of grid points in our three-dimensional grid for momenta with magnitudes below Λ vol .
In Sec. III below, we will discuss results from the unimproved treatment of the system in the box as well as the UV improved one and systematically compare both approaches.

Inclusion of the zero mode
In order to properly include the PBC zero modes of the external momentum p and the loop momentum q into our numerical setup, we have to make some minor adjustments.
(i) Our kernels are not well-defined at q = 0, but smooth in the limit |q| → 0. Therefore, we set the magnitude of the zero mode to a small but finite value, |q zero | = ε, introducing an effective infrared cutoff. We have checked explicitly that variations of ε over several orders of magnitude from ε = 1 MeV down to ε = 10 −7 MeV lead to no noticeable differences in our results.
(ii) The DSE kernels depend not only on the internal loop momentum q but also on the external momentum p and their dot product, which contains information on directions. The evaluation of this expression has to be modified if either p or q corresponds to a zero mode. In case of an internal zero mode, we set p · q = 0. On the other hand, in case of an external zero mode, the angular information contained in p · q is important for the spherical sum over the loop momentum q with multiplicity index m in Eq. (9). Therefore, in this case, we use |p zero | = ε and employ the same directions/angles for the zero-mode momenta as for the first finite momentum shell.
As already mentioned above, we refer to our setup with periodic boundary conditions including the zero mode as PBC and compare later with results with discarded zero mode denoted by PBC * . Of course, the numerical effort for PBC is increased compared to PBC * . Thus, while we will show results for PBC * in both setups-with and without UV improvement-we will provide results for PBC only in the higher-quality UV-improved setup.

B. Dyson-Schwinger equations
In order to obtain the dressed quark and gluon propagators, Eqs. (1) and (2), at nonzero temperature and quark chemical potential, we solve a coupled set of truncated DSEs in Landau gauge. Our setup takes the back reaction of the quarks onto the Yang-Mills sector explicitly into account; for recent works (in infinite volume) see, e.g., Refs. [11,19,20] as well as the comprehensive review [21] and references therein.
Considering the introduction of finite volume as described in the previous section, the dressed quark and gluon propagators are solutions of the DSEs where p = (ω n , p il ); i and l denote the sphere and multiplicity, respectively, of the external momentum, which is discrete, too. Here, D YM νσ denotes the sum of the inverse bare gluon propagator and all diagrams with no explicit quark content while is the bare quark propagator with wave function renormalization constant Z f 2 , mass renormalization constant Z m f , and renormalized bare quark mass m f . The quark self-energy Σ f and the quark-loop contribution Π νσ to the gluon self-energy read and where q = (ω k , q jm ) denotes the loop momentum. Furthermore, g is the strong coupling constant,Z 3 the ghost renormalization constant, and Γ f σ the dressed quark-gluon vertex. The prefactors 4/3 in Eq. (14) and 1/2 in Eq. (15) result from the color traces.
To get a closed system of equations, we proceed as follows. First, we replace D YM νσ by quenched temperaturedependent lattice data [22,23]. The explicit form of these fits can be found in the Appendix of Ref. [24]. Second, we use an ansatz for the dressed quark-gluon vertex: the leading Dirac tensor structure γ σ of the Ball-Chiu vertex [25] is supplemented by a phenomenological dressing function that accounts for non-Abelian effects and guarantees the correct perturbative running of the propagators in the UV. The explicit form of the vertex ansatz is detailed, e.g., in Ref. [10]. Since this truncation is identical to the one used in previous works, we refer the reader to Refs. [10,11,21,24] and do not show more detailed expressions in order to keep this section concise.
The set of truncated DSEs obtained in this way is solved self-consistently. This gives us access to the nonperturbative quark and gluon propagators at arbitrary temperature and chemical potential, i.e., across the phase diagram of QCD. Finally, by varying L, we can investigate how the structure of the latter changes in a finite volume. We use 2 + 1 quark flavors that are nontrivially coupled through the quark loop, which is evaluated explicitly. As a consequence, the gluon becomes sensitive to the chiral dynamics of the quarks.
The phase structure is obtained by monitoring the behavior of the quark condensate. It is the order parameter for chiral symmetry breaking and obtained from the quark propagator according to where the factor three stems from the color trace, and the remaining trace is evaluated in Dirac space. The quark condensate is divergent for all flavors with a nonzero bare mass because it contains a term proportional to m f Λ 2 . In our (2 + 1)-flavor setup, a regularized expression can be obtained by considering the difference which defines the subtracted quark condensate. The divergent part of the strange quark condensate cancels the corresponding one of the up quark condensate when the former is multiplied with the up-to-strange mass ratio. We use the inflection point of ∆ us with temperature to define the pseudocritical chiral transition temperature, We work in the isospin-symmetric limit of two degenerate light quarks (m u = m d , µ u = µ d ) and choose µ s = 0. The baryon chemical potential is then given by µ B = 3µ u . We use an up/down quark mass of m u = m d = 0.8 MeV at a renormalization point of 80 GeV (see Ref. [11] for details), and the mass ratio m s /m u = 25.7 is obtained by applying the Bethe-Salpeter framework of Ref. [26] to our truncation scheme and demanding physical pion and kaon masses in vacuum.

III. RESULTS AND DISCUSSION
Before we start with a discussion of our findings, let us remark how we determine the pseudocritical chiral transition temperature. To damp numerical instabilities in the condensates at finite volumes in temperature direction, we employ a hyperbolic tangent fit that represents the condensate very well up to chemical potentials around the CEP, i.e., in the crossover region. The inflection point of the fit function determines T c . For the sake of consistency, we also apply this fit procedure to the infinite-volume analysis. This, in turn, leads to (very) small changes in the transition temperatures as compared to previous works. For example, within the same truncation scheme we find T c (L → ∞) = (155 ± 1) MeV at zero chemical potential in this work compared to T c (L → ∞) = (156 ± 1) MeV in the infinite-volume calculation of Ref. [11]. We emphasize that the difference is purely technical and very small. In the following, we discuss and compare our findings for finite-volume effects without UV improvement with the ones including the UV improvement as discussed in Sec. II A 2. We studied systems in boxes with edge lengths of L = 3, 4, 5, 6 and 8 fm and compare to the infinite-volume limit that is explicitly calculated using the framework of Ref. [11].

A. Pseudocritical chiral transition temperature at vanishing chemical potential
In Fig. 2, we display the pseudocritical chiral transition temperature T c at vanishing chemical potential for antiperiodic and periodic spatial boundary conditions with and without zero mode for the quarks. For comparison, the infinite-volume result is indicated by a black, dotted line. Comparing both figures, we clearly see the effect of the UV improvement at large volumes. In the unimproved case, our results for both ABC and PBC * suffer from cubic artifacts and overshoot the infinite-volume line. In contrast, the improved results approach the infinitevolume results smoothly. 3 In this case, we also performed calculations at even larger box sizes, but since the results 3 Technically, the lack of a consistent infinite-volume limit in the unimproved calculations can be attributed to two reasons that are hard to disentangle. First, the necessarily rather small UV cutoff, typically around Λ = 10 GeV, in the case of the unimproved framework leads to non-negligible cubic artifacts as already discussed above. In the improved framework, the cutoff can be arbitrarily large and we work with Λ = O(100 GeV). Second, the renormalization point in the unimproved framework is necessarily located at even smaller momenta (8 GeV in our calculations) than the already small cutoff and thus too close to the infrared-momentum are similar to the infinite-volume ones on the per mill level, we did not include them in the plot. At smaller box sizes, T c decreases monotonously. While the decrease is rather moderate down to L ≈ 5 fm, volume effects become much more pronounced for even smaller volumes. For example, at L = 3 fm and PBC * we find that T c is almost halved as compared to infinite volume.
In general, one can observe that quarks with PBC * are much more sensitive to finite-volume effects across all investigated box sizes. This is caused by the missing zero mode and the associated larger infrared cutoff introduced by the discrete momentum grid. From Eq. (8), the ratio of the smallest possible momentum magnitude p min between antiperiodic and periodic boundary conditions is given by p ABC min / p PBC * min = 3/4 ≈ 0.866. From the right diagram of Fig. 2, it becomes apparent that the full PBC results, i.e., with zero mode, resemble closely the ones for ABC. Down to L = 4 fm, the ABC results lie on top of the PBC results with zero mode, and they differ only by around three percent at our smallest investigated box size of L = 3 fm.
For volumes as small as a box size of L = 3 fm, the system begins to enter what is called the epsilon regime in chiral perturbation theory [27]. In this region, the product α = m f V ψ ψ L=∞ f of quark mass, four-volume V , and infinite-volume quark condensate becomes of order one and smaller and chiral symmetry starts to get restored already in the vacuum theory. While the full effect (including critical scaling with α) sets in only at much smaller volumes (L 2 fm), first effects are already seen at our smallest box size of L = 3 fm. region, where medium and finite-volume effects become important. As a consequence, the renormalization constants of the unimproved calculation are contaminated by medium and finitevolume artifacts. For the improved torus, this problem disappears because we are able to renormalize in exactly the same fashion as in infinite-volume calculations, i.e., with a renormalization point located deep inside the perturbative region.

B. Crossover line and critical endpoint
The finite-volume modifications on the phase structure are summarized in Fig. 3. We show the phase diagrams for ABC (upper row), PBC * (center row), and PBC (lower row) of the quarks. For comparison, the infinite-volume crossover line including the CEP is added, too. The left diagrams correspond to results without UV improvement while the diagrams on the right are obtained with UV improvement.
Similar to the results at zero chemical potential, we note the drastic effects of the UV improvement. Whereas the CEPs of the series of larger and larger box sizes do not approach the infinite-volume CEP without improvement (diagrams on the left), they do so after the improvement has been implemented (diagrams on the right). Both the crossover line and the CEP at L = 8 fm are very close to the infinite-volume limit. The remaining discrepancy is within the numerical error of the infinite-volume calculation. 4 Furthermore, the volume-dependent shift of the crossover line and the CEP for increasing box size approach the infinite-volume result uniformly. Thus, while the overall qualitative behavior with and without improvement is the same, quantitative aspects can only be discussed in the improved framework.
For ABC and PBC * , both phase diagrams show a similar trend when the box size is decreased: the CEP moves toward smaller temperatures and larger chemical potentials. In more detail, we find that the increase of its location in µ B direction is larger than the decrease in T direction leading to a flattening of the chiral crossover line; see Sec. III C below where we discuss the volume dependence of the curvature of the crossover line for different boundary conditions. Down to L = 4 fm, the PBC results are generally very similar to the corresponding ABC results. However, the volume-dependent movement of the CEP is slightly slower and the endpoint values are closer to infinite volume. The differences to ABC are around/below the ten-percent level, though. On the other hand, at small volumes we observe a qualitative difference for full PBC including the zero mode: the shift of the CEP in µ B direction inverts. In particular, the CEP in our smallest volume of L 3 = (3 fm) 3 is located at smaller chemical potential than the infinite-volume result. Nonetheless, the crossover lines again become flatter for decreasing system sizes, which is in agreement with ABC and PBC * .
Overall, the temperature dependence of the crossover line and CEP is analogous to that of the pseudocritical chiral transition temperature at vanishing chemical potential discussed above in Sec. III A. Results for box sizes L 5 fm are rather similar to one another for the different boundary conditions, while the phase diagrams for ABC, PBC * , and PBC display markedly different structures for small system sizes L 4 fm. Again, finite-volume effects are much more pronounced for PBC * . Especially, the result for L = 3 fm stands out with a very flat crossover line and a CEP at (µ B , T ) ≈ (1070, 40) MeV with UV improvement.
The volume dependence of the location of the CEP has been investigated also in an FRG treatment within a two-flavor quark-meson model truncation [8]. In this work, the position of the CEP has been extracted from the maximum of the scalar susceptibility and its shift with L has been calculated for PBC between L = 4 fm and L = 10 fm. Compared to the DSE calculation presented here, the infinite-volume CEP of the FRG analysis is generally located at (much) higher chemical potential and lower temperature since its precise location depends on the chosen infrared input parameters, in particular on the value of the sigma-meson mass. However, a qualitative comparison of the results of Ref. [8] with the ones presented here yields a satisfying agreement between the present DSE and the FRG findings. Specifically, both the L-dependent relative shift of the CEP as well as the onset of finite-volume effects below L = 8 fm coincide well. Below L = 4 fm, the CEP disappeared completely in the FRG framework.

C. Curvature of the chiral crossover line
Finally, we discuss the volume dependence of the curvature of the crossover line. At small baryon chemical potential µ B , the crossover line can be parameterized as Here, T c (µ B ) and T c = T c (0) are the pseudocritical chiral transition temperatures at nonzero and vanishing chemical potential, respectively, while the coefficient κ 2 is the curvature of the transition line. We obtain κ 2 by fitting the values of the crossover line at small baryon chemical potentials, µ B ≤ 240 MeV, to the parametrization given in Eq. (19). For comparison, lattice calculations yield a curvature in the range 0.0120 ≤ κ 2 ≤ 0.0153 [28][29][30][31].
Our results obtained with the UV improved framework are shown in Fig. 4. As already apparent from the phase diagrams (see Fig. 3), the curvature is consistently smaller than the infinite-volume result and decreases for smaller box sizes. Overall, this flattening resembles the volume dependence of the pseudocritical chiral transition temperature. That is, the results for L = 8 fm are closest to the infinite-volume value and drop monotonously with decreasing L for all boundary conditions. Compared to the pseudocritical chiral transition temperatures shown in Fig. 2, we find that the curvature displays a somewhat stronger reaction to finite volume. Whereas the ABC temperatures are already close to the infinite-volume result for L = 6 fm, the curvature parameter κ 2 for both ABC and PBC * is still off by more than ten percent. Only for very large box sizes of L 8 fm, we observe agreement with the infinite-volume limit within errors. Here, it is important to note that the fit is quite sensitive to details in the input data and choices of fit intervals, such that κ 2 can only be extracted within a margin of several percent.
In contrast to ABC and PBC * , the PBC curvature is very close to the infinite-volume result already above L 6 fm. For L = 5 fm and below, the PBC results are again very similar to the ABC ones with slightly less pronounced finite-volume effects for L = 3 fm.
The curvature for PBC between L = 2 fm and L = 5 fm was studied with FRG techniques in Ref. [3]. Above L ≈ 3 fm, there is qualitative agreement with our results: the curvature increases with L. However, for smaller box sizes an interesting discrepancy occurs. In our case, we find a monotonic decrease for smaller and smaller box sizes, whereas the FRG results show an increase of the curvature when L gets smaller than L ≈ 3.5 fm, resulting in an overall nonmonotonic behavior. Even though this increase occurs for L lower than we have investigated here, the sharp drop of κ 2 for L = 3 fm appears to contradict such a scenario in our calculations. The reason for this deviation is unknown. While the PBC zero mode is attributed to be the driving force in the small-volume limit of Ref. [3], its inclusion does not qualitatively change the behavior of the curvature at small L in our case. However, there are indications that these deviations might be rooted in truncation/model artifacts, and we comment on this in the Appendix.

IV. SUMMARY AND CONCLUSIONS
In this work, we studied the effects of a finite, uniform, three-dimensional cubic volume with equal edge lengths L and (anti)periodic boundary conditions on the phase diagram of QCD. To this end, we employed a well-explored truncation scheme for the coupled set of DSEs for the Landau gauge quark and gluon propagators and extracted the volume dependence of the chiral order parameter. In order to extract the volume behavior of the system, we found two technical procedures to be mandatory: first, the removal of cubic artifacts due to a UV improvement of the setup and, second, the explicit inclusion of the zero mode for periodic boundary conditions. For both types of boundary conditions, periodic and antiperiodic, we then find similar and only moderate volume effects of the order of ten MeV and smaller for box sizes L 5 fm. Only for very small volumes sizable shifts of the CEP and the associated crossover line occur. These shifts are almost monotonous: smaller volumes correspond to smaller transition temperatures and the CEP shifts toward larger chemical potential. The only deviation from this general behavior occurs for periodic boundary conditions at very small box sizes. Our findings are consistent with corresponding results from an FRG treatment of the quark-meson model.
In order to bridge the gap toward direct applications in the context of heavy-ion collisions, (at least) two extensions of the current framework are necessary. First, one needs to adapt the implementation of the finite-volume boundary conditions to the actual physics case. There are no (anti)periodic boundary conditions when two nuclei collide. The situation rather resembles the one of a sphere (at large centrality) or an almond-shaped volume (at small centrality) with a fuzzy boundary. This situation may be better represented in a calculation in a sphere with MIT boundary conditions, as has been advocated in Ref. [18] within a DSE model framework. Second, one needs to address the volume dependence of quantities such as fluctuations of conserved charges [4,32]. Work in this direction within our framework is in progress.

Appendix: Zero-mode contribution in different truncations
There have been discussions about the influence of the zero mode. In model calculations, one can find that its inclusion leads to an increase of chiral symmetry breaking for small box sizes L 3 fm compared to large volumes, see, e.g., Refs. [3,9,32,33]. In particular, the amount of chiral symmetry breaking for small box sizes is found to be significantly larger than in infinite volume.
In this work, however, we find that PBC with zero mode behaves qualitatively and quantitatively almost identically to ABC. In order to investigate this discrepancy, we also performed vacuum calculations for PBC with zero mode in two simpler truncations: the quenched version of the truncation scheme used in this work and one with a model gluon propagator. Both take no backcoupling of quarks onto the gluon into account, i.e., the gluon can be seen as a static input for the quark DSE, and unquenching effects are either absent (quenched gluon) or modeled (model gluon). In more detail, the former truncation is the one discussed in Sec. II B, but the gluon in the quark DSE is solely given by our fits D YM νσ to quenched lattice data, while for the model gluon we resort to the well-established Maris-Tandy model [34]. Furthermore, within both truncations, we also varied the vertex ansatz by employing both the Ball-Chiu-inspired vertex used in this work and the bare vertex Γ f (bare) σ = Z f 2 γ σ . In Fig. 5, we show the L dependence of the subtracted quark condensate [Eq. (17)] in vacuum normalized to its infinite-volume value for both truncations with the two vertex ansätze, respectively. We investigated system sizes in the range of L = 2-8 fm since visible volume effects generally occur at smaller L for these truncations. In the large-L limit, we notice that all results tend to the infinite-volume value regardless of truncation and vertex.
For small L, however, especially below L 3 fm, the qualitative behavior depends crucially on the vertex ansatz. In case of the bare vertex, we are able to reproduce the effects seen in model calculations: an increasing condensate with decreasing system size in both truncations. In case of the Ball-Chiu-inspired vertex, though, the L dependence of the condensate is in line with our findings in this work.
Therefore, we suppose that the zero-mode effects in the small-L limit seen in Refs. [3,9,32,33] are either model artifacts or artifacts induced by approximations (e.g., mean field) within the models. [1] X. Luo and N. Xu, Search for the QCD critical point with fluctuations of conserved quantities in relativistic heavy-ion collisions at RHIC: an overview, Nucl. Sci. Tech.