On the QCD critical point, Lee-Yang edge singularities and Padé resummations

Mapping the phase diagram of Quantum Chromodynamics (QCD) for different temperatures and baryon densities, probed by the baryon chemical potential, $\mu_{B}$, is currently one of the major open problems in nuclear physics. State-of-the-art lattice simulations offer extensive information about the thermodynamic properties at vanishing chemical potential. For instance, it is well established that the transition between the baryonic phase and quark-gluon plasma is a smooth crossover Aoki:2006we . However, the applicability of lattice QCD at finite chemical potential is severely limited due to the fermion sign problem (for reviews see e.g. Refs. Philipsen:2007aa ; deForcrand:2009zkb ; Ding:2015ona ). A central question regarding the QCD phase diagram at finite temperature and chemical potential is whether there is a second order critical point, and a first-order transition curve emanating from it Bzdak:2019pkr . From the experimental perspective, this question has been the main motivation behind the Beam Energy Scan program at the Relativistic Heavy Ion Collider as well as the future Compressed Baryonic Matter experiment at the Facility for Antiproton and Ion Research Almaalol:2022xwv .

From the theoretical standpoint, the first principle computations of the QCD phase diagram at finite chemical potential requires tackling the fermion sign problem in some way. The two main approaches that bypass the sign problem are Taylor expanding around $\mu_{B}=0$ Allton:2002zi , and performing lattice simulations at pure imaginary chemical potential and analytically continuing from it deForcrand:2002hgr ; DElia:2002tig ; Bellwied:2015rza ; Ratti:2018ksb ; Borsanyi:2021sxv . In the former approach the equation of state is expanded in $\mu_{B}$¹¹1More precisely the expansion is in $\mu_{B}^{2}$ due to charge conjugation symmetry. as a Taylor series whose coefficients are observables evaluated at $\mu_{B}=0$ which can be computed without a sign problem. At the same time, even in the absence of the fermion sign problem, the computation of the higher order Taylor coefficients on the lattice becomes progressively more difficult due to noise. The state-of-the-art computations go up to ${\cal O}({\mu_{B}^{8}})$ (for recent results see HotQCD Bollweg:2022rps ; Bollweg:2022fqq and Wuppertal-Budapest Borsanyi:2020fev Collaborations). Currently, based on the available data, the existence of a critical point for $\mu_{B}/T\lesssim 3$ is highly unlikely.

Even with access to a modest number Taylor coefficients, it is possible extract valuable information regarding the equation of state that goes beyond approximating it as a truncated Taylor expansion. The starting point is that, for any temperature, the equation of state possess singularities for complex values of $\mu_{B}^{2}$, the Lee-Yang (LY) edge singularities. With certain assumptions regarding location of the nearest singularity to origin, $\mu_{B}=0$, one could estimate the location of the singularity from the Taylor coefficients. Furthermore, in the vicinity of a critical point, the LY singularities obey a particular scaling form. The key idea is to map the extrapolated singularity from the Taylor coefficients at different temperatures, referred to as the “Lee-Yang trajectory”, and compare it with the scaling form predicted by critical scaling. This approach allows the information extracted from the Taylor coefficients to extend beyond the region constrained by the truncated expansion. The idea of extracting the location of LY singularities from the truncated Taylor expansion goes back to work of Fisher Fisher:1974series . The significance of the LY trajectory in the context of the QCD phase diagram and the critical point has been pointed out in Halasz_1997 ; Ejiri:2005ts ; Stephanov:2006dn . Recently, based on these observations, quantitative estimations of the LY trajectory were made by extracting the radius of convergence or using Padé type resummations Mukherjee:2019eou ; Connelly:2020pno ; Basar:2021hdf ; Basar:2021gyi ; Dimopoulos:2021vrk ; Nicotra:2021ijp ; Singh:2021pog ; Schmidt:2022ogw ; Bollweg:2022rps ; Clarke:2023noy ; Bollweg:2022fqq .

In this paper, I employ an efficient method that improves Padé resummation in order to construct the LY trajectory. I utilize the recent HotQCD results Bollweg:2022rps ; Bollweg:2022fqq for the Taylor coefficients. From the improved resummation, assuming that the trajectory, in part, can be captured by the scaling behavior of the critical point, I extrapolate the location of the critical point and constrain the critical contribution to the equation of state in its vicinity. This is the main result of this paper and is summarized in Table 1 in Section IV. My hope is that these results can be incorporated into the state-of-the-art parametrizations of the equation of state such as in Refs. Parotto:2018pwx ; An:2021wof in to order to assist the experimental effort for the search for the critical point. In order to take into account lattice systematics, I also utilize data from the Wuppertal-Budapest Borsanyi:2018grb .

The rest of the paper is organized as follows. In Sec. II I summarize the physics of LY singularities in the context of the QCD critical point and establish the notation that I use for the rest of the paper. In Sec. III, I build the necessary mathematical machinery that I use to extract the singularities from the Taylor coefficients. My results are presented in Sec. IV. The same analysis is then repeated by using the Wuppertal-Budapest data in Sec. V and compared with the main analysis that uses the HotQCD results. I extensively discuss these results in the final section, Sec. VI.

In their seminal work in 1952, Lee and Yang showed that phase transitions of a thermodynamic system can be understood in terms of the complex singularities of the grand canonical partition function Yang:1952be ; Lee:1952ig . In general, the partition function, $Z(\zeta)$, of a system with finitely many degrees of freedom is a polynomial in fugacity, $\zeta=e^{\mu/T}$, and is nonnegative for $\zeta>0$. Naturally, being a polynomial, it has zeros for complex values of $\zeta$, which, in the thermodynamic limit, coalesce into branch points that emanate from the so called Lee-Yang edge singularities. At a second order phase transition the LY edge singularities pinch the real axis.

I will focus on the LY singularities associated with the three-dimensional Ising model which belongs to the same universality class as QCD. In the Ising model, there are two relevant operators, energy density and spin, whose couplings are the reduced temperature, $r$, and the magnetic field, $h$. The critical point in the phase diagram sits at $r=h=0$. The transition between positive and negative magnetization phases (probed by $h>0$ and $h<0$ respectively) is a smooth crossover for $r>0$ and a first order transition for $r<0$. It is convenient to use the scaling variable $x=hr^{-{\beta\delta}}$ to express the location the LY singularity which is simply along the pure imaginary axis:

$\displaystyle x=\pm ix_{LY}\,.$

(1)

Here $\beta$ and $\delta$ are the usual Ising critical exponents ZinnJustin:2002ru . In this work I use the value computed from conformal bootstrap PhysRevD.86.025022

$\displaystyle{\beta\delta}\approx 1.5631\,.$

(2)

The value of $x_{LY}$ was recently computed via functional renormalization group Rennecke:2022ohx ; Connelly:2020gwa ; Johnson:2022cqv as well as using the Schofield representation Karsch:2023rfb . Based on these works I will use the value:

$\displaystyle x_{LY}=|z_{c}|^{-{\beta\delta}}\approx 0.246$

(3)

The LY singularity is a critical point as well and the equation of state in its vicinity belongs to the same universality class as the $\phi^{3}$ theory with a pure imaginary coupling Fisher:1978pf . The universal, singular contribution to the magnetization around the LY singularity is given by

$\displaystyle m-m_{c}\sim(x\pm ix_{LY})^{\sigma_{LY}}$

(4)

with the critical exponent $\sigma_{LY}\approx 0.074-0.085$ An:2016lni . For a more detailed analysis of the analytical structure of the LY singularities in three and two-dimensional Ising models I refer the reader to Refs. An:2016lni ; An:2017brc and Fonseca:2001dc respectively. In the latter case, the physics of the LY singularity is captured by a non-unitary conformal field theory. In the former case, an analogous analytical construction does not exist at the moment, and one has to rely on numerical results such as the epsilon expansion.

The above observations about the three-dimensional Ising model can now be translated to QCD, given that these two theories are in the same universality class. In the vicinity of the critical point, $(T_{c},\mu_{c})$, the relevant directions in the QCD phase diagram, $T$ and $\mu_{B}$, can be mapped to those of the Ising model, $h$ and $r$, via a linear map Parotto:2018pwx ; Pradeep:2019ccv

$\displaystyle\begin{pmatrix}r\\ h\end{pmatrix}:=\begin{pmatrix}r_{T}&&r_{\mu}\\ h_{T}&&h_{\mu}\end{pmatrix}\begin{pmatrix}T-T_{c}\\ \mu-\mu_{c}\end{pmatrix}.\quad$

(5)

This mapping then leads to the following expression for the trajectory of the LY singularities of QCD in the vicinity of the critical point Stephanov:2006dn :

	$\displaystyle\mu_{LY}(T)\approx\mu_{c}-c_{1}(T-T_{c})\pm ix_{LY}c_{2}(T-T_{c})^{\beta\delta},$
	$\displaystyle\text{where }c_{1}:={{h_{T}}\over{h_{\mu}}}:=\tan\alpha_{1}\quad c_{2}:={{r_{\mu}}^{\beta\delta}\over{h_{\mu}}}\left({{r_{T}}\over{r_{\mu}}}-{{h_{T}}\over{h_{\mu}}}\right)^{{\beta\delta}}.\quad$		(6)

Notice that $c_{1}$ is the slope of the crossover line, whereas $c_{2}$ depends on the relative angle, $\alpha_{2}$, between the $h$ and $r$ axes Parotto:2018pwx ; Pradeep:2019ccv . Remarkably, the trajectory in Eq. (II) depends not only on the location of the critical point, but also on the non-universal mapping parameters. My goal is to reconstruct Eq. (II) and constrain these non-universal values from the Taylor series coefficients of the equation of state.

The main objective is to extract the singular behavior of the equation of state $p(T,\mu_{B})$ from its Taylor series expansion around $\mu_{B}=0$,

$\displaystyle\frac{p(T,\mu_{B})-p(T,0)}{T^{4}}\approx\sum_{n=0}^{N}\frac{\chi_{2n}(T)}{(2n)!}\left(\mu_{B}\over T\right)^{2n}\,.$

(7)

In practice there is only access to a finitely many terms, and the truncated Taylor series is a polynomial which obviously has no singularities. However even with finitely many coefficients, one could extract the singular behavior of the equation of state. Assuming that the nearest singularity to the origin is the LY singularity, following the extended analyticity conjecture of Fonseca and Zamalodchikov Fonseca:2001dc , along with Darboux’s theorem, it is possible to estimate both $|\mu_{LY}|$ via the radius of convergence by using the standard root/ratio tests Bazavov:2017dus ; GIORDANO2021121986 ; Mukherjee:2019eou , as well as $\arg\mu_{LY}$ Basar:2021hdf ; Basar:2021gyi . However, here, the leading singularities are in fact a complex conjugate pair [see Eq. (II)] and consequently the roots/ratios of the coefficients exhibit oscillatory behavior due to interference between the phases of the two singularities which render these estimations numerically challenging.

A better approach is to use Padé resummation, which approximates the original function, $p(\mu_{B}^{2})$²²2To keep the notation compact, I suppress the temperature argument, and switch the argument of $p$ to $\mu_{B}^{2}$., by a rational function,

$\displaystyle{\rm P}[p](\mu_{B}^{2}):=\frac{q_{1}(\mu_{B}^{2})}{q_{2}(\mu_{B}^{2})}$

(8)

where $p_{1}$ and $p_{2}$ are polynomials of order $[N/2]$³³3 In this work I only consider diagonal Padé approximants where the polynomials $p_{1}$ and $p_{2}$ are of the same order., whose coefficients are determined by Taylor expanding ${\rm P_{N/2}}[p](\mu_{B}^{2})$ and matching the coefficients with the original ones in Eq. (7).⁴⁴4Note that even though the coefficients of these polynomials depend on $T$, they are not necessarily smooth functions of $T$. The underlying singularities of $p(\mu_{B}^{2})$ are then represented by accumulation points the poles and zeros of ${\rm P}[p](\mu_{B}^{2})$. Furthermore, if the underlying singularity is a branch point, the poles/zeros lie on curved arcs. Remarkably, the shapes of these arcs can be understood in terms of a two-dimensional electrostatic problem of finding a configuration of charges (poles and zeros) which minimizes an effective capacitance STAHL1997139 ; saff ; Costin:2020pcj . At the same time, when there is a complex conjugate pair of singularities, as in QCD, the arcs that emerge from each singularity coalesce along the real axis. These singularities along the real axis are unphysical, but they are not numerical artifacts; their existence is unavoidable with Padé resummation. To make things worse, their number grows as $N$ increases, limiting the applicability of the Padé resummation to $\mu\lesssim|\mu_{LY}|$.

In order to overcome this inherent shortcoming of Padé resummation, the next step is to pair it with a conformal map. In a way, Padé resummation already introduces its own conformal map by representing the branches with the curved arcs that minimize the effective capacitance. It is possible to further improve convergence properties of Padé, and extend its domain of applicability by using an additional conformal map. The idea is to map the complex $\mu_{B}^{2}$ plane, where the equation of state is expressed in, to a different (preferably compact) region, such as the unit disk, denoted by $\zeta$ via a conformal map,

$\displaystyle\mu_{B}^{2}:=\phi(\zeta)\,.$

(9)

The next step is to expand the equation of state as a series expansion in $\zeta$ and perform the Padé resummation for this expansion:

$\displaystyle{\rm CP}[p](\zeta):=\frac{\tilde{q}_{1}(\zeta)}{\tilde{q}_{2}(\zeta)}\,.$

(10)

where $\tilde{p}$ and $\tilde{q}$ are order $N/2$ polynomials whose coefficients are determined by the Taylor coefficients of $p(\phi(\zeta))$. This resummation will be referred to simply as “conformal Padé”. Different choices of the conformal map, $\phi$, leads to significant improvements over the usual Padé resummation. Moreover, extra information such as the singular behavior of the original function can be “baked in” into the resummation via choosing an appropriate conformal map. Note that, in contrast, the only input of Padé resummation is the set of Taylor coefficients. After performing the Padé resummation in the $\zeta$ plane, the original function is then represented as

$\displaystyle f(\mu_{B}^{2})\underset{\text{conf. Pad\'{e}}}{=}{\rm CP}[p](\zeta)\big{|}_{\zeta=\phi^{-1}(\mu_{B}^{2})}\,.$

(11)

For certain conformal maps, an analytical expression for the inverse function exists, and for others it does not and the inversion has to be computed numerically. Similar to Padé, the original singularities of the function are represented by accumulation of zeros and poles of conformal Padé in $\zeta$ plane. They can be mapped back to $\mu_{B}^{2}$ plane via $\phi(\zeta)$. Note that the computation of the singularities does not require the inverse function. The key point is that for the conformal maps that are considered here, the spurious poles appear outside the unit disk and therefore are not present in the $\mu_{B}^{2}$ plane Costin:2020hwg ; Costin:2021bay . The conformal maps that are used in this work are discussed below.

The first map, named as the“two-cut map”, is defined as

$\displaystyle\phi_{2}(\zeta)=4|\mu^{2}_{LY}|\zeta\left[\frac{\theta/\pi}{(1-\zeta)^{2}}\right]^{\theta/\pi}\left[\frac{1-\theta/\pi}{(1+\zeta)^{2}}\right]^{1-\theta/\pi}\,.$

(12)

It maps the complex plane $\mu_{B}^{2}$ plane with two radial cuts into the unit disk as shown in Fig. 1. The branch points located at $\mu_{B}^{2}:=|\mu_{LY}^{2}|e^{\pm i\theta}$ are mapped to the edge of the unit disk with the angle $\psi_{2}:=2\arcsin(\sqrt{\theta/\pi})$:

$\displaystyle|\mu_{LY}^{2}|e^{\pm i\theta}\rightarrow\zeta_{LY}=e^{\pm i\psi_{2}}$

(13)

and the branch cuts are mapped to the edge of the unit disk as show in Fig. 1. Notably, this map has been used in other physical applications, mostly within the context of extracting the Borel singularities associated with asymptotic series Guida:1998bx ; Rossi:2018 ; Serone:2019szm ; Costin:2020hwg ; Costin:2021bay . As opposed to Padé, conformal Padé with the two-cut map does not generate unphysical poles along the real axis, allowing one to reconstruct the original function beyond the radius of convergence. Furthermore it generally gives a better estimate for the location of the branch singularities compared to Padé with the same number of Taylor coefficients Costin:2020pcj . One shortcoming of the two-cut map, however, is that its applicability is limited to the first Riemann sheet. Namely, even though it gives a good estimate for the location of the branch singularities, it does not allow one to actually go through the branch cuts. This can be seen from the fact that it maps the first Riemann sheet, bounded with the radial branch cuts, inside the whole unit disk.

An improvement over the two-cut map is the “uniformizing map” defined as

$\displaystyle\phi_{u}(\tau)=|\mu^{2}_{LY}|\left(e^{-i\theta}+2i\sin\theta\lambda\left(\tau\right)\right)\,.$

(14)

Here $\lambda(\tau)=\theta_{2}^{4}(\tau)/\theta_{3}^{4}(\tau)$ is the modular lambda function with $\theta_{2}(\tau)=\sum_{n=\infty}^{\infty}e^{2\pi i\tau(n+1/2)^{2}}$ and $\theta_{3}(\tau)=\sum_{n=\infty}^{\infty}e^{2\pi i\tau n^{2}}$ being the usual Jacobi elliptic functions defined in the upper half plane ${\rm Im}\tau>0$. I further map the upper half-plane into the unit disk via the following Mobius transformation:

$\displaystyle\tau\rightarrow\tau(\zeta)=i\frac{\mathbb{K}\left(\frac{1}{2}-\frac{i}{2}\cot\theta\right)+\mathbb{K}\left(\frac{1}{2}+\frac{i}{2}\cot\theta\right)i\zeta}{\mathbb{K}\left(\frac{1}{2}+\frac{i}{2}\cot\theta\right)-\mathbb{K}\left(\frac{1}{2}-\frac{i}{2}\cot\theta\right)i\zeta}$

(15)

where $\mathbb{K}(m)$ is the complete elliptic integral of the first kind:

$\displaystyle\mathbb{K}(m)=\int_{0}^{\pi/2}\frac{d\theta}{\sqrt{1-m\sin^{2}\theta}}\,.$

(16)

This map “uniformizes” the multi-sheeted $\mu_{B}^{2}$ plane of the original function by mapping the entire multi-sheeted domain into a simply connected domain. A path that crosses the branch cuts in the $\mu_{B}^{2}$ plane is represented by a smooth curve in the mapped plane. How this works can be briefly explained in two steps.

The first step is to map the $\mu_{B}^{2}$ plane with vertical branches into the fundamental domain via Eq. (14). This map is shown in Fig. 2 (center). Notice that the entire $\mu_{B}^{2}$ is mapped into a subset of the upper half-plane, the fundamental domain. The remaining half-circular “gaps” are filled by Schwartz reflection of the fundamental domain. These reflections are transformations built out of the elementary modular transformations, $S:\tau-\rightarrow-1/\tau$ and $T:\tau\rightarrow\tau+1$. Each Schwartz reflection fills in a portion of the semi-circular gaps in the upper half-plane, a process that can in principle continued ad-infinitum, asymptotically filling the whole gap. These regions are the images of the higher Riemann sheets with each Schwartz reflection corresponding to a particular sheet. Therefore, moving across different sheets in the $\mu_{B}^{2}$ plane is simply represented by moving smoothly in the $\tau$ plane. Each Schwarz reflection generates a progressively smaller region in the $\tau$ plane. Therefore, in order to resolve the higher sheets one needs progressively more precise information, meaning more Taylor coefficients.

The second step is to further map the fundamental domain into the unit disk via the Mobius transformation given in Eq. (15). As seen in Fig. 2 (right), the branch points are mapped to the edge of the unit disk. Similar to the two-cut case, the angle is given by

$\displaystyle|\mu_{LY}^{2}|e^{\pm i\theta}\rightarrow\zeta_{LY}=i\frac{\mathbb{K}\left(\frac{1}{2}\mp\frac{i}{2}\cot\theta\right)}{\mathbb{K}\left(\frac{1}{2}\pm\frac{i}{2}\cot\theta\right)}:=e^{\pm i\psi_{u}}$

(17)

The branches, however are mapped into certain curves that for a boundary of a compact domain within the unit circle. This compact, skewed diamond shaped region bounded by these curves is image of the first sheet. The remaining gaps in the unit circle are filled with Mobius transformations of the Schwartz reflected regions in the modular $\tau$ plane. These are the images of higher sheets. Consequently the entire multi sheeted $\mu_{B}^{2}$ space is mapped into the unit circle and crossing between is represented by smooth paths within the unit circle.

This way, one can reconstruct the original function not only beyond its radius of convergence, but also in the higher Riemann sheets, using only the Taylor coefficients of the expansion around the origin in the first Riemann sheet. In the context of this work, the higher Riemann sheets encode the scaling equation of state in the first order phase transition region, $T<T_{c}$. In Ref. Basar:2021gyi I demonstrated this by reconstructing the equation of state of the mean field Ising model in the first-order transition region ($T<T_{c}$) by using the Taylor expansion coefficients of the high-temperature ($T>T_{c}$) equation of state obtained in the first Riemann sheet. In addition to reconstructing the equation of state in the higher sheets, the uniformizing map also provides a better approximation to the function compared to Padé and two-cut conformal Padé.

This work will focus on the first sheet and the uniformizing map will be uses to get a better estimate for the LY singularities and the equation of state. This is because having access to only four Taylor coefficients with sizable statistical uncertainties makes the analytical continuation to higher sheets extremely challenging. As opposed to the two-cut map, there is an analytical expression for the inverse of the uniformizing map, given as,

$\displaystyle\phi_{uni}^{-1}(\mu_{B}^{2})=-i\frac{\mathbb{K}\left(\kappa_{+}\right)\mathbb{K}\left(\kappa_{-}+\frac{i}{2\sin\theta}\frac{\mu^{2}_{B}}{|\mu^{2}_{LY}|}\right)-\mathbb{K}\left(\kappa_{-}\right)\mathbb{K}\left(\kappa_{+}-\frac{i}{2\sin\theta}\frac{\mu^{2}_{B}}{|\mu^{2}_{LY}|}\right)}{\mathbb{K}\left(\kappa_{-}\right)\mathbb{K}\left(\kappa_{-}+\frac{i}{2\sin\theta}\frac{\mu^{2}_{B}}{|\mu^{2}_{LY}|}\right)+\mathbb{K}\left(\kappa_{+}\right)\mathbb{K}\left(\kappa_{+}-\frac{i}{2\sin\theta}\frac{\mu^{2}_{B}}{|\mu^{2}_{LY}|}\right)}$

(18)

where $\kappa_{\pm}:=(1\pm i\cot\theta)/2$ and $\theta=|\arg\mu^{2}_{LY}|$.

Notice that both the two-cut map and the uniformizing map depend on the location of the singularity, $\mu_{LY}$ which, in fact, is what one wants to extract from these resummation techniques. This seemingly paradoxical situation can be overcome by a simple procedure: 1) guess the location from ordinary Padé and use this initial guess as the value of $\mu_{LY}$ in the conformal map. 2) Extract the poles from conformal Padé in $\zeta$ plane whose images in the $\mu_{B}^{2}$ plane gives a refined estimate for $\mu_{LY}$. 3) Iterate the same procedure by using this refined estimate for the value of $\mu_{LY}$ in the conformal map. I observed that this iteration converges to a value which constitutes the final estimate for $\mu_{LY}$. This iterative process is used for different temperatures to construct the LY trajectory.

Before presenting the results it is worth to comment on an alternative resummation technique. The equation of state in the vicinity of the LY singularities has a singular contribution, but this singular contribution actually vanishes at the singularity for the pressure and generates a cusp for the density. However, a true divergence occurs for the susceptibility. For real values of $\mu_{B}$ and $T\gtrsim T_{c}$, the susceptibility does not diverge but sharply peaks around ${\rm Re}\mu_{LY}$. For this reason, sometimes performing a Padé resummation for the susceptibility instead of the pressure can give a better estimate of the underlying singularities as well as the function itself. Therefore I also performed resummation for the susceptibility:

$\displaystyle\chi(T,\mu_{B})={\partial^{2}p\over\partial\mu_{B}^{2}}\approx\sum_{n=0}^{N-1}\frac{\chi_{2n+2}(T)}{(2n)!}\left(\mu_{B}\over T\right)^{2n}\,.$

(19)

Using pressure or susceptibility in the resummation just amounts to reshuffling the same Taylor coefficients. However this reshuffling, especially when $N$ is small, does make a difference. For example in Refs. Basar:2021hdf ; Basar:2021gyi which studied the Gross-Neveu and the Chiral Random Matrix models, using the susceptibility led to much more accurate results than the pressure. For QCD, this situation is more complicated due to the statistical uncertainties in the coefficients as will be discussed further in the next section.

The results for estimates for the Lee-Yang singularities as a function of temperature, obtained by the resummation methods are presented in this section. The calculations for the locations of the critical point and the non-universal mapping parameters in Eq. (5) based on these estimates are presented as well, along with results for the susceptibility as a function of chemical potential for a few different temperatures.

I use the recent results for the Taylor coefficients up to ${\cal O}(\mu_{B}^{8})$ computed by the HotQCD collaboration Bollweg:2022rps for $\mu_{S}=\mu_{Q}=0$. The continuum extrapolation for the first two terms, $\chi_{2}$, $\chi_{4}$ are used. Unfortunately, continuum extrapolations for the remaining two terms, $\chi_{6}$ and $\chi_{8}$, were not available at the time this work was completed. I use the spline extrapolation for the $N_{\tau}=8$ data given in Ref. Bollweg:2022rps . For sake of completeness the HotQCD data are plotted in Fig. 3.

With four terms in the Taylor expansion, the diagonal Padé resummation is a ratio of two quadratic polynomials in $\mu_{B}^{2}$. Analytical expressions for the poles and zeros of Padé Bollweg:2022rps as well as conformal Padé can be obtained; however, their functional forms do not play a central role in this analysis and to keep the discussion concise they will not be included here.

In order to take into account the statistical uncertainties, I sampled an ensemble of coefficients from a Gaussian distribution and used diagonal (2,2) Padé resummation which has a complex conjugate pair of poles. Since there are no other poles to form an accumulation point, these poles are used as estimators for the LY singularities, $|\mu_{LY}^{2}|e^{\pm i\theta}$. The LY trajectory constructed by repeating these steps for different temperatures is shown in Fig. 7. I also found that the zeros did not follow any meaningful pattern. This is likely due to $N$ being relatively small.

With the conformal maps, there is one more step in extracting the singularity. As mentioned in the previous section, the conformal maps explicitly depend on the location of the singularity. I first used the Padé estimate for the singularity in the conformal map and then followed the iterative procedure described in Section III to refine the estimate for $|\mu_{LY}^{2}|e^{\pm i\theta}$. The results of this procedure are shown in Figs. 4 and Fig. 5.

Figure 4 shows the trajectory of the iteration for different set of Taylor coefficients sampled from a Gaussian ensemble. Each line with color represents the iteration obtained with fixed set of Taylor coefficients sampled from a Gaussian ensemble, in the $\zeta$ (left) and $\mu_{B}$ (right) planes. The final estimate for each trajectory is denoted by a small solid disk. Each disk and the iteration curve are color coded. The pale blue disk in the center figure represents $1\sigma$ uncertainty region. For comparison, the $1\sigma$ uncertainty region for the ordinary Padé resummation with the same temperature and ensemble of Taylor coefficients is also included. The right figure shows the evolution of the LY trajectory constructed by averaging over the Gaussian ensemble with the iteration. Different opacities denote different steps in the iteration, darker being later. Recall that the image of the true singularity of the equation of state in the $\zeta$ plane is along the unit disk (see Figs. 1 and 2). Remarkably the iteration indeed converges to the edge of the unit disk. Fig. 5 shows the same trajectories for the uniformizing map. The same ensemble of Taylor coefficients is used both for ordinary Padé and the two conformal maps. In Fig. 6 the result of the iteration for a single trajectory is shown. For both conformal maps, the initial point is the same Padé pole obtained from the Gaussian ensemble. The left figure show the convergence towards the edge of the unit disk.

The LY trajectories constructed from different resummations are shown in Fig. 7. Each point for conformal Padé is obtained after 100 iterations as described above. The error bars represent the $1\sigma$ uncertainties inherited from the statistical uncertainties in the Taylor coefficients.

The next goal is to extract the location of the critical point from the LY trajectory. The critical chemical potential, $\mu_{c}$ is the point where ${\rm Im}\mu_{LY}$ vanishes which is clearly beyond the available data as ${\rm Im}\mu_{LY}>0$ for the temperature range in hand. At the same time, the fact that ${\rm Im}\mu_{LY}$ is decreasing with decreasing temperature is suggestive that if there is a critical point, it lies at $T_{c}<135$ MeV. By extrapolating the trajectory to the point where ${\rm Im}\mu_{LY}=0$ I estimate $\mu_{c}$ and $T_{c}$ using the following fits for the extrapolation

	$\displaystyle{\rm Re}\mu_{LY}(T)$	$\displaystyle=$	$\displaystyle a_{0}+a_{1}(T-T_{c})+a_{2}(T-T_{c})^{2}$		(20)
	$\displaystyle{\rm Im}\mu_{LY}(T)$	$\displaystyle=$	$\displaystyle a(T-T_{c})^{\beta\delta}$		(21)

whose form is motivated by the scaling form given in Eq. (II). The results for the best-fits for different resummations are shown in Fig. 7 as solid lines. In these fits, I used the first 20 terms in the trajectory with $134.5\leq T\leq 144$. Finally, the location of the critical point, as well as the non-universal mapping parameters, the slope of the crossover line at the critical point, $\alpha_{1}$, and $c_{2}$ as given in. Eq. (II), extrapolated from these fits, are listed in Table 1.

The final set of results concerns the Padé and conformal Padé resummations by using the Taylor coefficients of the susceptibility given in Eq. (19). A particularly interesting result obtained from this resummation is the susceptibility as a function of chemical potential, shown in Fig. 8. The band denotes the $1\sigma$ uncertainty as before. Finally Fig. 9 shows the singularities obtained this way. For comparison, uniformizing map results for the singularities obtained from the expansion of pressure are also shown in the same figure (in purple).

Notice that using the coefficients of susceptibility reproduces qualitatively the expected behavior of $\chi$ as a function $\mu_{B}$, namely exhibiting a peak at some nonzero value of $\mu_{B}/T$, albeit with sizable statistical uncertainty. At the same time, the singularities extracted this way are less reliable compared to using the coefficients of pressure. I found that the conformal Padé singularities in the $\zeta$ plane were much further away from the edge of the unit disk than those obtained from the coefficient of pressure. Furthermore the statistical errors are very large, especially for $T\lesssim 145$ MeV. This is ultimately because the Padé and conformal Padé singularities depend rather sensitively on $\chi_{6}$, which changes sign around $140$ MeV, making the relative error quite large. Due to the difference in the functional dependence of the Padé singularities in the coefficients, this is not the case for the expansion of pressure. Therefore I conclude that using coefficients of pressure in the estimates for the location of singularities, using pressure as opposed to susceptibility is more reliable. However, using the susceptibility seems to reproduce the qualitatively expected form of the equation of state better. A more detailed analysis of this observation is left for future work.

In addition to the statistical uncertainties that are present in any stochastic computation, and the uncertainties in extracting the singularities due to finite number of Taylor coefficients, there are also systematic uncertainties inherited from the lattice computations.The statistical uncertainties are already quantified in my analysis. The uncertainties due to the finite number of Taylor coefficients are more difficult to quantify but they can be gauged by analyzing how close the images of the extrapolated singularities are to the edge of the unit disk in the conformal plane, as explained in the previous section. At the same time, the systematic uncertainties are much more challenging to quantify. In order to shed light on these systematic uncertainties, in this section the same analysis described above is repeated by using a different data set of Taylor coefficients computed by the Wuppertal-Budapest (WB) Collaboration Borsanyi:2018grb . The WB simulations are done at imaginary chemical potential with the same physical volume as HotQCD, in particular on $48^{3}\times 12$ (WB) and $32^{3}\times 8$ (HotQCD) lattices.

Even though there are four coefficients in the Taylor expansion, via rescaling $p$ and $\mu_{B}^{2}$, one can show that the Padé poles can be expressed as a function of two variables, $\chi_{2}\chi_{6}/\chi_{4}^{2}$ and $\chi_{8}\chi_{2}^{2}/\chi_{4}^{3}$, up to an overall factor of $\chi_{2}/\chi_{4}$. By using the same normalization in Ref. Bollweg:2022rps , these parameters computed by the two lattice collaborations are shown in Fig. 10 for the temperature range that this work focuses on. Notably the coefficients are quantitatively different (whereas the overall factor $\chi_{2}/\chi_{4}$ varies by a few percent between WB and HotQCD data), signaling the importance of the systematics in any kind of resummation framework that is based on Taylor coefficients.

Based on the observation that the original Taylor coefficients computed by different collaborations are different, one might expect that the LY singularities obtained from Padé resummations would also differ substantially. However, I found that this is not entirely the case. Before doing so, it is worth commenting on the statistical uncertainties first. Even though the overall statistical errors in the WB results are smaller in magnitude, the overall signal to noise ratio for the LY singularities is actually higher compared to HotQCD results. This is because $\chi_{8}$ is closer to zero. As an illustration, the trajectories of the iteration procedure described in Section IV for two different conformal maps are shown in Fig. 11. As explained in Section IV, each line represents a trajectory whose initial point is the Padé singularity associated with the Taylor coefficients sampled from a Gaussian ensemble, and the final point is obtained after 100 iterations of the conformal map. For comparison, the results obtained from the HotQCD Collaboration are shown in red. The qualitative behavior of the trajectories is similar; however, the WB results have larger statistical uncertainties, most likely due to the small mean value of $\chi_{8}$ as mentioned.

The LY trajectory for different temperatures is shown in Fig. 12. Qualitatively, it looks similar to the one obtained from the HotQCD data (see Fig. 7). In particular, its imaginary part decreases with decreasing temperature. This behavior is consistent with what one would expect if there is a critical point, since at a critical point the imaginary part vanishes. Although the error bars are too large to determine the exact functional form of the trajectory with good accuracy, using the functional form assuming the ${\mathbb{Z}}_{2}$ scaling form given in Eq. (II) leads to reasonable values for the critical temperature listed in Table 2. These results are consistent with those obtained from the HotQCD data (see Table 1) albeit with larger statistical uncertainties. At the same time, even though the real part of the LY trajectory is qualitatively similar to the one obtained from HotQCD data, the systematic difference is more pronounced. Furthermore, due to the larger statistical uncertainties combined with those that already enter the estimation of $T_{c}$, it is not possible to meaningfully extract $\mu_{c}$ with the available data.

In conclusion, the most robust feature of the Lee-Yang trajectory constructed from two separate lattice computations, HotQCD and Wuppertal-Budapest, is that its imaginary part decreases temperature, consistent with the existence of a critical point. However, it does not vanish within the available temperature range (within 1$\sigma$ uncertainty); signaling that, if exists, the critical point is likely at $T_{c}<130$ MeV. Determining the quantitative features of the LY trajectory is more challenging with the current data, however, the extrapolation of using the ${\mathbb{Z}}_{2}$ scaling ansatz leads to a robust value of $T_{c}\sim 100$ MeV, both for the WB and HotQCD data. Due to the statistical and systematic uncertainties, estimating the value of the critical chemical potential is even more challenging. At the same time, the current lattice data is consistent with the estimation $\mu_{c}\sim 600$ MeV.

The task of extracting the location of the Lee-Yang singularities from the presently available lattice QCD data is a challenging one. Currently, there is only access to four coefficients with sizable statistical uncertainties. In this paper, I utilized various conformal maps to improve the accuracy of the usual Padé resummation. By choosing appropriately designed conformal maps, I incorporated further analytical information regarding the equation of state, in addition to the Taylor coefficients; namely the closest singularities must be complex conjugate pair, i.e. the equation of state must be defined on a two-cut Riemann surface. I then conformally mapped this surface into the unit disk and extracted the Padé singularities there. Performing Padé resummation on a compact space gave a better handle in pinning down the location of the true singularity.

In order to take into account the statistical uncertainties in the Taylor coefficients, they are sampled from a Gaussian ensemble whose variance is matched by the lattice data. The error due to the small number of Taylor coefficients, on the other hand, is harder to quantify. For large number of Taylor coefficients, there is a scaling relation between the magnitude of noise in the Taylor coefficients and the expansion order before (conformal) Padé resummation breaks down Costin:2022hgc . However it is an asymptotic result whose applicability to four terms is unclear.

In order to refine the estimation for the location of the singularity a novel iterative tool is used. Remarkably for both conformal maps, the iteration brought the images of the singularities extracted from conformal Padé closer to the edge of the unit disk where the real singularity lies, as seen in Figs. 4 and 5. This observation increases the confidence in the results of conformal Padé. Note that the uniformizing map is exponentially sensitive to the location of the singularity. Therefore it is not surprising that the poles in $\zeta$ plane for the uniformizing map are further away from the edge of the unit disk since neither the number Taylor coefficients is large enough nor their precision high enough to resolve singularity with exponential accuracy. It is also worth pointing out that the final results for the LY trajectory obtained from both conformal maps agree with each other and systematically differ from those from ordinary Padé, with the conformal Padé results being slightly less sensitive to the statistical uncertainties compared to Padé which can be seen in Figs. 4, 5 and II. Furthermore, as explained in e.g. Ref. Costin:2020pcj conformal maps provide a better estimate for the singularity compared to ordinary Padé. The various test cases where the ordinary and conformal Padé predictions can be compared with the exact results also confirm this observation Basar:2021hdf ; Basar:2021gyi . For these reasons, I think that conformal Padé results are more trustworthy than those of Padé.

Then, from the LY trajectory constructed via the resummations, I extrapolated the location of the critical point, as well as constrained values of the non-universal mapping parameters in its vicinity. These results are given in Table 1. This is the central result of this paper. Notice that several significant figures are used in these extrapolations. This is not to claim such precision in the final estimate, but to illustrate the quantitative differences between different resummations. The fact that they all lie in the same overall region is encouraging. Furthermore, the results obtained from two independent lattice computations indicate a robust trend in ${\rm Im}\mu_{L}{LY}$ consistent with $T_{c}\approx 100$ MeV. My results also agree with other results computed via similar Padé type resummations Clarke:2023noy , as well as other methods such as functional renormalization group, ($T_{c},\mu_{C}$) = (107, 635) MeV Fu:2019hdw and truncated Dyson Schwinger equations, ($T_{c},\mu_{C}$) = (117, 600) MeV Bernhardt:2021iql ; Fischer21 . The estimate based on the best fit, $\mu_{c}/T_{c}\approx 6$, is also consistent with the constraints from the recent lattice QCD data which strongly disfavor the existence of a critical point for $\mu_{c}/T_{c}\lesssim 3$ Borsanyi:2020fev .

Because the extrapolated $T_{c}$ ($\sim$ 100 MeV) is fairly lower than the minimum available temperature (135 MeV) I had to rely on best fits for the estimates for the critical point. Especially ${\rm Re}\mu_{LY}$, from which $\mu_{c}$ is extrapolated, depends sensitively on the quadratic fit. Combined with the current statistical and systematic uncertainties in the Taylor coefficients, such an extrapolation introduces a sizable uncertainty in the estimate for $\mu_{c}$. Having access to lattice data for lower temperatures in the future would reduce the sensitive reliance on the extrapolation and therefore significantly improve the accuracy of the estimation of the location of the critical point.

Another important point worth discussing is that I assumed the equation of state obeys the scaling obtained from the universality class of the “$Z_{2}$ critical point” for part of the available temperature range for the LY trajectory in order to extrapolate the data to $\mu_{c}$ and $T_{c}$. In this scaling, the relevant axes mapped to the Ising parameters $h$ and $r$, are identified as $\mu_{B}$ and $T$. At the moment, it is unclear whether this scaling is valid at these temperatures. It is widely believed that around the pseudo-critical temperature, the closest singularity is associated with the “$O(4)$ critical point” (more precisely a line of second order points) that is located at $m_{u,d}=0$ and belongs to the $O(4)$ universality class. Therefore the LY trajectory should obey the “$O(4)$ scaling” which has a different form than the $Z_{2}$ scaling used in Eq. (II). This is because the $O(4)$ scaling identifies one of relevant axes with the quark mass. Even though the agreement between the data and the form predicted by the $Z_{2}$ scaling seems suggestive that, at least for $T\in(130-145)$ MeV range, the $Z_{2}$ scaling is valid, this point requires further investigation. This issue will be addressed in a forthcoming publication.

Finally, a related issue concerns the shape of the LY trajectory. In various models that exhibit similar critical phenomena to the conjectured QCD phase diagram, the LY trajectory can be exactly calculated. Among them are the Gross-Neveu model Basar:2021hdf , quark meson model Mukherjee:2021tyg and the Chiral Random Matrix model Stephanov:2006dn ; Basar:2021gyi where ${\rm Re}\mu_{LY}(T)$ and ${\rm Im}\mu_{LY}(T)$ are monotonically decreasing and increasing functions of $T$ respectively, for $T>T_{c}$. As seen in Fig. 7, the QCD data indicates that ${\rm Re}\mu_{LY}$ is not monotonically decreasing in $T$. For $T\gtrsim 145$ MeV, it is actually increasing. If this behavior is indeed representative of the physical LY trajectory and not an artifact of the resummations due to a combination of noise and small number of Taylor coefficients, it asks for further investigation. In fact, at higher temperatures, the shape of the LY trajectory is expected to be controlled by yet another singular point; the Roberge-Weiss singularity Roberge:1986mm and for some $T=T_{RW}$, the trajectory should pass through the point(s) ${\rm Re}\mu_{LY}(T_{RW})=0$, ${\rm Im}\mu_{LY}(T_{RW})=\pm\pi$. A recent estimate for the Roberge-Weiss temperature based on multi-point Padé approximants is $T_{RW}\approx 170$ MeV Schmidt:2022ogw . This means that, if the non-monotonic behavior is correct, then for some $T\in(160-170)$ MeV , ${\rm Re}\mu_{LY}(T)$ has to peak and then go down to zero. At the same time, because $\chi_{8}$ crosses zero around $T\approx 170$ MeV, the Padé singularities are too noisy around these temperatures to lead to any meaningful estimation of $\mu_{LY}$. Note that the Roberge-Weiss point is related to confinement/deconfinement transition and therefore is not present in any of the aforementioned models. This discussion is left for future work as well.

I thank M. Stephanov and G. Dunne for fruitful discussions. The author is supported by the National Science Foundation CAREER Award No. PHY-2143149.