Towards the chiral phase transition in the Roberge-Weiss plane

We discuss the interplay between chiral and center sector phase transitions that occur in QCD with an imaginary quark chemical potential $\mu=i(2n+1) \pi T/3$. Based on a finite size scaling analysis in (2+1)-flavor QCD using HISQ fermions with a physical strange quark mass and a range of light quark masses, we show that the endpoint of the line of first-order Roberge-Weiss (RW) transitions between center sectors is second order for light quark masses $m_l\ge m_s/320$, and that it belongs to the $3$-d, $Z(2)$ universality class. The operator for the chiral condensate behaves like an energy-like operator in an effective spin model for the RW phase transition. As a consequence, for any non-zero value of the quark mass, the chiral condensate will have an infinite slope at the RW phase transition temperature, $T_{RW}$. Its fluctuation, the disconnected chiral susceptibility, behaves like the specific heat in $Z(2)$ symmetric models and diverges in the infinite volume limit at the RW phase transition temperature $T_{RW}$ for any non-zero value of the light quark masses. Our analysis suggests the critical temperatures for the RW phase transition and the chiral phase transition coincide in the RW plane. On lattices with temporal extent $N_\tau=4$, we find in the chiral limit $T_{\chi}=T_{RW}=195(1)~$MeV.


I. INTRODUCTION
Strong interaction matter is described by Quantum Chromo Dynamics (QCD). One of the central problems in studies of QCD at non-zero temperature (T ) and nonzero values of the chemical potentials (µ f ) for the different quark flavors (f ) is to establish the phase diagram of the theory. Aside from the external control parameters (T, µ f ) this also includes an exploration of the dependence of the phase diagram on the quark masses m f . In particular, it is a long standing open question whether in the limit of vanishing light quark masses 1 , m u = m d = 0, the chiral phase transition is of first or second order [1]. In the former case a critical light quark mass, m crit l , would exist, where the first-order transition terminates. Evidence for such a scenario has been found in calculations on coarse lattices that used an unimproved staggered fermion discretization scheme, both for N f = 3 [2][3][4] and N f = 2 [5] degenerate light quarks. This firstorder phase transition, however, turns out to be strongly cut-off dependent [6], with similar behavior found for N f = 3 O(a)-improved Wilson fermions [7][8][9]. In fact, recent chirally extrapolated results seem to rule out a first-order phase transition in the continuum for number of light flavors being smaller than six [10]. Furthermore, calculations performed with improved staggered 1 Here and in the following we will always consider degenerate values for the light up and down quarks, m l ≡ mu = m d .
fermions, using the Highly Improved Staggered Quark (HISQ) [11] or stout [12] discretization schemes, so far did not find any evidence for the existence of a first order transition region [13][14][15][16][17]. In cases where a region of first order transitions has been found, it has also been observed that the upper bound, m crit l , for such a transition region increases when simulations are performed with a non-zero, purely imaginary chemical potential, both for staggered [4,5,18] and Wilson [19] discretizations. When searching for a possible region of first order transitions in simulations with an improved staggered fermion action, it thus is meaningful to try to establish the existence of such a region in simulations with an imaginary chemical potential, µ ≡ iµ i . QCD thermodynamics at imaginary chemical potential has a rich phase structure on its own, based on two exact symmetries, which hold for any quark mass configuration: The Z(3) periodicity of the QCD partition function with imaginary quark chemical potential in Eq. 2 corresponds to the global center subgroup of the SU (3) gauge symmetry [20]. While thermodynamics is invariant under the shifts in Eq. 2, the phase of the Polyakov loop distinguishes between the different center sectors, as indicated in Fig. 1. At the boundaries between the center sectors, (µ i /T ) RW ≡ (2n + 1)π/3, the system undergoes a first-order phase transition at high tempera-End or meeting points First order lines Figure 1. The QCD phase diagram with imaginary chemical potential. Vertical lines mark first-order transitions between different center sectors, the dotted lines represent the analytic continuation of the thermal transition at zero and real µ, whose nature depends on the quark masses. The black dots, where the RW transitions terminate, can then be first-order triple points, tricritical points or second-order endpoints.
The question whether a non-zero value m crit l of the light quark masses exists, below which QCD undergoes a first order chiral phase transition, may then be rephrased somewhat differently at (µ i /T ) RW , where it is related to the nature of the endpoint of the line of first order phase transitions between center sectors at T RW [23,24]. The dotted line in Fig. 1 represents the analytic continuation of the chiral transition line T χ (µ i ) to imaginary chemical potentials. For intermediate quark mass values, where this is a crossover, the endpoint of the RW transition is of second order in the 3-d, Z(2) universality class. On the other hand, if the chiral transition is of first order, as is the case for unimproved actions on coarse lattices, the endpoint of the RW-transition represents a first-order triple point. The boundary between these scenarios, corresponding to a specific quark mass value m tric RW , is marked by a tricritical RW endpoint. The nature of the RW endpoint thus depends on the quark mass configurations (m l , m s ), constituting the so-called Roberge-Weiss plane, which is analogous to a Columbia plot for (µ i /T ) RW . All three situations have been observed explicitly as a function of quark mass for N f = 2 unimproved Wilson [25,26] and staggered [27,28] fermions on N τ = 4, 6, again with a strong cutoff dependence of the first-order region.
Here we address the nature of the RW endpoint using simulations of HISQ fermions on N τ = 4 lattices for a physical strange quark mass and a sequence of decreasing light quark masses. As the chiral limit is approached, it is conceivable that the chiral and RW transitions split, so that the cusps of the dotted line in Fig. 1 would not be connected to the first-order RW-lines. The answer to the question whether or not T χ = T RW , when the chiral transition intersects (µ/T ) RW , does not seem to be obvious and we want to address it here. More generally, we will analyze the influence of the RW phase transition on chiral observables at non-zero values of the quark masses, determine their quark mass dependence and explore the interplay between the RW and the chiral phase transition. Our findings based on simulations with the HISQ action are fully compatible with those of a similar previous study using stout-smeared staggered fermions [14,29].
This paper is organized as follows. In the next section we describe the lattice setup for our calculations and introduce the basic observables we are going to analyze. In Sec. III we present results on the endpoint in the second (n = 1) Roberge-Weiss plane. Sec. IV is devoted to an analysis of chiral observables and the chiral phase transition in the limit of vanishing light quark masses. We present our conclusions in Sec. V. Two appendices are devoted to a summary of all our simulation parameters (Appendix A) and a new analysis and parametrization of finite size scaling functions for the order parameter, its susceptibility and the Binder cumulant in the 3-d, Z(2) universality class (Appendix B).

II. LATTICE SETUP AND OBSERVABLES
All our calculations are done on lattices of size N 3 σ ×N τ with N τ = 4 and a non-zero imaginary value of the chemical potential, µ/T = i(µ i /T ) RW = i(2n + 1)πT /3, corresponding to the (n + 1) th Roberge-Weiss plane [20]. Although thermodynamics is equivalent in all RW planes, we will work in the second RW plane (n = 1). We use the tree-level HISQ action [11] and a tree-level improved Symanzik action in the gauge sector. This setup eliminates O(a 2 ) discretization errors at the tree-level. At vanishing values of the chemical potential the HISQ action has been used extensively in studies of the thermodynamics of QCD with two degenerate light up and down quarks (m l ) and a heavier strange quark mass (m s ) [30,31]. In particular, it recently has been used to study thermodynamics close to the chiral limit of QCD with a strange quark mass tuned to its physical value and light quark masses corresponding to a light Goldstone pion mass in the range 55 MeV < ∼ m π < ∼ 140 MeV [15,32,33].
In our simulations we follow a line of constant physics determined in calculations with the same action at vanishing chemical potential and for physical values of the light and strange quark mass [31]. We use the parametrization given in [34] that uses the kaon decay constant f K = 155.7(2)/ √ 2 MeV [35] to set the scale. The line of constant physics has been determined by tuning the strange quark mass to its physical value. The physical value of the light quark mass is then chosen as 1/27  135  16, 24, 32  1/160  55  16,24, 32  1/320  40  24, 32   Table I. Lattice sizes N 3 σ × 4 used for calculations with quark mass ratios H = m l /ms. The second column gives the corresponding pseudo-scalar Goldstone masses mπ. m l = m s /27. As the crossover temperature at non-zero imaginary chemical potential is shifted towards larger values than those at µ = 0, our calculations are performed in a range of gauge couplings, 5.8 ≤ β ≤ 6.1, which corresponds to a parameter range typically studied in finite temperature calculations at µ = 0 on lattices with temporal extent N τ = 6 [30]. For this range of couplings the line of constant physics, defined by a strange quark mass tuned to its physical value, thus is known from these studies.
When approaching the chiral limit, we vary the light quark mass, while keeping the strange quark mass fixed to its physical value. The partition function for (2+1)flavor QCD with two degenerate light quark masses, a strange quark mass, and identical quark chemical potentials, (µ i /T ) RW = π, for all flavors may be written as, Here M f = D(iµ i )+m f ·1 is the staggered fermion matrix for quarks of mass m f that is obtained after integrating out the fermion fields of the HISQ action [11] and S G is the tree level improved Symanzik gauge action. Further details on the HISQ action as used by us can be found in [30,31]. To perform calculations at non-zero values of an imaginary chemical potential we only needed to replace all temporal link variables, U0(n 0 , n) by e iμi U0(n 0 , n), with iμ i denoting the imaginary quark chemical potential in lattice units, i.e.μ i N τ ≡ (µ i /T ) RW . Here Uν(n 0 , n) ∈ SU (3) are the standard gauge link variables defined on a 4-dimensional lattice and pointing into directionν at a lattice point (n 0 , n).
The light quark masses have been varied towards the chiral limit, starting at the physical value, m l = m s /27. The smallest quark mass value used in our calculations, m l = m s /320, corresponds to a Goldstone pion mass of about 40 MeV. A summary of the quark mass values and lattice sizes used in our simulations is given in Tab. I. More details on the parametrization of the line of constant physics are given in [30,31,34].
For each value of the light to strange quark mass ratio, H = m l /m s , we typically perform calculations at 8-10 values of the gauge coupling that have been chosen such as to cover temperatures in the range, 0.9 < ∼ T /T RW (H) < ∼ 1.1, with T RW (H) denoting the RW phase transition temperature for given H. At temperatures close to T RW (H) we generated about 200,000 Hybrid Monte Carlo trajectories of unit length for H = 1/27 and half-unit length for smaller H. Away from the critical region fewer trajectories have been generated. Further details on the statistics collected for different lattice sizes and run parameters are given in Appendix A. All basic observables needed for the analysis presented here have been calculated at the end of each Rational Hybrid Monte Carlo (RHMC) trajectory.

A. The Polyakov loop and its susceptibility
In the Roberge-Weiss plane [20] a first-order phase transition occurs at all temperatures above a critical temperature T ≥ T RW . An order parameter for the spontaneous breaking of the Z(2) symmetry between different center sectors is the expectation value of the imaginary part of the Polyakov loop, P , While the Euclidean action of QCD is invariant under center transformations of the gauge links, Uν(n) → e 2πk/3 Uν(n), an additional phase remains in the Polyakov loop, whose imaginary part thus changes its sign. The RW transition can thus be mapped straightforwardly to an Ising-like effective Hamiltonian, with energy-like and magnetization-like operators In the following, we will use the imaginary part of the Polyakov loop P for the magnetization, as it couples to the magnetic field-like coupling 2 h = µ i /T − (µ i /T ) RW . An alternative order parameter would, e.g., be the quark number density n f = (T /V )∂ ln Z/∂µ f , which is odd under the Z(2)-reflections about the center boundary. It has been used to study the scaling in the vicinity of the RW transition in [29,[36][37][38]. The energy-like operator will be a superposition of lattice QCD operators which is even under the Z(2)-reflections about the center boundary.
As we perform calculations on finite lattices at vanishing external field h, the expectation value of the imaginary part of the Polyakov loop, ImP , vanishes exactly at any value of the temperature. In this case an approximate order parameter for spontaneous symmetry breaking, commonly used in finite-size scaling studies, is the absolute value of ImP . For the finite-size scaling analysis of the RW phase transition we thus will use as an order parameter M and its susceptibility χ M the observables, We also calculate two ratios of moments of the order parameter, involving first, second and fourth order moments, respectively. We consider the order parameter ratio introduced by Kiskis [39,40], and the kurtosis, related to the Binder cumulant [41], In the scaling regime of a 2 nd order phase transition, where contributions from regular terms are negligible, these ratios are scaling functions with a fixed absolute normalization. Ratios calculated on different size lattices have unique crossing points at a phase transition temperature, T c . The values at these crossing points are universal, i.e. characteristic for a given universality class. From a new determination of finite-size scaling functions of the 3-d, Z(2) universality class, presented in Appendix B, we obtain for these universal crossing points K 2 (T c , ∞) = 0.240(3) and B 4 (T c , ∞) = 1.606 (2). In order to determine the RW phase transition in the RW-plane for different values of the quark mass ratio H we perform a finite-size scaling analysis and compare our results, obtained at fixed H and for several spatial lattice sizes N σ , with scaling functions in the 3-d, Z(2) universality class. If the transition turns out to be second order for the quark mass value H under consideration, these scaling functions will describe the behavior of M and χ M in the vicinity of the RW critical point, (T → T RW , N σ → ∞).
For vanishing external field, h = 0, finite-size scaling functions are functions of a single scaling variable, Here the non-universal parameters l 0 , t 0 and T RW are functions of H, and ν is a critical exponent of the 3-d, Z(2) universality class, which is the relevant universality class for any H = 0, if the RW transition turns out to be 2 nd order.
In the limit (T → T RW , N σ → ∞) two finite-size scaling functions, f G,L and f χ,L , control the behavior of M and χ M , respectively, with critical exponents β, γ, ν which are related to each other through the hyper-scaling relation, γ + 2β = dν, where d = 3 in our case. The relation between the amplitudes of the singular parts follows from Eq. 8.
The scaling functions f G,L and f χ,L for the 3-d, Z(2) universality class have been determined by us from a new finite-size scaling analysis for the 3-d, λφ 4 model, where the coupling λ has been tuned to suppress contributions from corrections-to-scaling in the scaling functions [42,43]. We summarize results on the finite-size scaling functions f G,L and f χ,L as well as the scaling functions for the ratios K 2 and B 4 in Appendix B.

B. Chiral condensate and chiral susceptibility
In the limit of vanishing light quark masses, m l , and for any value of the strange quark mass (which also may be infinite), a chiral phase transition occurs in (2 + 1)-flavor QCD for any value of the imaginary chemical potential µ i . In any Roberge-Weiss plane at µ i /T = (2n + 1)π/3, however, the standard chiral observables, i.e. the chiral condensates and the related chiral susceptibility, also are sensitive to the Roberge-Weiss phase transition, which persists also for any non-zero value of the light quark masses. In the vicinity of this phase transition chiral observables thus show "unconventional" behavior. In fact, we will show in the next section, that chiral observables behave like energy-like observables [32,44] in an effective theory for the Roberge-Weiss phase transition, Eq. 5.
For the study of the properties of the chiral phase transition we use the additively and multiplicatively renormalized order parameter, introduced in [45], which, for light up (u) and down (d) quark masses, is given in terms of the light quark condensate ψ ψ l = ( ψ ψ u + ψ ψ d )/2 and the strange quark condensate ψ ψ s . The condensates are obtained as derivatives of the partition function, Z(T, V, m u , m d , m s ), with respect to one of the quark masses, m f , where M f denotes the staggered fermion matrix for quark flavors, f = u, d or s. For degenerate light quark masses we obviously have ... u = ... d ≡ ... l and we also use the short-hand notation, M l ≡ M u = M d . Furthermore, we used in Eq. 12 the kaon decay constant, f K , for normalization and to make the order parameter dimensionless.
For the determination of a pseudo-critical temperature for chiral symmetry restoration we use the fluctuations of the light quark chiral condensate, i.e. the so-called disconnected part of the chiral susceptibility, which for two degenerate light quark flavors is given by, The fluctuations of the chiral order parameter, χ dis , only give rise to a part of the total chiral susceptibility, Both susceptibilities have pronounced peaks as function of temperature which commonly are used to define a pseudo-critical temperature T χ . The location of pseudo-critical temperatures determined from χ dis and χ m , respectively, coincide in the chiral limit. As long as the susceptibilities are not influenced by the presence of another critical point, their quark mass dependence at low temperature, in the chiral symmetry broken phase, is dominated by contributions from the light Goldstone modes with mass m G , which become massless in the chiral limit, m G ∼ √ m l . A consequence of this is that the leading quark mass dependence of the chiral order parameter is proportional to √ m l . The chiral susceptibility consequently diverges for all T < T χ , The location of a peak in either χ m or χ dis defines a pseudo-critical temperature for a chiral transition. In the absence of any influence from other singularities the peak height of chiral susceptibilities stays finite at any nonzero value of the light quark masses (H > 0). In fact, it is a characteristic feature of phase transitions in the 3-d, O(N ) universality classes that for H > 0 the maxima of the order parameter susceptibility decrease with increasing volume. In the thermodynamic limit they thus approach the asymptotic value from above [46]. In the chiral limit, H → 0, the chiral susceptibility diverges like χ dis ∼ H 1/δ−1 . We note that this divergence is stronger than that induced by the Goldstone modes. Thus also the peak in the rescaled disconnected chiral susceptibility, H 1/2 χ dis will diverge in the 3-d, O(N ) universality class.
In the vicinity of the Roberge-Weiss phase transition the behavior of the chiral observables also is strongly influenced by this transition. In particular, they will show a volume dependence that is quite different from what one would expect for a chiral observables in the O(N ) universality classes. The chiral condensate as well as the chiral susceptibility are even operators under the Z(2) symmetry transformation that controls the Ising-like Roberge-Weiss transition for n = 1, Uν(n) → U † ν (n). In the vicinity of the RW-transition temperature, the chiral condensate thus may be considered as an energy-like operator in an effective Hamiltonian for the RW-transition. Since non-vanishing quark masses, which appear as couplings in the QCD Lagrangian, do not break the center or reflection symmetry that controls the RW-transition, these couplings (masses) will appear, to leading order, only in the energy-like coupling (reduced temperature) of the effective Hamiltonian describing the RW-transition. As a consequence the chiral condensate will behave like an energy and the chiral susceptibility will behave like a specific heat in the scaling regime of the RW phase transition. For any non-zero value of the light quark masses, H > 0, this will lead to a volume dependence of the disconnected chiral susceptibility that is quite different from that for O(N ) symmetric theories. Rather than staying finite in the infinite volume limit the disconnected chiral susceptibility is expected to diverge at T RW for any H > 0, with α being the specific heat critical exponent in the 3-d, Z(2) universality class (see Appendix B).
As we perform calculations on rather coarse lattices with temporal extent N τ = 4, and as we use a fermion discretization scheme (staggered fermions) which only preserves a U (1) subgroup of the full SU (2) L × SU (2) R flavor symmetry group, it is expected that the chiral phase transition for µ i /T = (µ i /T ) RW belongs to the O(2) universality class, if this transition turns out to be a 2 nd order phase transition. Whether this will also be the relevant universality class to describe the chiral observables in the limit |h| = |µ i /T − (µ i /T ) RW | → 0 crucially depends on the answer to the question whether or not in the RW-plane the chiral transition temperature and the critical temperature for the RW transition coincide. We will address this question in Sec. IV.

III. THE ROBERGE-WEISS ENDPOINT
The endpoint of the line of first-order transitions in the RW plane may be a first or second order transition, depending on the magnitude of the light quark mass m l [24,25,28]. If this transition is of 2 nd order and H ≡ m l /m s > 0, critical behavior in the vicinity of this endpoint will be controlled by the 3-d, Z(2) universality class. The volume dependence of e.g. the maximum of the order parameter susceptibility at the pseudo-critical temperature T pc (N σ ) is distinctively different from that expected in the case of a first-order transition, that may occur at the RW endpoint below some critical value of the quark masses, In order to determine the order of the transition at the endpoint of the RW transition as a function of the light quark mass we performed simulations with H ∈ [1/320, 1/27], which, in the continuum limit, corresponds to pseudo-scalar Goldstone masses (pion masses) simulation parameters and the statistics gathered in our simulations are summarized in Appendix A. In Fig. 2 we show simulation results for the order parameter M and its susceptibility χ M for three quark mass ratios, H = 1/27, 1/160 and 1/320.
It is obvious that for all values of H the peak heights of the order parameter susceptibility (Fig. 2 (bottom)) grow with increasing volume. The growth rate is substantially smaller than N 3 σ , which one would expect to find in the case of a 1 st order phase transition. We also note that the peak height of the susceptibility shows only a weak dependence on the quark mass ratio H.
In Fig. 3 we show results for the order parameter susceptibility, rescaled by the volume factor expected to be relevant for a first-order transition (top) and a Z(2) second-order phase transition (bottom), respectively. The rescaled results, obtained for H = 1/160 on lattices with spatial extent N σ = 16, 24 and 32, clearly show that at the pseudo-critical temperature, T RW (N σ ), a first-order volume scaling of the peak height is ruled out. The observed volume dependence is in good agreement with the expected 3-d, Z(2) scaling behavior. We find similar results also at our smallest quark mass value, H = 1/320. We thus proceed in the following with an analysis of our data in terms of 3-d, Z(2) scaling functions, taking into account deviations from universal scaling that may arise from corrections-to-scaling or regular contributions to the free energy density, wherever needed.
In order to minimize contributions from such terms we only fitted results obtained on lattices with N σ ≥ 24. We first fitted only data for the order parameter using the ansatz M = A M N −β/ν σ f G,L (z f ), which depends on three non-universal parameters, A M , T RW , and z f,0 . As regular contributions to the order parameter would be linear in the external field h, one expects that regular contributions do not contribute to M in the case of h = 0. Indeed we find that a fit using the universal scaling ansatz only already gives a good description of the order parameter calculated on smaller lattices with extent N σ = 16. This suggests that contributions from corrections-to-scaling or regular terms are indeed negligible for the order parameter. We then used the parameters A M , z f,0 and T RW obtained from fits to the order parameter to compare data for χ M as well as the cumulant ratios K 2 and B 4 with the corresponding scaling ansätze. While we find that these ansätze describe the structure of χ M and the ratios K 2 , B 4 , we observe deviations which hint at the importance of contributions from regular terms or correctionsto-scaling. In order to take, on the one hand, information from these observables into account but, on the other hand, keep the number of fit parameters small we then decided to perform our final data analysis using a simultaneous fit to the order parameter and its susceptibility, including a regular contribution to the latter observable, We note that in the QCD case, the symmetry broken phase corresponds to the high temperature region whereas in the 3-d Ising model this is the low temperature region. The sign of the scaling variable z f thus needs to be interchanged when comparing both models, i.e. the non-universal scaling variable t 0 is negative in the case of QCD with an imaginary chemical potential. Consequently the infinite volume limit also is approached differently. From Fig.10 (right), shown in Appendix B, it is obvious that the pseudo-critical temperatures increase with increasing volume in the 3-d Ising model. We thus expect the pseudo-critical temperatures to decrease with increasing volume in the case of QCD. This is apparent from the volume dependence of the peak locations in the susceptibilities shown in Fig. 2.
In the fit ansatz for χ M we add a regular term χ reg (t) = a 0 + a 1 t. This is consistent with the fact that regular terms do not break the Z(2) symmetry and thus are even functions in the external field h. The coefficients a 0 and a 1 will depend on the quark mass ratio H. Furthermore, an implicit quark mass dependence arises through the dependence of T RW on H. Our final fit, per-formed simultaneously to M and χ M , thus depends on five parameters, the non-universal parameters A M , T RW , and z f,0 as well as the leading regular coefficients a 0 and a 1 . The results for the fit parameters in the universal part of the fit ansatz, z f,0 and T RW , are then used for translating results for the cumulant ratios K 2 and B 4 into a scaling plot. These cumulant ratios are not fitted, but serve as a consistency check using the fit parameters we obtain from fits to M and χ M , respectively.
For our fits we use data in a range of temperatures that deviate maximally about 5% from the critical temperature T RW . Details on the fit ranges and the resulting chi-squares of our fits are summarized in Tab. II. The resulting scaling fits are shown in Fig. 2. The solid lines indicate the temperature range and data sets used in the fits. Data outside this region and those obtained on lattices with spatial extent N σ = 16 have not been included in the fits. This is indicated by dashed lines in Fig. 2. Results for the three fit parameters that enter the scaling functions and the two parameters of the polynomial ansatz for the regular part are summarized in Tab. II. In Fig. 4 we show the rescaled order parameter M and its susceptibility χ M where the fitted regular part has been subtracted from the data.
The cumulant ratio K 2 , introduced in Eq. 8, is constructed from M and χ M . It approaches zero at high and (π/2 − 1) at low temperature, respectively [39]. In the infinite volume limit results for different lattice sizes approach a unique crossing point. From our analysis of the 3-d, Z(2) finite size scaling functions, presented in Appendix B, we find K 2 (T RW , ∞) = 0.240(3). This behavior is reflected in the data for K 2 (T, N σ ), shown in Fig. 5 (top). Deviations from the unique crossing point result from the regular contribution to χ M . In the vicinity of T RW the leading correction to scaling arising from the regular term is obtained from Eqs. 8 and 18, The value of K 2 (T RW , N σ ) thus approaches the unique infinite volume value from above, which is consistent with the data shown in Fig. 5 (top). In Fig. 5 (bottom) we show the rescaled data for K 2 (T, N σ ), which only amounts to a change of the temperature axis from T to z f , i.e. it only involves the non-universal parameters, T RW and z f,0 , which are taken from the joint fit to M and χ M .
The qualitative behavior of the Binder cumulant ratio B 4 introduced in Eq. 9 is quite similar to that of the ratio K 2 . It too has a unique crossing point at T RW , but depends on a scaling function f B (z f ) not directly related to those of the order parameter and its susceptibility. While a comparison of our simulation data with the scaling function f B (z f ) is parameter free, taking care of correction to this universal behavior arising from regular terms does require additional parameters entering a fit to the data for B 4 (T, N σ ). We have not done this here. The scaling function f B (z f ), calculated for the 3-d,    Table II. Fit ranges, number of degrees of freedom for the five parameter fits and χ 2 /dof of the joint fits to M and χM for three values of the quark mass ratio H = m l /ms. Results for the non-universal parameters z f,0 and TRW of the scaling functions, the overall amplitude AM controlling the strength of the singular contribution to the order parameter M , and the parameters (a0, a1) appearing in the polynomial ansatz for the regular contribution to χM are given in column 5 to 9 of the table.
Z(2) universality class, is given in Appendix B. It yields for the crossing point B 4 (T RW , ∞) ≡ f B (0) = 1.606(2). In Fig. 6 we give a parameter free comparison of data obtained from our QCD calculations at two different quark mass values to the scaling function f B (z f ). The results shown for the order parameter as well as its susceptibility and the various ratios for order parameter fluctuations clearly show that we do not have any evidence for the occurrence of a first order phase transition for H ≥ 1/320, i.e. for Goldstone pion masses larger than about 3 40 MeV. In fact, these observables show good 3 We stress that we do not perform a continuum extrapolation here. All mass and temperature values expressed in physical units thus should only be used for orientation. agreement with the expected finite size scaling behavior in the 3-d, Z(2) universality class. As can be seen from the fit parameters listed in Tab. II the amplitude A M and the scale parameter z f,0 , which control the peak height of the susceptibility χ M as well as its width show only a mild dependence on the quark mass, i.e. on H ≡ m l /m s . We will come back to a discussion of the quark mass dependence of the Roberge-Weiss transition endpoint temperature after the analysis of the chiral observable in the vicinity of T RW , which is done in the next section. We note in passing that the value for T RW which we obtain at H = 1/27 is consistent with the scaling of the quark number density n f (an alternative order parameter), as well as with the scaling of the Lee-Yang edge singularities extracted from n f by analytic continuation to real µ f [36][37][38].

IV. CHIRAL SYMMETRY RESTORATION IN THE ROBERGE-WEISS PLANE
As discussed in Sec. II.B the chiral condensate and its susceptibility play a particular role when studying the phase structure of QCD with an imaginary chemical potential. On the one hand the chiral condensate is the order parameter for the chiral phase transition. In the chiral limit, it will vanish above the chiral phase transition temperature T χ (µ) for any value of the imaginary chemical potential. In the vicinity of T χ and for nonzero values of the light quark masses, the properties of the order parameter ∆ ls and the chiral susceptibility χ m (or χ dis ) will be controlled by universal scaling relations in the 3-d, O(4) universality class 4 , if the critical region is not influenced by the presence of other singularities. The latter scenario becomes of relevance in the Roberge-Weiss plane, where an additional Z(2) symmetry emerges that leads to a second order phase transition at non-zero values of the quark mass, as discussed in the previous section. It thus needs to be clarified to what extent this transition influences the chiral phase transition, which in turn is reflected in the properties of ∆ ls and χ dis . In Fig. 7 (top) and (bottom) we show results for ∆ ls and χ dis , respectively. The finite size dependence of the chiral condensate follows a pattern, which is qualitatively consistent with the finite size scaling behavior of O(N ) symmetric spin models. However, in particular when looking at the volume dependence of the chiral susceptibility in the vicinity of the Roberge-Weiss phase transition temperature, it is clear that the behavior is distinctively different from that in O(N ) models. In the latter case the peak height of the susceptibilities (slightly) decreases with increasing volume for any non-zero, fixed values of H [15,46]. However, the results for χ dis , presented in Fig. 7 (bottom) for three values of the light quark masses, clearly show an increase of the peak height with increasing volume.
The volume dependence of the peak height is weaker than that of the order parameter susceptibility discussed in the previous section. This is expected for the susceptibility of an energy-like observable. The finite-size scaling of the peak height of a specific-heat-like susceptibility is given by Eq. 16, i.e. χ peak dis ∼ N α/ν σ , with α/ν = 0.1726(65). This should be compared to an exponent γ/ν 1.966, which is the relevant scaling exponent for a magnetization-like susceptibility.
On the lattices analyzed so far the peak height of χ dis increases faster than expected for a specific-heatlike observable in the 3-d, Z(2) universality class. When comparing the peak heights on subsequent lattice sizes, N σ = 16, 24 and 32, for our largest quark mass ratio, H = 1/27, we find for the ratios of peak heights (24/16) 1.38 and (32/24) 0.97 , respectively. Already in this case, and even more so for the smaller quark mass ratios, it is evident from Fig. 7 (bottom) that the disconnected chiral susceptibility suffer from large finite volume effects, which become more severe for smaller quark masses. It thus is conceivable that the lattice sizes used in our current analysis are not sufficient to determine the correct scaling exponents for chiral observables.
A similar behavior is found for the temperature derivative of ∆ ls . Also its peak height increases with increasing volume, suggesting that ∆ ls will have an infinite slope at T RW . The observed volume dependence of ∆ ls and χ dis thus is qualitatively consistent with the finite-size scaling behavior expected for energy-like observables in the vicinity of a phase transition controlled by Z(2) symmetry breaking, although a detailed quantitative analysis will require calculations on larger lattices and smaller quark masses.
It is obvious from Fig. 7 that, unlike the order parameter susceptibility, the disconnected chiral susceptibility shows a strong quark mass dependence. As discussed in Sec. II, in a phase with broken chiral symmetry, i.e. for temperatures T < T χ , the disconnected chiral susceptibility will diverge as H −1/2 for any value of the temperature, due to the presence of a light Goldstone mode. In Fig. 8 (top) we show the disconnected chiral susceptibility on the largest lattice available (N σ = 32) and in Fig. 8 (bottom) we show the rescaled susceptibility, H 1/2 χ dis (T ). At T < T peak χ dis this rescaled susceptibility seems to approach a constant value, while for T > T peak χ dis it approaches zero.
In the vicinity of a critical point belonging to the 3-d, O(N ) universality classes also the peak of the rescaled susceptibility is expected to diverge in the limit H → 0, T → T χ . As can be seen in Fig. 8 the peak height of χ dis clearly rises with decreasing H. However, from the rescaled data shown in that figure one can deduce that at least for fixed finite volume the apparent exponent c controlling the quark mass dependence, is smaller than 1/2. In fact, the ratio of the two subse- While it is thus evident that the chiral susceptibility will diverge in the chiral limit at low temperatures, quantifying the quark mass dependence of the peak height and its interplay with the volume dependence induced  Table III. Pseudo-critical temperatures deduced from the inflection point of the chiral condensate (T infl. ∆ ls ) and the peak of the disconnected chiral susceptibility (T peak χ dis ). Results are from calculations on a 32 3 × 4 lattice.
by the RW transition is beyond the scope of our current analysis. Disentangling the subtle interplay between infinite volume and chiral limits requires further analyses on larger lattices.
Finally we want to discuss whether or not the chiral phase transition temperature T χ and the RW transition temperature T RW coincide. The temperature at the peak of χ dis , i.e. T peak χ dis (H), decreases with decreasing mass. In fact, this closely follows the quark mass dependence found for the critical temperatures determined from the scaling fits to the order parameter M and its susceptibility, χ M . In Tab. IV we summarize results for the location of peaks of the chiral susceptibility as well as the inflec- tion point of the chiral order parameter. These results have been obtained using the largest lattices available, 32 3 × 4. In Fig. 9 we compare the pseudo-critical temperatures obtained from chiral observables on this largest available lattice with the infinite volume critical temperature of the RW phase transition. Obviously the pseudo-critical temperatures agree with each other within errors and also are in good agreement with T RW (H) given in Tab. II. Differences are about 1 MeV and slightly decrease with decreasing H.
Our current data suggest that the chiral and RW phase transitions may coincide in the chiral limit. In that case the RW endpoint will be a multi-critical point [47,48] and one does not expect to find Z(2) or O(N ) critical behavior. The critical temperature generally is expected to scale like In order to extrapolate the RW phase transition temperature and the chiral pseudo-critical temperatures to the chiral limit, H → 0, we performed two parameter fits, using Eq. 21, and keeping the combination of critical exponent βδ fixed. Since we are not aware of a settled set of exponents for this specific multi-critical case, we used the range 0.599 ≤ 1/βδ ≤ 0.8, which covers the values for the three dimensional O(2) (0.599) and Z(2) (0.639) universality classes, as well as mean-field exponents for critical (2/3) and tri-critical (4/5) behavior. Fits with these ansätze indeed give a viable description of the data for T RW (H) and the pseudo-critical temperatures for the chiral phase transition. The extracted critical temperatures grow slightly with increasing 1/βδ. Averaging over the fit results for the four sets of critical exponents given above, we find for the RW phase transition temperature Similarly we find as an estimate for the chiral phase transition temperature from the results for pseudo-critical temperatures obtained on a 32 3 × 4 lattice, The corresponding fits are shown in Fig. 9. In view of the uncertainties regarding our extrapolations and taking into account that our estimates for the pseudo-critical temperatures T χ (H) are not yet extrapolated to the thermodynamic limit, we quote T RW = 195(1) MeV as our first estimate for the phase transition temperature in the chiral limit.

V. CONCLUSIONS
We analyzed scaling behavior in the vicinity of the endpoint of the line of first order phase transitions in the Roberge-Weiss plane. Our calculations, have been performed in (2+1)-flavor QCD using the Highly Improved Staggered Quark action on coarse lattices of size N 3 σ × 4. The strange quark mass has been tuned to its physical value. By varying the spatial extent of the lattice in the range 16 ≤ N σ ≤ 32 we performed a finite size scaling analysis for different values of the light quark mass corresponding to the ratios, H ≡ m l /m s = 1/27, 1/160, 1/320, which in the continuum limit corresponds to pion masses in the range 40 MeV < ∼ m π < ∼ 135 MeV.
We find that the RW phase transition at the endpoint of a line of first order phase transition is second order in the 3-d, Z(2) universality class for the entire quark mass regime explored by us. The chiral order parameter and the related disconnected chiral susceptibility show scaling behavior that is consistent with that of energy-and specific-heat like observables in this universality class. The pseudo-critical temperatures obtained from peaks of the disconnected chiral susceptibility and the susceptibility of the order parameter for the RW transition differ by at most 1 MeV at the quark mass ratios considered by us and show similar quark mass dependence. This suggests that the RW and chiral phase transition temperatures coincide or differ by not more than 1 MeV in the chiral limit. On the coarse lattice used for these calculations we estimate for the RW transition temperature in the chiral limit T RW = 195(1) MeV.

ACKNOWLEDGMENTS
This work was supported in part by the Deutsche Forschungsgemeinschaft (DFG) through the grant 315477589-TRR 211 and the grant 05P18PBCA1 of the German Bundesministerium für Bildung und Forschung and the grant 283286 of the European Union. F.C. and O.P. also acknowledge support by the state of Hesse within the Research Cluster ELEMENTS (Project ID 500/10.006). Numerical calculations have been made possible through Partnership for Advanced Computing in Europe (PRACE) grants at the Swiss National Supercomputing Centre (CSCS), Switzerland, and grants at the Gauss Centre for Supercomputing at Jülich supercomputer center, Germany. These grants provided access to resources on Piz Daint at CSCS as well as on JUQUEEN and JUWELS in Jülich. Additional calculations have been performed on the GPU clusters at Bielefeld University, Germany, the LOEWE-CSC cluster at Goethe University Frankfurt, Germany, and the P C 2 at Paderborn University, Germany. We thank the HPC.NRW team in Bielefeld for its support.   Table IV. Parameters used in simulations with the HISQ action on lattices of size N 3 σ × 4 with Nσ = 16, 24 and 32. Calculations are done on a line of constant physics defined by m l /ms = 1/27 [31]. The strange quark masses, ms in units of the lattice spacing a are given in column 3. In columns 4 to 6 the number of configurations (Uτ ) is given. They have been stored and analyzed after RHMC trajectories of unit length (τ = 1) for m l = ms/27 and half-length (τ = 1/2) for m l = ms/160 and ms/320.  Table V. Parametrization of the finite size scaling functions for the order parameter (fG,L(z f )), its susceptibility (fχ,L(z f )) and the Binder cumulant (fB(z f )). Given are the expansion coefficients a ± X,n used for expansions around z f = 0. The asymptotic coefficients used for z f → ±∞ are given in Eq. B13. For the specific λφ 4 model, with λ = 1.1, used in our simulations, the critical temperature is well-known [49] T c = 2.665980(3) , and the scale parameter t 0 = 0.302(1), which is the only relevant scale parameter for our analysis, is taken from Ref. [43]. Also the critical exponents of the 3-d, Z(2) universality class are known quite accurately. For consistency with [43] we use critical exponents taken from [50] β = 0.3258 (15), ν = 0.6304 (13) , and determine other exponents from hyper-scaling relations, e.g. δ = dν/β − 1,γ = dν − 2β and α = 2 − dν.
Close to the critical point (t, h, N −1 σ ) = (0, 0, 0) the universal critical behavior is controlled by the scaling function f s . Choosing the scale factor b = N σ , taking derivatives with respect to the external field h and finally choosing h = 0 we may write the order parameter M and its susceptibility χ M as, with z f = N 1/ν σ t/t 0 . The scaling functions f G,L and f χ,L are defined such that the amplitude A M introduced in Eqs. 10 and 11 equals unity. The singular part of the order parameter ratio, introduced in Eq. 8, is then just a ratio of these scaling functions, i.e., M,∞ = 1 arises from the normalization of the order parameter scaling function. In the polynomial ansätze for positive and negative z f we also fix the first three expansion coefficients to be identical in both regimes, i.e. a 0− X,n = a 0+ X,n , for n = 0, 1, 2 and X = M, χ, B. The fit parameters for the polynomial ansätze in the central fit interval around z f = 0 are summarized in Tab. V.
We also note that the fit results for a + M,∞ and a + χ,∞ yield a + χ,∞ /(a + M,∞ ) 2 0.58 in good agreement with the infinite volume value of the ratio K 2 in the symmetric phase, i.e., K 2 (T, ∞) = π/2 − 1 for T > T c .