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





philipsen@itp.uni-frankfurt.de

    Institut für Theoretische Physik - Goethe-Universität, Germany 

    Max-von-Laue-Str. 1, 60438 Frankfurt am Main


    John von Neumann Institute for Computing (NIC)
    GSI, Planckstr. 1, 64291 Darmstadt, Germany



sciarra@itp.uni-frankfurt.de

    Institut für Theoretische Physik - Goethe-Universität, Germany 

    Max-von-Laue-Str. 1, 60438 Frankfurt am Main





    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, μ=π T/3, which separate mass regions with chiral and deconfinement phase transitions from the
    intermediate region, for QCD with =2 unimproved staggered quarks on =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 μ=0.
    Comparing with previous results from =4 
    we find a ∼50% reduction in the light tricritical pion mass. The heavy tricritical pion mass stays roughly
    the same, but is too heavy to be resolved on =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 ∼100%, implying a drastic shrinking of the chiral first-order region to possibly zero.


12.38.Gc, 05.70.Fh, 11.15.Ha


    Finite size and cut-off effects on the Roberge-Weiss transition in Nf=2 QCD with Staggered fermions

    Alessandro Sciarra
    September 27, 2019
=========================================================================================================




§ 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 <cit.>. 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 <cit.>, 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 <ref>. 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 <cit.>. Note that this qualitative picture is the same for unimproved staggered <cit.>
and unimproved <cit.> as well as improved <cit.> 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 =4 staggered fermions without rooting <cit.>.
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 <cit.> nor imaginary chemical potential <cit.> 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 <ref>, we investigate how the (red) tricritical points on the =2 line in the 
Roberge-Weiss-plane (bottom plane at (μ/T)^2=-(π/3)^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 chiral 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 by-product 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 <ref>. We then proceed to describe our numerical methodology in 
<ref> and our novel analysis method in <ref>. Our numerical results are given in
<ref> before we conclude in <ref>.




§ 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, 
Z(μ)=Z(-μ). For purely imaginary chemical potential, μ=, ∈ℝ,
it is furthermore periodic <cit.>,

    Z(/T)=Z(/T+2π k/N_c),    k=0,… N_c-1 ,

and we use N_c=3 colors for the QCD gauge group. 
These symmetries imply the phase structure shown in <ref>, with three
different ℤ(3) center sectors, which are periodically repeated for higher . Physical observables,
and in particular the thermodynamic functions, are invariant under a change of sectors, which are 
characterized by different phases of the Polyakov loop

    L(𝐱)=1/3Tr∏_τ=1^N_τ-1U_4(τ,𝐱)≡ |L(𝐱)|e^-iφ ,

with ⟨φ⟩ =2kπ/3, k∈{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 <ref>.
The purpose 
of the present work is to locate these tricritical masses on N_τ=6 lattices with =2 and compare their values with previous determinations
on a coarser N_τ=4 lattice <cit.>, as well as with those of other discretization schemes.






§ NUMERICAL SETUP


We consider the QCD partition function of =2 mass-degenerate quarks with a purely imaginary chemical potential. After integration over the fermionic fields it can be written as

    (T,)=∫𝒟U  ( D[U,])^1/2 e^-[U] ,

where  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

    [U]= β∑_P{ 1 - [_C P]} .

The fermion matrix reads

    D_i,j   =  δ_i,j + 
       + 1/2∑_ν=1^4 η_i,ν(Ũ_i,ν δ_i,j-ν̂ - Ũ^†_i-ν̂,ν δ_i,j+ν̂) ,

where =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 =a in the temporal direction,

    Ũ_i,ν={
            U_i,ν       ν∈{1,2,3}
    
            e^iU_i,ν       ν=4
        .  .

The temperature is specified by the inverse euclidean time extent of the lattice,

    T=1/a(β)  .






In order to locate a phase transition and study its nature, we calculate standardized cumulants

    B_n(X,β,,)
             ≡⟨(X-⟨ X⟩)^n⟩/⟨(X-⟨ X⟩)^2⟩^n/2 ,

constructed from an (exact or approximate) order parameter X. 
In particular, a non-trivial zero of the skewness of the X-distribution,

    ()=0 ,

determines at which value of β= a thermal transition takes place, while the value of the kurtosis  
in the thermodynamic limit, evaluated at the critical coupling, will determine the order of the phase transition 
(refer to <ref> for common kurtosis values).
Note that the so-called Binder cumulant <cit.>,

    𝒰_4(X)≡ 1-⟨(X-⟨ X⟩)^4⟩/3⟨(X-⟨ X⟩)^2⟩^2 
                         =1-1/3 B_4 ,

is trivially related to the kurtosis (of the same observable) and contains the same information.

We fix /T=π since in this case the imaginary part of the  is an exact order parameter,

    X=≡1/^3∑_𝐱Im(L(𝐱)) .

Referring to <ref>, it distinguishes between the low T disordered phase and the high T ordered phase with two-state coexistence,

    {⟨φ⟩     =0  ⇒ ⟨⟩     =0       low T
    ⟨φ⟩≠ 0  ⇒ ⟨⟩≠ 0       high T
    .  ,

with the advantage of knowing its mean value exactly.

On a =6 lattice, we thus set =π/6. Since this is the boundary between two Roberge-Weiss sectors for all temperatures, 
()=0 for any value of β and we cannot use <ref> to locate the  end-point.
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 non-analytic step function in the V→∞ 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 β=, provided that the spatial lattice extent is large enough.
This crossing provides the location of the  end-point.
In the neighborhood of the critical point , the kurtosis shows a well-defined finite size scaling behavior
as a function of the scaling variable

    x≡(β-)^1/ν .

Its Taylor expansion around the critical point x=0 is 

    (β, x) = (,∞) + a_1 x + a_2 x^2 + (x^3)  .

Sufficiently close to the thermodynamic limit, the coefficient (,∞) and the critical exponent ν take their universal values depending on the type of transition.

In order to locate the two tricritical points in the  plane, we performed simulations at different values of  and  different values of 
β around the critical temperature.
Evaluating the kurtosis in the critical region and fitting it to <ref>, considering the linear term only, gives 
(,∞), a_1,  and ν for every value of .
The change of ν as a function of  then permits to locate the light and heavy  values.

Although our main quantitative analysis is based on the kurtosis of the order parameter, we also calculated the susceptibility of |L|,

    χ(|L|)≡^3 ⟨(|L|-⟨ |L|⟩)^2⟩ .

which is expected to scale around  according to

    χ=^γ/ν f(t ^1/ν)    .

Here t≡(T-)/ 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 χ/^γ/ν evaluated on different lattice sizes against t^1/ν 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 of  in the intervals [0.004, 0.011] and [0.15,0.85].
For each value of , three to five spacial lattice sizes have been used, keeping =6 and =π/6 fixed.
This corresponds to aspect ratios /∈[2,7].
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 <ref>).
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) <cit.> to increase resolution (see the discussion at the end of this and in the next section).

Configurations were generated by a standard RHMC algorithm <cit.>, 
producing four different Monte Carlo chains per β with unit-length trajectories.
Where advantageous (<0.007), the multiple-pseudofermions technique <cit.> 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  and the chiral condensate for small masses) have been computed for every trajectory.
We increased statistics until the standard deviation of the kurtosis (X) decreased below ∼0.2, and (X) 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 <ref>.
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 , especially when entering the first-order regions (cf. <ref> in <ref>).
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 3216 lattices for each value of .
Using the publicly available code described in Borsanyi:2012zs, the scale was set by the Wilson flow parameter w_0.
Pion masses, instead, were measured with standard spectroscopy techniques <cit.>.

All our numerical simulations (except those for scale setting purposes) have been run using the publicly available <cit.>  based code  <cit.>, which is optimized for GPUs.
The  <cit.> 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  <cit.>.

Our quite intricate fitting procedure used to extract the critical exponent ν is completely analogous to the one previously 
described in Appendix B of Cuteri:2015qkq.
For each value of the quark mass, nearly all possible fits of the data to the linear part of <ref> are performed, and a filtering procedure is applied afterwards in order to pick the best fits.
This is needed because the range in which () can be considered linear for each  value is not known a priori.
However, here we differ in a technical detail from the previous study in Cuteri:2015qkq.
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  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 , making the information coming from the smallest volume systematically less important, see <ref> for detailed overview.
This aspect is even more explicit looking at how many points per  have been included in the fit.






§ 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 () as function of the scaling variable x defined in <ref>.
An example is shown in <ref>.

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 <ref>) 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 <ref>).
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 <cit.> for the thermal random-field Ising model (see also a similar analysis in Appendix A of Houdayer_2004), we associate a quantitative quality to the collapse of () by estimating the average variance of
the data as

    Q(,)≡1/Δ x∫_x_min^x_max{   N_V∑_i=1^N_V[(x(β̅_c,,V_i))]^2 
       - [∑_i=1^N_V(x(β̅_c,,V_i))]^2 }dx .

Here, Δ x≡ x_max-x_min is the considered range in the scaling variable, V≡, N_V is the number of considered lattice volumes, while  and  are fixed values for the critical temperature and for the critical exponent ν, respectively.
A normalisation factor N_V^-2 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  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 <ref> must be done numerically, since the exact functional form of (x) 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  is measured at fixed values of β, and the mapping (<ref>) between β and x depends on the unknown parameters  and ν.
Therefore, it is not possible to have simulated data uniformly spaced in x for any  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 (,) 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 x(β), 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 <ref>, the kurtosis of the imaginary part of the  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 physics – as 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  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 Q(,ν) as function of two variables.

Next, we need to estimate the statistical error on  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  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[In , for example, it is possible to use the  function, but a user implemented minimization based on a scan in  and in ν is more efficient, though less precise.].

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 <cit.>, we chose

    | x_min|=| x_max| .

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 <cit.>, consists of repeated analyses 
for successively decreasing Δ x followed by an extrapolation of the resulting parameters to Δ x→0.
An example is reported in <ref>.






§ NUMERICAL RESULTS AND DISCUSSION


As mentioned in <ref>, 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.6301(4) 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 (, ∞), since it is known to suffer less from finite volume corrections <cit.>.
The main result of our investigation is reported in <ref> (more detailed information about the displayed data is available in <ref> in <ref>).
For each value of the quark mass, the critical exponent ν is extracted, both with the fit analysis used already in Cuteri:2015qkq and with the new collapse strategy introduced in <ref>.
The agreement between the two methods to extract the critical exponent ν is evident in <ref>.
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 / needed to extract ν increases significantly from 2 to 6 compared to studies at μ=0.
This has been remarked already in previous studies <cit.>.
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 <ref>, the mass region ≤ 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 for =0.006 in <ref> where a clear decrease in ν is visible adding =30 and removing =12 
in the analysis.
Another aspect that made us confident to be entering the first-order region for ≲0.007 is the typical “Binder bump” behavior discussed in detail in Cuteri:2015qkq.
At =0.007 the kurtosis of the order parameter starts to overshoot the value 3 for β≲, which is due
to the coexistence of three states and thus clearly signals entering the first-order region.

A few points in <ref> 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 ∈{12,18,24} that =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.
About =0.75, instead, we considered the outcome of the analysis with ∈{36,42} satisfactory, since the kurtosis of the order parameter reaches values larger than 3.5 for β≲ and the bump shrinks and get larger increasing , behavior typical of the first-order region.

After these considerations, our estimates of the tricritical bare quark masses are

    [light]   =0.007^+0.002_-0.003
    [heavy]   =0.55(10)  ,

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  for all simulated bare quark masses, by running T==0 simulations at the  found in the  plane.
The outcome is reported in <ref>.
Having fixed the scale, it is possible to express <ref> in terms of pion masses in physical units,

    [π, light]   =328^+44_-81 
    [π, heavy]   =2813^+235_-261   .

Our results in physical units are given in <ref>, where the critical exponent ν obtained with our new analysis strategy 
is plotted as a function of .
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 a>1, implying sizeable cut-off effects on this value. 

We now compare with the previous results obtained on =4 lattices <cit.>.
However, there, only the tricritical bare quark masses and a rough estimate of [π, light] were reported.
We therefore improve the determination of the latter, by performing additional scale setting simulations, whose outcome can be found in <ref>.
In particular, we measured  and a for three values of the quark bare mass, corresponding to the light tricritical point quoted in Bonati:2010gi (the central value and at one standard deviation apart from it).
The value of β has been chosen using a polynomial interpolation of the  obtained by the authors at the simulated masses.
Taking as error on the tricritical pion mass the difference between its value and  resulting from the neighboring bare masses, 
we obtain

    m_π, light^tric, =4=473^+29_-28  .

We thus conclude that, in the light region, a shift of around 44% is found when moving from a =4 to a finer 
=6 lattice.
In the heavy mass region, only a rough comparison is possible, since no pion mass is reported in Bonati:2010gi and in any case 
on =4 the pion is resolved even less.
However, it is possible to compare the dimensionless ratio /T at the tricritical point,

    m_π, heavy^tric, =4/T   =2.9(3)
    
            m_π, heavy^tric, =6/T   =3.3(6)
    
           ,

which turn out to be compatible.






§ 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 Cuteri:2015qkq with Wilson fermions on =6 lattices, i.e. with similar lattice spacing,
appear to have considerably larger cut-off effects.
For example, comparing  from  Cuteri:2015qkq with our , the pion-resolution problem is milder in the present study.
It is also interesting to compare the position of the tricritical points in physical units,

    m_π, light^tric, Wilson   =669^+95_-81 
    
            m_π, light^tric, Staggered   =328^+44_-81

and

    m_π, heavy^tric,Wilson   =3659^+589_-619 
    
            m_π, heavy^tric, Staggered   =2813^+235_-261   .

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 <ref>, 
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 =4 to =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 <ref> 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.



    We thank Francesca Cuteri for useful discussions and input for <ref>.
    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  for their support. 




apsrev4-1


*



§ SIMULATION DETAILS


It is known that different powers of the same observable have different integrated autocorrelation times , which can be estimated using the Wolff algorithm <cit.>.
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 , 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.
Clearly, the larger this number is the more accurate the result will be.
However, simulations in full  are costly and a compromise is needed.
We always had at least 100 independent events for () 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 <ref>.

<Ref> contain, instead, the detailed outcome of our analysis, whose data were plotted in <ref>.


    
    
    
    
    
    
    
    
    cyan!20
    magenta!20
    yellow!20
    orange!20
    
    

    
    
    
    
    
    
    L>c

    


/Rotate 90


    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
     
    

    L>r
    G>[gray]0.2[.5][]c
    I>[gray]0.2[][.5]c

    
        
        
        table-number-alignment = center
        
                       points     δβ·10^3          ν    
            (,∞)     a_1     χ^2__d.o.f.     Q(%)     Ω_min     Ξ_min 

            0.004             18  2̊4̊             24  8̊                                       0.3  0̊.̊2̊         5.43241(10)     0.52  (8)      1.794(14)     -0.44(40)      0.864     67.18     52.21     72.17 

            0.005             18  2̊4̊             11  7̊                                       0.5  0̊.̊5̊         5.43649(11)     0.49  (4)      1.699(14)     -0.27(14)      1.004     44.56     96.06     99.33 

            0.006             18  2̊4̊  3̧0̧     8  5̊  1̧3̧                     1  1̊  0̧.̧3̧         5.43918(7)      0.554(26)      1.823(11)     -0.55(14)      1.005     45.27     84.11     78.89 

            0.007             18  2̊4̊  3̧0̧     5  4̊  3̧                 1  0̊.̊5̊  0̧.̧5̧             5.44247(10)     0.48  (6)      1.780(14)     -0.21(17)      0.768     63.12     81.44     62.67 

            0.008             1̱2̱  18  2̊4̊     7̱  6  6̊                 2̱  1  0̊.̊5̊         5.44513(11)     0.472 (26)     1.878(10)     -0.23(7)       0.995     45.70     94.98     85.20 

            0.009             1̱2̱  18  2̊4̊     6̱  8  9̊                 2̱  1  0̊.̊5̊         5.44828(10)     0.540 (25)     1.871(10)     -0.44(11)      1.001     45.51     80.52     85.60 

            0.010             1̱2̱  18  2̊4̊     6̱  5  8̊                 3̱  2  1̊     5.45154(14)     0.623 (26)     1.894(12)     -0.88(16)      1.005     44.65     80.74     72.80 

            0.011             1̱2̱  18  2̊4̊     6̱  8  1̊1̊                     4̱  2  1̊     5.45350(16)     0.650 (27)     1.999(18)     -1.07(20)      0.994     46.71     81.53     64.29 

            [0.2pt]
            0.150             1̱2̱  18  2̊4̊     6̱  16  1̊0̊                         4̱  1  1̊     5.64787(23)     0.66  (4)      1.854(13)     -0.65(18)      1.007     45.43     83.79     68.40 

            0.250             1̱2̱  18  2̊4̊     8̱  23  1̊1̊                         4̱  1  1̊     5.71176(22)     0.652 (26)     1.914(10)     -0.57(10)      0.990     48.67     81.19     87.43 

            0.350             1̱2̱  18  2̊4̊     5̱  8  1̊0̊                     5̱  2  1̊     5.75551(16)     0.596 (18)     1.900(7)      -0.39(6)       0.997     46.08     81.96     78.00 

            0.400             18  2̊4̊  3̧0̧     5  7̊  1̧0̧                     2  1̊  0̧.̧5̧         5.77364(14)     0.59  (4)      1.777(11)     -0.36(12)      0.992     46.57     85.59     72.89 

            0.450     18  2̊4̊  3̧0̧  36     4  3̊  5̧  5     2  2̊  0̧.̧5̧  0.5      5.78783(7)      0.54  (3)      1.777(8)      -0.23(8)       0.997     45.11     83.15     83.00 

            0.500             2̊4̊  3̧0̧  36     4̊  7̧  7                 2̊  0̧.̧5̧  0.5             5.80037(18)     0.60  (6)      1.735(23)     -0.44(23)      0.994     45.61     83.99     58.00 

            0.550             2̊4̊  3̧0̧  36     4̊  1̧2̧  8                     2̊  0̧.̧5̧  0.5             5.81091(13)     0.51  (4)      1.783(20)     -0.17(10)      0.997     46.22     82.53     90.86 

            0.600             2̊4̊  3̧0̧  36     6̊  8̧  6                 1̊  0̧.̧5̧  0.5             5.82012(6)      0.457 (24)     1.758(9)      -0.09(4)       0.991     46.26     90.35     70.40 

            0.650             2̊4̊  3̧0̧  36     7̊  8̧  8                 1̊  0̧.̧5̧  0.5             5.82785(9)      0.481 (28)     1.813(12)     -0.13(5)       0.991     46.80     81.26     94.29 

            0.750                     36  ø42                    15  ø16                        0.3  ø0.2             5.84112(10)     0.39  (4)      1.86(3)       -0.031(29)     0.992     47.54     93.13     58.67 

            0.850             3̧0̧  36  ø42     5̧  8  ø8                 1̧  0.3  ø0.2             5.85124(5)      0.358 (25)     1.807(17)     -0.014(9)      0.996     45.88     83.39     51.43 

            
        Overview of the selected fits to extract the final value of ν (results on grey background are preliminary).
                 The fits have been performed according to <ref>, considering the linear term only.
                 The  column contains the spatial lattice extents that have been included in the fits.
                 Ω_min and Ξ_min are respectively the minimum overlap percentage and the minimum symmetry percentage as defined in Eqs. (B3) and (B4) of Cuteri:2015qkq.
                 In the third and the fourth column, the number of fitted points per  and the reweighting resolution in β have been reported, respectively.
                 The colors should help as guideline to distinguish the information among the different volumes.
                 Observe how the smallest volume has systematically a smaller weight in the fit (apart from the two smallest masses, which should anyway investigated more).
        
    

    
/Rotate 0


