Universality, Lee-Yang singularities and series expansions

We introduce a new way of reconstructing the equation of state of a thermodynamic system near a second order critical point from a finite set of Taylor coefficients computed away from the critical point. We focus on the Ising universality class (${\mathbb Z}_2$ symmetry) and show that in the crossover region of the phase diagram it is possible to efficiently extract the location of the nearest thermodynamic singularity, the Lee-Yang edge singularity, from which one can (i) determine the location of the critical point, (ii) constrain the non-universal parameters that maps the equation of state to that of the Ising model in the scaling regime, and (iii) numerically evaluate the equation of state in the vicinity of the critical point. This is done by using a combination of Pad\'e resummation and conformal maps. We explicitly demonstrate these ideas in the celebrated Gross-Neveu model.


INTRODUCTION
In the vicinity of a critical point, the correlation length of a thermodynamic system grows, the underlying microscopic properties become irrelevant, and consequently many different substances exhibit the same behavior. This powerful notion of universality allows one to classify second order phase transitions based on the underlying symmetries of the system without detailed knowledge of the microscopic dynamics. Some famous examples are the liquid-gas and ferromagnetic transitions, both of which possess critical points that are in the (static) universality class of the 3d Ising model. Universality identifies the singular behavior of the equation of state, say pressure as a function of temperature and chemical potential p(T, µ), with the Ising equation of state as a function of the Ising variables, i.e. the reduced temperature, r, and the magnetic field, h, near the critical point. However neither the location of the critical point, nor the relation between (T, µ) and (r, h) are universal, and for a quantitative description of the critical phenomena they have to be computed from the microscopic dynamics of the theory of interest. This is in many cases a mathematically intractable task, and often the solution can only be obtained as a series expansion at a point away from the critical point with a finite number of terms [1].
In this work we consider the following general problem: given a truncated local series expansion of the equation of state in some parameter such as µ, obtained away from the critical point, what can we say about the critical behavior of the system? More precisely, how can we determine whether a critical point exists, and if it does how can we reconstruct the singular behavior of the equation of state near the critical point from the truncated expansion obtained away from it? We show that it is possible to obtain a surprisingly large amount of information about the critical behavior of the system from the series coefficients, even if we have access to a modest number of them.
A major motivation for this work is the search for the conjectured critical point in the phase diagram of Quantum Chromo Dynamics (QCD) [2] which is one of the major outstanding problems in nuclear physics both theoretically and experimentally [3]. The theoretical approaches are severely limited by the sign problem that prevents first-principle lattice computations at nonzero baryon chemical potential, µ B and one of the methods to deal with this obstacle is to expand the equation of state around µ B = 0 and compute the Taylor coefficients on the lattice without a sign problem (see [4] for a recent review). Another motivation is to understand the properties of strongly interacting fermions, such as unitary Fermi gases. where our theoretical knowledge of the equation of state is typically limited to the first few terms in the virial expansion [5].

LEE-YANG EDGE SINGULARITIES
Before detailing our method we briefly summarize the Lee-Yang edge singularities that play a crucial role in our analysis. In their seminal work phase transitions [6,7] Lee and Yang showed that the thermodynamic properties of a system is encoded in the distribution of the zeroes of the partition function Z(ζ) as a function of fugacity, ζ = e µ/T . For a finite system the partition function is a positive polynomial for ζ ≥ 0. However in general it has zeroes for complex values of µ and T . In the thermodynamic limit the zeroes coalesce into branch cuts emanating from the so called Lee-Yang (LY) edge singularities. When the LY singularities pinch the real axis, the system exhibits a second order phase transition. Likewise the branch cut associated with a LY singu-arXiv:2105.08080v2 [hep-th] 15 Jun 2021 larity crosses the real line when there is a first order phase transition.
The LY singularities are critical points themselves and have their own universality class. For example, for the Ising model, and O(N ) models in general, they are characterized by the φ 3 theory with a pure imaginary coupling [8]. For these models the analytical structure of the LY singularities has been studied [9,10]. In the complex plane of the scaling variable, x = hr −βδ 1 , the singularities are located in the pure imaginary axis, x = ±ix LY where x LY ∈ R has recently been calculated via the functional renormalization group method [12]. The magnetization around the LY singularity behaves as m − m c ∼ (x ± ix LY ) σ with the critical exponent σ ≈ 0.074 − 0.085 for d = 3 [13]. The LY singularities in the context of QCD critical point has been discussed in, for example, [14][15][16][17][18][19][20].
Consider the equation of state of a theory near a critical point, (T c , µ c ), in the Ising universality class. From universality we can relate it to the Ising model via a linear map [21,22] This relation then leads to the following expression for the trajectory of the LY singularities [16]: Notice that c 1 is the slope of the crossover line, whereas c 2 depends on the relative angle between the h and r axes [21,22]. Therefore the trajectory in Eq. (2) depends on not only the location of the critical point, but also on the non-universal mapping parameters. We now explain how to optimally construct Eq. (2) from a series expansion.

THE METHOD
Consider a thermodynamic function, f (T, µ), (pressure, density, susceptibility etc.) given a series expansion with finite number of terms: f (T, µ) ∼ N n=0 f n (T )µ 2n . Our goal is to extract as much information from it as possible, especially about its singular behavior near a critical point if there is one. Obviously, in this form what we have is a polynomial which does not exhibit any singular behavior. 1 β and δ are the standard Ising critical exponents [11].
In principle from a ratio test it is possible determine the radius of convergence which would indicate the location of the closest singularity to origin. However when the nearest singularities are a complex conjugate pair of LY singularities (which is the case for a smooth crossover) the ratios of series coefficients do not converge monotonically but rather have an oscillating envelope as a result of Darboux's theorem 2 . This makes numerically extracting the singularity from the ratios challenging. Alternatively, the singular behavior of f (T, µ) can be approximately constructed by a Padé resummation, P N/2 [f ](µ 2 ) := p(µ 2 )/q(µ 2 ) where p and q are polynomials of order N/2 whose the coefficients are determined by expanding P N/2 [f ](µ 2 ) and identifying the coefficients with the Taylor coefficients. The singularities of f , typically branch points, are approximated by the poles and zeroes of the Padé approximant. However as we will demonstrate later, Padé resummation has a known shortcoming; it creates spurious singularities which limits its range of applicability. In addition, it is also not the most optimal approximation scheme and can be dramatically improved by pairing with a conformal map [25]. The key idea is to apply an appropriately chosen conformal map, µ 2 := φ(z), and do the Padé reummation in z: CP N/2 [f ](z) :=p(z)/q(z) wherep andq are order N/2 polynomials whose coeffcients are determined by the Taylor coefficients of f (φ(z)). The singularities of CP N/2 [f ](z) are mapped to the complex µ 2 plane as via φ(z). For the remainder of the letter we shall refer to this method of extracting singularities simply as "conformal Padé". Unlike Padé, Conformal Padé does not suffer from the spurious pole problems. It also provides an optimal approximation to the original function for a wide range of functions as proven in [25]. Our strategy is to first construct the trajectory (2) by extracting the µ LY (T ) via conformal Padé for a sequence of temperatures and then to obtain the location of the critical point as well as the coefficients c 1 , c 2 from it. This is possible since the critical exponents are fixed by universality and, as mentioned, x LY known [12].
Conformal Padé methods are typically used to reconstruct Borel plane singularities in resumming asymptotic series, such as the expansion [26] or perturbation series for relativistic [27] and nonrelativistic [28] systems. Here we take a different approach and apply it to a convergent series to di-rectly extract its singular behavior near the critical point. Our input, the series expansion of the equation of state, does not have to come from perturbation theory or the expansion. Our approach is similar to that of [1] but with a key differences: we focus on the complex LY singularities and conformal maps play a crucial role in reconstructing the equation of state.

THE GROSS-NEVEU MODEL
To concretely demonstrate these ideas we focus on the celebrated Gross-Neveu (GN) model [29] which is a four-fermion theory with the action where ψ is a Dirac fermion with N f flavors. It exhibits some of the key features of QCD, such as asymptotic freedom, chiral symmetry breaking and dimensional transmutation. Notably it was also shown that the GN model gives a reasonably good description of the first order phase transition in doped polyacetlyne [30]. The theory has a discrete (Z 2 ) chiral symmetry, ψ → γ 5 ψ, for m q = 0. We will work in the large N f limit N f → ∞ with g 2 N f =fixed, where the fluctuations are suppressed and the mean field solution is exact. The exact phase diagram of the GN model is known [31,32] and has a rich structure such as spatially inhomogeneous kink/anti-kink crystals at high densities [32] and an exactly soluble, all-orders Ginzburg Landau expansion [33]. However, in order to keep the discussion simple, we shall assume that the translational symmetry remains unbroken in this work. Furthermore we focus on the crossover region of the phase diagram which is not affected by the existence of inhomogeneous phases. Thermodynamics of the model follows from minimizing the grand potential with respect to φ which determines the fermion mass, M as a function of T, µ and γ. The parameter γ ∝ m q is a measure of explicit chiral symmetry breaking and it vanishes in the chiral limit [31]. In our analysis we will work with a fixed, nonzero value of γ.  3 . When the chiral symmetry is explicitly broken, for a fixed γ = 0, the transition is a smooth crossover for hight T that ends at a critical point, (T c , µ c ), and turns into a first order transition. Since the theory has a Z 2 chiral symmetry the second order transition is in the same static universality class as the mean field Ising model. The trajectory of the LY singularities, µ LY (T ), is determined by the condition In the crossover region, T > T c , this condition leads to a pair of complex solutions µ LY (T ) = Reµ LY (T ) ± iImµ LY (T ) which coalesce and pinch the real axis at the ordinary critical point as expected: µ LY (T c ) = µ c . In the vicinity of the critical point µ LY (T ) takes scaling form given in Eq. (2) with the mean field exponents, β = 1/2, δ = 3: where in the mean field limit x LY = 2/3 √ 3. In the next section we compute T c , µ c , c 1 , c 2 directly from a truncated series expansion of the equation of state and compare these results with the exact solution obtained by numerically solving Eq. (5).

RESULTS
We computed the equation of state perturbatively in µ 2 by first solving ∂ φ Ω(φ) = 0 order-by-order for a range of temperatures with γ = 0.1. By plugging this solution into Eq. (4) and expanding in µ 2 we obtained the Taylor series expansion for the pressure p(T, µ) ≈ N n=0 p 2n (T )µ 2n . To illustrate the numerical evaluation of the equation of state we focus on the susceptibility, (7) because its singular part in the vicinity of the critical point it grows as χ(µ) ∼ Re(µ 2 − µ 2 LY ) σ−1 where σ = 1/2 in the mean field limit. Of course, in many cases it is very difficult, to generate such large number of terms. Therefore we also show results obtained by 11 terms for comparison. We computed the singularities both from Padé and conformal Padé which are shown in Fig. 2 for two different temperatures very close to and away from the critical point. We used a simple conformal map, defined over one-cut complex plane with a singularity located at µ 2 LY to resolve f near the singularity µ 2 LY . Since the other singularity is the complex conjugate pair, µ * 2 LY , including its contribution is trivial. Of course, a priori, we don't know what µ LY . Therefore we first obtained a crude estimate for it from regular Padé and we used it in φ 1 (z), and refined this estimate via conformal Padé.
In order to reconstruct the trajectory, µ LY (T ), we repeated this procedure for different temperatures. In order to smooth out the T dependence of µ LY we used fits whose forms are determined from Eq. (6); namely a linear fit for Reµ LY (T ) and a y = ax 3/2 fit for Imµ LY (T ). The results are shown in Fig. 3. From these fits we obtained the values of T c , µ c , c 1 and c 2 shown in Table I.
Finally we computed the susceptibility as a function µ via Padé and conformal Padé. In order to capture the global behavior of the equation of state we used a different conformal map, defined on a two-cut complex plane with two branch points located at |µ LY | 2 e ±iπθ [25,27,28]. The results for two representative temperatures near and away from T c are shown in  solve Imµ LY and creates a sequence of poles/zeroes along the real axis. Even when it does, away from the critical point, T = 1.60T c , there are still spurious poles along the real axis which makes Padé useless for µ Reµ LY as seen in Fig. 4 (gray curves). As seen from the same figures, conformal Padé does not suffer from such a problem and even for 11 terms, it does capture the qualitative behavior of the peak. For 21 terms the agreement with the exact result up to µ ≈ 0.8 is quite remarkable. As mentioned above, in order to refine the location of a given singularity, µ LY , we found the one-cut conformal map given in Eq.(8) to be very accurate and easy to execute numerically. However globally it generates singularities between µ 2 LY and µ * 2 LY as seen in Fig. 2. That is why we used a two-cut map given in Eq. (9) to evaluate χ(T, µ) shown in Fig. 4.
From the singularities extracted from conformal Padé we constructed the LY singularity trajectory Eq. (6). The critical temperature, T c , is obtained from the point Imµ LY vanishes (the vertical line in Fig. 3). Therefore it is crucial to be able to accurately compute Imµ LY where conformal Padé has a significant advantage over other methods, even with an input of 11 Taylor coefficients. The remaining parameters µ c , c 1 , c 2 are extracted from fits explained above. Notably the location of the critical point can be determined with less than %5 error even with 11 terms. The error is larger for the crossover slope, c 1 and c 2 , however it is below %10 for all cases except for c 2 with 11 terms (≈ %21).

SUMMARY AND CONCLUSIONS
In this letter we tackled a general problem: given a truncated series expansion of the equation of state of a theory obtained away from any critical point, what can we say about a the existence of a critical point in the phase diagram and the singular behavior of the equation of state around it? We considered a fairly common case of a first order transition that ends at a second order critical point (which is in the Ising universality class) and turns into smooth crossover. We introduced a new way of extracting the leading singular behavior of the equation of state around the Lee-Yang edge singularity in the complex µ plane by using a combination of a conformal map and Padé resummation. Equipped with the knowledge of this trajectory, we showed how to extract the location of the critical point and to constrain the non-universal Ising mapping parameters in the scaling region. Finally we demonstrated that this method significantly improves the numerical evaluation of the equation of state even with a modest number of Taylor coefficients.
There are various extensions and refinements of these ideas that we left for future work. For example, one can extend this analysis to the first order side, T < T c where one expects a the scaling part of the equation of state to jump to a different Riemann sheet [9,34]. This can be achieved by using the so-called "uniformization map" as a part of conformal Padé which gives access to higher sheets. Another possible application is be to study the Fonseca-Zamolodchikov conjecture this way [34]. How to handle various uncertainties in the original Taylor coefficients is another non-trivial problem that is currently under investigation.
An immediate application of the ideas we developed in this letter would be to assist the search for the critical point of QCD both by constraining its location and the Ising mapping parameters, as well as improving the numerical implementation of the equation of state in the hydrodynamic simulations to smoothly interpolate the small µ behavior to the scaling region. Another possible extension of this work is to incorporate it with analytical continuation from imaginary µ, especially with the recently introduced resummation method [35].