Deconfinement critical point of lattice QCD with Nf “ 2 Wilson fermions

Francesca Cuteri,1, ̊ Owe Philipsen,1, 2, : Alena Schön,1, ; and Alessandro Sciarra1, § Institut für Theoretische Physik, Goethe-Universität Frankfurt Max-von-Laue-Str. 1, 60438 Frankfurt am Main, Germany John von Neumann Institute for Computing (NIC) GSI, Planckstr. 1, 64291 Darmstadt, Germany The SUp3q pure gauge theory exhibits a first-order thermal deconfinement transition due to spontaneous breaking of its global Z3 center symmetry. When heavy dynamical quarks are added, this symmetry is broken explicitly and the transition weakens with decreasing quark mass until it disappears at a critical point. We compute the critical hopping parameter and the associated pion mass for lattice QCD with Nf “ 2 degenerate standard Wilson fermions on Nτ P {6, 8, 10} lattices, corresponding to lattice spacings a “ 0.12 fm, a “ 0.09 fm, a “ 0.07 fm, respectively. Significant cut-off effects are observed, with the first-order region growing as the lattice gets finer. While current lattices are still too coarse for a continuum extrapolation, we estimate mπ « 4GeV with a remaining systematic error of „ 20%. Our results allow us to assess the accuracy of the leadingorder and next-to-leading-order hopping expanded fermion determinant used in the literature for various purposes. We also provide a detailed investigation of the statistics required for this type of calculation, which is useful for similar investigations of the chiral transition.


I. INTRODUCTION
For physical quark mass values, the thermal QCD transition is known to be an analytic crossover [1]. Its continuation to small baryo-chemical potentials has been studied in detail by means of Taylor expansion or analytic continuation from imaginary chemical potential, with so far no hints of a non-analytic phase transition [2]. Since a severe sign problem precludes lattice QCD simulations at finite baryon density, one way to constrain the phase diagram is to study the thermal transition with the QCD parameters varied away from their physical values. Non-analytic phase transitions related to the spontaneous breaking of the global center and chiral symmetries, respectively, are explicitly seen in simulations employing unimproved fermion discretizations on coarse lattices, in the heavy and light quark mass regime, as indicated in Figure 1. One can then study how these critical structures evolve when a chemical potential is switched on, for an overview see Ref. 3.
In this work we focus on the heavy mass corner, whose thermodynamics can be addressed also at finite baryon chemical potentials, either by means of effective lattice theories obtained from hopping expansions [4][5][6], or by effective Polyakov loop theories in the continuum [7,8]. In order to assess the accuracy of these approaches and their possible extensions to light quarks, reliable benchmarks are warranted. Moreover, the phase transitions in this parameter regime are interesting in their own right. In the infinite quark mass limit, the theory reduces to SU p3q cuteri@itp.uni-frankfurt.de : philipsen@itp.uni-frankfurt.de ; schoen@itp.uni-frankfurt.de § sciarra@itp.uni-frankfurt.de pure gauge theory with its first-order transition [9], which is caused by the spontaneous breaking of the Z 3 center symmetry. Recently, the latent heat associated with this transition has also been determined [10]. When dynamical quarks are added to the theory, the center symmetry is broken explicitly and the phase transition weakens, i.e, the latent heat decreases until it vanishes at a critical quark mass. The value of the critical quark mass and the latent heat are thus intimately related by the dynamics of the deconfinement transition, so that valuable non-perturbative insights are possible in this parameter region.
An early lattice investigation for N f " 1 Wilson fermions was restricted to N τ " 4 and, for large volumes, required mapping to an auxiliary model [11]. The first systematic studies in the standard Wilson discretization [12,13] on N τ P {4, 6, 8} lattices are based on an effective potential for the order parameter determined by a histogram method [5], which is reweighted with a hopping expanded quark determinant from quenched configurations. This allows for a flexible investigation and interpolation between the N f P {1, 2, 3} cases. In this work we focus on N f " 2 only, but simulate the full fermion determinant with no approximation beyond the lattice discretization. Preliminary results for N τ P {6, 8} have been reported in Refs. 14 and 15. Here we significantly increase statistics and add a third lattice spacing with a set of simulations on N τ " 10 lattices. In agreement with earlier work, a shift of the deconfinement critical point towards smaller bare quark masses is observed as the lattice is made finer, Figure 1 (right). On the other hand, on the finer lattices we observe a quantitative deviation regarding the location of the critical point with respect to the hopping expanded results [13]. It is interesting to observe that the direction of the cut-off effect in the bare parameter space is the same for the critical de- confinement boundary in the heavy quark mass regime as for the critical chiral boundary in the light quark mass regime with both Wilson [16][17][18] and staggered [19][20][21][22][23] discretizations (in some other studies employing improved staggered fermion discretizations [24,25], no chiral critical line is seen, thus bounding a potential chiral first-order region). The shift of the chiral Z 2 boundary implies an increase in the simulation cost while the continuum is approached, which is particularly drastic in the chiral transition region. This further motivates to attempt a continuum limit in the heavy quark mass regime first, where it should be more feasible.
After devoting Section II to the description of our lattice setup and of the relevant symmetry at work in the infinite mass limit, we discuss our finite size scaling analysis in Section III. Details on the simulations and the analysis strategy are provided in Section IV, followed by a critical appraisal of the growing statistics requirement for decreasing lattice spacings in Section V. Finally, our results for the deconfinement critical point are reported in Section VI and conclusions are drawn in Section VII.

II. LATTICE ACTION AND CENTER SYMMETRY
We work with the standard Wilson gauge action with the lattice coupling β " 2N c {g 2 , the plaquette P µν pnq and n labeling the lattice sites. The standard Wilson fermion action for N f mass-degenerate quarks is defined as with the fermion matrix where γ´µ "´γ µ . The bare fermion mass m is adjusted via the hopping parameter The lattice coupling β controls the lattice spacing apβq, and temperature is defined as We do not use any improvements on the Wilson discretization to make sure the phase structure does not get modified by any unknown effects of improvement terms. E.g., the twisted mass formulation introduces additional unphysical phases, which then require care to be avoided [26]. Quite generally, approaching the continuum limit in the right way requires knowledge of the lattice phase diagram in the bare parameter space, for which we stick to the simplest action. The Polyakov loop Lpnq is defined by the product of all temporal links at space point n, closing through the periodic boundary, Physically, it describes the (Euclidean) time evolution of a static quark. At finite temperature, the periodic boundary in the temporal direction permits topologically non-trivial gauge transformations, which are twisted by a global center element of the group when crossing the boundary, gpτ`N τ , nq " z gpτ, nq , g P SU p3q , z " e i 2πk 3 1 , with k P {0, 1, 2}. The pure gauge action is invariant under such transformations, S G rU g s " S G rU s, while the Polyakov loop picks up the twist factor, In the pure gauge theory and in the thermodynamic limit, the expectation value L is therefore an order parameter for center symmetry breaking, with non-zero values signaling deconfinement [27]. In the presence of dynamical quarks, S F rU g s ‰ S F rU s and the center symmetry is explicitly broken, i.e., L ‰ 0 always, with 1{m playing the role of the symmetry-breaking field. Nevertheless, a rapid change accompanied by large fluctuations of L still signals the deconfinement transition, which weakens with decreasing quark mass until a critical point is reached.
Being the endpoint of a first-order transition, this critical point will be in the same universality class as 3d liquid-gas transitions, i.e., the 3d-Ising model. This is fully confirmed by our data. In the effective 3d-Ising Hamiltonian governing the vicinity of the critical point, the Polyakov loop will then be the dominant contribution to the magnetization-like variable, coupling to the symmetry-breaking magnetic field 1{m. It is interesting to contrast this situation with the chiral transition in the light quark regime, where the Polyakov loop contributes to the energy-like variable in the corresponding effective Hamiltonian, since it is invariant under chiral transformations [28].

III. FINITE SIZE SCALING ANALYSIS
To determine the order of the phase transition as a function of the quark mass, we compute standardized moments of X " |L|, the absolute value of the Polyakov loop, and study their behavior as a function of the physical volume. Here n P {3, 4} and {α n } is a set of relevant physical parameters, in our case β, κ, N τ , N s (in the following, the dependence of B n on X will be understood).
The phase boundary at a pseudocritical coupling β c is defined by the vanishing of the skewness B 3 pβ c q " 0. We use finite size scaling of the fourth moment, the kurtosis B 4 , trivially linked to the Binder cumulant [29], to locate the critical point κ Z2 . For V Ñ 8, B 4 evaluated on the phase boundary takes the form of a step function, assuming values that are characteristic for the respective nature of the phase transition (c.f. Table I). On finite volumes, B 4 is an analytic curve that gradually approaches the step function with increasing volume. In the vicinity of the critical point κ Z2 , universality implies it to be a function of the scaling variable x " pκ´κ Z2 qN 1{ν s only, Table I. Values of the kurtosis at the transition and of the relevant critical exponents [30].
which can be expanded around κ " κ Z2 in a Taylor series. For sufficiently large volumes, the series can be truncated after the linear term, and fitted to extract the critical parameters. This simple picture holds for asymptotically large volumes, and previous studies show that these can be prohibitively expensive to attain, cf. [31,32]. Moreover, we find the required aspect ratios in the heavy mass region to be even significantly larger than those in our light quark studies. The reason cannot be related to the lightest state in the spectrum, the glueball in this case, whose correlation length is shorter than that associated with the pion in the light quark mass regime. Instead, a possible explanation is that regular (non-divergent) contributions to the fluctuations grow towards the heavy mass region, so that it takes larger volumes for the diverging terms to dominate.
It is therefore expedient to also consider the leading finite volume corrections to Eq. (10). Since for finite quark masses the center symmetry is explicitly broken, the Polyakov loop is no longer a true order parameter, but a mixture of the magnetization-and energy-like observables entering the effective 3d-Ising Hamiltonian which governs the vicinity of the critical point. Magnetizationand energy-like observables scale with different exponents, and only for asymptotically large volumes the larger of them dominates to produce the simple expression Eq. (10). With the leading mixing correction taken into account, the kurtosis takes instead the form where B is another non-universal parameter to be extracted from a fit. For a detailed derivation and explanation in the context of the light quark mass regime, see [17]. As a consequence of this correction, the B 4 curves are no longer fitted by straight lines, and different volumes intersect at different points, with the pairwise intersection approaching the universal value with growing volumes. This is shown schematically in Figure 2.

IV. SIMULATION AND ANALYSIS DETAILS
We study cut-off effects by simulating three different temporal lattice extents N τ P {6, 8, 10}, corresponding Qualitative behavior of the kurtosis B4 for spatial volumes too small to be described by Eq. (10). The effect of the leading finite size correction, Eq. (11), is to shift B4 at κZ 2 (i.e. for x " 0) to larger values and approach the universal value with growing volumes. The enlargement shows that different volumes do not cross in one point. Pairwise intersections converge to B4pβc, κZ 2 , 8q in the thermodynamic limit. β Raw data from multiple chains Raw (merged) and reweighted data Reweighted Raw merged Figure 3. Example of the data analysis for B3, B4 at κ " 0.1100 on a 6ˆ36 3 lattice. In the left plots different Markov chains are displayed, with data points slightly shifted horizontally for better visibility. The label nσ denotes by how many standard deviations the maximally incompatible pair is apart. On the right, the merged raw and reweighted data for B3 (top) and B4 (bottom) are shown. The determination of βc and B4pβcq is depicted in red.
to different lattice spacings at the respective critical couplings. For every value of N τ , the critical κ Z2 in the heavy quark mass region was determined (i.e. the position of red circle in Figure 1) for five to six values of κ. At each value of κ up to four spatial volumes have been simulated, corresponding to aspect ratios of N s {N τ P r4, 7s. For N τ " 8, also simulations at N s " 80 have been done to give better insights into the position of κ Z2 . The decon-finement transition has been located simulating at two to four different values of β around the pseudocritical coupling. Configurations have been produced using a standard Hybrid Monte Carlo algorithm [33] with unit trajectories and the acceptance rate tuned to stay between 75% and 85%. At least 5000 thermalization trajectories have always been discarded. Observables (the plaquette and the Polyakov loop) have then been measured on the fly for all Monte Carlo steps. At each value of β, between 56k and 800k trajectories have been accumulated, making sure to always have over « 50 independent events per β. In total 28, 39.5 and 21.8 millions trajectories have been produced for N τ " 6, N τ " 8 and N τ " 10, respectively. Details about simulation parameters and statistics are listed in Tables IV to VI. All the finite temperature simulations, except from those on the 8ˆ80 3 lattices, have been performed using v1.0 of our publicly available [34] OpenCL based CL 2 QCD code [35], which is optimized to run on GPUs. The L-CSC supercomputer [36] at GSI in Darmstadt has been used for this set of runs. The few simulations on aspect ratio 10 lattices have been performed with openQCD-FASTSUM [37], a highly MPI-optimized software that has been run on the Goethe-HLR supercomputer. In all cases, monitoring and handling of thousands of jobs has been efficiently automatized using the software package BaHaMAS [38,39], whose new version v0.3.1 can also be used in conjunction with the openQCD-FASTSUM code.
In order to faster accumulate statistics and to gain better control over statistical errors, for each parameter set we simulated four different Markov chains, except at the outermost β values for aspect ratio 10 where only three Markov chains were run. In this way a criterion to decide when the gathered statistics is large enough can be established: All replicas included in the final analysis were run until B 3 was compatible within at most 3 standard deviations between all of them. An example of one dataset can be found in Figure 3 (left). At the beginning of the simulations, for very low statistics, this criterion can be trivially fulfilled due to very large statistical errors. Therefore, to ensure a proper bracketing of β c and to avoid stopping the simulations too early, it has been almost always demanded that B 3 is incompatible with zero at 1 standard deviation at the smallest and largest β on every chain. Once both these guidelines are satisfied, the chains were merged to the raw data points shown in Figure 3 (right).
In order to have precision on β c , the multi-histogram method [40], also known as reweigthing, was used to interpolate between our measurements. This interpolation is repeated for the B 4 data, thus pinning down B 4 pβ c q. Carrying out the same procedure for each simulated value of κ, the resulting B 4 pβ c , κ, N s q data points can be plotted as a function of κ and fitted according to Eqs. (10) and (11), as it will be discussed in Section VI.
Finally, to set the physical scale for our simulation points we use the Wilson flow parameter ω 0 {a, which we determined and converted using the publicly available software described in Ref. 41. To this end, we produced 800 independent zero-temperature configurations on 32ˆ16 3 lattices for each κ at the value of the critical coupling β c . As a physical quantity to parametrize the determined critical point κ Z2 , we then computed the pseudoscalar meson mass m π corresponding to those bare parameter values. Note that for all lattice spacings considered here, am π ą 1 at and around κ Z2 . Therefore, our  pion mass estimates are afflicted by large cut-off effects. This problem naturally reduces as N τ is increased. All critical couplings β c , the corresponding lattice spacings a, m π and critical temperatures T c are summarized in Table III for each value of κ and N τ .

V. STATISTICS REQUIREMENTS TOWARDS THE CONTINUUM
A crucial parameter to judge the statistical quality of the analysed ensembles is the integrated autocorrelation time τ int of the observables. Qualitatively speaking, τ int parametrizes the memory that the simulated system has of its dynamics, which in our case strongly depends on the order of the phase transition. Consider a first-order transition to start with and let us briefly summarize what happens to the behaviour of the order parameter. In a finite volume, far away from the phase transition, the system stays in a given phase. The order parameter will fluctuate around its mean value and its probability distribution will be approximately Gaussian, with some associated τ int . As soon as β gets sufficiently close to β c , the system will explore both phases and the order parameter will jump from time to time from fluctuating around its mean value in one phase to fluctuating around its mean value in the other one. These fluctuations between the phases are slower than the fluctuations within one phase, since tunneling between different phases is exponentially suppressed by a potential barrier growing with volume [42]. As a result, the system has a much longer-term memory since its dynamics now occurs on a larger timescale. τ int is related to the average tunneling rate between the two phases and thus increases significantly. This in turn implies that sufficiently many tunneling events in a simulation are needed to reliably estimate τ int . At a crossover transition, instead, the distribution of the order parameter does not show a twopeak structure for any β value, even around β c , where only its variance slightly increases. Therefore, τ int is expected to be relatively small in this case. Finally, at a second-order phase transition the critical slowing down phenomenon will play an important role and the diverging correlation length of the system will be reflected in an increase of τ int , too. In a finite volume this effect will be only partially felt, though, because the correlation length cannot literally diverge.
We estimate τ int for both B 3 and B 4 using a PYTHON implementation of the Γ-method [43], and distribute our data in bins of size 2τ int to remove autocorrelations. Tables V and VI give an overview of τ int for the skewness and the kurtosis (averaged over the simulated β values) and the number of statistically independent events obtained from the binning procedure. The qualitative expectation outlined above is completely met in Figure 4, where for each pκ, N s q pair, τ int is displayed at the simulated β closest to β c . In Figure 4(a) the autocorrelation time is shown as a function of κ and spatial volume for N τ " 8. For each volume, a maximum is seen around κ Z2 , showing critical slowing down near a second-order transition (that maximum would be sharper if τ int had been evaluated exactly at β c ). The more drastic effect, however, is the increase of τ int with increasing volume. In Figure 4(b) we compare different N τ values at fixed aspect ratio. As expected when approaching the continuum limit, we observe another increase of τ int around the critical as well as the first-order region. The statistics for N τ " 10 is effectively smaller compared to N τ P {6, 8}, which explains the larger error bars and possibly an overestimate in the first-order region.
Note also, that the observed rise in the autocorrelation time feeds back into the practical organization of the simulation. The choice of β-values to be simulated is an op-timization of having a sufficiently narrow spacing needed for reweighting, and covering a large enough interval to bracket β c , with as few values as possible. This requires monitoring and frequently analysing running simulations in order to optimally adjust the parameters. Since a given statistics which turned out to be sufficient on a coarser lattice, will have to be increased on a finer one, even this process of parameter tuning requires more and more trajectories as the lattices get finer. Altogether, our study of the autocorrelation time illustrates the crucial importance and formidable difficulties of obtaining sufficient statistics for a reliable determination of the order of the phase transition as the continuum is approached.

VI. RESULTS AND DISCUSSION
We are now ready to extract our desired results from fits of our data to Eqs. (10) and (11), measuring the fit quality by the reduced chi-square χ 2 d.o.f and the Qparameter. The latter amounts to the probability of getting a value of chi-square larger than χ 2 d.o.f , assuming that the data have Gaussian noise, and gives a measure of the quality of the fit (its optimal value is 50%). In the ideal situation of negligible finite size effects, we expect to obtain a good fit to Eq. (10) and, simultaneously, a good fit to Eq. (11) with a consistent κ Z2 and B compatible with zero. The general strategy then is to compare the fits performed with and without the correction term described in Eq. (11), with the goal to isolate the leading terms. For this purpose, it is in principle enough to successively exclude points from the smallest volumes in the fit to Eq. (11), until the coefficient B is compatible with zero and check for consistency of the fit parameters.
However, we generally observe that the inclusion of all volumes corresponding to the smallest κ-values, i.e., the ones deepest in the first-order region, significantly deteriorates the fit quality, irrespective of the ansatz. Because tunneling between the phases is exponentially suppressed by volume in the first-order region, these points are increasingly difficult to determine accurately, as discussed in the previous section. Moreover, the entire B 4 curve is asymmetric about its inflection point, because its two asymptotes are not symmetrically displaced with respect to B 4 pβ c , κ Z2 , 8q. As the entries in Table III illustrate, the masses associated with the smaller κ-values rise very quickly and are thus quite far from the critical quark mass in which we are interested. Including these smallest κ-values would therefore require higher order terms in the Taylor expansion of the kurtosis, in both brackets of Eq. (11), in order to account for the larger distance from the critical point as well as for the larger finite size corrections. This would increase the number of fit parameters and weaken the conclusions from the fits, so we excluded these points. On the other hand, linear fits over ranges in κ that only cover the crossover region can result in biased estimates for κ Z2 (cf. the shift in κ Z2 comparing the two e.8.# and i.8.# fits at N τ " 8 represented in   Table II. Shaded points have not been included in the fits. Table II). Consequently, we made sure that for each N τ we have included at least one κ in the fits, which reliably belongs to the first-order region.
Distinctive features for κ values in the first-order region are the reversed N s ordering of the corresponding kurtosis (central) values with respect to what happens in the crossover region (the value of the kurtosis decreases with N s in the first-order region), as well as a more pronounced two-peak structure in the distribution of the order parameter, which we explicitly checked in every case. In fact, all data points with κ ă κ Z2 fulfill these criteria, including the ones omitted from our final fits. For N τ " 8 and κ " 0.115, a parameter set in doubt, we simulated aspect ratio 10 specifically to check these features. As can be seen in Figure 5(b), no hints for a first-order phase transition were found which is consistent with the fit indicating κ Z2 ă 0.115. In our general strategy it is important that the fit gives a κ Z2 larger than the small-est simulated κ value, and that the κ Z2 extracted from the fit is cross-checked against evidence for a first-order phase transition for all simulated κ ă κ Z2 .
Some representative fits are shown in Table II, where the best ones chosen as our final result are highlighted. Note that for N τ " 6, 8 we are able to find good fits with B " 0, as well as several consistent ones including additional data points and B ‰ 0. This indicates that all simulation data are indeed described by Eq. (11) and the coefficient B is fully controlled. Our data for B 4 together with the best fits highlighted in the table are also shown graphically in Figure 5.
Our final results for the critical hopping parameter κ Z2 as a function of lattice spacing are collected in Figure 6 (top). Strong cut-off effects are apparent. To compare with other approaches and get a feeling for the physical scales involved, Figure 6 (bottom) shows the critical couplings converted to pseudoscalar meson masses. As   Table II. An overview of the outcome of the final fit analysis is presented here. For each value of Nτ the effect of excluding some data points is shown. The fits are labeled by x.y.z, where x refers to the fit ansatz, y is the value of Nτ and z simply a counter. For x " e the fit has been done according to Eq. (10) (excluding the correction term, B " "´"), while for x " i Eq. (11) has been used and the correction term has been included. z " 1 always represents the fit with the least excluded data points, for z ą 1, more and more data points get excluded. The rows with blue background contain the best fit as a compromise between all parameters. Subtables in the second column show which κ (columns) have been included () or excluded () for which Ns (rows). Gray cells denote simulations that were not done. already indicated in Section IV, these numbers have to be taken with care because am π ą 1 in all cases. This problem could in principle be circumvented by heavy quark effective theory methods [44]. However, it is apparent that at least two finer lattice spacings are needed before a continuum extrapolation can be attempted. For these the mass in lattice units might well be small enough, so we postpone this issue until such data are available. Nevertheless, the shift in the critical pion mass is significantly reduced between N τ " 8, 10 compared to N τ " 6, 8.
A linear extrapolation using the last two points (unimproved Wilson fermions have Opaq effects) would then predict m π « 4 GeV, where twice the shift to the extrapolated value should amount to a conservative estimate of the remaining systematic error of " 20%.
It is now interesting to compare our results with those of Ref. 13, where the critical points for N f P {1, 2, 3} Wilson fermions were determined on N τ P {4, 6} by means of reweighting with a next-to-leading-order (NLO) hopping expanded fermion determinant and a histogram technique instead of the cumulant analysis. For N f " 2 on N τ " 6, N s " 24, Ref. 13 reports κ Z2 " 0.1202p19q at NLO hopping expansion, which translates to m c π {T c « 11.2. Compared with our κ Z2 " 0.0877p9q or m c π {T c « 18.1 (which is obtained from our scale setting summarized in Table III), one observes a large discrepancy of « 50% in the critical pion mass. Note that, since the same lattice action is employed, this discrepancy is not related to the cut-off error on the determination of m c π . As discussed in Ref. 13, the histogram method is applied at a fixed spatial volume and does not include an extrapolation to the thermodynamic limit. Indeed, the authors report a reduction of κ Z2 by « 6% when increasing volume to N s " 32 and by « 13% on N s " 24 when going from LO to NLO in the hopping expansion, i.e., both truncations have a systematic error in the same direction. The comparison to our result highlights the necessity for either refined approximations or a full calculation in order to achieve a reliable continuum limit.

VII. CONCLUSIONS
While qualitatively well understood, the heavy mass (top right) corner of the Columbia plot in Figure 1 (left) still lacks a quantitative determination, in the continuum limit, of the location of the deconfinement critical boundary. In this work we focused our attention on the N f " 2 deconfinement critical point and studied its location on progressively finer lattices simulating at N τ P {6, 8, 10}.
Our results show that the continuum limit is not yet within reach given the extent of the observed cut-off effects. It is, indeed, apparent that for a continuum extrapolation to become feasible at least two finer lattice spacings will be needed. Nevertheless, this work documents the progress that has been made both in refining the fitting strategy for the finite size scaling analysis and in appraising the growing statistics requirements towards the continuum limit. Concerning the former, we showed how a correction term can be used as a probe for finite size effects, which allowed us to isolate the leading terms in the linear kurtosis expansion around κ Z2 and identify the required aspect ratios for the linear regime. With regards to the latter, our detailed analysis of the integrated autocorrelation times of the relevant observable illustrates the alarming prospects in terms of the statistics needed to reliably establish, as the continuum limit is approached, the order of the phase transition and the location of κ Z2 . Furthermore, our results allow for quantitative comparisons with those obtained using hopping expansions and thus to assess their systematic error.

ACKNOWLEDGMENTS
We thank Christopher Czaban for collaboration during the early stages of this project [14,15] and Reinhold Kaiser for providing a very handy GUI for the fitting analysis. The authors acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the CRC-TR 211 "Stronginteraction matter under extreme conditions" -project number 315477589 -TRR 211. We also thank the computing staff of the Goethe-HLR and L-CSC clusters for their support.  Table III. Outcome of the pion mass and the scale setting simulations, performed on 32ˆ16 3 lattices, accumulating 800 independent configurations. The values of βc and the critical temperature Tc of the deconfinment phase transition are also included to provide a more complete overview. a and mπ have also been measured at the κZ 2 -values obtained from the fits and are reported on gray background. The value of βc at κZ 2 , at which the scale is set and the pion mass is measured, has been obtained via interpolation of βc-values at simulated κ nearby.