Finite size and cut-off effects on the Roberge-Weiss transition in $N_\text{f}=2$ QCD with Staggered fermions

In the absence of a genuine solution to the sign problem, lattice studies at imaginary quark chemical potential are an important tool to constrain the QCD phase diagram. We calculate the values of the tricritical quark masses in the Roberge-Weiss plane, $\mu=\imath\pi T/3$, which separate mass regions with chiral and deconfinement phase transitions from the intermediate region, for QCD with $N_\text{f}=2$ unimproved staggered quarks on $N_\tau=6$ lattices. A quantitative measure for the quality of finite size scaling plots is developed, which significantly reduces the subjective judgement required for fitting. We observe that larger aspect ratios are necessary to unambiguously determine the order of the transition than at $\mu=0$. Comparing with previous results from $N_\tau=4$ we find a $\sim50$% reduction in the light tricritical pion mass. The heavy tricritical pion mass stays roughly the same, but is too heavy to be resolved on $N_\tau=6$ lattices and thus equally afflicted with cut-off effects. Further comparison with other discretizations suggests that current cut-off effects on the light critical masses are likely to be larger than $\sim100$%, implying a drastic shrinking of the chiral first-order region to possibly zero.


I. INTRODUCTION
The theoretical prediction of the QCD phase diagram as a function of temperature T and baryon chemical potential µ B has proved to be a difficult challenge for several decades. Because of the non-perturbative nature of the strong interactions on hadronic scales, a first principles approach such as lattice QCD is required. On the other hand, because of the severe sign problem of lattice QCD at finite µ B , standard Monte Carlo simulations are limited to addressing small densities, µ B ă 3T , only [1,2]. Even at zero baryon density, there remain open questions. While the thermal transition from a hadron gas to a quark gluon plasma is well established to be an analytic crossover for physical quark masses [3], the universality class of the transition in the chiral limit of the u, d-quarks is still not settled since it cannot be simulated directly.
For these reasons, it is useful to study the dependence of the thermal transition on QCD parameters like quark masses, numbers of flavors, imaginary chemical potential, for which there is no sign problem, as well as on the lattice spacing. The current knowledge of the nature of the QCD thermal transition as a function of the three light quark masses and imaginary chemical potential, as obtained on coarse lattices with unimproved actions, is sketched in philipsen@itp.uni-frankfurt.de : sciarra@itp.uni-frankfurt.de Figure 1. For large and small quark masses, there are regions with first-order deconfinement and chiral phase transitions, which in the infinite and zero mass limits are associated with the breaking and restoration of the center and chiral symmetries, respectively. These are separated by surfaces of second order transitions from a large region where the transition is merely an analytic crossover, to which also QCD with physical parameters belongs [3,4]. Note that this qualitative picture is the same for unimproved staggered [5,6] and unimproved [7] as well as improved [8] Wilson discretizations, whereas the precise location of the boundary at µ " 0 differs significantly between them, indicating large cut-off effects. These are also observed for N f " 4 staggered fermions without rooting [9]. By contrast, simulations with improved staggered actions do not see any region of first-order chiral transitions within the available mass range, neither at zero [10] nor imaginary chemical potential [4,11] thus providing upper bounds on the critical mass values.
In the present work we continue earlier studies using the unimproved staggered discretization at imaginary chemical potential on finer lattices. In particular, referring to Figure 1, we investigate how the (red) tricritical points on the N f " 2 line in the Roberge-Weiss-plane (bottom plane at pµ{T q 2 "´pπ{3q 2 ) move as the lattice spacing is reduced to " 2{3 of its previous values. Together with similar investigations at µ " 0, this establishes the behavior of the critical surfaces when approaching the continuum. Such studies are complementary to ones with improved actions, where no non-analytic chi-T ri cr it ic al Tric ri ti ca l Figure 1.
Qualitative sketch of the three-dimensional Columbia plot realized on coarse lattices. Whether the chiral, first-order triple region in the Roberge-Weiss plane shrinks on finer lattice enough to make an Op4q-region appear in the m u,d " 0 plane remains unclear.
ral transition is seen, and necessary, if all discretizations are to be understood in the same manner, with expected agreement in an eventual continuum limit. As a byproduct of our study, we develop a new analysis of the finite size scaling of cumulants, which significantly reduces the amount of subjective judgement required for fitting.
In order to render the paper self-contained, we briefly summarize the main features of QCD at imaginary chemical potential in Section II. We then proceed to describe our numerical methodology in Section III and our novel analysis method in Section IV. Our numerical results are given in Section V before we conclude in Section VI.

II. QCD AT IMAGINARY CHEMICAL POTENTIAL
Because of charge conjugation symmetry and its explicit breaking by a non-vanishing baryon density, the QCD partition function is an even function of quark chemical potential, Zpµq " Zp´µq. For purely imaginary chemical potential, µ " ıµ i , µ i P R, it is furthermore periodic [12], and we use N c " 3 colors for the QCD gauge group. These symmetries imply the phase structure shown in Figure 2, with three different Zp3q center sectors, which with ϕ " 2kπ{3, k P {0, 1, 2}. At high temperatures, there are first-order phase transitions between the center sectors, whereas at low temperatures they are analytically connected. The dotted line represents the analytic continuation of the thermal transition, whose order depends on the quark masses. For large and small quark masses, these lines represent first-order deconfinement and chiral transitions, respectively, whereas for intermediate quark masses they correspond to an analytical crossover. Consequently, there are three possibilities for the end-point of the Roberge-Weiss transition: for large and small quark mass it is a first-order triple point, where the thermal first-order transition lines meet that of the center transition. For intermediate quark masses, the thermal transition is only a crossover and the center transition ends in a critical end-point in the 3D Ising universality class. At the boundaries between these situations, corresponding to specific quark mass values, the end-point is tricritical and corresponds to the red boundary points in the Roberge-Weiss plane of Figure 1. The purpose of the present work is to locate these tricritical masses on N τ " 6 lattices with N f " 2 and compare their values with previous determinations on a coarser N τ " 4 lattice [13], as well as with those of other discretization schemes.  Table I. Critical values of ν, γ and B4 " B4pX, . . .q for some universality classes [14].

III. NUMERICAL SETUP
We consider the QCD partition function of N f " 2 mass-degenerate quarks with a purely imaginary chemical potential. After integration over the fermionic fields it can be written as where S g is the gauge part of the action and D is the fermion matrix. For our investigation we used the standard Wilson gauge action and the standard staggered discretization of dynamical fermions. Denoting the lattice gauge coupling by β " 6{g 2 , with the continuum gauge coupling g, and an elementary plaquette by P , we have The fermion matrix reads wherem u,d " am u,d is the quark bare mass in lattice units, a is the lattice spacing, i, j refer to lattice sites, η i,ν are the staggered phases,ν is a unit vector on the lattice andŨ i,ν are the gauge links, which include the purely imaginary chemical potentialμ i " aµ i in the temporal direction,Ũ i,ν " The temperature is specified by the inverse euclidean time extent of the lattice, In order to locate a phase transition and study its nature, we calculate standardized cumulants B n pX, β,m u,d ,μ i q " pX´ X q n pX´ X q 2 constructed from an (exact or approximate) order parameter X. In particular, a non-trivial zero of the skewness of the X-distribution, determines at which value of β " β c a thermal transition takes place, while the value of the kurtosis B 4 in the thermodynamic limit, evaluated at the critical coupling, will determine the order of the phase transition (refer to Table I for common kurtosis values). Note that the so-called Binder cumulant [15], is trivially related to the kurtosis (of the same observable) and contains the same information. We fix µ i {T " π since in this case the imaginary part of the Polyakov loop is an exact order parameter, Referring to Eq. (2), it distinguishes between the low T disordered phase and the high T ordered phase with twostate coexistence, with the advantage of knowing its mean value exactly. On a N τ " 6 lattice, we thus setμ i " π{6. Since this is the boundary between two Roberge-Weiss sectors for all temperatures, B 3 pL Im q " 0 for any value of β and we cannot use Eq. (9) to locate the Roberge-Weiss endpoint. However, the kurtosis B 4 is expected to vary from values close to 3 (crossover) at low T to values close to 1 (first order) at high T . Although it becomes a nonanalytic step function in the V Ñ 8 limit, it is a smooth function on finite volumes, with the curves for different volumes crossing at a universal value for B 4 at the critical point β " β c , provided that the spatial lattice extent is large enough. This crossing provides the location of the Roberge-Weiss end-point. In the neighborhood of the critical point β c , the kurtosis shows a well-defined finite size scaling behavior as a function of the scaling variable Its Taylor expansion around the critical point x " 0 is Sufficiently close to the thermodynamic limit, the coefficient B 4 pβ c , 8q and the critical exponent ν take their universal values depending on the type of transition. In order to locate the two tricritical points in the Roberge-Weiss plane, we performed simulations at different values ofm u,d and different values of β around  the critical temperature. Evaluating the kurtosis in the critical region and fitting it to Eq. (14), considering the linear term only, gives B 4 pβ c , 8q, a 1 , β c and ν for every value ofm u,d . The change of ν as a function ofm u,d then permits to locate the light and heavym tric u,d values. Although our main quantitative analysis is based on the kurtosis of the order parameter, we also calculated the susceptibility of |L|, which is expected to scale around β c according to Here t " pT´T c q{T c is the reduced temperature and f is a universal scaling function. Comparing the collapse plots obtained by fixing the critical exponents γ and ν to the first-order or second-order values, and by plotting χ{N γ{ν s evaluated on different lattice sizes against tN 1{ν s also provides information about the nature of the thermal transition and serves as a cross-check of the kurtosis analysis. A similar cross-check using the chiral condensate ψ ψ was occasionally performed in the small-mass region, leading to consistent conclusions. We investigated 19 values ofm u,d in the intervals r0.004, 0.011s and r0.15, 0.85s. For each value ofm u,d , three to five spacial lattice sizes have been used, keeping N τ " 6 andμ i " π{6 fixed. This corresponds to aspect ratios N s {N τ P r2, 7s. Larger spatial volumes than initially chosen were added whenever the kurtosis of the order parameter on different volumes was not crossing at the same point (an example is reported in Figure 3). For every lattice size, between three and seven values of β around the critical temperature have been simulated. In between those, the observables have been evaluated at additional β-values using the Ferrenberg-Swendsen multiple histogram method (also known as reweighting technique) [16] to increase resolution (see the discussion at the end of this and in the next section).
Configurations were generated by a standard RHMC algorithm [17], producing four different Monte Carlo chains per β with unit-length trajectories. Where advantageous (m u,d ă 0.007), the multiple-pseudofermions technique [18] has been used. The algorithm acceptance has been tuned to be not lower than 80%. At least 5k trajectories were always discarded as thermalization and afterwards observables of interests (i.e. the plaquette, the Polyakov loop and the chiral condensate for small masses) have been computed for every trajectory. We increased statistics until the standard deviation of the kurtosis B 4 pXq decreased below " 0.2, and B 4 pXq was the same on all the four chains at the same β within two (three) standard deviations in the large (small) mass region. For this reason the collected statistics per β is not uniform, detailed information is given in the Appendix. In order to satisfy these strict requirements, millions of trajectories per volume and almost half a billion in total have been produced. This large statistics is necessary due to the large autocorrelation times shown by L Im , especially when entering the first-order regions (cf. Table III in the Appendix). We always insisted on having at least 100 independent events per β in the analysis.
In order to determine the lattice spacing and the pion mass, also zero-temperature simulations have been performed. We produced 800 independent configurations on 32ˆ16 3 lattices for each value ofm u,d . Using the publicly available code described in Ref. 19, the scale was set by the Wilson flow parameter w 0 . Pion masses, instead, were measured with standard spectroscopy techniques [20,21].
All our numerical simulations (except those for scale setting purposes) have been run using the publicly available [22] OpenCL based code CL 2 QCD [23], which is optimized for GPUs. The L-CSC [24] supercomputer at GSI in Darmstadt has been used, and the thousands of jobs needed in the study have been efficiently handled using the simulation monitoring package BaHaMAS [25]. Our quite intricate fitting procedure used to extract the critical exponent ν is completely analogous to the one previously described in Appendix B of Ref. 26. For each value of the quark mass, nearly all possible fits of the data to the linear part of Eq. (14) are performed, and a filtering procedure is applied afterwards in order to pick the best fits. This is needed because the range in which B 4 pL Im q can be considered linear for each N s value is not known a priori. However, here we differ in a technical detail from the previous study in Ref. 26. As already mentioned, reweighting was used not only to smoothen the signal, but also to supply additional β-points for the fit. This approach allows to reduce the number of required simulations, provided there are clear criteria according to which such points are added. Reweighted points introduce a correlation with the others, and too many of them would render the fits unreliable. Hence, we added more reweighted points between simulated points only if with lower resolution it was not possible to obtain a good fit. Moreover, another important aspect should be considered in choosing the reweighting resolution. The β-region where the kurtosis is linear shrinks on larger volumes. Thus, choosing the same resolution in β on different N s would imply to include fewer points from larger volumes in the fit and, consequently, to enhance finite size effects. Therefore, we increased the reweighting resolution in β on larger N s , making the information coming from the smallest volume systematically less important, see Table V for detailed overview. This aspect is even more explicit looking at how many points per N s have been included in the fit.

IV. AN ALTERNATIVE TO FITTING: QUANTITATIVE COLLAPSE PLOTS
Any fitting procedure, however careful, relies on a few subjective decisions like the number of reweighted points and the filtering parameters to judge a good fit. We now propose an alternative procedure to independently determine the critical exponent from scaling/collapse plots, which can then be compared with the results of the fitting procedure.
A collapse plot is obtained if an observable, which displays universal finite size scaling, is plotted as a function of its scaling variables, such that the curves for different volumes fall on top of each other provided the volumes are sufficiently large to represent the thermodynamic limit. There are several common observables used for this purpose, here we focus on B 4 pL Im q as function of the scaling variable x defined in Eq. (13). An example is shown in Figure 5.  Whenever the lattice volume is not large enough, scaling is violated and no good collapse is obtained, even when the known critical values (listed in Table I) are used. Finite size corrections are responsible for that and, in principle, a better collapse can be obtained using different (non-universal) values of the exponents. The quality of the collapse is usually judged by eye, which is mostly sufficient to distinguish between a first and a second order phase transition with known exponents (like in Figure 4). However, a more rigorous method is clearly needed in a situation where the scaling exponents change from first to second-order. On any finite volume, this will lead to intermediate values of the exponents, which should be determined unambiguously together with an associated error. For this purpose we now construct a quantitative measure of the collapse of our data.  Considering how we judge a collapse plot by eye, the measure of the quality has to be related to the distance between different points at the same value of the universal scaling variable. Inspired by the method used by Barkema and Newman [27] for the thermal random-field Ising model (see also a similar analysis in Appendix A of Ref. 28), we associate a quantitative quality to the collapse of B 4 pL Im q by estimating the average variance of the data as Here, ∆x " x max´xmin is the considered range in the scaling variable, V " N s , N V is the number of considered lattice volumes, whileβ c andν are fixed values for the critical temperature and for the critical exponent ν, respectively. A normalisation factor N´2 V was neglected in front of the expression, since it is irrelevant for the estimate of the critical exponent. It is now possible to obtain an estimate for β c and ν by minimizing Q as function of these two variables. Nevertheless, this is a non-trivial task and there are some problems to be addressed. The integration in Eq. (17) must be done numerically, since the exact functional form of B 4 pxq is unknown. This, in principle, would not be a problem, if only we had the kurtosis at the same x on different volumes. However, in lattice simulations the kurtosis B 4 is measured at fixed values of β, and the mapping (13) between β and x depends on the unknown parameters β c and ν. Therefore, it is not possible to have simulated data uniformly spaced in x for any β c and ν. On the other hand, the measured data can be interpolated in β using the multiple histogram method. Hence, it is possible to reweight the kurtosis in β in such a way that its values at the same x are available for all the volumes. After this step the calculation of Q is trivial. In practice, this implies an interpolation for each pair pβ c ,νq at which Q has to be evaluated, which is too costly if a precise determination of the final value of the critical exponent is desired. A cheaper alternative is to use the reweighting technique to obtain the kurtosis as an approximately continuous function of β, i.e., to add a large number of points between two simulated temperatures. The numerical integration to obtain Q can then be performed with negligible additional error. However, due to the particular form of the map xpβq, sometimes, especially for small values of ν, the number of interpolated points needed to have a sufficiently precise numerical integration can become very large and, therefore, the reweighting very costly. A smarter approach is then required.
As can be seen in Figure 3, the kurtosis of the imaginary part of the Polyakov loop is a quite regular function of β, in the sense that no sudden variations are present. This means that a numerical interpolation of the data which does not take into account the physicsas the multiple histogram method does -will probably still find the correct value of the kurtosis. Clearly, this is true under the assumption that the resolution of the data to be interpolated is high enough. For example, the simulated data are usually too distant in temperature to be correctly interpolated without the reweighting technique. But after having applied the multiple histogram method to the data, a second interpolation can be done very cheaply. In practice, we used the software Mathematica to obtain an interpolated function out of a set of points and perform numeric operations on it. The advantage of having a kurtosis as a function makes the calculation of Q straightforward. Furthermore, it is then possible to automatically minimize Qpβ c , νq as function of two variables.
Next, we need to estimate the statistical error on β c and on ν, which has to contain the error on the reweighted points and the statistical error of our simulations. An error on reweighted points is often obtained using the bootstrap method. This means that, in the reweighting procedure, N boot sets of reweighted kurtosis values are calculated, and the bootstrap errors are extracted from them. Now, instead of using these sets to compute errors on the kurtosis, they can be used to minimize Q, obtaining N boot different estimates of β c and of ν, which will give the desired final error. Since, typically, the number of bootstrap samples is of the order of some hundreds, it is clear that the minimization of Q should not take too much time 1 .
Finally, let us discuss how x min and x max should be chosen. Clearly, no extrapolation outside the simulated interval in β should be done. Thus, the largest ∆x is the interval in x around 0 where data from all volumes are available. Since x c " 0, we have x min ă 0 and x max ą 0 and, in order to have a symmetric Taylor expansion interval [26], we chose Using too large an interval of integration is, in general, not correct, since it assumes data collapse possibly outside the critical region. On the other hand, the width of the scaling region is not known a priori. A clever solution to this problem, successfully applied in [27], consists of repeated analyses for successively decreasing ∆x followed by an extrapolation of the resulting parameters to ∆x Ñ 0. An example is reported in Figure 5.

V. NUMERICAL RESULTS AND DISCUSSION
As mentioned in Section III, our strategy to locate the tricritcal points is to measure the critical exponent ν for different quark masses and see where it changes from its first-order value 1{3 for small and large masses to the 3D Ising value 0 6301p4q for intermediate masses. The changes approach a step function in the thermodynamic limit but remain smooth as far as finite lattice volumes are used to extract ν. The critical exponent ν is preferable over B 4 pβ c , 8q, since it is known to suffer less from finite volume corrections [26,29,30]. The main result of our investigation is reported in Figure 6 (more detailed information about the displayed data is available in Tables IV and V in the Appendix). For each value of the quark mass, the critical exponent ν is extracted, both with the fit analysis used already in Ref. 26 and with the new collapse strategy introduced in Section IV. The agreement between the two methods to extract the critical exponent ν is evident in Figure 6. The quantitative collapse analysis has systematically smaller errors on ν, though. This, together with the fact that no arbitrary decision in the analysis may affect the outcome, should make this method preferable.
In the large mass region, the signal is quite smooth and ν changes monotonically from the second-order to the first-order value. However, approaching and entering the first-order region, the minimal aspect ratio N s {N τ needed to extract ν increases significantly from 2 to 6 compared to studies at µ " 0. This has been remarked already in previous studies [13,26]. Presumably this is due to the fact that in the Roberge-Weiss plane we are dealing with the more complex three-state coexistence and its coalescence in a tricritical point.
The same behavior is expected also in the small-mass first-order region, where simulations with an aspect ratio larger than 4 are too costly. Therefore, in Figure 6, the mass regionm u,d ď 0.007 has been marked with a gray background to stress that larger volumes are required to polish the result. However, we have reasons to believe that the tricritical point is already located. It is always possible to compare the critical exponent ν extracted using only part of the available volumes, while leaving the smallest or largest out of the analysis. In this way finite size effects are made visible by checking whether ν drifts towards first-order or second-order values upon inclusion of larger volumes. This is shown form u,d " 0.006 in Figure 6 where a clear decrease in ν is visible adding N s " 30 and removing N s " 12 in the analysis. Another aspect that made us confident to be entering the first-order region form u,d À 0.007 is the typical "Binder bump" behavior discussed in detail in Ref. 26. Atm u,d " 0.007 the kurtosis of the order parameter starts to overshoot the value 3 for β À β c , which is due to the coexistence of three states and thus clearly signals entering the firstorder region.
A few points in Figure 6 were accepted to be obtained from two spatial volumes only. For the two smallest quark masses which were simulated, it was clear from the crossing point of the kurtosis on N s P {12, 18, 24} that N s " 12 was too far away from the thermodynamic limit. On the other hand, to add a larger spatial extent would have been very costly without the guarantee to be sufficient for a conclusive statement. Aboutm u,d " 0.75, instead, we considered the outcome of the analysis with N s P {36, 42} satisfactory, since the kurtosis of the order parameter reaches values larger than 3.5 for β À β c and the bump shrinks and get larger increasing N s , behavior typical of the first-order region.
After these considerations, our estimates of the tricritical bare quark masses arê m tric light " 0.007`0 .002 where the conservative choice of having an asymmetric error in the chiral region is to stress that further investigation would be needed in the chiral limit to polish the measurement.
In order to asses how much the results are affected by cut-off effects, we measured both the lattice spacing a and the pion mass m π for all simulated bare quark masses, by running T " µ i " 0 simulations at the β c found in the Roberge-Weiss plane. The outcome is reported in Table II. Having fixed the scale, it is possible to express Eq. (19) in terms of pion masses in physical units,  Tables IV and V for the the data). Atm u,d " 0.006 the outcome of the data analysis without the largest spatial value (shaded points) has been included to guide the discussion in the text. For each mass, a badge containing the aspect ratios Ns{Nτ used in the analysis has been drawn. The horizontal colored lines are the critical values of ν for some universality classes. The mass axis has been broken and two different scales have been used in order to improve readability. Points in the grey-background region have to be taken with a pinch of salt, since it is reasonable to believe that finite size effects are for them dominant.
Our results in physical units are given in Figure 7, where the critical exponent ν obtained with our new analysis strategy is plotted as a function of m π . It is important to stress that in the large-mass region the lattices used are still too coarse to correctly resolve the pion and we have am π ą 1, implying sizeable cut-off effects on this value. We now compare with the previous results obtained on N τ " 4 lattices [13]. However, there, only the tricritical bare quark masses and a rough estimate of m tric π, light were reported. We therefore improve the determination of the latter, by performing additional scale setting simulations, whose outcome can be found in Table II. In particular, we measured m π and a for three values of the quark bare mass, corresponding to the light tricritical point quoted in Ref. 13 (the central value and at one standard deviation apart from it). The value of β has been chosen using a polynomial interpolation of the β c obtained by the authors at the simulated masses. Taking as error on the tricritical pion mass the difference between its value and m π resulting from the neighboring bare masses, we obtain m tric, Nτ "4 π, light " 473`2 9 28 MeV .
We thus conclude that, in the light region, a shift of around 44% is found when moving from a N τ " 4 to a finer N τ " 6 lattice. In the heavy mass region, only a rough comparison is possible, since no pion mass is reported in Ref. 13 and in any case on N τ " 4 the pion is resolved even less. However, it is possible to compare the dimensionless ratiom u,d {T at the tricritical point, which turn out to be compatible.

VI. DISCUSSION AND CONCLUSIONS
Since we are still far from being able to perform a continuum extrapolation, it is instructive to compare with other discretizations. The results obtained in Ref. 26 with Wilson fermions on N τ " 6 lattices, i.e. with similar  Table II. Results of the scale setting. T " 0 simulations have been performed on 32ˆ16 3 lattices always collecting 800 independent configurations. w0{a has been determined and converted to physical scales using the publicly available code described in [19]. [26] [13] and this study This study [4] [11] Ref. lattice spacing, appear to have considerably larger cut-off effects. For example, comparing am tric π, heavy " 2.2302p2q from Ref. 26 with our am tric π, heavy " 1.7260p3q, the pionresolution problem is milder in the present study. It is also interesting to compare the position of the tricritical points in physical units, The large differences between discretizations again imply being far from the continuum limit, where results from all discretizations have to merge. The observed trend is consistent with the findings of simulations with improved staggered actions, where the tricritical points can only be bounded to be at much smaller masses, as indicated in Figure 8, as well as with the analogous findings at zero chemical potential (see discussion in the introduction). In particular the comparison across discretizations implies enormous cut-off effects in the critical masses, which could end up being over "100% of an eventual continuum limit. We remark that cut-off effects in the critical temperatures are much milder. At present, there is no theoretical explanation as to why the discretization effects on critical quark masses in the Columbia plot are so strong.
In conclusion, we have determined the shift of the tricritical points in the Roberge-Weiss plane of unimproved staggered fermions by changing from N τ " 4 to N τ " 6 lattices. The aspect ratios and statistics required to extract the correct order of the phase transition are found to be larger in the Roberge-Weiss plane than at µ " 0. We find the cut-off effect on the tricritical masses to be smaller but qualitatively the same as that observed with Wilson fermions, and consistent with results for both discretizations at zero chemical potential. This implies in particular, that the entire chiral critical surface depicted in Figure 1 is shifted significantly towards smaller (and possibly zero) light quark masses, as the lattice spacing decreases, which is also consistent with results from improved staggered actions. Unfortunately, our study also implies that much finer lattices at inevitably smaller quark masses are necessary, before one can hope the results of the light tricritical mass to stabilize in a continuum limit.

ACKNOWLEDGMENTS
We thank Francesca Cuteri for useful discussions and input for Figure 8. The authors acknowledge support by the Deutsche Forschungsgemeinschaft (DFG) through the grant CRC-TR 211 "Strong-interaction matter under extreme conditions" and by the Helmholtz International Center for FAIR within the LOEWE program of the State of Hesse. We also thank the computing staff of L-CSC for their support.
It is known that different powers of the same observable have different integrated autocorrelation times τ int , which can be estimated using the Wolff algorithm [31]. This is important to be taken into account when it comes to measure standardized cumulants, like the kurtosis of a given observable. Binning, i.e. substituting a block of data with its average, allows to obtain uncorrelated data from the correlated ones. This is true if the size of a block is at least twice τ int , though. It is then possible to understand how many independent measurement of the quantity of interest are available in a Monte Carlo simulation, just by dividing the number of the trajectories produced by 2τ int . Clearly, the larger this number is the more accurate the result will be. However, simulations in full QCD are costly and a compromise is needed. We always had at least 100 independent events for B 4 pL Im q in the merged chain obtained by putting together the four independent Markov chain that we produced for each β value. A detailed overview of the collected statistics is presented in Table III.
Tables IV and V contain, instead, the detailed outcome of our analysis, whose data were plotted in Figure 6.  Table III. Overview of the statistics accumulated in all the simulations (red entries are preliminary). Since the resolution in β is not the same at differentm u,d , the number of simulated β has been reported per each range. The accumulated statistics per β varies because of the criterion adopted to stop to increase the statistics on the 4 chains. Therefore we reported here the total number of trajectories produced per given Ns. For each Ns, the number of simulated β, the average integrated autocorrelation time and the smallest number of independent events per chain of B4pLImq can be found in the brackets next to the total statistics. Observe thatτint and n events min are not connected. The former is an average among all the different chains run at one fixed spatial lattice extent, while the latter is the effective length of the shorter chain for that given Ns. The number of independent events is obtained as ratio between the number of produced trajectories and the bin size, which is roughly 2 τint. and ν extr. are the outcome of a linear extrapolation for ∆x Ñ 0. Note that the reweighting resolution in β used to add new points between simulated ones varied between 0.004 and 0.0002 and it has been chosen in order to have around 20 values of the kurtosis to be later interpolated.