Radius of convergence in lattice QCD at finite $\mu_B$ with rooted staggered fermions

In typical statistical mechanical systems the grand canonical partition function at finite volume is proportional to a polynomial of the fugacity $e^{\mu/T}$. The zero of this Lee-Yang polynomial closest to the origin determines the radius of convergence of the Taylor expansion of the pressure around $\mu=0$. The computationally cheapest formulation of lattice QCD, rooted staggered fermions, with the usual definition of the rooted determinant, does not admit such a Lee-Yang polynomial. We argue that the radius of convergence is then bounded by the spectral gap of the reduced matrix of the unrooted staggered operator. This is a cutoff effect that potentially affects all estimates of the radius of convergence with the standard staggered rooting. We suggest a new definition of the rooted staggered determinant at finite chemical potential that allows for a definition of a Lee-Yang polynomial, and, therefore of the numerical study of Lee-Yang zeros. We also describe an algorithm to determine the Lee-Yang zeros and apply it to configurations generated with the 2-stout improved staggered action at $N_t = 4$. We perform a finite-volume scaling study of the leading Lee-Yang zeros and estimate the radius of convergence of the Taylor expansion extrapolated to an infinite volume. We show that the limiting singularity is not on the real line, thus giving a lower bound on the location of any possible phase transitions at this lattice spacing. 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. Our simulations are performed at strange quark chemical potential $\mu_s=0$, but the method can be straightforwardly extended to strangeness chemical potential $\mu_S=0$ or strangeness neutrality.


I. INTRODUCTION
One of the open problems in the study of QCD at finite temperature and density is determining the phase diagram of the theory in the temperature (T )-baryon chemical potential (µ B = 3µ q ) plane.It is by now established that at µ B = 0 there is an analytic crossover [1,2] at a temperature [3][4][5][6] of T c ≈ 150-160 MeV.It is further conjectured 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 crossovers and first-order transitions is known as the critical end point (CEP), and the transition is expected to be of second order there.
The conjectural phase diagram discussed above is mostly based on effective models of QCD [7].To settle the issue, one needs a first-principles study of QCD at finite temperature and density, which requires nonperturbative tools like the lattice formulation of the theory.Unfortunately, the introduction of a finite chemical potential µ B makes the direct application of traditional importance-sampling methods impossible due to the notorious sign problem.For this reason, lattice QCD could so far give very limited information about the phase diagram away from µ B = 0.
Among these methods, reweighting has the advantage of having no other systematic error besides the overlap problem and in principle would lead to the correct results in the limit of infinite statistics, so that it can be at least used as a "brute-force" approach to the sign problem.On the other hand, the Taylor expansion method is also affected by systematic errors from the truncation of the Taylor series, as well as the existence of a finite radius of convergence, while extrapolation from imaginary chemical potential involves a rather uncontrolled analytic continuation in µ B .
The aim of this paper is to obtain as much information as possible about the analytic structure of the pressure in the complex baryon chemical potential plane.In particular, we will estimate the position of the singularity closest to µ B = 0, which provides both a lower bound on the location of possible phase transitions on the phase diagram as well as the limit of reliability of the equation of state coming from a Taylor expansion, which is an important input for the phenomenology of heavy ion collisions.
In a finite volume the analytic properties of the pressure of a typical statistical mechanical system are governed by the zeros of the partition function in the complex-µ B plane, the so-called Lee-Yang zeros [40].In general, the grand-canonical partition function of a relativistic lattice system at finite μ = µ q T = µ B 3T in a finite spatial volume V is a polynomial in the fugacity z = e μ, which we may call the Lee-Yang polynomial, times a nonvanishing factor e −kV μ for some model-dependent constant k.The Lee-Yang zeros are the singular points of the pressure, p = − T V log Z, as a function of complex fugacity or chemical potential.The accumulation of such zeros near the real µ B axis in the thermodynamic limit V → ∞ signals the presence of a genuine phase transition.In the case of a crossover, no nonanalyticity develops on the real line and the distance of the Lee-Yang zeros from the real axis provides a measure of the strength of the transition.The radius of convergence of the Taylor expansion of the pressure around µ B = 0 is equal to the distance from the origin of the Lee-Yang zero closest to it, which we will refer to as the leading Lee-Yang zero.Depending on the behavior of the leading zero in the infinite-volume limit, the radius of convergence could correspond to the chemical potential at which an actual phase transition takes place (in case that the imaginary part of the leading zero extrapolates to zero) or just give a lower bound on the location of a phase transition (in case that the imaginary part extrapolates to a nonzero value).
Conversely, the position of the closest Lee-Yang zero can be inferred from the high-order behavior of the Taylor coefficients.In fact, at a fixed lattice spacing, as long as one uses a discretization where the partition function is an entire function of µ, knowing the position of the leading Lee-Yang zero is completely equivalent to knowing the asymptotically high-order behavior of the Taylor coefficients of the pressure.This was shown in [41], where explicit formulas are given for the conversion [42].The determination of the leading Lee-Yang zero from the Taylor coefficients however involves an extrapolation to high orders of the Taylor expansion, which is technically challenging.It turns out that a direct determination of the Lee-Yang zeros using reweighting techniques is instead more straightforward.One might wonder how this is possible: if it is so difficult to estimate reliably the high-order Taylor coefficients in order to determine the radius of convergence, determining the latter directly seems hopeless.The answer is that the large (above 100%) errors on the high-order coefficients are strongly correlated, and cancel out in the particular combinations that give the leading Lee-Yang zero, and therefore the radius of convergence.This surprising conclusion was discussed in Ref. [41] and demonstrated explicitly in a numerical study of unrooted staggered fermions on a small lattice.This suggests that it is more efficient to calculate the radius of convergence first at a finite lattice spacing, where the strong correlations between the Taylor coefficients are present, than taking the continuum limit of the coefficients first and calculating the radius of convergence of the continuum expansion, in which case part of the correlations is lost.
A convenient way to do reweighting, which allows for a straightforward determination of the Lee-Yang zeros, is to compute the spectrum of the so-called reduced matrix P [32,35,37] on an ensemble of gauge configurations at µ q = 0, which then allows one to reweight to any finite µ q using the relation det M (μ) = e −kV μ det(P − e μ), where M (μ) is a lattice discretization of the QCD Dirac operator at finite µ q and k is the same model-dependent constant appearing in the partition function.The reduced matrix of Ref. [32] is instrumental to the approach of this paper.In fact, expressed in terms of the reduced matrix, the fermionic determinant is a polynomial in fugacity on each configuration up to a nonvanishing prefactor, which allows for a straightforward reconstruction of the polynomial part of the grand-canonical partition function at finite µ q and the subsequent determination of the Lee-Yang zeros by means of standard numerical techniques.
The discussion above is quite general, and it applies whenever a reduced matrix formulation is available [32,43,44].Unfortunately, this does not include the computationally most convenient formulation of lattice fermions, namely the rooted staggered discretization.Near the continuum the spectrum of the staggered Dirac operator shows quartets of near-degenerate eigenvalues, with relative splittings of order O(a), corresponding to the so-called staggered tastes [45][46][47].Taking roots of the staggered determinant should then fully solve the doubling problem of lattice fermions, by reducing the N f = 4 flavor theory of unrooted staggered fermions down to the desired number of degenerate flavors.In the N f = 2 case, for example, the square root of the determinant is taken.This procedure has become standard, and although there is no rigorous proof that it ultimately provides us with a genuine local continuum quantum field theory, the results obtained are in good agreement with experiments and with lattice results obtained using other fermion discretizations [2,[48][49][50][51].At finite chemical potential one can still recast the determinant of the staggered Dirac operator in terms of a reduced matrix, but the rooted determinant is not a polynomial in fugacity anymore, and so the essentially polynomial character of the grandcanonical partition function is lost.
The lack of a Lee-Yang polynomial is not the only problem afflicting rooted staggered fermions at finite µ B .Unlike at µ B = 0, where one can simply take the real positive root of the real positive determinant, at finite µ B where the determinant is complex there is not such a natural choice, and some other criterion is needed to resolve the intrinsic ambiguity of rooting.The standard choice for both reweighting and the Taylor-expansion method is to take the root that, on a given gauge field configuration and as a function of µ B , continuously connects to the real positive root at µ B = 0.While this choice is perfectly fine in the continuum limit, where the formation of eigenvalue quartets is expected to make the ambiguities related to rooting go away, at any finite spacing it leads to serious analyticity problems.As we will argue the resulting partition function as a function of complex fugacity will be nonanalytic everywhere on the support of the spectrum of the reduced matrix P .The radius of convergence of the Taylor expansion of the pressure will then be given by the spectral gap of P , i.e. the closest distance where the spectral density of P is nonzero, which on a finite ensemble corresponds to the distance of the closest eigenvalue on the entire ensemble.This spectral gap a priori has nothing to do with Lee-Yang zeros or phase transitions.This means the radius of convergence will not be given by a partition function zero, unless it happens to be inside the gap.Even if this lucky coincidence happens, one would not be able to tell from the Taylor expansion.
We will circumvent this problem by defining the "rooted" staggered determinant not by choosing a particular branch of the root function, but by setting it equal to a polynomial in the fugacity that is expected to converge to the same continuum limit as the standard procedures discussed above.Such a polynomial will differ from the standard definitions of the rooted determinant by terms that are nonanalytic in µ B but vanish in the continuum.This way the radius of convergence of the expansion of the pressure will certainly not be related anymore to the spectral gap (which depends on the support of the spectrum on the whole ensemble and is thus determined by some "extreme" configuration) and will be given by an actual Lee-Yang zero (which depends only on averages over the ensemble).
Our construction of a rooted determinant in polynomial form is motivated by the very idea behind rooting of staggered fermions, i.e., the formation of taste quartets in the continuum.In Ref. [52] it was suggested that the optimal way to do rooting is to identify the taste multiplets in the spectrum of the staggered Dirac operator and replace them with their average to define the rooted determinant.This was argued to considerably reduce the finite-spacing effects compared to other procedures.As we will argue, the formation of quartets in the spectrum of the staggered Dirac operator leads to the formation of quartets in the spectrum of the corresponding reduced matrix P .In the spirit of Ref. [52], we can therefore define a rooted determinant in the reduced matrix approach by judiciously grouping the eigenvalues of P and replacing them with a properly defined average.In this way we automatically obtain a definition of the rooted determinant which is a polynomial in fugacity and so analytic in µ B .One can then define a Lee-Yang polynomial and carry out the study of its zeros by standard methods.The purpose of this paper is thus twofold: i) show how to conveniently group eigenvalues and how to average them in order to provide a more convenient definition of the theory at finite µ B with two light flavours of rooted staggered fermions; ii) obtain a direct and complete determination of the Lee-Yang zeros of the so-defined partition function.
The plan of the paper is the following.In Sec.II we review the reduced-matrix approach, focusing in particular on its problems of analyticity at finite µ B when one takes roots of the determinant, and how to solve these problems by grouping eigenvalues.We also briefly review Lee-Yang zeros and discuss how to numerically compute them.In Sec.III we perform a numerical study with 2stout improved N t = 4 staggered lattices.Finally, in Sec.IV we draw our conclusions and discuss future prospects.

II. REDUCED MATRIX FOR STAGGERED FERMIONS
A. Generalities The introduction of a finite quark chemical potential µ = µ q in the staggered Dirac operator is usually done by coupling it to the temporal links as follows: where a is the lattice spacing, (T α ) xy = δ x+ α,y are the translation operators and η α are the usual staggered phases.It is easy to show that det M (μ) ≡ det [2(D stag (aµ) + am)] depends only on μ = µ T .It has been shown in Ref. [32] that [53] det where V = N 3 s is the spatial volume and N s is the spatial linear size of the lattice in lattice units (which must be an even number), and ξ i are eigenvalues of the reduced matrix P .In the temporal gauge [U 4 (t, ⃗ x) = 1 for 0 ≤ t < N t − 1], this reads i.e., η 4 B i is the sum of the spatial derivatives and mass parts of the staggered matrix on the ith time-slice, and L is the block-diagonal matrix of temporal links on the last time slice (i.e., the untraced Polyakov loops).Since P is µ independent, knowledge of the ξ i for a given gauge configuration allows to compute the corresponding unrooted quark determinant for arbitrary µ.From a Monte Carlo simulation at µ = 0 one can then obtain the grandcanonical partition function at any µ via where the subscript 0 indicates that the expectation value is computed at µ = 0.The quantities P ξ and P, defined in Eq. ( 4), are polynomials of degree 6V of the fugacity z = e μ.The coefficients of the Lee-Yang polynomial P are the average of those of P ξ , and coincide with the (normalized) canonical partition functions, and as such they are positive quantities.Notice that Roberge-Weiss symmetry [54] imposes that only coefficients of order 3n can be nonzero.The canonical partition functions are usually obtained as the coefficients of a Fourier expansion of the grand-canonical partition function at imaginary chemical potential [32,[55][56][57][58][59][60][61][62][63][64][65][66][67][68][69].Our direct determination from the eigenvalues of P is free from the systematic uncertainty associated with the extraction of Fourier coefficients from a discrete set of imaginary chemical potentials.
The matrix P has a few nice properties, which we list here: 1. det P = 1.

The product
′ i ξ i of eigenvalues inside the unit circle is real positive.
A proof of these properties can be found in the Appendix.
The Lee-Yang zeros are the roots of the polynomial P(e μ).Due to the symmetries of the partition function, the Lee-Yang zeros will be symmetric under the reflection µ → −µ due to CP symmetry and under the reflection µ → µ * due to the fact that Z(μ) is real analytic.Furthermore, Roberge-Weiss symmetry implies that we can restrict ourselves to the strip Im μ ∈ [− π 3 , π 3 ], since the zeros are then repeated with a period of 2π 3 in the imaginary μ direction.There can be no zeros on the real axis due to positivity of Z(µ) and on the imaginary axis since the determinant is real positive there.Note that these properties of Z originate in properties 1-3) of P .

B. Problems with rooting
When using rooted staggered fermions, one must replace the ratio of determinants in Eq. ( 4) with its appropriate root.The intrinsic ambiguity associated with rooting is easily solved at µ = 0, where it is natural to take the real positive root of the real positive ratio of determinants.Such a natural possibility is not available at nonzero µ where the fermionic determinant is generally complex.The choice in previous reweighting studies [35] was to use the root that on a fixed gauge configuration continuously connects to the real positive root as a function of µ.We focus here on the N f = 2 case.In terms of the eigenvalues of the reduced matrix one has with the branch cut of the square root chosen to lie along the negative real axis.Note that this equation actually provides a definition of the square root on the lefthand side.On a single configuration, this quantity has 6V branch cuts parallel to the real axis.Restricting to real chemical potentials one always moves parallel to the branch cuts and never crosses them, making the function continuous at real µ.The Taylor-expansion method in turn is formulated by taking the square root of the entire ratio of determinants det M (μ) det M (0) , and expanding around µ = 0, where the argument of the square-root function is one.This leads to an expansion ) On a single configuration the radius of convergence of this Taylor expansion is given by the point where the unrooted determinant first becomes zero, i.e., by the closest eigenvalue of the reduced matrix.Notice that one would get the same formula if one expanded log det M (μ) det M (0) rew around zero, meaning that on a single configuration the two methods define the same rooted determinant within the radius of convergence.
The discussion above extends naturally to the Taylor expansion of log Z around µ = 0, but now the radius of convergence is determined by the eigenvalue of the reduced matrix whose logarithm is closest to 0 on the entire ensemble.The issue is easily understood with a very simple example, in which we have only two configurations and two eigenvalues on each configuration.In this case the toy partition function is defined as where a 1,2 , b 1,2 are complex numbers, the square roots are all defined with the branch cut on the negative real axis and ζ is the complex fugacity parameter ζ = e μ − 1 = μ + O(μ 2 ).No matter how close a 1 , a 2 and b 1 , b 2 are, the radius of convergence of the Taylor expansion of the "pressure" log Z toy around ζ = 0 cannot be larger than min( a 1 , b 1 , a 2 , b 2 ), leading to our statement that the radius of convergence is bounded from above by the spectral gap.This fact remains true for any a 1 ≠ a 2 , b 1 ≠ b 2 , and only changes in the case of exact degeneracy, which in this toy example mimics the "continuum limit".The partition function Z toy can still have a zero, that in general can be either closer or farther away from ζ = 0 than the closest square-root branch point.In the limit of exact degeneracy, a 1 = a 2 ∶= a and b 1 = b 2 ∶= b the radius of convergence of log Z toy is given by the partition function zero at (a + b) 2.
The situation is not expected to improve if one has a larger number of branch points in the square root and if instead of the sum of two terms one has an average over the position of the branch points with some probability distribution, which is the general case to which the rooted fermion determinant in the reduced-matrix approach belongs.In this case one expects that after averaging/integration over gauge fields, the partition function will be nonanalytic over the support of the spectrum ×4 and β = 3.35 superimposed into a single plot.A small gap is also present near 0, but it is not visible on the plot. of P .In fact, the only way to avoid this conclusion is the existence of some cancellation mechanism between the branch points, which seems unlikely.This is problematic not only for the reduced-matrix approach, but for any approach involving the rooting of the fermionic determinant, like the Taylor-expansion method or analytic continuation from imaginary chemical potential.In fact, this argument suggests that since the eigenvalues in the full gauge ensemble are expected to fill densely some region inside the unit circle (see Fig. 1), at any finite spacing there will be an analytically inaccessible region, whose boundary is determined by the extreme edges of the spectrum.Numerical results seem to indicate the existence of a gap around the unit circle and thus of a finite domain of analyticity [70][71][72].
From a practical point of view, with the usual definitions of the rooted determinant, the analyticity domain of the partition function in fugacity on any finite ensemble of configurations will be determined by the position of the eigenvalue of P closest to 1 in the ensemble.
A consequence of this discussion is that with the usual definition of the rooted staggered determinant at finite µ, the radius of convergence at a finite lattice spacing may not be related to the physics of phase transitions.This also means that the radius of convergence of the continuum Taylor series may not be equal to the continuum limit of the finite-spacing radius of convergence.
Before discussing the solution to these problems in the next subsection, it is worth remarking that the possible easing of problems near the continuum limit is based on the expected formation of quartets of eigenvalues.So far, the formation of quartets near the continuum limit has been investigated only for the usual staggered operator [45][46][47].It is then worth checking that they do indeed form also in the spectrum of the reduced matrix.A first check is provided by solving analytically the eigenvalue problem in the free case, in which quartets of eigenval- ues are explicitly shown to appear (see the Appendix).Since in the continuum limit the relevant configurations fluctuate around the free one, quartets of eigenvalues are expected to show up for sufficiently small lattice spacing.We also conducted a second, numerical check.While a direct study is currently out of reach, here we mimic the approach to the free case by applying a large number of stout smearing steps to a fixed gauge configuration on a small 6 3 × 4 lattice.After 50 steps of stout smearing [73] with ρ = 0.05 doublets appear (see Fig. 2), while a much larger number of smearing steps (∼ 150) is needed for the appearance of quartets.This behavior matches what has been observed with the usual staggered operator [47].Ultimately, one will have to go to fine lattices at a fixed temperature, and look at the spectrum of P to make sure there really is a quartet structure.

C. Geometric matching
Besides the analyticity problems discussed in the previous subsection, rooting of the fermion determinant in the reduced-matrix approach is expected to introduce systematic finite-spacing effects analogous to those discussed in Ref. [52] for the usual fermionic matrix.This is because of the existence of configurations where the taste multiplets are cut through by the branch cut on the negative real axis.While a sensible rooting procedure should correspond in practice to replacing a multiplet with a single effective eigenvalue at the multiplet position, taking roots eigenvalue by eigenvalue as in Eq. ( 5) can lead to the effective eigenvalue being far away, having the opposite sign.
The solution to both kind of problems is to identify taste quartets of eigenvalues and replace them by a properly defined average.In principle this could be done also on not so fine lattices by, e.g., tracking the eigenmodes as the gauge configuration undergoes a number of smearing steps, seeing which ones end up forming quartets, and grouping them accordingly.This procedure is computationally very expensive, and so one would rather opt for the next best thing: find the group of closest eigenvalues and treat them as if they were taste quartets.This will eventually become equivalent to the optimal procedure in the continuum limit, since the near-continuum taste quartets are well separated from each other and will automatically be identified as the groups of closest eigenvalues.After identifying the quartets, they are replaced by the fourth power of an appropriate average, and the rooted determinant is obtained by retaining a single power.
Before detailing the procedure for actual lattice QCD let us come back briefly to our toy example of the previous subsection and see how this approach improves the analyticity of the partition function.Using the geometric mean as the average within each pair, we define the rooted toy partition function as The partition function is now a polynomial, and the radius of convergence of log Z toy,poly around 0 is determined by its Lee-Yang zero , which continuously tends to a+b 2 as the splitting is diminished.Extending the toy example to the case of the lattice QCD partition function, estimated using N conf gauge field configurations and computing the 6V eigenvalues of the reduced matrix, with the usual rooting procedure we expect to find 6V N conf branch points of square-root type, and no easy way to count the zeros of the partition function.This also means the presence of 6V N conf square-root type branch points in log Z, besides the logarithmic ones originating from the zeros.Using a matching procedure to replace pairs of eigenvalues with their geometric mean one finds instead no branch points in Z, since the resulting partition function is a polynomial of order 3V , thus with exactly 3V zeros.This also means that log Z has only logarithmic singularities, and that the radius of convergence of its expansion around µ = 0 is determined by the closest one of them.
We now give details on how the matching procedure is actually implemented in practice.The guidelines are the following: • The correct quartets must be automatically selected in the continuum.
• The new eigenvalues should retain the properties 1-3 of the unrooted reduced matrix discussed at the end of Sec.II A, in order to retain properties of the partition function Z itself.
In this way we expect that in the continuum limit one finds a well-defined continuum theory of a single fermion flavor with the desired properties.Since we are interested in the case N f = 2 of two light flavors, we can simplify the procedure and content ourselves with identifying doublets of eigenvalues instead of quartets.To make our proposal concrete, we still have to specify how to identify doublets and how to average them.The identification of doublets is achieved by minimizing the total sum of the distances within each pair, which we dub "geometric matching".As already remarked above, this will provide the correct identification in the continuum limit, where nearly degenerate taste quartets are well separated from each other (except possibly for UV related modes that do not affect the longdistance physics).The relevant optimization problem can be solved in polynomial time by means of the so-called Blossom algorithm [74].A publicly available implementation of this algorithm exists [75], which we used in our analysis.
Once doublets (ξ 1 , ξ 2 ) of eigenvalues of P are identified, their product is replaced by their geometric mean, leading to the desired rooted determinant.This procedure is illustrated in Fig. 3.More precisely, we replace ξ 1 ξ 2 → ξ obtained as follows: • We first compute ξ′ = √ ξ 1 ξ 2 , choosing among the two roots the one which lies closer to ξ 1 and ξ 2 .
• We then compute the product of the 3V 2 such defined ξ′ inside the unit circle and determine its phase e iδ .
• we finally multiply each of the ξ′ , both inside and outside the unit circle, by the same phase e −i 2δ 3V , i.e., The choice of the geometric mean automatically leads to satisfying properties 1 and 2 of our list.At this stage property 3 will in general not be satisfied, but it will after our global phase correction.The phase of the averaged doublet is actually not constrained by properties 1 and 2 if we choose it in the same way for (ξ 1 , ξ 2 ) and its symmetric partner . For sufficiently large volumes the phase correction is tiny, and it also disappears in the continuum limit when the doublets become exactly degenerate.
Having obtained the averaged eigenvalues { ξi } i=1,...,3V , we define the rooted determinant as where "P" stands for "paired".Equation ( 4) for the rooted determinant becomes with P a polynomial in the fugacity.We can now solve for the roots of this polynomial to determine the Lee-Yang zeros of the partition function.
Finally, we note that both the usual definition and our new definition of the rooted determinant are nonanalytic in the gauge fields at finite fixed µ.The source of the nonanalyticity is however slightly different.The eigenvalues of the reduced matrix are of course analytic functions of the link variables.When rooting eigenvalue by eigenvalue as in Eq. ( 5), the nonanalyticity comes from the eigenvalues of the reduced matrix crossing the branch cut of the square root on the negative real axis.In the case of geometric matching, nonanalyticity comes from the minimization procedure, with eigenvalues sometimes changing pairs as the gauge fields vary.A further source of nonanalyticity is the small phase correction to ensure property 3.The crucial difference is that in the case of geometric matching we can be sure that after integration over the gauge fields the resulting partition function -being a polynomial -is analytic in µ on the whole complex plane.This very nice feature of our approach is a consequence of the simple fact that the sum of polynomials is again a polynomial.With our approach, one can therefore take the continuum limit of the radius of convergence directly, without encountering analyticity issues, which was one of our goals stated in the Introduction.
Finally, we note that the method can be straightforwardly generalized to the N f = 1 case, simply by performing the geometric matching of pairs twice or by defining an objective function which searches for quartets directly instead of pairs.Close to the continuum limit, this is expected to lead to a correct identification of the taste quartets.

III. NUMERICAL RESULTS
We have performed numerical simulations at µ = 0 using a tree-level Symanzik improved gauge action for the gauge fields and 2 + 1 flavors of rooted 2-stout improved staggered fermions [3] at physical quark masses.We used lattices of temporal size N t = 4 and spatial size N s = {6, 8, 10, 12} at β = {3.32,3.33, 3.34, 3.35}, corresponding to temperatures below and up to T c .A chemical potential is then introduced only for the light quarks, with µ u = µ d = µ B 3, µ s = 0.For each simulation point we gather on the order of 20000 configurations, separated by 10 HMC trajectories each.For each gauge configuration we performed a full diagonalization of the reduced matrix P , using the publicly available MAGMA linear algebra library for GPUs [76].This is the most computationally intensive step of the analysis, as it roughly scales with V 3 = N 9 s .We then followed with the geometric matching of the eigenvalues to compute the partition function at finite µ via the reweighting formula Eq. (11).
In order to determine Z(µ) Z(0), for each configuration we calculate the coefficients of the polynomial P ξ (e μ) from the geometrically matched eigenvalues ξi using arbitrary precision arithmetic.Using Roberge-Weiss symmetry we can set to zero all coefficients of the polynomial not of order 3n.Using charge conjugation symmetry instead we can set all coefficients to be real.After averaging the coefficients over gauge configurations, we determine all the roots of P(e μ) using the Aberth method [77], which we also implemented in arbitrary precision.In Fig. 4 we show all the Lee-Yang zeros in the strip Im µ ∈ [− 2π 3 , 2π 3 ] of the complex µ plane for a 12 3 ×4 system at β = 3.35.(Due to rooting there are now half as many zeros in this strip than in the unrooted case.)Error bars are not shown in Fig. 4, and while they are typically quite large, due to the sign problem, the ones closest to µ = 0 zero have reasonably small errors.See Figs.7,8,9 and 11 (shown later).
As a first check of our method, we have studied the correlation between the phase θ 1 of the rooted determinant obtained via the geometric matching procedure used in this paper and the phase θ 2 of the rooted determinant obtained by rooting each term of the product separately, Eq. ( 5).The two methods are expected to give µ q a=0.20 µ q a=0.15 µ q a=0.10 FIG. 5. Correlation of the phase of the determinant as obtained via geometric matching or by rooting separately the contribution of each eigenvalue, i.e.Eq. ( 5), on a 12 3 × 4 lattice at β = 3.34.Here θ1 is the phase defined by the standard rooting procedure, while θ2 is the phase defined by our novel definition of the N f = 2 determinant.Imµ q /T Reµ q /T full truncated FIG. 6.Comparison between Lee-Yang zeros obtained with the full and the truncated polynomial, as described in the text.Here the lattice size is 12 3  × 4 in lattice units and β = 3.35.The range of the Reµ axis is chosen such that all roots of the truncated polynomial can be seen.
the same continuum limit, but on a rather coarse lattice they might be very different, leading to huge systematic uncertainties.In particular, a frequent relative change of sign would make the root of the fermionic determinant a particularly ill-behaved quantity.In Fig. 5 we plot the correlation between cos θ 1 and cos θ 2 for three values of the chemical potential on a 12 3 × 4 lattice at β = 3.34.The two quantities are nicely positively correlated, especially at small real µ, indicating that the two methods will give similar results on the real axis, even at a finite lattice spacing.
A possible systematic effect to take into account comes from the relatively poor determination of the high-order coefficients of P. While all the exact coefficients of the polynomial must be positive, they can turn out to be negative on a finite sample due to limited statistics and the sign problem.In such cases, while the average is negative, the statistical error on the coefficient is above 100%, making them consistent with zero.We have then checked the zeros in Fig. 4 against those obtained truncating the Lee-Yang (LY) polynomial, removing those terms for which we obtained a numerical estimate of the coefficient that is compatible with zero.The comparison is shown in Fig. 6: the physically relevant LY zeros near µ = 0 are unaffected by the truncation.This is true also for the zeros at Im µ = ± π 3 , corresponding to the thermal cut of a fermion gas [78].Notice that the unphysical zero at real µ visible in Fig. 4 is not in the "safe" region, and turns out to be a numerical artifact.
The LY zero closest to the origin determines the radius of convergence of a Taylor expansion of log Z(μ) around µ = 0.In Fig. 7 we show the LY zeros closest to the origin, including their error bars, on a 10 3 × 4 lattice at β = 3.34, from which the radius of convergence is easily determined.
We also show how the Lee-Yang zero determined with our method compares to the spectral gap of the reduced matrix in Fig. 8.The radius of convergence of the pressure defined with our method reaches inside the region inaccessible by the traditional definition.
We have then repeated the procedure for all the available volumes.The imaginary part of the leading Lee-Yang zero in general is expected to scale as where for a first-order phase transition a = 0 and c = 1, while for a second-order phase transition a = 0 and c < 1 is given by the critical exponents of the theory [79,80].In the absence of a phase transition a > 0 and c is in general not known.Our data are not consistent with a first-order transition, neither it is consistent with a second-order transition in the 3D Ising or O(4) universality classes.At the moment our data is not precise enough to fit for all of a, b and c.Empirically we find that the infinite volume extrapolations with c = 1 lead to good fits.In this first study, we hence use a linear function in 1 V with the free parameters a and b.We will also extrapolate the radius of convergence itself with the same ansatz.In Fig. 9 we illustrate this for the case β = 3.35.The results for the extrapolated radius of convergence in the thermodynamic limit for the various temperatures investigated in this work are collected in Fig. 10.The radius of convergence is µ B T ≈ 2 and almost constant in the range of temperatures investigated in this paper.Note that this radius is larger than the radius of 3 2 mπ T ≈ 1.4 where reweighting from the phase quenched theory is expected to break down due to the onset of pion condensation.
The existence of a genuine phase transition in the ther- modynamic limit is signaled by the vanishing of the imaginary part of the Lee-Yang zero closest to the real axis.In Fig. 11 we show the imaginary part of the LY zero closest to the origin as a function of the volume at β = 3.35.This extrapolates to a finite value, indicating that the radius of convergence is not determined by a phase transition.Also this value turns out to be almost temperature independent in the range of temperatures considered here.

IV. SUMMARY AND OUTLOOK
The radius of convergence of the Taylor expansion in µ B is one of the most sought-after quantities in the finitetemperature QCD community [26,41,[81][82][83][84].The main reasons are that it gives a lower bound on the location of the critical end point, and also that it gives us insight in how far one can trust the equation of state calculated from a Taylor expansion around µ B = 0.Even if the lower bound on the location of the CEP happens to be not very stringent, this validity region is still important, since viscous hydrodynamic simulations of heavy ion collisions usually explore a very wide range of temperatures and baryochemical potentials [85].
In this paper we have made two suggestions about how to obtain a reliable estimate of this quantity in the framework of reweighting methods using the reduced-matrix approach [32-35, 37, 38].The first one concerns the definition of the rooted staggered determinant itself at finite µ B .We have proposed an alternative definition that takes care of the analyticity issues of the partition function by making it -up to an exponential factor -a polynomial in the fugacity.This involves a procedure we dubbed geometric matching, where to define the N f = 2 determinant at finite µ B we judiciously group eigenvalues of the reduced matrix in pairs, and substitute each pair with the geometric mean of its members.This follows the spirit of the proposal of Ref. [52].While on the coarse N t = 4 lattices used in our study there are no clear taste quartets or pairs yet, in the continuum our procedure is expected to give a reasonable definition of rooted determinant on each gauge configuration.
Our definition allows one to calculate the radius of convergence at any finite lattice spacing and later take the continuum limit, instead of having to take the continuum limit of the Taylor coefficients first.We believe this is beneficial, since the strong correlation between the statistical errors of the high-order Taylor coefficients present on a single ensemble are what makes possible the determination of the radius of convergence in the first place [41], and taking the continuum limit of the individual Taylor coefficients first can potentially wash out these strong correlations.
The second suggestion concerns the way the radius of convergence is calculated numerically: we have demonstrated that a brute-force calculation of the Lee-Yang polynomial via reweighting, and a brute-force calculation of its roots is possible, at least for small lattices.This provides a direct determination of the radius of convergence without relying on a finite-order Taylor expansion.This second suggestion could also be used with other fermion discretizations, like Wilson fermions, where the rooting ambiguity discussed in this paper is not present.
Our numerical results on N t = 4 2-stout improved staggered lattices suggest a radius of convergence of µ B T ≈ 2 at and slightly below T c , for µ s = 0.While extending such a study to finer lattices is certainly a challenge, we believe it is a challenge worth pursuing, based on the conceptual advantage of having an actual Lee-Yang polynomial at a finite lattice spacing, instead of a function that is nonanalytic in dense regions of the complex chemical potential plane, as is the case with the standard definition of staggered rooting.
While the most important shortcoming of the present work is the use of a rather coarse lattice, there are further improvements that could be made in the future.For one, for a precision determination of the radius of convergence our assumption of c = 1 in Eq. ( 12) should be relaxed.This requires more statistics on the currently used volumes and also data on larger volumes.Another possible improvement is the implementation of the strangeness neutrality condition, so that the choice of our µ S more closely resembles that of heavy ion collision experiments.This requires generalizing our procedure to identify taste quartets of eigenvalues.The most straightforward way to achieve this is to simply repeat the pairing procedure twice.This would allow for an arbitrary choice of µ s , in particular ⟨S⟩ = 0 could also be implemented.
Our discussion on the radius of convergence with the usual staggered determinant suggests that high-order cumulants of the baryon number are quite sensitive to taste symmetry breaking.It is therefore also an important question for the future to what extent the Taylor coefficients in µ B obtained with our new definition of the N f = 2 determinant differ from those obtained with the standard definition, and which of the two has a better continuum scaling.Once the issue of the continuum limit is under control, the finite-volume scaling of the Taylor coefficients in the continuum is also an important future question.

FIG. 1 .
FIG.1.Eigenvalues of the reduced matrix P close to the unit circle on 3000 configurations at a lattice volume of 123  ×4 and β = 3.35 superimposed into a single plot.A small gap is also present near 0, but it is not visible on the plot.

FIG. 2 .
FIG. 2. Eigenvalues of P after 50 stout smearing steps with ρ = 0.05 on a single gauge configuration on a 6 3 × 4 lattice.The formation of doublets of eigenvalues is apparent.

FIG. 3 .
FIG.3.The original eigenvalues of P (red circles) together with the new eigenvalues after the rooting via geometric matching (black boxes) for a 123  × 4 lattice at β = 3.35.

FIG. 4 .
FIG.4.Lee-Yang zeros in the complex µ T plane.Here the lattice size is 123  × 4 in lattice units and β = 3.35.The statistical errors are not shown.

FIG. 7 .
FIG. 7. Radius of convergence on a 10 3 × 4 lattice at β = 3.34, together with the first few Lee-Yang zeros with error bars.

FIG. 8 .FIG. 9 .
FIG.8.The logarithms of the eigenvalues of the reduced matrix P for 3000 configurations superimposed on a single plot, as well as the leading Lee-Yang zeros determined with the new method for β = 3.35 on a 123  ×4 lattice.The radius of convergence of the pressure defined with our method reaches inside the region inaccessible by the traditional definition.