Roberge-Weiss endpoint and chiral symmetry restoration in $N_f = 2+1$ QCD

We investigate the fate of the Roberge-Weiss endpoint transition and its connection with the restoration of chiral symmetry as the chiral limit of $N_f = 2+1$ QCD is approached. We adopt a stout staggered discretization on lattices with $N_t = 4$ sites in the temporal direction; the chiral limit is approached maintaining a constant physical value of the strange-to-light mass ratio and exploring three different light quark masses, corresponding to pseudo-Goldstone pion masses $m_\pi \simeq 100, 70$ and 50 MeV around the transition. A finite size scaling analysis provides evidence that the transition remains second order, in the 3D Ising universality class, in all the explored mass range. The residual chiral symmetry of the staggered action also allows us to investigate the relation between the Roberge-Weiss endpoint transition and the chiral restoration transition as the chiral limit is approached: our results, including the critical scaling of the chiral condensate, are consistent with a coincidence of the two transitions in the chiral limit; however we are not able to discern the symmetry controlling the critical behavior, because the critical indexes relevant to the scaling of the chiral condensate are very close to each other for the two possible universality classes (3D Ising or O(2)).

We investigate the fate of the Roberge-Weiss endpoint transition and its connection with the restoration of chiral symmetry as the chiral limit of N f = 2 + 1 QCD is approached. We adopt a stout staggered discretization on lattices with Nt = 4 sites in the temporal direction; the chiral limit is approached maintaining a constant physical value of the strange-to-light mass ratio and exploring three different light quark masses, corresponding to pseudo-Goldstone pion masses mπ ≃ 100, 70 and 50 MeV around the transition. A finite size scaling analysis provides evidence that the transition remains second order, in the 3D Ising universality class, in all the explored mass range. The residual chiral symmetry of the staggered action also allows us to investigate the relation between the Roberge-Weiss endpoint transition and the chiral restoration transition as the chiral limit is approached: our results, including the critical scaling of the chiral condensate, are consistent with a coincidence of the two transitions in the chiral limit; however we are not able to discern the symmetry controlling the critical behavior, because the critical indexes relevant to the scaling of the chiral condensate are very close to each other for the two possible universality classes (3D Ising or O(2)).

I. INTRODUCTION
Numerical investigation of QCD or QCD-like theories in the presence of imaginary chemical potentials coupled to quark number operators has been the subject of various lattice studies . The main source of interest is the possibility of obtaining information about QCD at finite baryon density via analytic continuation, thus partially avoiding the sign problem. Moreover, numerical results at imaginary µ are also a relevant test bed for effective models trying to reproduce the properties of QCD at finite density [27][28][29]. Furthermore, imaginary chemical potentials are an interesting extension of the QCD phase diagram per se, as, for particular choices of the chemical potentials, one recovers exact symmetries even in the presence of finite quark masses, leading to the presence of interesting phase transitions and critical points which, in principle, could be relevant also for the physical region of the phase diagram.
A well known example is QCD with an imaginary baryon chemical potential µ B , which, for particular values, known as Roberge-Weiss (RW) points [30] (µ B ≡ iµ B,I = ikπT where k is an odd integer), has an exact Z 2 symmetry, which is a remnant of the original Z 3 symmetry present in the pure gauge case; this simmetry gets spontaneusouly broken at a critical temperature T RW which fixes the endpoint (RW endpoint) of first order transition lines which are present (at fixed µ B ) in the high-T region of the phase diagram, as sketched in Fig. 1. More exotic combinations have been also considered, like those in which an exact Z Nc center symmetry is recovered (N c being the number of colors) by locking it to flavor symmetry in the presence of N f = N c degenerate flavors [31][32][33][34][35][36][37].
The RW transition lines and their endpoints have been thoroughly investigated by lattice simulations [3, 4, 12-15, 20, 26, 38-49] and effective models [50][51][52][53][54][55][56][57][58][59][60][61][62]. Early studies, performed on lattices with N t = 4 sites in the temporal direction and using unimproved staggered fermions, have shown interesting features for the RW endpoint transition for both N f = 2 and N f = 3 degenerate flavors: the transition is first order for small quark masses, likely down to the chiral limit, second order for intermediate masses, and first order again for large quark masses; the three regions are separated by tricritical points [38][39][40]; for N f = 2 the tricritical point delimiting the first order chiral region takes place for m π ≃ 400 MeV [40]. These results, which suggest a strict relation of the RW endpoint transition to the chiral properties of the theory, have been confirmed by simulations employing standard Wilson fermions, even if with indications of a strong cut-off dependence for the location of the tricritical points: indeed, the chiral tricritical light pion mass has been located at m π ≃ 910 MeV for N t = 4 and at m π ≃ 670 MeV for N t = 6 [49]. A systematic study adopting stout improved staggered fermions has been reported in Ref. [46] for N f = 2 + 1 QCD with physical quark masses, employing lattices with N t = 4, 6, 8 and 10, i.e. going down to lattice spacings of the order of 0.1 fm. This has permitted to obtain a reliable continuum extrapolation for the endpoint transition temperature, T RW ≃ 208(5) MeV, corresponding to T RW /T c ≃ 1.34 (7) where T c is the pseudo-critical chiral crossover temperature at zero baryon chemical potential. A finite size scaling (FSS) analysis has been performed only for N t = 4 and 6 and has provided evidence for a second order transition, meaning that the chiral tricritical pion mass, if any, is lower than the physical pion mass, m π ≃ 135 MeV.
In this study we extend the analysis of the RW endpoint adopting the same improved discretization already used in Ref. [46], exploring lower than physical quark masses, going down to a pseudo-Goldstone pion mass of the order of 50 MeV. Our purpose is twofold. First, in view of the apparent strong reduction of the chiral first order region, we would like to understand if a chiral tri-critical pion mass can still be located. The exploration of the QCD phase diagram at zero chemical potential, usually summarized in the so-called Columbia plot, has provided evidence for a general shrinking of the first order regions as the continuum limit is approached, and presently it is not even clear if a first order survives in the chiral limit for the N f = 4 or N f = 3 case [63], where standard universality arguments would predict it [64]. Therefore, the fact that a first order RW endpoint transition is still found in the chiral and continuum limit of N f = 2 + 1 QCD is not guaranteed.
Let us say right from the beginning that the task itself is highly non-trivial. Indeed, due to well known problems in the lattice discretization of fermion degrees of freedom, a reliable approach to the chiral limit is only possible if the continuum limit is approached first [65], while keeping finite size effects under control. In other words, for a reliable investigation one should guarantee at the same time that: i) one gets close enough to the chiral limit; ii) one stays close enough to the continuum limit so that the chiral properties of dynamical fermions are effective (in the present context of staggered fermions, that means that taste symmetry breaking is negligible and all pions becomes effectively light); iii) the physical volume of the system is still large enough, in particular Lm π ≫ 1 as m π → 0. Satisfying all these criteria is still an unbearable task, even for present computational resources.
The present study is limited to N t = 4 lattices and therefore represents just a small step. We anticipate that we have not been able to detect any signal of a first order transition down to m π ≃ 50 MeV. On one hand, one might consider this result as inconclusive for the reasons exposed above: as we will show, our approach to the chiral limit on N t = 4 actually means that just one pion mass goes to zero, while the others stay finite and quite heavy (larger than 400 MeV), so that it is not clear which kind of "chiral" theory one is really approaching. Nevertheless, on the other hand, it is a striking fact that results change so drastically with respect to earlier results on N t = 4 lattices [38][39][40] (a tricritical pion mass decreasing by at least one order of magnitude or vanishing at all), by just improving the discretization of the theory.
The second purpose of our investigation is an improved understanding of the relation between chiral and center symmetry and the RW transition. In the massless limit, one could expect two different transitions temperatures, T χ and T RW , along the RW chemical potentials (µ B = ikπT with k odd), one corresponding to the restoration of chiral symmetry and the other to the breaking of the remnant center symmetry. Examples where the chiral restoration transition is well decoupled from the center symmetry breaking transition are well known in the literature, like for instance QCD with fermions in the adjoint representation [66]. In this case, results obtained at finite quark mass show that the two transitions are generically close to each other; however what happens in the chiral limit, where both symmetries are exact, is unknown. In this case the task is more feasible. Indeed, even at finite lattice spacing, the staggered discretization provides a remnant of the chiral symmetry which becomes exact as the bare quark mass is extrapolated to zero: it corresponds to a single generator of the original chiral group, it breaks spontaneously at low temperature, leading to a single massless pion, and it gets restored at the chiral transition temperature T χ . Therefore, it makes sense to investigate the relation between T RW and T χ in the chiral limit also for finite values of N t , even if of course the answer itself could be N t -dependent. Whether the two transitions coincides and, in this case, which symmetry controls the critical behavior, is a clear-cut question which can and should be answered.
The paper is organized as follows. In Section II we review the general properties of QCD in the presence of imaginary chemical potentials, illustrate the lattice discretization adopted for our study and give details on our numerical setup and analysis. In Section III we report our numerical results regarding the order of the transition for different values of the bare quark mass, discussing also the corresponding values of the pion masses (pseudo-Goldstone and not) and the quality of our approach to the chiral limit. In Section IV we investigate the relation between T χ and T RW as the chiral limit is approached. Finally in Section V we present our concluding remarks.

II. NUMERICAL SETUP
We consider a rooted stout staggered discretization of N f = 2 + 1 QCD in the presence of imaginary quark chemical potentials µ f,I , its partition function reads: S YM is the tree level Symanzik improved gauge action [67,68] constructed in terms of the original link variables, W n×m i; µν being the trace of a n × m rectangular loop, while the staggered fermion matrix (M f st ) i, j is built up in terms of the two times stout-smeared [69] links U (2) i; ν , with an isotropic smearing parameter ρ = 0.15.
Adopting thermal boundary conditions (periodic/antiperiodic in Euclidean time for boson/fermion fields), the temperature is given by T = 1/(N t a); we have fixed N t = 4 in all simulations, while the lattice spacing a is a function of β and of the bare quark masses. In this study,  the value of the bare gauge coupling β, it is possible to make use of standard reweighting methods [70] in order to optimize the numerical effort; that was not possible in Ref. [46], where also the weight of the fermion determinant changed from one simulation to the other (because of the tuning of the quark masses), making reweighting not feasible in practice.
In Ref. [46], the critical β reported for N t = 4 is β RW (N t = 4) ≃ 3.45, which corresponds to am l ≃ 0.00558 according to the LCP determined in Refs. [71][72][73]. Based on that, we have decided to run simulations for three different values of the quark masses, namely am l = 0.003, am l = 0.0015 and am l = 0.00075: for each value we have located the pseudo-critical coupling β RW (am l , N t ) and performed a series of run at different values of β around β RW which have then been used for reweighting. In each case, simulations have been performed on lattices L 3 s × 4, where different values of the spatial extent L s (in the range 16 → 32) have been considered to perform a FSS analysis. For some selected values of β values for each mass we have performed numerical simulations also on T ∼ 0 lattices, which have been used for renormalization and scale setting purposes. Table I shows a complete list of our finite T simulation parameters; statistics reach up to 50K Rational Hybrid Monte-Carlo unit length trajectories for simulation points around the transition.
The bare quark masses have been chosen in order to reach, for the lowest mass, a pion mass approximately equal to m π = m (phys) π × 0.00075/0.00558 ≃ 50 MeV around the transition point. This estimate is only qualitative, as also the critical bare coupling moves as we change am l . For this reason, simulations at T ≃ 0 have been performed in order to obtain a direct determination of m π at the different simulation points. Zero temperature runs have been used also to determine the lattice spacing, exploiting a technique based on the gradient flow [74] and in particular the so-called w 0 parameter [75]. All scale setting and pion mass determinations are shown in Table II:  trajectories for each simulation point. Pion masses have been obtained from standard Euclidean time correlators of appropriate staggered quark operators (see, e.g., Refs. [76,77]). In this case, in addition to the lowest (pseudo-Goldstone) pion state, we have also determined other pion masses, which are expected to be higher, at finite lattice spacing, because of the taste violations of the staggered discretization. With the purpose of estimating the magnitude of such taste violations, which fix the quality of our actual approach to the chiral limit, we report in Table II also the value of the mass of the first excited pion, m (1) π . In order to better visualize the quality of our approach to the chiral limit, in Fig. 2 we show, for a fixed value of the bare gauge coupling β = 3.39, the values obtained for the pseudo-Goldstone pion and for m (1) π as a function of the square root of the light bare quark mass (in physical units). It is quite striking that, while m π approaches zero as m l → 0 following quite closely the prediction of chiral perturbation theory, m π ∝ √ m l , the first excited pion is instead much less affected by the change of am l . Therefore, in our approach to the chiral limit and for what concerns the critical behavior around the transition, we are effectively considering a theory with no more than one light pion: that is quite different from the physical theory and, eventually, one would like to understand how this fact may bias the results obtained for the order of the phase transition.
In order to implement a purely baryonic chemical potential (i.e. µ Q = µ S = 0) we have set µ u = µ d = µ s ≡ µ q = µ B /3. An imaginary µ q is equivalent to a rotation of fermionic temporal boundary conditions by an angle θ q = Im (µ q )/T , there is therefore a periodicity in θ q , which however is 2π/N c (instead of 2π) because this ro-tation can be exactly canceled by a center transformation on gauge fields. This periodicity is smoothly realized at low T , while at high T the value of θ q selects among the three different minima of the Polyakov loop effective potential, leading to first order phase transitions which occur when θ q crosses the boundary between two adiacent center sectors. These transitions form first order lines (RW lines) located at θ q = (2k + 1)π/N c and k integer: there the average Polyakov loop L jumps from one center sector to the other and serves as an order parameter for such transitions. A sketch of the phase diagram is reported in Fig. 1: each RW line terminates with an endpoint located at a temperature T RW , where an exact Z 2 symmetry breaks spontaneously. Therefore, moving in temperature along these lines, one can meet either a second order critical point in the 3D-Ising universality class, or a first transition; in the latter case the endpoint is actually a triple point.
In the following we shall consider one particular RW line, θ q = π, for which the imaginary part of the Polyakov loop can serve as an order parameter. In order to identify the universality class of the endpoint, a FSS analysis will be performed for the susceptibility of the order parameter As an alternative order parameter, one could take any of the quark number densities (where q = u, d, s) which should vanish for θ q = (2k + 1)π/N c (because of the mentioned periodicity and because they are odd in θ q ) unless the Z 2 symmetry (which is equivalent to charge conjugation) is spontaneously broken. However, our analysis will be based exclusively on the Polyakov loop.
Numerical simulations have been performed on the COKA cluster, using 5 computing nodes, each with 8 NVIDIA K80 dual-GPU boards and two 56 Gb/s FDR InfiniBand network interfaces. Our parallel code (Open-StaPLE) is a single [78] and multi [79] GPU implementation of a standard Rational Hybrid Monte-Carlo algorithm. It is an evolution of a previous CUDA code [80], developed using the OpenACC and OpenMPI frameworks to manage respectively parallelism on the GPUs and among the nodes. The multi-GPU implementation [79] has been essential in order to perform some of the zero temperature runs, which otherwise would have not fitted on a single GPU for memory reasons.
Of course, the most expensive simulations have been those regarding the lowest explored quark mass, am l = 0.00075. On the whole, a rough estimate of the total computational cost of our investigation is 3 × 10 5 equivalent run-hours on a K80 GPU.

III. FINITE SIZE SCALING ANALYSIS AND ORDER OF THE TRANSITION
The susceptibility χ L , defined in Eq. (4), is expected to scale as where t = (T − T RW )/T RW is the reduced temperature and one has t ∝ β − β RW close enough to T RW . This means that χ L /L γ/ν s , measured on different spatial sizes,  should lie on a universal scaling curve when plotted as a function of (β − β RW )L 1/ν s . The critical exponents which are relevant to our analysis are reported in Table III. Apart from first order and 3D-Ising exponents, we also report tricritical indexes: they are expected to describe the critical behavior exactly at the separation point between the first order and the second order region, however, before the thermodynamic limit is really approached, they could describe the critical behavior in a finite neighborhood of the tricritical point [84].
for the three different masses is reported in Figs. 3, 4 and 5, respectively for first order, 3D-Ising and tricritical indexes. It clearly appears that a first order transition is excluded for all masses, while a reasonable scaling is obtained when considering both the 3D-Ising and the tricritical critical behavior.
As a further confirmation of the absence of a first order transition for all explored masses, in Fig. 6 we report, just for the lowest quark mass, am l = 0.00075, the probability distribution of the plaquette and of the unrenormalized quark condensate at the critical point for the different lattice sizes. A vague double peak structure is visible only in the distribution of the chiral condensate and for small L s , however it tends to disappear as the thermodynamic limit is approached. Therefore, our results suggest that a chiral first order region, if any, is limited to a region of pion masses below 50 MeV. There are of course many systematics that should be considered before drawing a definite conclusions. First of all, as we have already discussed, our approach to the chiral limit actually means that just one pion becomes massless, while all other pion masses stay above 400 MeV. Therefore one should repeat this study with significantly larger values of N t (smaller lattice spacings), so that also the other pions become lighter. In principle, additional chiral degrees of freedom could change the scenario and make the first order region larger, even if this is at odds with the common experience of shrinking of first order regions as the continuum limit is approached. Unfortunately, going to significantly larger values of N t is not feasible with our present computational resources, so this is left for future work.
A second remark regards the lattice sizes that we have adopted in our study, in particular the maximum values of aL s m π that we have reached are 2, 3, and 4 respectively for am l = 0.00075, am l = 0.0015 and am l = 0.003. The values are not particularly large, especially for the lowest explored quark mass. However, we have seen no significant deviation from a second order scaling, and no signal for the development of a double peak structure as the volume is increased; on the contrary, some weak double peak signals visible in the chiral condensate distribution for small L s have shown a tendency to disappear when going to larger volumes.

IV. CHIRAL SYMMETRY RESTORATION AND THE ROBERGE-WEISS ENDPOINT TRANSITION
The existence, for the staggered fermion discretization, of an unbroken remnant of the full continuum chiral symmetry group, permits to consider a well posed question, regarding the connection between chiral symmetry restoration and the Roberge-Weiss transition, even on the coarse lattices explored in our investigation. In short, the question is the following: in the chiral limit and for µ B = ikπT , with k odd, the theory enjoys both chiral symmetry and the Z 2 RW symmetry, which are both expected to undergo spontaneous symmetry breaking (or restoration) at two temperatures, T χ and T RW . While results obtained for finite quark masses indicate a generic closeness of the two phenomena, one would like to know if actually T χ = T RW or not. Moreover, if the two temperatures coincide, which of the two symmetries dominates the transition and fixes its universality class? The latter question is important to understand what are the relevant degrees of freedom around the transition in a non-trivial theory like QCD, where chiral and gauge degrees of freedom are strictly entangled 1 .
In order to answer the questions above, we consider the behavior of the (light) chiral condensate, which is the order parameter for chiral symmetry breaking and is defined as follows: where V = L 3 s is the spatial volume and the contribution from each light flavor f is expressed in terms of the following lattice observablē Tr 1 which has been evaluated by means of noisy estimators (in particular up to 16 Z 2 random vectors have been used for each measurement). The light quark condensate is affected by additive and multiplicative renormalizations, which can be taken care of by, respectively, appropriate  subtractions and ratios. In particular, in this study we consider the following prescription [87]: where the leading additive renormalization, which is linear in the quark mass, cancels in the difference with the strange condensate, while the multiplicative renormalization, being independent of T , drops out by normalizing with respect to quantities measured at T = 0 and at the same UV cutoff. This prescription neglects contributions to the additive renormalization which are of higher order in the light quark mass; it is therefore particularly well suited for the present study, in which we consider the approach to the m l = 0 limit. In order to determine the relevant quantities at T = 0, we have exploited the same set of runs already used for the determination of the pion masses and of the physical scale; determinations at intermediate values of the inverse gauge coupling have been obtained by spline interpolations. An example of the renormalized chiral condensate obtained for am l = 0.0015, and expressed as a function of T , is shown in Fig. 7, where determinations corresponding to different spatial extents L s are present. It is quite clear from the figure that the dependence on L s is not negligible and larger in the region around and below the critical temperature. For this reason, before performing an analysis of the approach to the chiral limit, we have extrapolated the chiral condensate to the infinite volume limit at each value of the temperature. The chiral condensate is not an order parameter for the Roberge-Weiss transition, therefore, for finite quark mass, it is expected to have a smooth approach to the thermodynamic limit; however the degree up to which chiral degrees of freedom are entangled in the Roberge-Weiss transition is not known, moreover the behavior could be already affected by the closeness of the chiral transition. For these reasons, the extrapolation to the thermodynamic limit has been performed, for each temperature, trying different fitting functions which assume either an exponential (in L s ) suppression of finite size effects or power law correc-tions in the spatial size: the error on the final extrapolation, which is reported in Fig. 7 as well, takes into account the spread among the different fitting ansätze as a source of systematic uncertainty, and is particularly more pronounced around and below the critical temperature. The infinite volume extrapolations obtained for the different quark masses are reported in Fig. 8, where data are also fitted to an arctangent function, A = P 1 + P 2 arctan (P 3 (T − T χ (am l ))), obtaining the values of the chiral transition temperature T χ (am l ) reported in Table IV and diplayed, together with the Roberge-Weiss critical temperatures T RW , in Fig. 9. One aspect which is already clearly visibile is that T χ (am l ) is very close to T RW (am l ) and, even if it the two temperatures are actually always compatible within errors, they seem to approach each other more closely as the chiral limit is approached.
This is already a good piece of evidence for the coincidence of T RW and T χ in the chiral limit. However, in order to complete the picture, one would like to know if the drop of the condensate at T χ (am l ) is actually associated to a critical behavior around T χ (am l = 0), corresponding to the vanishing of the condensate and the restoration of chiral symmetry at that point. Trying to answer this question, one can also obtain information about the universality class.
The critical temperatures themselves do not provide much information. Around the chiral transition the pseudo-critical temperatures obtained for finite quark mass, are expected to scale like T χ (am l ) = T χ (0) + C · (am l ) 1/(βδ) (10) where β and δ are the critical indexes of the relevant universality class. Two possibilities that we have taken into account are the 3D O(2) and Z 2 critical behaviors: the first one is naturally associated with a second order chiral transition in the presence of just one Goldstone pion (it would be O(4) in the continuum case, which however has practically indistiguishable critical indexes); the second is the relevant universality class for the RW transition and would also be associated with a critical endpoint of a first order line present at very small quark masses.  The corresponding critical indexes are reported in Table V. In Fig. 9 we report best fits of T RW (am l ) according both to the critical behavior in Eq. (10) and to a regular behavior T RW (am l ) = T RW (0)+C(am l )+O((am l ) 2 ). As one can easily appreciate, even if the two ansätze lead to different chiral extrapolations, they are not distinguishable in the quark mass range which has been actually explored and both fits yield acceptable values of the chisquared test. Thus, we cannot state, just according to T RW (am l ), if the explored transition is entangled with a chiral critical behavior as am l → 0; this is similar to what happens at µ B = 0, where the analysis of T c (am l ) alone is not enough to fix the universality class of the chiral transition [89].
Therefore, we turn our attention to the order parameter for chiral symmetry, i.e. the chiral condensate, which around a chiral transition and in the chiral restored phase is expected to scale like [90,91] ψ ψ r (T, am l ) = (am l ) 1/δ φ (T − T χ )(am l ) −1/(βδ) (11) where φ is an appropriate scaling function. In Fig. 10 we show a plot of ψ ψ r (T, am l )(am l ) −1/δ (extrapolated to the infinite volume limit) versus (T − T χ )(am l ) −1/(βδ) for both choices of critical indexes (O(2) or Z 2 ). The value of T χ has been chosen in both cases so as to maximize the collapse of the condensates obtained at different values of am l , obtaining T χ = 169.2 MeV and T χ = 169.6 MeV respectively for O(2) or Z 2 . Such values are compatible with those obtained by fitting directly T χ (am l ) and the observed scaling is pretty good for both universality classes: this is also due to the fact that the critical indexes δ and β are quite similar for O(2) and Z 2 .
Therefore, our present results are consistent with a scenario in which chiral symmetry is restored exactly at T RW in the chiral limit, i.e. T RW = T χ . In order to distinguish the correct universality class one should explore quantities characterized by different critical indexes, like the specific heat, which however are less trivial to determine; in this respect, the situation is quite similar to the present status of the determination of the universality class of the chiral transition at zero chemical potential.
Of course, our conclusions are still an extrapolation of results obtained at finite, even if small, quark masses, i.e. one cannot completely exclude a priori that going to even lighter quark masses T RW and T χ separate. Moreover, results obtained at finer lattice spacing (i.e. at larger values of N t ) could in principle be different.

V. DISCUSSION AND CONCLUSIONS
We have investigated the fate of the Roberge-Weiss endpoint transition and its relation with the restoration of chiral symmetry as the chiral limit of N f = 2 + 1 QCD is approached. The study has been performed on lattices with N t = 4 sites in the temporal direction, a stout staggered discretization for the fermion sector and the tree level Symanzik improved action for the pure gauge sector. We have worked at fixed values of the bare quark masses around the transition points, in order to easily exploit multi-histogram methods, maintaining a physical strange-to-light mass ratio (m s /m l = 28.15) and exploring three different light quark masses, am l = 0.003, 0.0015 and 0.00075, corresponding respectively to pseudo-Goldstone pion masses m π ≃ 100, 70 and 50 MeV around the transition.
The imaginary quark chemical potential has been fixed to µ f,I /T = π for all flavors, so that the imaginary part of the Polyakov loop has been taken as an order parameter for the RW transition. An analysis of the finite size scaling of its susceptibility has excluded the presence of a first order transition for all values of the quark mass; this fact has been confirmed by an inspection of the probability distribution of the plaquette and of the chiral condensate at the transition points, which have revelead no double peak structures as the thermodynamic limit is approached. On the contrary, a good scaling has been observed according to the predicted second order critical behavior, i.e. that of the three dimensional Z 2 (Ising) universality class.
Therefore, our results still provide no evidence of a first order region around the chiral point for the RW transition, which for N f = 2 unimproved staggered fermions was located below m π ≃ 400 MeV [40] for N t = 4. A strong cutoff dependence of the tricritical pion mass has been found also in studies with Wilson fermions [42,49], however it is striking that the tricritical pion mass can go down at least one order of magnitude (or disappear at all) by just improving the discretization at fixed N t .
Our results clearly need further refinement in some respects. Indeed, because of the taste symmetry breaking of staggered fermions, the chiral limit is approached only by the pseudo-Goldstone pion directly linked to the residual staggered chiral symmetry, while all the others stay above 400 MeV and do not seem to be much affected by the am l → 0 limit (see Fig. 2). Therefore, even if not looking quite natural, it cannot be excluded apriori that, as the continuum limit is approached and the full set of chiral degrees of freedom come into play, the critical behavior changes and the (possible) first order region around the chiral point enlarges again. Unfortunately, exploring larger values of N t while approaching the chiral limit would require computational resources which are presently not available to us.
In spite of these caveats, thanks to the exact residual chiral symmetry of staggered fermions, we have been able to answer a different but related question regarding the relation between the RW transition and the chiral restoration transition. In the chiral limit both symmetries are exact and predict the existence of a phase transition with well defined critical temperatures, T RW and T χ : whether the two transitions coincides and, in this case, which symmetry controls the critical behavior, is a clearcut question. Our results have shown that, for all explored masses, the renormalized chiral condensate drops sharply with an inflection point in coincidence (within error bars) with the location of RW endpoint transition; morover, the behavior of the condensate around the transition for the different masses and temperatures scales consistently with a critical behavior corresponding to chiral symmetry restoration in the chiral limit (see Eq. (11) and Fig. 10). Therefore, our results are consistent with T RW = T χ . Regarding the critical behavior, we have not been able to distinguish between an O(2) (chiral) or Z 2 (Roberge-Weiss) universality class in the chiral limit, mostly because the critical indexes associated with the chiral order parameter (δ and β) are almost coincident (in comparison with our numerical precision) in the two cases. Of course, our present results do not exclude that the situation might be different as the continuum limit is approached.