{'fraction_non_alphanumeric': 0.04527067431144254, 'fraction_numerical': 0.016813816877132504, 'mean_word_length': 4.113309352517986, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A QCD chiral critical point at small chemical potential: is it there or not?

For a QCD chiral critical point to exist, the parameter region of small quark masses for which the finite temperature transition is first-order must expand when the chemical potential is turned on. This can be tested by a Taylor expansion of the critical surface (m_{u,d},m_s)_c(mu). We present a new method to perform this Taylor expansion numerically, which we first test on an effective model of QCD with static, dense quarks. We then present the results for QCD with 3 degenerate flavors. For a lattice with N_t=4 time-slices, the first-order region shrinks as the chemical potential is turned on. This implies that, for physical quark masses, the analytic crossover which occurs at mu=0 between the hadronic and the plasma regimes remains crossover in the mu-region where a Taylor expansion is reliable, i.e. mu less than or similar to T. We present preliminary results from finer lattices indicating that this situation persists, as does the discrepancy between the curvature of T_c(m_c(mu=0),mu) and the experimentally observed freeze-out curve.


Introduction
The experimental determination of the QCD phase diagram is underway via extensive heavyion collision programs. At the same time, much effort is being devoted to its theoretical determination via numerical lattice simulations. In the latter case, one can ask the more general theoretical question: what is the phase diagram of 2 + 1-flavor QCD in the (m u,d , m s , µ u,d , µ s , T ) parameter space?
For µ u,d = µ s = 0, theoretical expectations are summarized in Fig. 1 (left). In the chiral (m u,d = m s = 0) and the quenched (m u,d = m s = ∞) corners, order parameters probing the breaking of the chiral and the center symmetries exist, and the symmetry breaking or restoring transitions are firstorder. For intermediate quark masses, simulations show a crossover. This leads to the existence of two lines of critical points, both in the 3d Ising universality class, delimiting the first-order regions. In the quenched corner, the critical line has been studied in [2,3]. In the chiral corner, Fig. 1 (right) shows the result of [4]. Good agreement with expectations was found, including consistency with a tricritical point (m u,d = 0, m s = m tric s ) with a rather heavy mass m tric s ∼ 2.8T c . The physical point, marked by an X, lies in the crossover region. These results were obtained with standard staggered fermions on an N t = 4 (a ∼ 0.3 fm) lattice. One task is now to quantify cutoff effects and extrapolate to the continuum limit. This extrapolation has been performed for the physical point [5], confirming that it remains in the crossover region. Our preliminary results on an N t = 6 lattice, consistent with those of [6] and with earlier indications from improved actions on N t = 4 [7], are presented Sec. 3.5. Cutoff effects are large, and the transition becomes much weaker on a finer lattice.
The next issue is to determine the effect of a baryonic chemical potential. For physical quark masses, it is expected that the transition becomes stronger with µ, so that it turns from a crossover at µ = 0 into a first-order transition, thus defining the QCD critical point (µ E , T E ) where the transition is second-order. The fermion determinant becomes complex when µ = 0, and the ensuing "sign problem" forbids standard Monte Carlo sampling. Nevertheless, a determination of this critical point has been performed, using the same staggered action and N t = 4 lattice spacing as us [8]. While the cutoff error is likely to be comparable to that at µ = 0, the numerical method used,  Figure 2: For physical (m u,d , m s ) quark masses, the finite-temperature "transition" at µ = 0 is really a crossover. As µ is turned on, this crossover becomes a genuine phase transition at the QCD chiral critical point (left), provided the region of first-order transitions at small quark masses expands with µ. If not (right), there is no QCD chiral critical point. The curvature of the critical surface dm c /dµ 2 at µ = 0 distinguishes between these two possibilities for small µ. Note that, for heavy quarks, the first-order region shrinks when µ is turned on [9]. namely reweighting, introduces new possible systematic errors. Contrary to the well-known multihistogram reweighting where one interpolates between Monte Carlo results obtained at several values of the coupling, here one extrapolates results obtained at a single value of the coupling, µ = 0, to µ = 0. While this procedure is in principle exact for infinite statistics, the practical question is whether the finite Monte Carlo sample contains the relevant information at the extrapolated coupling. This difficulty is known as the overlap problem. Thus, a crosscheck of the results of [8] using a different approach seems worthwhile. To avoid all difficulties caused by a potential lack of overlap, we set a more modest goal, and study the effect of an infinitesimal chemical potential. More precisely, we consider the critical surface, swept by the chiral critical line of Fig. 1 (right) as a function of µ, in terms of a Taylor expansion in (µ/T ) 2 about µ = 0, (for the N f = 3 case to keep the notation simple) The sign of the first Taylor coefficient c 1 is of crucial relevance to the QCD critical point. Since for µ = 0 the physical point is in the crossover region, the first-order region must expand with µ ( Fig. 2 left). If instead the first-order region shrinks ( Fig. 2 right), then the crossover persists and there is no QCD critical point unless another critical surface, topologically unrelated to the chiral critical surface we consider, is present in the phase diagram. The simplest corresponding (µ, T ) phase diagram in the two cases is depicted in Fig. 3, for a quark mass below and above m c (0). In the "exotic" scenario (right), the first-order line present for small quark mass turns into a crossover as the mass increases, starting with large µ. If m > m c (0) as in the physical situation, the "transition" is simply a crossover, for all µ's, until different physics take over. Distinguishing between these two scenarios requires full knowledge of the critical surface. However, by measuring the first Taylor coefficient, we can determine the behaviour of the critical surface for µ/T 1, i.e. in the region µ B 500 MeV where experimental searches are considered. In our first study [4], we performed simulations at imaginary µ and determined m c (µ = iµ I ), then fitted our results by a truncated Taylor expansion which can be trivially continued to real µ. While this method appears to work very well for the pseudo-critical line T c (µ) [10,11], we faced two related difficulties when computing m c (µ): the signal is very noisy, and this makes the determination of the systematic error due to truncation of the Taylor series difficult. For the case of three degenerate flavors, m s = m u,d , µ s = µ u,d , we performed quadratic and quartic fits in µ/T . The quartic term was statistically insignificant, so we set it to zero and obtained for the leading term which favors the "exotic" scenario Fig. 2 (right), but not even at the 2σ level. Note that including the quartic term in the fit ( Table 2, line 2 of Ref. [4]) changes the leading coefficient to -2.6(1.2). The lack of compelling numerical evidence motivated us to improve our statistics and our numerical methods.
Here, we compare three methods to obtain the derivative dm c dµ 2 | µ=0 : A. Analytic continuation from imaginary µ B. Direct measurement of the derivative at µ = 0 C. Noisy reweighting to small imaginary µ In all three methods, the same criterion is used to determine the pseudo-critical couplings in a finite volume: the Binder cumulant B 4 ≡ (δ X) 4 (δ X) 2 2 , with δ X = X − X and X =ψψ, takes the value 1.604 characteristic of the 3d Ising universality class at the critical point. In a finite volume, B 4 is an analytic function of (m, T, µ), and T is fixed to the pseudo-critical temperature T c (m, µ) by requiring, for instance, that the susceptibility (δ X ) 2 be maximum. We use the equivalent prescription (δ X ) 3 = 0. The two constraints {B 4 = 1.604, (δψ ψ) 3 = 0} define a line in the (T, m, µ) space. We want to extract the curvature dm c dµ 2 of this line at µ = 0. Method A is the one used in [4]. The results of different simulations at various imaginary values of µ are fitted with a truncated Taylor expansion about µ = 0. One difficulty comes from the systematic error associated with the order of the truncation. Another comes from the final statistical error, which in [4] amounted to about half the signal.
Method B is reminiscent of the Bielefeld-Swansea approach [11], where the derivatives of the free energy with respect to µ/T are expressed as expectation values of non-local operators, which involve traces of inverse powers of the Dirac operator. In our case, the derivatives of B 4 and (δψψ) 3 with respect to m, T and µ can be expressed as expectation values of complicated operators, which can be measured in a single simulation at µ = 0.
Method C, to our knowledge, is a new application of reweighting, where the reweighting factors are not evaluated exactly but estimated by a stochastic method. This noise does not prevent reweighting provided it is unbiased. We can apply this strategy to the results of µ = 0 simulations and reweight them to imaginary µ, monitoring the change in B 4 and (δψψ) 3 . Keeping this imaginary µ very small has two advantages: (i) the reweighting factor remains close to 1 and the overlap of the µ = 0 Monte Carlo ensemble with the target µ = 0 ensemble is guaranteed; (ii) fluctuations in the two ensembles are strongly correlated, and cancel in the observables.
In Section 2, we first test methods B and C on a Potts model representative of the heavy-quark limit of QCD, for which we have previously determined the critical line for both imaginary and real µ [9]. All methods agree, with a clear advantage to method C for its simplicity. Then in Section 3, we compare methods A and C for N f = 3 QCD, on a coarse 8 3 × 4 lattice. Again, methods A and C agree, leaving no doubt that for this system the first-order region shrinks when the chemical potential is turned on. Finally, in Section 4 we present some preliminary results on a finer 18 3 × 6 lattice, and discuss the implications and limitations of our findings.

Potts model 2.1 An effective description of dense static quarks
When the quark mass m becomes infinite, quarks become static and a quark source is a Polyakov loop. The canonical partition function of QCD with n static quarks andñ static antiquarks is [12] is the Polyakov loop integrated over space and S g the gauge action. The grandcanonical partition function is then simply Simplifying the gauge action to a Potts interaction between neighbouring Polyakov loops, one obtains the partition function of a 3d q = 3 Potts model in an external field: For heavy masses, the transition is first order. It turns into a crossover for light masses. The critical line can be determined for imaginary µ (red points), showing the expected Z 3 transition at µ/T = i π 3 , or directly for real µ, because the sign problem is mild. The line is analytic at µ = 0, but shows some curvature. Methods B and C measure its Contrary to [12] who consider the limitm,μ → ∞, (m −μ) finite, and drop the last term corresponding to antiquarks, we keep the complete expression in order to study the (m,μ) phase diagram. Because h ′ = h * when µ = 0, the action becomes complex and a sign problem appears. However, it is very mild, so that simulating with the real part of the action and reweighting with the imaginary part works efficiently and reliably. Two phases are present: confined/disordered at small κ, deconfined/ordered at large κ. The phase transition is first-order for the ordinary Potts model, i.e., in the absence of any external field h = h ′ = 0. This case corresponds tom → ∞.
Decreasingm turns on the magnetic field. This weakens the phase transition, so that it becomes a crossover for sufficiently large h m , i.e., sufficiently smallm = m/T . The critical value ofm corresponding to the end of the first-order region depends on the chemical potentialμ. The critical linē m c (μ) has been determined in [9], for both real and imaginary chemical potentials. Fig. 4 summarizes the results of [9]. It represents the qualitative behaviour of QCD with heavy quarks near the critical line in the upper right corner of Fig. 1 (left), because the symmetries and the infrared degrees of freedom are the same. Given T c ≈ 270 MeV, the value of m c (0) even lies in the 1-2 GeV range estimated from full QCD simulations [2,3]. Fig. 4 shows that, as the chemical potential is turned on, the first order region shrinks as indicated in the quenched corner of Fig. 2 (right).
As a function ofμ 2 , the critical line is analytic atμ = 0, so that the real-μ dependence can be reconstructed by analytic continuation of imaginary-μ simulation results. However, some curvature is visible, particularly as the imaginaryμ approaches the Z 3 transition pointμ = i π 3 . A global fit of all our data requires a fourth-order polynomial inμ 2 : where the first coefficient is 0.585(3) including its fitting error. This is essentially the result of method A (although we have actually augmented the imaginary-µ data with real-µ data). One may be concerned that a similar curvature occurs in QCD as well, which would make it difficult to analytically continue data from a few discrete values of imaginary µ. Therefore, the object of this study of the Potts model is to see how well the first coefficient of the polynomial, 0.585(3), can be reproduced by methods B and C which involve simulations at µ = 0 only.

Method B
The order parameter probing the Z 3 symmetry is the magnetization M ≡ ∑ x Φ(x). Therefore, we determine the critical line by requiring Staying on the critical line implies that the changes in B 4 and in (δ M) 3 are both zero. From the resulting two linear equations, one can eliminate dκ, and obtain which is the desired slope of the critical line at µ = 0. The coefficients A, B,C, A ′ , B ′ ,C ′ are the following expectation values Thus, all 6 coefficients can be measured during a Monte Carlo simulation at the µ = 0 critical point, i.e., where κ andm have been tuned to satisfy eq.(2.5). We have performed such measurements on a 72 3 lattice, of the same size used to obtain Fig. 4. Substituting into eq.(2.7), and performing a jackknife bin analysis, we obtain dm/dμ 2 = 0.593 (8), which agrees very well with the earlier result from method A. Note again that method B is insensitive to the curvature of the critical line, unlike method A.

Method C
The previous method enforced algebraically that the Binder cumulant B 4 stay constant under an infinitesimal change dm, dμ 2 . Instead, one can measure the change in B 4 , ∆B 4 , under a small variation ∆m or ∆μ 2 , thus estimating the finite differences ∆B 4 /∆m and ∆B 4 /∆μ 2 . For sufficiently small variations, these discrete differences will approach the derivatives ∂ B 4 /∂m and ∂ B 4 /∂μ 2 . Finally, once the gradient of B 4 is known, the direction which keeps B 4 constant is given by The second equation in (2.5) is not needed here. The measurement of B 4 at the shifted couplings (m c + ∆m, µ = 0) or (m c , ∆μ) is performed by the usual analysis, which automatically tunes κ to satisfy (δ M) 3 = 0. To measure B 4 at the shifted couplings, one could perform new Monte Carlo simulations. But this is not necessary. Because the shift in the couplings is very small, it is adequate and safe to use the original Monte Carlo ensemble and simply reweight the results in the standard way [13]. Moreover, to avoid complex weights one reweights to imaginaryμ, which simply introduces a sign flip in eq.(2.9). The reweighting factors remain real positive, and close to 1.
With reweighting, the fluctuations in the original and the reweighted ensembles are strongly correlated. This can be turned into a virtue, as these correlated fluctuations drop out of our observable, which is the change in B 4 , rather than B 4 itself. This procedure is illustrated Fig. 5. The top row shows the change in B 4 under a change in h m = exp(−m) (left), and under a change in imaginaryμ (right). One observes the expected dependence, respectively linear and quadratic. Note that the changes in B 4 are small, O(10 −2 ), and measured very accurately -much more accurately than B 4 itself. The bottom row of Fig. 5 shows the change in B 4 , divided by the change in h m (left) or by (μ) 2 (right). If the Taylor expansion could be truncated to leading order, the data would be constant. Instead, one sees the small influence of the next Taylor order. A fit to (constant + linear) gives the desired partial derivatives, marked by a black circle, which can be substituted in eq.(2.9). The resulting slope is 0.589 (7), again consistent with the other two methods.
Note that the statistical errors on the various points Fig. 5 are extremely correlated with each other, so that a jackknife bin analysis is required to obtain reliable errors on the final slope. Note also that there is a broad optimum for the shift in the couplings: too small a shift produces too small a change in B 4 , and the estimates of the derivatives approach 0/0; too large a shift introduces a systematic error from higher-order Taylor terms and a potential overlap problem.
The errors from methods B and C are similar. This is normal, since these two methods make use of the same Monte Carlo data and extract the same information. So the preference for method C in this case comes from its simplicity: there is no need to measure the relatively complicated observables eqs. (2.8). This will become a more serious issue in the case of QCD.

N f = 3 QCD
In [4], for N f = 3 QCD, with standard staggered fermions on an N t = 4 lattice, we determined the critical quark mass at µ = 0: am c 0 = 0.0263 (3), and proceeded to Taylor expand the pseudocritical coupling β c and the Binder cumulant ofψψ: from which one can extract the variation of the pseudo-critical temperature and of the critical quark mass with the chemical potential.

Recall: method A
In [4], we performed independent simulations at different values of the quark mass and imaginary chemical potential µ = iµ I , which we then fitted with the Taylor expansions above.
For the pseudo-critical temperature, a leading order fit was satisfactory, yielding for the µ 2 dependence c 10 = 0.781 (7). At subleading order, there was no evidence for a cross-term (am − am c 0 )(aµ) 2 , and the (aµ) 4 term was barely statistically significant. Including it in the fit yielded c 10 = 0.759(22) (see Table 1 of [4]). Both fits are shown in Fig. 6 (left), by the narrow blue and broader green error band, respectively.
For the curvature of the critical surface, the coefficient c ′ 1 eq.(3.2) encodes the relevant information: to keep B 4 constant, one must satisfy d(am) d(aµ) 2 = c ′ 1 . Again, no cross-term C was visible, and fits including µ 2 only, or µ 2 and µ 4 terms were performed, resulting in the narrow blue and broader green error bands Fig. 6 (right), corresponding to Table 2, lines 2 and 3 of [4]. They suggest a negative value for c ′ 1 , but do hardly more than that.

Method B
The Taylor coefficients in eqs. (3.1,3.2) can be expressed as expectation values to be measured at µ = 0. We wrote down these expressions, analogous to eqs.(2.8). But we did not implement a program to measure these operators, for several reasons. First, the programming effort is nontrivial. For instance, the trace of the 5 th inverse power of the Dirac operator must be evaluated to obtain the derivative of (δψψ) 4 with respect to the quark mass. Moreover, important cancellations will take place among the various contributions, leading to further difficulties with optimizing the number of noise vectors to be used as stochastic estimators of the various traces. Finally, we realized, from the Potts test case, that method C makes use of the same information contained in the µ = 0 Monte Carlo ensemble, and gives the same output as method B with less effort.

Method C
As in the Potts case, one can shift very slightly the quark mass and the chemical potential, and reweight the µ = 0 Monte Carlo ensembles to these shifted couplings. The effect of a shift in the quark mass was already measured in [4] with sufficient accuracy, so we were interested in the effect of a small chemical potential, taken as imaginary to preserve positivity of the weights. The difference with the Potts case is that the reweighting factors, for each configuration {U } were not computed exactly, but only estimated as where η is a Gaussian random vector. Since N f = 3 in our case, the fractional powers of the Dirac operator were approximated to high precision by a ratio of polynomials. To reduce the variance, we actually formed uncorrelated estimators for ρ(U, µ 2 , µ 1 ) and multiplied them together. Further, ρ(U, µ 2 , 0), µ 2 > µ 1 was constructed as ρ(U, µ 2 , µ 1 ) × ρ(U, µ 1 , 0). In this way, highly correlated reweighting factors were obtained for 6 values of µ = iµ I , with aµ I ranging from 0.01 to 0.1. Crucially, the noise in the reweighting factors does not prevent the standard application of reweighting, because these factors enter linearly in the numerator and denominator of the reweighted expectation value W = ∑ i ρ i W i ∑ i ρ i . Difficulties arise only when one has to form W k , because the same stochastic estimator for configuration i is used k times and the statistical error, when taken to an even power, introduces a bias. This is a relevant concern, since we want to measure the 4 th cumulant ofψψ. However, this bias is of order O(1/N k−1 ) and disappears as the Monte Carlo sample size N grows. To test for such bias, we magnify it by multiplying every estimated weight by a random number drawn uniformly in ]0, 2[. The results of the analysis, for the derivative c 10 of the pseudo-critical coupling eq.(3.1) and for the coefficient c ′ 1 in the Binder cumulant expansion eq.(3.2), are presented Fig. 7. They are to be compared with Fig. 6, where the original weights were used. No bias is visible in Fig. 7, leading us to conclude that no bias is present in Fig. 6 either. Note that our sample size is very large: we analyzed about 5 million 8 3 × 4 configurations at 21 β -values. This large statistics was made possible by using the Grid at CERN. The actual running time was less than two weeks.

Comparing methods A and C
The final results from method C (Fig. 6) are c 10 = 0.746 (1) or for the pseudo-critical temperature, using the two-loop β -function to convert to physical units, and (15) or m c (µ) m c (0) = 1 − 3.3 (5) µ πT 2 + · · · (3.6)  for the curvature of the critical surface, where the value of b 10 eq.(3.2) was taken from [4]: b 10 ≈ 13.6. A subleading term would be visible as a slope in the fits Fig. 6. Indeed, a small effect is visible on the right, corresponding to c ′ 2 = −2.5 ± 1.2. These results are consistent with those of method A, provided subleading terms are included in fitting the latter (green error bands in Fig. 6), even if they are statistically almost unconstrained. This illustrates the difficulty of estimating the systematic error of truncating the Taylor expansion used in fitting. This difficulty is eliminated in our new method, which in addition is about 100 times more efficient.
However, method A gathers statistics over a broader range of chemical potentials, and with updated statistics exceeding 15 million configurations, now allows us to clearly identify subleading Taylor terms. They are visible from the curvature of the fits of the Binder cumulant as a function of the quark mass Fig. 8 (left). The complete ansatz eq.(3.2) was used, now yielding c ′ 1 = −0.074(28) and c ′ 2 = −1.0(5) with a χ 2 of 23 for 21 d.o.f. Note the consistency of methods A and C, also for the subleading order whose sign contributes to further shrinking the region of first-order transitions. Fig. 8 (right) shows the same results, after subtraction of the fitted mass dependence, as a function of (aµ I ) 2 . The fit is shown by the lower parabola. Now, the results of Fig. 6 (right) (leading and subleading terms) are shown in the same figure as the upper parabola. The agreement between the two independent methods is remarkable, given that method C only probes the region (aµ I ) 2 ≤ 0.01 where agreement is near-perfect.

Towards the continuum limit
Of course, cutoff errors on our N t = 4 (a ∼ 0.3fm) lattice can be large, and it is essential to perform a continuum extrapolation. To this end, we are pursuing our project on N t = 6 lattices, and present some preliminary results Fig. 9.
The left figure illustrates cutoff effects on the critical bare quark mass, m c 0 , corresponding to a second-order transition at µ = 0 in N f = 3 QCD. One can see that the quark mass, expressed in units of the temperature, must be reduced by a factor ∼ 5 on N t = 6 lattices. A similarly large effect is present in the resulting pion mass, m c π , measured at zero temperature for quark mass  m c 0 . The ratio m c π /T c decreases from 1.680(4) (N t = 4) to 0.954 (12) (N t = 6), so that a naive a 2 extrapolation would give ∼ 0.4 in the continuum! This very large cutoff effect is consistent with earlier indications [7,14] and with a new study [6], all suggesting that the transition becomes much weaker in the continuum limit. Note that the cutoff effect on the hadron spectrum is comparatively mild, so that the net effect of a finer lattice is to dramatically push the critical surface Fig. 2 toward the origin, while leaving the physical point untouched. Thus the gap between the critical surface and the physical point widens, pushing the critical point in Fig. 2 (left) to larger values of µ E .
Our second preliminary result, Fig. 9 (right), shows the curvature of the pseudo-critical coupling dβ c /d(aµ) 2 | µ=0 for the critical quark mass m c 0 . The error band corresponds to the N t = 4 study. The trend is for dβ c /d(aµ) 2 to be smaller for N t = 6, while if we use the two-loop β -function to convert to physical units, one should observe dβ c /d(aµ) 2 ∝ N 2 t . Instead of increasing by (6/4) 2 , our measured value seems to decrease. Now, for N t = 4 already, the estimated curvature of the pseudo-critical line T c (m c 0 , µ) was about 3 times less than that of the experimental freeze-out curve [15]. These two curves appear to become more clearly separated as a → 0, which also reduces m c 0 .

Discussion
Our findings including next-to-leading terms on N t = 4 predict the "exotic" scenario of Fig. 2  (right): the region of first-order transitions shrinks as the chemical potential is turned on, so that the chiral critical surface does not intersect the physical line, and there is no chiral critical point in QCD. This statement, which goes against conventional wisdom, gets qualified by a number of systematic errors. Like [4][5][6][7][8]10,11,14], we use staggered fermions with the rooting trick, which is potentially unsafe for very light quark masses. Next, the curvature of the critical surface varies with the cutoff, and presently the rate of this change is unknown. It also changes from the N f = 3 case presented here to the N f = 2 + 1 theory, although we found in [4] that for N t = 4 the sign of the curvature remains unchanged. Finally, once |µ| ∼ T c higher order terms may become relevant.
Despite these caveats, we believe the qualitative picture to be robust. If the continuum critical quark mass can be Taylor expanded as per eq.(1.1) with coefficients O(1), the critical surface in Fig. 2 rises "almost vertically" no matter the sign of its curvature, and a critical point at small |µ/T | implies a fine-tuning of the physical quark masses, so as to be very close to the critical line at µ = 0. Such a fine-tuning seems unnatural, and indeed the µ = 0 critical line seems to recede considerably in the continuum limit, now requiring a large curvature of the opposite sign to what we observe in order to accomodate a critical point at |µ/T | < 1.
Finally, the object of our study is the chiral critical surface. Our findings do not exclude additional critical structure due to non-chiral physics causing the phase transitions Fig. 3, bottom right.