Contribution to understanding the phase structure of strong interaction matter: Lee-Yang edge singularities from lattice QCD

We present a calculation of the net baryon number density as a function of imaginary baryon number chemical potential, obtained with highly improved staggered quarks (HISQ) at temporal lattice extent of $N_\tau=4,6$. We construct various rational function approximations of the lattice data and discuss how poles in the complex plane can be determined from them. We compare our results of the singularities in the chemical potential plane to the theoretically expected positions of the Lee-Yang edge singularity in the vicinity of the Roberge-Weiss and chiral phase transitions. We find a temperature scaling that is in accordance with the expected power law behavior.


I. INTRODUCTION
The phase diagram of quantum chromodynamics (QCD) belongs to the most pressing open issues in high energy physics. With large scale experimental programs at RHIC and LHC the phase diagram is scanned for hints of a critical point or a first order phase transition. In addition, many ab initio calculations of lattice QCD are performed to infer on the QCD phase diagram.
Unfortunately, the notorious sign problem hampers numerical studies of the QCD phase diagram. At vanishing baryon chemical potential (µ B ≡ 0), lattice QCD calculations rely on Monte Carlo methods for an efficient sampling of the QCD partition sum. At nonvanishing baryon chemical potentials (µ B > 0), standard Monte Carlo methods cease working as the fermion determinant becomes genuinely complex. Hence, the kernel of the QCD partition sum is strongly oscillating with increasing lattice volumes.
With this study we systematically investigate singularities of the grand canonical potential in the complex chemical potential plane, which we identify from an (analytically continued) rational approximation of lattice data obtained at purely imaginary chemical potentials. The rational approximation of the net baryon number density is done in consistency with the second, third and fourth order cumulants of the baryon number density. In this sense our method could be seen as a combination of the Taylor expansion approach and the imaginary chemical potential method. The position of those singularities provides very valuable information on the QCD phase diagram. We find that they are in agreement with the critical scaling of the Lee-Yang edge singularities in the vicinity of the Roberge-Weiss transition and the chiral transition. We also discuss the scaling of the Lee-Yang edge singularity in the vicinity of a hypothetical critical end point. Finally, we point out that the position of the singularities can be used to estimate the radius of convergence of any analytic expansion and to extract nonuniversal parameters that map QCD to the universal scaling function. Among the latter might also be the position of the QCD critical end point, which was demonstrated recently in the case of the Gross-Neveu model [27]. Genuine Lee-Yang zeros have been recently also studied in other works [28,29]. This paper is organized as follows. In Sec. II we introduce the scaling theory of the Lee-Yang edge singularities and apply them to the cases of the Roberge-Weiss transition and the chiral transition of QCD. We also discuss the case of the QCD critical endpoint. In Sec. III, we provide details of our lattice QCD calculations and in Sec. IV we discuss our strategy for determining rational approximations to our lattice data. We present our findings of the extracted singularities in Sec. V. We summarize and conclude in Sec. VI. In order to consolidate our findings, we compiled a number of Appendixes (B-F), where we discuss numerical issues related to our rational approximations, which are based on a multipoint Padé method.

II. EXPECTED SINGULARITIES
The grand canonical partition function in lattice QCD, Z GC = Z(V, T, µ B ) has the form of a (high dimensional) polynomial at any finite volume V and is positive for arXiv:2110.15933v2 [hep-lat] 7 Mar 2022 real values of the temperature T and baryon chemical potential µ B . However, following Lee, Yang [30,31] and Fischer [32], we point out that Z GC exhibits many zeros in the complex chemical potential plane, which can be used to extract valuable information on the phase transitions that may occur in the system. A physical phase transition in the thermodynamic limit can be identified when for V → ∞ one of the complex zeros approaches a point with real parameters T, µ B .
The grand canonical potential log(Z GC ) diverges when one approaches a zero of Z GC . Singularities of the grand canonical potential and its derivatives will limit any analytic expansion performed at zero or purely imaginary chemical potentials. The positions of the zeros of Z GC can thus also be used to estimate the radius of convergence of the Taylor expansion method.
In different (T, µ B ) regions, the QCD partition function can be approximately described by effective theories. At high temperatures, far above the QCD transition, we might be able to use a free Fermi gas to describe the thermodynamic behavior of the quarks. Below the QCD crossover temperature T pc , the Hadron resonance gas is known to describe the bulk thermodynamics of QCD matter quite well [33][34][35]. However, we can also consider universal behavior in the vicinity of the Roberge-Weiss, the chiral transition or even in the vicinity of the QCD critical end point (if existing). In particular, we can predict the positions of singularities in the complex µ B plane, by applying suitable mappings from the parameter space of QCD to the relevant scaling fields in the vicinity of a critical point. In the following, we make extensive use of the fact that the scaling function of the order parameter f G (z) exhibits a branch cut singularity at z = z c , where the scaling variable z is expressed in terms of the reduced temperature t and symmetry breaking field h as z = t/|h| 1/βδ . The universal position z c of the universal singularity, known as the Lee-Yang edge singularity, has been recently determined for different universality classes [36]. The three distinct scaling approaches are visualized in Fig. 1. At this point not all of the nonuniversal normalization constants that fix the above mentioned mappings are known. We thus vary some of the parameters to give an impression of the functional dependence, which is discussed in more detail below. Also shown are identified Lee-Yang-edge singularities from our lattice QCD calculations on N τ = 4, 6 lattices. The determination of the data points is discussed in Sec. V.

A. Singularities in the vicinity of the Roberge-Weiss(RW) transition
The RW critical point (µ B /T = iπ) is a remnant of the Z(3) symmetry in the quenched (m q → ∞) limit of QCD, where m q specifies the quark masses, and the QCD partition function has a reflection symmetry around µ B /T = iπ [37]. The nature of the RW end point could be either a first order triple point or a second order Z(2) critical point depending on the value of the quark masses [38]. In (2 + 1)-flavor QCD with a physical value of the light quark masses with improved lattice discretization, which is discussed here, one finds a second order Z(2) critical point [39][40][41]. The order parameter in the vicinity of a second order transition can be written as the sum of a universal and a regular part, where, t, h are scaling fields and β, δ are critical exponents. Also, f G is the universal scaling function for the order parameter, whereas M reg accounts for regular contributions.
The relevant scaling fields for a Z(2) symmetric second order RW transition can be defined as whereμ B = µ B /T and t 0 , h 0 , and T RW are nonuniversal parameters. T RW is the RW transition temperature, which is known for our particular lattice setup [41]. We /t 0 . This solution has also been used in [48] to derive an estimate of the radius of convergence. Equation (8) is visualized in Fig. 1 as green band, where we chose m l /m phys s = 1/27 (physical mass ratio) and T c = 147 MeV which is our best estimate for the chiral transition temperature for N τ = 6. It is however obvious that T c does not alter the line of constant z = z c much; it mainly alters the normalization of the temperature behavior. The curvature κ B 2 is chosen as κ B 2 = 0.012. We vary z 0 from 1.5 to 2.5 which generates the width of the green band. Our best estimate for N τ = 6 is z 0 = 2.35, which stems from scaling fits to the magnetic equation of state.

C. The QCD critical point
The same kind of scaling is expected to hold close to the QCD critical point. Unfortunately, the mapping to the universal theory with Z(2) symmetry is unknown. A frequently used Ansatz for the scaling fields is based on a linear mapping where the critical point is located at (T cep , µ cep ). This Ansatz leads to [49] where c 1 is given by the slope of the transition line at the critical point and c 2 is related to the angle between the (t = 0) and (h = 0) lines. For the red band in Fig. 1, we chose c 1 = −2κ B 2 = −0.024, assuming that the transition line is a quadratic function all the way down to the critical point. We vary µ cep from 500 to 630 MeV, which generates the width of the red band. For T cep , we chose in accordance with the Ansatz for the transition line: , with T c = 156.5 MeV. The prefactor c 2 |z c | −βδ was chosen to be 0.5. In principle, we need at least four data points in the scaling region of the critical point to determine all the unknown nonuniversal parameters, including the location of the QCD critical point T cep , µ cep . For the Gross-Neveu model, it has been recently shown that T cep , µ cep , c 1 , c 2 can be determined from a fit to the above Ansatz [27]. It will be very interesting to apply such scaling fits also to lattice QCD data, which is however currently beyond the scope of this exploratory study.

D. Thermal singularities
The Lee-Yang edge singularities are not the only singularities in the complex µ B /T plane. At temperatures above the QCD crossover (T > T pc ), we expect that quasifree quarks are the relevant degrees of freedom in the system. Quarks are distributed according to the Fermi-Dirac distribution f p (T, µ) = 1/(exp (ε p − µ)/T + 1). The singularities of this function are located at ±iπT ±ε p . In particular, the singularities which are closest to the origin are located at ±iπT ± 0 , where ε 0 = m is the rest mass of the particle. The thermal singularities of quasifree quarks are modified by residual inter-actions as long as we are not considering the Stefan-Boltzmann limit (T → ∞). We expect that to leading order these modifications are expressed through a substantially larger thermal massm(T ) m. The analytic structure of the Fermi-Dirac distribution function interferes with the scaling of the Lee-Yang edge singularity of the Roberge-Weiss transition, as given in Eqs. (4) and (5). It is a priori not clear which type of singularities are closer to the origin/imaginary axis and can thus be found by a Padé/rational approximation of the data. As a result of this study, we find that the leading singularities at Im[µ B /T ] = ±π follow the RW scaling, see Sec. V A.

III. LATTICE SETUP AND OBSERVABLES
The partition function of (2+1)-flavor of highly improved staggered quarks (HISQ) [50] with imaginary chemical potential can be written as where M (m, iµ I ) represents the fermion matrix of a HISQ flavor with mass m and chemical potential µ = iµ I . The first determinant represents the two degenerate light flavors (up and down quarks). For the gauge part S G (U ), we are using the Symanzik improved Wilson action, which is correct to O(a 2 ) in the lattice spacing. For the gauge field generation, we were using the SIMULATe-QCD package [51] with and implementation of the rational hybrid Monte Carlo algorithm (RHMC) [52]. The lattice bare parameters are used from various publications of HotQCD. The lattice bare quark masses are varied with the lattice coupling such that for each coupling physical meson masses are obtained; i.e., we stay on the line of constant physics (LCP). Here, we make use of the parametrization of the LCP (for the physical value of the pion mass, m l /m phys s = 1/27) obtained and refined in previous works [53][54][55]. The same holds true for the scale setting, where we used the parametrization of the β function based on the kaon decay constant. For simplicity, we fix the ratio of the explored chemical potential in this study to µ l /µ s = 1.
The observables we calculate are the cumulants of the net baryon number density, given as tives generate traces of the type Tr[(M −1 ∂M/∂µ X ) n ], which we evaluate with the random noise method, using O(500) random vectors. For more details on the required traces and the method of evaluation, see, e.g., [3].
For obvious reasons (sign problem), we perform our calculations at purely imaginary baryon chemical potential iμ I B , withμ I B ∈ R. Exploiting all symmetries, we restrict values forμ I B to half the period, i.e.,μ I B ∈ [0, π]. The symmetries of the partition function generate specific properties of the observables χ B n . At imaginary chemical potential, they are imaginary and odd functions ofμ I B for odd n and real and even functions ofμ I B for even n. These properties have been verified by us and can be seen from Fig. 2, where we show results for the first three cumulants. Preliminary results were presented in [56]. The data are tabulated in Appendix A.
Since this is an exploratory study, we leave the continuum extrapolation for later publications and perform the calculations on rather course lattices, 24 3 ×4 and 36 3 ×6. We note however that we have deliberately chosen rather large spatial volumes (we have aspect ratio N σ /N τ = 6) in order to minimize finite size effects which are expected to become large in the vicinity of a phase transition.

A. Padé approximants
Padé approximants [57] are a popular subject in approximation theory. The main idea is to approximate a given function f (x) with a rational function whose derivatives agree with those of f (x) up to a given order. This can be easily rephrased in terms of power series. To set up our notations, we first consider a so-called single point [m/n] Padé. Suppose the Taylor expansion of f (x) about a single point (in what follows, it is useful to take this point as x = 0) is known up to a certain order (say L), We denote R m n (x) the [m/n] Padé approximant we are looking for, which is the ratio of two polynomials (P m andQ n ) of order m and n, respectively. By discarding a nontrivial b 0 and writingQ n (x) = 1 + Q n (x), in Eq. (15) we have made a definite choice for the coefficients {a i , b j } our approximant essentially depends on. Given our knowledge of the Taylor expansion for f (x), in principle, our choice can be for any order [m/n] such that m + n + 1 = L + 1.
Strictly speaking, not all such [m/n] approximants exist. Rational functions are not the only viable solutions in approximation theory; polynomial approximants are very popular as well. We stress from the very beginning a main virtue of rational functions we are interested in; they provide a natural handle to probing the singularities structure of f (x). Quite trivially, the singularities of R m n (x) are coming from the zeros of 1 + Q n (x). To make the latter statement more precise, we need to more precisely state what we mean by a [m/n] approximant. Our attitude is quite pragmatic; we consider a R m n (x) for a given choice of [m/n] and solve for the {a i , b j } coefficients given that the power series of f (x) is known to a given order L such that m + n + 1 = L + 1. Given our input, we are not guaranteed that P m and 1 + Q n (x) are coprime polynomials [58]. If this is not the case (i.e., there is a nontrivial greatest common divisor), singularities of R m n (x) are coming from those zeros of 1 + Q n (x) which either are not zeros of P m (x) or are zeros of 1 + Q n (x) of a higher order than they are of P m (x). While trivial, this very last statement will be of some relevance in the following. We are also quite concerned with yet another property of R m n (x). Poles are the only singularities we can find in a rational function; still, there are signatures of other singularities of f (x) (e.g., branch cuts) which we can recognize in R m n (x) once the coefficients of the latter have been determined to reconstruct the power series of f (x). We discuss this topic (which indeed is relevant in the our analysis) in Appendix E. The most direct way of solving for the unknown coefficients is by approximation through order; we match coefficients of different powers of x between the (unknown) rational function and the (known) Taylor series by demanding functional independence (up to order x L ). That means we rewrite and match Equation.
(16) defines a set of simultaneous, linear equations, which can be solved by a convenient linear solver. We stress that most often the linear system is not singular but ill-conditioned, which is a warning that determining the solution can be hard. Most importantly, we should keep in mind that our knowledge of the derivatives of f (x) is coming from stochastic evaluation via Monte Carlo simulations. Despite this, we see that we can manage to find the information we are aiming at. Approximation through order somehow hides that the information we are making use of is coming from derivatives of f (x). The solution encoded in Eq. (16) can of course also be obtained by evaluating in x = 0 the tower of relationships Yet another way of obtaining the unknown a i , b j is to solve the set of [59], While Eqs. (16), (17), and (18) are equivalent, the latter is somehow less convenient, not being linear and typically requires the use of computer algebra tools like Mathematica. In our particular problem, we explicitly showed that the three return the same results (to a very good approximation); Eq. (17) has been to a large extent our preferred choice.
There exists a significant amount of literature on single point Padé approximants (about existence, uniqueness, and convergence). This is not true for the so-called multipoint Padé to the same extent. The construction of a multipoint Padé can be extremely useful in situations when Taylor coefficients for a function about a single point are not known to higher orders, but instead either the function values are known or a few Taylor coefficients are known about (possibly) many points. Since this is precisely the situation we face in our lattice studies of QCD at finite chemical potential, it is useful to understand how to build rational approximations from these multipoints.
Extending what we saw above to multipoints is straightforward; in particular, we extend the formalism encoded in Eq. (17). A few Taylor coefficients (i.e., derivatives) of a function f (x) known at a collection of points {x i | i = 1 . . . N } are consistent with the approximation to f provided by Eq. (15) if they satisfy the set of equations which is once again a linear system in n + m + 1 unknowns where now n + m + 1 = N i=1 (L i + 1). In the previous formula, the highest order of derivative which we know (i.e., L i ) can be different for different points. Multipoints Padé will be our choice for the analysis of this work. We now proceed to discuss a few technical details of our implementation, referring the reader to Appendixes B-F for extra remarks and comments.
B. Padé approximants for the QCD net baryon number density from imaginary chemical potential This is not the first time Padé approximants are applied to the study of finite density (lattice) QCD [60][61][62]. A recent paper has, in particular, proposed to join Padé analysis and Bayesian methods, with applications to the study of the crossover line [63]. In a way that is close to the spirit of this work, in recent times, Padé approximants have been successfully used to probe the singularity structure of simple theories in the context of the Lefschetz thimble method [64]. To our knowledge, this work is in a sense the first attempt at a systematic study of the QCD phase diagram and, in particular, of Lee-Yang edge singularities, building on Padé analysis. The function we want to approximate by rational functions is the net baryon number density (19) as derivatives. Our main goal is to get signatures of singularities of ) plane (at fixed values of T and V ). Ultimately we aim to understand the phase diagram of the theory. In particular, we find clear evidence of the Roberge-Weiss transition in the µ I B − T plane. Most importantly, if at some point we could find evidence of singularities eventually pinching the (real) µ R B axis, then we would be in the presence of a QCD critical point candidate.
We have already seen that Eq. (19) is not the only way to solve for the coefficients entering the Padé approximants (15). Not only, e.g., does the multipoints version of Eq. (18) works as well, but also other formalisms could be (and actually were) used in our analysis. In the construction of these other formalisms, a key point is that the values of the χ B n , i.e., the function f and its derivatives in Eq. (19), are known to a limited precision since they are evaluated by Monte Carlo. Our Padé analysis was performed following three different approaches, aiming at assessing their mutual consistency.
1. The solution of the linear system (19) has been worked out in two different ways, namely, • One can build the system by writing the most general form for R m n (x), i.e., that of (15). • One can instead impose the form with the coefficients {a i } and {b j } that turn out to be real. This form ensures the following: (a) The function χ B 1 (T, V, µ B ) has the right parity (it is an odd function). (b) As a consequence of the coefficients being real valued, for imaginary µ B = µ I B , the odd cumulants χ B 2n+1 (T, V, µ I B ) are imaginary valued, while the even χ B 2n (T, V, µ I B ) are real valued, as it must be. (c) When Eq. (20) is computed for real µ B = µ R B , the cumulants are real; i.e., the analytic continuation one is typically interested in is guaranteed to be meaningful. Notice that taking into account different functional forms for R m n is not the end of the story. Another alternative which can (and actually was) taken into account is whether one • performs the Padé analysis in the (original) complex-µ B plane or • goes through a conformal map µ B = φ(ν) and performs the Padé analysis in the complex-ν plane.
2. Because of the cumulants being known to finite precision, the minimization of a generalized χ 2 is an obvious alternative to the solution of (19). Suppose we want R m n (x) to be a Padé approximant for the function f (x) whose values and derivatives we know at given points n depends on can be fixed minimizing the generalized χ 2 , Of course, all the alternatives that we commented in 1 (namely, different functional forms for R m n , use of conformal maps) can be also implemented in this approach.
3. Both 1 and 2 make use of the knowledge of f (x) (and its derivatives) at given points; i.e., the only information on f (x) we have is at a finite (possibly small) number of points. One could instead compute a smooth interpolation of f (x) before entering the Padé analysis.

C. Results of Padé analysis of net baryon number density
The focus of our analysis is on singularities of the net baryon number density. Still, before proceeding to this, we make a short digression on a feature which is worth discussing. In investigating the phase diagram of QCD in the (imaginary chemical potential-temperature) µ I B -T plane, and, in particular, in the study of the Roberge-Weiss transition, a prominent role is played by the free energy as a function ofμ I B (at given values of the temperature T ); a cartoon for this quantity is often plotted. Since we have a function R m n (μ I B ) approximating the net baryon density, we can obtain the free energy F (μ I B ) by (numerical) integration. In Fig. 3, we display the free energy F (μ I B ) at three different temperatures; the profile clearly gets closer to a cusp as the temperature gets closer to T = T RW (to the extent that the transition can be detected on a finite volume).
We now inspect how well our rational approximants describe the data. On top of that, we are interested in the analytic continuation of results from imaginary to real values of the baryonic chemical potential (this is in the end a key issue in any imaginary-µ B study of finite density lattice QCD). Finally, we present the relevant singularity pattern which emerges from our analysis. In Fig. 4, we display what we get both for imaginary and for real baryonic chemical potential. For three of the temperatures we probed on N τ = 4, we plot the results we got from the solution of Eq. (19) for both functional forms (15) and (20). For imaginary values of the baryonic chemical potential, the two solutions are de facto indistinguishable. For real values (analytic continuation), the real parts are quite close to each other, a significant discrepancy between the two different Ansätze being there only at T = T RW and forμ R B > π. As for imaginary parts, Eq. (20) is guaranteed to return zero; it is interesting to notice that also the solution we got for the Ansatz (15) has a quite tiny imaginary part (at least up toμ R B ∼ π). All this can be taken as an indication of reasonably tiny systematic effects as far as the dependence on the precise form of the Padé approximants is concerned. All in all, the indeterminations we have to live with when we analytically continue our results to real baryonic chemical potential seem to be competitive when we compare to other methods. This is true despite the fact that, inspecting Fig. 4, a few spikes are clearly visible; we see in what sense they do not come as a surprise and are in fact harmless. In Fig. 5, we plot the singularity pattern we get at three of the temperatures we probed on N τ = 4. We once again present results we get for both functional forms (15) and (20). A few remarks are in order, which we invite the reader to consider taking into account the points that we make in Appendixes C, E, and, F.  (15) and (20); for the latter, the imaginary part is guaranteed to be zero. µ I B = π and indeed to a very good accuracy they do.
• The signature for a branch cut is clearly visible at T = T RW = 201.4 MeV, for both Ansätze (15) (see upper row) and (20) (lower row). Notice that the latter is by construction sensitive to all the four replicas of the same singularity, as expected for symmetry reasons: if we find a singularity in z, then also −z and −z must be singular points. This is a general feature, which is clear at all temperatures; (15) instead only captures singularities in the upper half plane (i.e., we see one of the two symmetries).
• At T = 186.3 MeV (this is the next to highest temperature that we probed), plots apparently allude to a branch cut as well, while at T = 167.4 MeV the singularity shows up as a simple pole. Much the same happens at the remaining temperatures (i.e., T = 176.6 MeV and T = 160.4 MeV are consistent with the appearance of simple poles).
• While the pattern of the relevant pieces of information (i.e., true poles and zeros) is the same for different functional forms, the pattern of other zeros and poles depends on the functional form. Notice that the mechanism of zero-pole cancellations is manifest; these cancellations are due to numerical noise. In a sense, we see fake information, which would not be there for exact data, but since noise is not that much, this fake information is close to disappearing.
• While the zero-pole cancellations seem almost perfect in Fig. 5, Fig. 4 is warning us that this is not really the case, and this is the reason for the spikes which we see there. We encourage the reader to spot which singular points are responsible for the spikes. Again, the almost perfect cancellations in Fig. 5 reveal that these spikes are harmless.
• As should be clear, the big spike atμ B = iπ for T = T RW is a different story: this is indeed the Roberge-Weiss transition showing up.
We repeated our analysis hunting for singular points in the complex-fugacity plane; after mapping our measurements to this plane, we performed Padé analysis in the z = eμ B variable. There are at least two reasons for such an (additional) analysis. First of all, we want to make sure that the information which we get is stable and does not disappear once we change the variable. Also, the linear systems which we have to solve are typically illconditioned. Due to the nature of the conformal map, this feature disappears (in a sense, we can trust results to a higher level of confidence). In Fig. 6, we present the singularity pattern in the complex-fugacity plane.
• Since our original data are taken on the imaginary axis (µ B = iµ I B ), in the complex fugacity plane, we end up on the unit circle |z| = 1. In view of the observations that we make in Appendix D, notice that this is a very convenient location with respect to the location of the singularities that we detect.
• Since singularities are expected atμ I B = π, in the fugacity plane, they should show up on the real axis, and indeed, they do. Due to the relative positions with respect to the input data (|z| = 1), we are not sensitive to any (symmetry) replica (that is, one single singularity shows up).
• All the other features (e.g., zero-pole cancellations) show up much the same as they do in the original µ B plane.
In Table I, we collect all the findings that we discussed so far. In particular, for each temperature that we probed at N τ = 4, we list the nearest singularities as obtained (a) from the solution of the linear system (19) in the µ B /T plane (method I), (b) from the minimisation of the generalised χ 2 (21) (method II), and (c) from the solution of the linear system (19) in fugacity plane, both mapping back results in the original plane and inspecting them in the fugacity plane (method III* and III). The errors are computed out of a bootstrap procedure in which we repeat our Padé analysis letting the input data (i.e., the results of our Monte Carlo measurements) vary within errors. As one can see, results are well consistent. The singularities which we have been discussing so far (and that are listed in Table I) are not the only ones on display in Fig. 1. Results obtained on N τ = 6 at T = 145 MeV apparently point at a singular point that could be consistent with a chiral singularity. While this result is intriguing, in this case extra care is in order.
• In this case, we have a (far) enhanced dependence on the interval our Padé analysis takes into account.
In particular, this singularity shows up if we limit our analysis toμ I B ∈ [0, π]. • The result which is shown in Fig. 1 comes from the minimization of the generalized χ 2 (21) taking (15) as an Ansatz, with m = n = 4. This choice returns the bestχ 2 value.
• A consistent result for the singularity is found from other methods if we limit the analysis to the same interval (μ I B ∈ [0, π]) (even changing a bit the degree, which thing is easier in this approach). This singularity is not that stable under the variation of the interval. While this is not a priori that surprising (our multipoint Padé analysis is interval sensi-  (12)  In Fig. 7, we display the net baryon number density as obtained on N τ = 6 at T = 145 MeV (notice that the signal is substantially tinier than at higher temperatures).
Here we have 2 × 2 options displayed: data are compared to rational approximants obtained (a) either from the generalized χ 2 (21) or from the solution of Eq. (19) for the basic functional form (15) and (b) taking into account data either forμ I B ∈ [0, π] or forμ I B ∈ [0, 2π]. As one can see, the rational approximants we get from the two methods are (always) substantially equivalent; indeed, there is a difference when it comes to taking into account a larger or smallerμ I B interval. Notice, however, that for µ I B ∈ [0, π] every solution is de facto indistinguishable from any other. We also show analytical continuations; taking into account data from an extendedμ I B interval results in an imaginary part staying very close to zero in a wider interval of real chemical potential. In Fig. 8, we display the analytic structure we get from the different 2×2 options. As one can see, results are again very much consistent whatever method we choose for computing the rational approximants. As anticipated, the singularities we found are indeed different for different inputμ I B interval. When we take input from the extendedμ I B ∈ [0, 2π] interval, we apparently get what one would regard as a thermal singularity. The outcome is pretty different when input is taken from the restrictedμ I B ∈ [0, π] interval. We see that what we get in this case is a chiral singularity candidate. We stress that this ambiguity shows up only in this case (i.e., for this lowest temperature, which is only probed on N τ = 6). We stress once again that our multipoint Padé analysis is interval sensitive and it could well be that in different intervals we have access to different pieces of information. The fact that this is possibly the only piece of information related to chiral symmetry breaking makes all this intriguing and definitely deserving further investigation. Indeed we are working on this, in particular, repeating our Padé analysis for the chiral condensate, which is the relevant order parameter for chiral symmetry breaking.

A. The Roberge-Weiss critical region
The nearest singularities which we identified in the temperature range 201 MeV < T < 160 MeV from our Padè approximations presented in the last section and which appeared to be stable are listed in Table I. Despite the fact that we could not observe indications for a branch cut connected to all these singularities, we now demonstrate that they can indeed be identified with Lee-Yang edge singularities of the Roberge-Weiss critical point, i.e., that they scale in accordance with our expectations presented in Sec. II A. In particular, it is obvious that we obtain for the imaginary partμ I LY = Im[ µ B T ] = π within errors for all temperatures and methods as demanded by Eq. (5). In order to show that the real part scales in accordance with Eq. (4) we perform fits to the data listed in Table I with the Ansatẑ with fit-parameter a, b. For the Roberge-Weiss critical temperature, we set T RW = 201.4 MeV in accordance with [40]. We fixed the critical exponents to that of the Ising universality class; i.e., we have βδ ≈ 1.5635. The parameter b is added to capture the leading order finite size effects. Since our calculations are done in a finite volume, we expect that the Lee-Yang edge singularities will not reach the real h axis, which is here theμ I B axis. Or with other words, there is no phase transition in a finite volume. A more elaborate finite size analysis will be left for future publications.
The fits work quite well [67] and are shown in Fig. 9. Results for the fit parameter and reduced χ 2 values are given in Table II.
Besides the demonstration for scaling, we can relate our results for the fit parameter a to the nonuniversal constant z 0 . From Eq. (4), we obtain Using the value |z c | = 2.452 for the 3d-Ising universality class [36], we obtain z 0 ≈ 9.2 -9.5. The specific values for our three data sets are given in Table II. To our knowledge, that is the first determination of z 0 for the Roberge-Weiss transition. Note, however, that the value is obtained on course lattices (N τ = 4) with no proper continuum extrapolation yet.
In principle, z 0 could also be determined from a fit to the magnetic equation of state (EoS). Here we anticipate more severe corrections from the finite size, as well as large contributions from regular terms. In particular, more data close to the RW transition are needed to obtain a reliable fit. A determination of z 0 from a fit to the magnetic EoS is thus beyond the scope of this work.

B. The chiral critical region
We have probed one additional temperature below the pseudocritical phase transition temperature of (2 + 1)flavor QCD, namely, T = 145 MeV. For this low temperature, the calculations have been done on 36 3 ×6 lattices. To put the temperature value into perspective, we recall the continuum extrapolated numbers for the pseudocritical temperature T pc = (156.5 ± 1.5) MeV [45] and the chiral critical temperature T c = 132 +3 −6 MeV [44]. We also note that the corresponding N τ = 6 results are 10-15 MeV higher. In conclusion, the probed temperature of T = 145 MeV is compatible with the chiral critical temperature, and we thus expected it to be sensitive to chiral scaling.
We now compare the position of the singularity we find for this temperature with the expected position of the Lee-Yang edge singularity, governed by O(2) critical behavior [68]. Hence, we fix the critical exponents to βδ = 1.6682. The chiral transition has been subject to various lattice QCD studies in the past, the nonuniversal parameters that appear in Eq. (8) are known to some extent, as discussed already in Sec. II B. In Fig. 10, we calculate the 68% and 95% confidence areas of the expected Lee-Yang edge singularity when we vary the nonuniversal parameter under the assumption of Gaussian distributed errors. In particular, we chose for the N τ = 6 specific values and errors [69], and in addition, we take |z c | = 2.032 [36]. As can be seen from Fig. 10, the results from the rational approximation to our data (method II), (μ R B ,μ I B ) = (3.03(28), 1.61(10)), lie within the 68% confidence area of this prediction.

VI. SUMMARY AND CONCLUSIONS
We computed cumulants of the net baryon number density as a function of the imaginary baryon number chemical potential in lattice QCD, the fermionic regularization being that of highly improved staggered quarks (HISQ). The results were the input for a multipoint Padé analysis by which rational approximations were calculated, with various choices of both the functional forms of the latter and of the methods by which the approximants were determined. The results have been shown to be stable, in particular, also if we repeat our analysis in the fugacity plane (i.e., after a conformal map). Our rational approximations not only describe very well the data but appear to be quite well under control when we analytically continue them to real values of the baryonic chemical potential. The main focus of our analysis has been on the singularity structure that we can infer from the complex poles of our rational approximations. By comparing the latter with the theoretically expected Lee-Yang edge singularities in the vicinity of the Roberge-Weiss phase transition, we found a quite good agreement. In particular, the temperature scaling of the singularities is consistent with the expected power law behavior. We also found a preliminary evidence of a singular point consistent with the phase transition which is expected in the chiral limit of (2+1)flavor QCD in a staggered regularization. All our findings (and, in particular, the last that we mentioned) deserve further investigation by getting more precise measurements and probing less coarse lattices, a target that we are aiming at in the near future. An interesting task that is also in front of us is that of comparing our methodology and results with other recent approaches to the study of the analytical structure of finite density QCD, e.g. that of [29,70]. All data from our calculations, presented in the figures of this papers can be found in [71].

ACKNOWLEDGEMENTS
Our work is dedicated to our late colleague and friend E. Laermann. J.G. and C.S. would like to thank F. Karsch for stimulating discussions. C.S. would also like to thank S. Mukerherjee and V. Skokov for discussions and A. Lahiri for providing data on the cutoff dependence of z 0 , T c . This work was supported by The gauge fields have been generated with a rational hybrid Monte Carlo algorithm (RHMC). In Tables III  and IV we list results from calculation on the 24 3 × 4 and 36 3 × 6 lattices, respectively. Also listed are the number of configurations on which we have measured the observables and which are separated by 10 RHMC trajectories of length 0.5-1.0.

Appendix B: MULTIPOINT VS SINGLE-POINT PADÉ
As mentioned before, most of the literature that exists on existence, uniqueness and convergence [57,72] of Padé sequences exists mainly for single-point Padé expansions wherein the rational approximation is constructed from a single Taylor expansion with arbitrarily many Taylor coefficients [73]. On the other hand, we may be presented with a situation in which we have low-order Taylor data but at arbitrarily many points. This is known in the literature as multipoint Padé -but most commonly only values at other points are used.
In our work we also use higher Taylor coefficients at other points. Since not a lot of literature on multipoint Padé exists, we will validate our findings with numerical experiments conducted on a number of test functions. Based on our numerical experiments it is also shown that there are situations in which a multipoint Padé does better while there are other situations in which a single point may be better.
Most of our numerical experiments are based on the 1D Thirring model (a model that was studied in [64]). This model is chosen because its partition function has a known analytical solution. For the purposes of these experiments we will simulate the number density of the 1D Thirring model: (B1) where because we know the exact location of poles of the number density -it is easy to validate/invalidate our approximation. Shown in Fig. 11 are the approximations and singularity structure of the number density simulated at β = 1, L = 8 and m = 2 (same parameters used in all figures depicting the 1D Thirring model).

Appendix C: MANIFESTATION OF SPURIOUS POLES
As will be shown below, "spurious" poles can enter our analysis in two ways: 1. Firstly, we can (will) get spurious poles in noisy data -we will discuss this in the last section. But the important message is that if our function has a genuine pole, we will find a quasi-stable pole from our approximation even in noisy data and the effect of decreasing (increasing) noise will be that the pole becomes more (less) stable -eventually converging to (diverging from) the correct value in the absence of noise.
2. Secondly, even in the absence of noise -when we simulate a test function with its clean data -we will find that after a certain (optimal) order, our Padé will start spitting out spurious poles which will be exactly cancelled by corresponding zeroes. This is a clear result of demanding a very high order of approximation. This happens both for single point Padé and multipoint Padé. This effect can be seen in Fig. 12.
A note on Froissart doublets [74]: these appear as zeropole doublets in a unit circle in a Padé approximation to a  series perturbed by noise. The separation between such pairs is proportional to the scale of the noise present. Fig. 13 is an example of Froissart doublets in case of simulating pure noise.

Appendix D: INTERVAL DEPENDENCE OF MULTIPOINT PADÉ
Based on our numerical experiments, it was observed that the approximation obtained from the multipoint Padé approach is sensitive to the interval sampled. By this we mean that for some functions the signature of the singularity might be missed if the interval is not chosen appropriately.
In Fig. 14, this sensitivity is demonstrated with the help of the 1D Thirring model as our test case, since we know the positions of its singularities. The singularities which we detect from rational approximations obtained in different intervals are shown and compared with the analytic (exact) positions.
Another example where the interval dependence of a function is manifest is when the function has a branch cut. We would ask the reader to keep in mind that range dependence may not be apparent from the functional form always, whereas it is can be manifest in the    structure of zeroes and poles. In Fig. 15 is shown the example of such a function ( 2µ+1 µ+6 ).  The message to be conveyed is the appearance of "spurious" poles and zeroes exactly canceling each other on right half plane when we go very high in the order.

Appendix E: MORE NUMERICAL EXPERIMENTS
Some more numerical experiments were performed to see how the Padé approximation treats different types of singularities (poles, branch points, essential singularities etc). While some literature exists on how the single point Padé treats these singularities with varying the order of the approximation [75] -not much exists on how multipoint Padé treats these singular points. This is also a nice way to test that the approximation works. The reason we are focusing on a "cusp" like singularity is motivated from the periodicity properties of the partition function, which at and above the Roberge-Weiss temperature behaves like a cusp at the RW point and multiples of π. The free energy, by definition, also has this structure and hence it's first derivative becomes discontinuous.
The first of the examples of known functions shown below mimics the above mentioned behaviour (corner function) while the second example is just to show how the multipoint Padé handles a genuine cusp [76]. The exponent function with argument of negative absolute values and it's derivative showing a discontinuity similar to our case are shown in Fig. 16.
Poles are strictly speaking the only singularity that a Padé approximation can have. When we demand a rational approximation of an irrational function, such as the square root function above, the only way rational function can mimic the branch cut is by placing a sequence of zeroes and poles alternately along the branch cut.

Function with a genuine cusp singularity:
Just as an extra example the Padé analysis of a genuine cusp singularity is also presented below (Fig. 17 Since our data from the lattice simulations comes with noise, it is important to study the effect of noise on functions containing genuine singularities. It is already well known that in the presence of noise, poles move about around the true singularity. Also, the mean distance from the true pole increases with increasing the magnitude of error. This can be seen (again) with the help of the Thirring model (but the reader is free to choose their own test function). The figures below (Figs. 18,19,20) are intended to mimic the lattice data at least where statistical errors are concerned. For instance, from our QCD data, we have around 1% error on the χ B 1 and 10% errors on χ B 2 . In the figures, we show the effect of adding this combination of errors to the Thirring model. The take away message is that even though the poles move around -the signature of the singularity is present and consistent within errors with the true singularity.

Systematic error:
We have already seen the interval dependence of poles. Padé theory dictates that the true poles of a function remain fixed when changing orders of the Padé. While it is very clear in building single point Padé approximations what increasing or decreasing the order of a Padé means, it is not the case for multipoint. We can change the order in at least two distinct ways or a combination of them, by either increasing the number of points sampled or keeping the points fixed and increasing the derivatives at those points. The systematic errors that we want to highlight in this section are those that cause the pole to move around even at a fixed instance of statistical error while varying the order of the approximant as mentioned above (see Figs. 21 and 22).