Towards a reliable lower bound on the location of the critical endpoint

We perform the first direct determination of the position of the leading singularity of the pressure in the complex chemical potential $\mu_B$ plane in lattice QCD using numerical simulations with 2-stout improved rooted staggered fermions. This provides a direct determination of the radius of convergence of the Taylor expansion of the pressure that does not rely on a finite-order truncation of the expansion. The analyticity issues in the complex $\mu_B$ plane of the grand canonical partition function of QCD with rooted staggered fermions are solved with a careful redefinition of the fermion determinant that makes it a polynomial in the fugacity on any finite lattice, without changing the continuum limit of the observables. By performing a finite volume scaling study at a single coarse lattice spacing, we show that the limiting singularity is not on the real line in the thermodynamic limit, thus showing that the radius of convergence of the Taylor expansion gives a lower bound on the location of a possible phase transition. In the vicinity of the crossover temperature at zero chemical potential, the radius of convergence turns out to be $\mu_B/T \approx 2$ and roughly temperature independent.


Introduction
The theoretical study of the phase diagram of QCD in the temperature (T )-baryon chemical potential (µ B = 3µ q = 3µ) plane is of considerable interest for the physics of high energy nuclear collisions and for the understanding of the early stages of the universe. It is by now established through first-principles lattice QCD calculations that at µ B = 0 there is an analytic crossover [1] at a temperature [2] of T c ≈ 150 − 160 MeV. It is further conjectured [3] that in the (T, µ B ) plane there is a line of crossovers, departing from (T c , 0), that eventually turns into a line of first-order phase transitions. The point (T CEP , µ CEP ) separating the two lines is known as the critical endpoint (CEP), and the transition is expected to be of second order there. Full confirmation of this picture is hindered by the notorious sign problem, which prevents direct simulation at µ B > 0, but several extrapolation techniques have been developed to study QCD at finite small densities. These include Taylor expansion at µ B = 0 [4,5], analytic continuation from imaginary chemical potential [6,7,8], and reweighting [9,10]. The basic idea of these methods is to reconstruct the behavior of the theory at real µ B > 0 by extrapolating from zero or purely imaginary µ B , where the sign problem is absent. In the context of these extrapolation methods, one of the most sought-after quantities in the finite temperature lattice QCD community is the radius of convergence of the Taylor expansion of the pressure around µ B = 0 [5,11,12]. The reason for this is twofold. First, the radius of convergence gives information about the location of the supposed CEP. More precisely, the radius of convergence gives a rigorous lower bound on the location of possible phase transitions. Moreover, if the leading singularity in the thermodynamic limit is on the real axis then it corresponds to a genuine phase transition. Second, one of the most important inputs of lattice QCD to the phenomenology of heavy ion collisions is the equation of state, including its extrapolation to finite µ B . Knowledge of the radius of convergence would tell us how far such an extrapolation is trustworthy. Even if the lower bound on the CEP location happens to be not very stringent, this is still very useful information, since viscous hydrodynamics simulations of heavy ion collisions at intermediate energies already explore a very wide range of temperatures and baryo-chemical potentials. Currently available lattice estimates of the radius of convergence are obtained using the first few Taylor coefficients, and with their statistical significance strongly deteriorating with order, it would seem rather unlikely that the convergence of the estimators of the radius of convergence can be demonstrated in the near future. This contribution introduces a novel method, with which this problem may be circumvented. More details can be found in Refs. [13] and [14].

Lee-Yang zeros and the exact asymptotics of the Taylor expansion in a finite volume
In a typical statistical mechanical system the grand canonical partition function is, up to a non-vanishing factor, a polynomial of the fugacity e µ/T (the Lee-Yang polynomial). In such a situation, the zeros of this polynomial, the so-called Lee-Yang zeros determine the analytic properties of the pressure ∝ log Z. Namely at each Lee-Yang zero, the pressure develops a logarithmic branch point singularity. The Lee-Yang zeros have some simple symmetry properties, namely, they have a µ → µ * symmetry due to the fact that the coefficients of the Lee-Yang polynomial are the canonical partition functions, i.e. they are real numbers Z n > 0. Second, in QCD-like theories, due to CP symmetry they also have a µ → −µ symmetry. This means in particular the in the Taylor expansion log Z = n c n µ 2n the high order terms are dominated by 4 logarithmic singularities (Lee-Yang zeros) at the same distance from µ = 0. This is enough information to calculate the exact asymptotic behavior of the Taylor expansion in any finite volume as [13]: as k → ∞, where r and θ are the polar coordinates of an arbitrary member of the closest lying quartet of Lee-Yang zeros. This formula has some interesting immediate consequences, namely: i) The ratio estimator, widely used in the literature to estimate the radius of convergence can never converge in any finite volume, where the lattice simulations are actually done; ii) Even if a given gauge ensemble is large enough to determine r and θ accurately, the errorbars of the high order coefficients will always blow up (as can be checked by applying linear error propagation applied to the above formula) iii) The fluctuations of asymptotically high order are strongly correlated on any given gauge ensemble, since they are all given by r and θ. The last point in turn implies that if the correlations between the c n are kept, even if they have > 100% errorbars, the leading Lee-Yang zero position may still be recovered. This was explicitly demonstrated on a small toy lattice in [13]. How the leading Lee-Yang zero can be recovered from the Taylor coefficients is also explained in [13], where improved estimators of the radius of convergence are constructed, that take into account the exact asymptotic form of the Taylor coefficients discussed above.
The main take-away of this discussion is that what allows one to reconstruct the leading Lee-Yang zero in the first place is the strong correlations between the large order Taylor coefficients, which make possible a very striking cancellation of statistical errors in certain combinations of them, corresponding to the position of the Lee-Yang zeros. This means that in order to eventually calculate the radius of convergence in continuum QCD, one is probably better served by calculating the radius of convergence on a fixed lattice spacing and volume first, and only later tackle the continuum and infinite volume limits at the level of the radius of convergence itself. Unfortunately, as we will soon discuss, this is not possible with the numerically cheapest formulation of lattice QCD, rooted staggered fermions, as the partition function with this discretization has some spurious singularities, discussed below.

Spurious non-analyticities with rooted staggered fermions and how to remove them
The staggered discretization corresponds to 4 degenerate flavors of fermions in the continuum limit. In order to describe N f = 1 or 2, the so-called rooting trick is used, where one writes the partition function as: , where det M is the determinant of the staggered Dirac matrix. When the baryon chemical potential µ B > 0 the determinant is a complex number. This is the famous sign problem of lattice QCD. In order to define the N f = 2 determinant, one then has to choose a branch of the square root function. The standard methods of analytic continuation and Taylor expansion use here a branch of the square root that continuously connects to the positive square root at µ B = 0. While this choice is expected to be perfectly fine near the continuum limit, at a finite lattice spacing it has the non desirable property that the partition function develops a square root type branch point whenever the determinant has a zero on a given gauge configuration. A natural consequence of this is that the partition function is non-analytic everywhere on the support of the probability density of the zeros of the staggered determinant. In such a case the radius of convergence of the pressure is bounded from above by the point of this support closest to µ = 0. This is a cut-off effect that at any finite lattice spacing potentially reduces the radius of convergence to a smaller, unphysical value, as the radius of convergence is no longer given by a partition function zero, but rather the closest zero of the staggered determinant on the entire ensemble. With the standard methods, the obvious solution to this problem is to calculate the continuum limit of the individual c n first, and only later study the high order behavior of the series. Unfortunately this destroys the strong correlations present between the c n on a single ensemble, which allows one to reconstruct the position of the leading Lee-Yang zero in the first place.
We now proceed to describe a possible resolution to the problem discussed above. This involves introducing a new definition of the N f = 2 determinant with staggered fermions, that does not involve square roots. It has been shown in Ref. [9] that the staggered determinant on a fixed gauge configuration can be calculated for any complex value of the quark chemical potential µ as det where V is the spatial volume in lattice units, and N s is the spatial linear size of the lattice in lattice units and the ξ i are eigenvalues of the so-called reduced matrix. The reduced matrix and therefore the eigenvalues ξ i only depend on the gauge configuration and not on the chemical potential. Near the continuum limit, the fact that the staggered operator describes 4 flavors of fermions manifests itself in a near degeneracy of the eigenvalues of the original Dirac operator D as well as those of the reduced matrix [14]. One can then define the two-flavor determinant by grouping closely lying pairs of eigenvalues of the reduced matrix and substituting each pair (ξ p 1 , ξ p 2 ) with a suitably defined "average eigenvalue"ξ p and defining the two-flavor determinant as det M(µ/T ) N f =2 = e −3Vµ/(2T ) 3V p=1 ξ p − e µ/T with the sum now running over the pairs. This way the definition of the 2 flavor determinant does not involve any square roots and therefore no branch point singularities. This then guarantees that the radius of convergence will again be given by the partition function zero closest to µ = 0. With this new definition in hand we can embark on the calculation of the leading Lee-Yang zero. To test the idea numerically we conducted simulations with N t = 4 2stout improved staggered fermions with lattice spatial volumes of 6 3 , 8 3 , 10 3 and 12 3 and calculated the position of the leading Lee-Yang zero for all cases. We also performed a simple infinite volume extrapolation with a linear ansatz in 1/V. The nonzero infinite volume limit of the imaginary part indicates that the radius of convergence of the pressure is not given by a genuine phase transition, but rather by a complex singularity. Our numerical results are summarized briefly in Fig. 1. On the left we show an infinite volume extrapolation of the imaginary part of the leading Lee-Yang zero at a specific temperature, demonstrating that the leading singularity is complex, and therefore that the radius of convergence is not given by a phase transition. On the right, we show the infinite volume extrapolation of the radius of convergence as a function of the temperature. For more details, we refer the reader to Ref. [14]. Extending such a study to finer lattices is challenging. However, our approach has the conceptual advantage over the standard definition of staggered rooting of dealing with an actual Lee-Yang polynomial at a finite lattice spacing, instead of a function that is non-analytic in dense regions of the complex chemical potential plane. This makes the challenge worth the pursuit.