Apparent convergence of Padé approximants for the crossover line in finite density QCD
AApparent convergence of Pad´e approximants for thecrossover line in finite density QCD
Attila P´asztor, ∗ Zsolt Sz´ep, † and Gergely Mark´o ‡ ELTE E¨otv¨os Lor´and University, Institute for Theoretical Physics,P´azm´any P. s. 1/A, H-1117, Budapest, Hungary. MTA-ELTE Theoretical Physics Research Group,P´azm´any P. s. 1/A, H-1117 Budapest, Hungary. Fakult¨at f¨ur Physik, Universit¨at Bielefeld, D-33615 Bielefeld, Germany.
We propose a novel Bayesian method to analytically continue observables to real baryochemicalpotential µ B in finite density QCD. Taylor coefficients at µ B = 0 and data at imaginary chemicalpotential µ IB are treated on equal footing. We consider two different constructions for the Pad´eapproximants, the classical multipoint Pad´e approximation and a mixed approximation that is aslight generalization of a recent idea in Pad´e approximation theory. Approximants with spuriouspoles are excluded from the analysis. As an application, we perform a joint analysis of the avail-able continuum extrapolated lattice data for both pseudocritical temperature T c at µ IB from theWuppertal-Budapest Collaboration and Taylor coefficients κ and κ from the HotQCD Collabo-ration. An apparent convergence of [ p/p ] and [ p/p + 1] sequences of rational functions is observedwith increasing p. We present our extrapolation up to µ B ≈
600 MeV.
I. INTRODUCTION
Despite considerable effort invested so far, the phasediagram of QCD in the temperature( T )-baryon chemicalpotential( µ B ) plane still awaits determination from firstprinciples. At the moment, the only solid informationavailable is the curvature of the crossover temperature [1–3], together with some upper bound on the absolute valueof the next Taylor coefficient of order O ( µ B ) [2, 3]. Theseresults come either from the evaluation of Taylor coeffi-cients with lattice simulations performed at µ B = 0 or viasimulations performed at imaginary µ B , where the signproblem is absent, with the Taylor coefficients obtainedfrom a subsequent fit.Whether the input data are the Taylor coefficients orthe values of a function at several values of the imagi-nary chemical potential, fact is that the numerical ana-lytic continuation needed to extrapolate the crossover toreal µ B is a mathematically ill-posed problem [4, 5]. Thismeans that although the analytic continuation of a func-tion sampled inside some domain D is uniquely deter-mined by the approximant used, the extention of a func-tion differing on D by no matter how small an amountcan lead to arbitrarily different values at points outside D . That is to say: analytic continuation is unique, butis not a continuous function of the data. For such ill-posed problems, the only way to achieve convergence inthe results is to use some kind of regularization. Thismakes sure that the noise in the data is not overempha-sized by the analytic continuation. As the noise is re-duced, the regularizing term is made weaker. This leadsto a kind of double limit when the regularization and ∗ [email protected] † [email protected] ‡ [email protected] the noise are taken to zero together. The simplest kindof regularization for analytic continuation is the use ofsome ansatz, which is assumed to describe the physicsboth in the range where data is available, and in therange where one tries to extrapolate. The conservativeview is to use for analytic continuation few-parameterapproximants, which all fit the data well, and performthe continuation only in a range where they do not devi-ate much from each other, assessing the systematic errorof the continuation from this deviation. Here we pursuea more adventurous approach, by considering a sequenceof approximants of increasing functional complexity, andtrying to observe whether they converge or not.In the absence of physically motivated ansatz, a goodguess is to to study the [ p/p ] (diagonal) and [ p/p + 1](subdiagonal) Pad´e sequences, as these are only slightlymore complicated to work with than polynomials, buthave far superior convergence properties. Ordinary Pad´eapproximants ( i.e. rational functions constructed usingapproximation-through-order conditions to match theTaylor expansion of a function at a given point) areknown to converge uniformly on the entire cut planefor functions of Stieltjes type [6] (which have a cut onthe negative real axis). For this class of functions thesubdiagonal sequence of multipoint (or N − point) Pad´eapproximants [7], also know as the Schlessinger pointmethod in the context of scattering theory [8], is alsoconvergent (see [9] and references therein). For a mero-morphic function, on the other hand, Pad´e approximantsare known to converge in measure [10, 11], i.e. almosteverywhere on the complex plane, in stark contrast topolynomial approximations, which stop converging at thefirst pole of such a function.While the convergence properties of Pad´e approxi-mants in exact arithmetic are often very good, even incases where the mathematical reason for the convergenceis not fully understood yet, these approximations tendto be very fragile in the presence of noise. This often a r X i v : . [ h e p - l a t ] O c t manifests itself in spurious poles, whose residue goes tozero as the noise level is decreased, as well as spuriouszero-pole pairs (called Froissart doublets [12]). The dis-tance between the zero and the pole goes to zero as thenoise decreases, eventually leading to the annihilationdisappearance of the doublet. There is a large body ofmathematical literature devoted to the removal of thesespurious poles. Procedures which do so typically involvesome further regularization, like in Ref. [13], where thisis based on singular value decomposition, or monitoringthe existence of Froissart doublets for later removal, likein Ref. [14]. In cases where the noise level on the datacannot be arbitrarily decreased , the exclusion of spuri-ous poles is mandatory if one wants to go to higher orderapproximants in the analysis.When dealing with numerical analytic continuation, weneed to select the approximant from a class of possiblefunctions ( i.e. a model) and a method to take into ac-count the data ( i.e. a fitting method). For the formerwe use two types of rational approximants, the classicalmultipoint Pad´e approximants recently used for analyticcontinuation in Refs. [15–17] and a slight generalizationof the the Pad´e-type approximant introduced and studiedrecently in [18]. The parameters of the multipoint Pad´eapproximant are determined solely in terms of the inter-polating points and information on Taylor coefficients, ifit exists, can be taken into account in the second, data fit-ting step. In contrast, the Pad´e-type approximant allowsfor a joint use of interpolating points and Taylor coeffi-cients in determining the parameters of the approximant.Although the focus in [18] was on the diagonal sequence[ p/p ] of Pad´e-type approximants, the method can be eas-ily generalized to construct the subdiagonal sequence aswell. For the data fitting step we use a Bayesian anal-ysis. The likelihood function ensures that approximantsare close to both the data on the Taylor coefficients at µ B = 0 and the data at purely imaginary µ B , while aBayesian prior makes sure that spurious poles are ex-cluded from the extrapolation. Considering two differenttypes of Pad´e approximants is a non-trivial consistencycheck, mainly because the exact form of the prior distri-bution will be different for the two cases, as the numberof interpolation point where the function values will berestricted is different.We note that while Bayesian methods - especially vari-ations of the Maximum Entropy method with differententropy functionals - for the analytical continuation toreal time are quite commonly used in lattice QCD [19–26], as far as we are aware, such methods have not beenapplied to the analytic continuation problem in µ B so far.The only related example we are aware of is Ref. [27],where a Bayesian method is used to extract high order E.g. if the noise is only coming from machine precision in floatingpoint arithmetic, the Froissart doublets can often be removedby simply using multiple precision arithmetic for the “naive”algorithm. derivatives of the pressure around µ B = 0 from data atimaginary µ B . One must note however, that the mathe-matical problem in that paper is that of numerical differ-entiation, which is distinct from the analytic continuationproblem discussed here. This paper, therefore, is the firstattempt of using this class of mathematical techniques tothe analytic continuation problem in finite density QCD.The paper is organized as follows. In Section II we in-troduce the mathematic tools used for our analysis. Firstwe treat the novel Pad´e-type approximants in the absenceof noise. Since the traditional multipoint Pad´e approx-imants are quite well known, they are relegated to Ap-pendix A. Next, we discuss the Bayesian analysis in thepresence of noise in a general manner that includes boththe multipoint Pad´e and the mixed Pad´e approximantcase. In Section III we turn to physical applications. Wefirst demonstrate the effectiveness of Pad´e approximantsin a chiral effective model. Finally, we perform a jointanalysis of the continuum extrapolated lattice data onthe Taylor coefficients at µ B = 0 and the crossover lineat imaginary µ B . Appendix B summarizes the formulasrelating our notational conventions on the Taylor coeffi-cients to those found elsewhere in the literature. II. NUMERICAL METHOD FOR ANALYTICCONTINUATIONA. Pad´e-type rational approximants in the absenceof noise
Using the notation of [18] the mathematical formula-tion of analytic continuation is as follows. Assuming theexistence of a continuous real function f : R → R , wewould like to know its value for t > interpolating points τ i < i =1 , . . . , l, the values f i := f ( τ i ) are known,2. a number of coefficients c i , i = 0 , . . . , k in the Tay-lor expansion f ( t ) = c + c t + · · · + c k t k (1)around t = 0 are known.A widely used method to tackle the problem is to fit arational fraction [ n/m ] ≡ R nm ( t ) ≡ N n ( t ) D m ( t ) := n (cid:88) i =0 a i t im (cid:88) i =0 b i t i , (2) One can choose b = 1 without loss of generality, having n + m +1independent coefficients. to the set of function values and/or available deriva-tives. When only derivatives at t = 0 are used, oneobtains the ordinary Pad´e approximant in which thecoefficients of both the denominator and the numera-tor are fully determined by the Taylor coefficients byimposing the approximation-through-order conditions R nm ( t ) = f ( t ) + O ( t m + n +1 ) . This condition implies thatthe Taylor expansion of the Pad´e approximant around t = 0 agrees with the Taylor expansion of the functionup to and including the order of the highest Taylor coef-ficient known.When only the set of function values at the interpo-lating points are used, one obtains the so-called multi-point Pad´e approximant , which is particularly useful innumerics in its continuous fraction formulation becausethe coefficients can be determined easily from recursionrelations. For odd number of points N = 2 k + 1 , k ≥ , one obtains the approximant [ k/k ], while for even num-ber of points N = 2 k, k ≥
1, one obtains the approximant[ k +1 /k ] . More details can be found in Appendix A, wherethe construction of the multipoint approximant C N issummarized.As mentioned in Ref. [28] (see p. 16), an obvious mod-ification of the multipoint Pad´e approximation can begiven if any number of successive derivatives exists atthe points where the value of the function is known. Re-cently such a modification, called Pad´e-type rational ap-proximant with n = m = k was constructed in [18], with k being the degree up to (and including) which the expan-sion of the approximant matches the Taylor expansion ofthe function. The denominator of this [ k/k ] approximantis fixed by function values at arbitrarily chosen interpo-lating points and the coefficients of the numerator areobtained by imposing the approximation-through-orderconditions.It is easy to generalize the construction used in [18] toobtain Pad´e-type approximants for which l (cid:54) = k. With k + 1 coefficients of the Taylor expansion (including thevalue of the function at zero) one can construct manyPad´e-type approximants of this type, one just has to sat-isfy the relation n + m = k + l . In this case k + 1 coef-ficients of the numerator N n ( t ) are determined from theapproximation-through-order conditions, meaning thatstrictly speaking R nm satisfies by construction n ≥ k, andthe remaining n + m − k = l coefficients are fixed byfunction values at l number of interpolating points via f i D m ( τ i ) = N n ( τ i ) , i = 1 , . . . , l. In what follows we shall use Pad´e-type approximants ofthe form R pp and R pp +1 , with p ≥
1, satisfying 2 p = k + l and 2 p + 1 = k + l, respectively. To construct for ex-ample R ( t ) using c , c , and c , one equates (1) with(2) and after cross multiplication one matches the coeffi-cients of t , t , and t . This gives a = c , a = c + b c , and a = c + b c + b c , which are common for all One equates the expansion of f ( t ) with R nm ( t ), cross multiply andthen equates the coefficients of t on both sides of the equation. approximants with n ≥ . Using these expressions for a , a , and a in N ( t ), one sees that the five conditions f i D ( τ i ) = N ( τ i ) , i = 1 , . . . , a , b , b , b , b , whichcan be easily solved numerically with some standard lin-ear algebra algorithm. As for the R ( t ) approximant,this is constructed very similarly to the original Pad´eapproximant, only the condition on the third derivative(unknown in our case) is replaced by R ( τ ) = f , where τ is an interpolating point. B. Bayesian approach
The Bayesian approach [29] that considers the datasample fixed and the model parameters as random vari-ables gives a perspective on the curve fitting problemwhich is particularly suited for a meta-analysis of datawith noise.We do not include Pad´e approximants of different orderin one large meta-analysis, rather we perform a separateBayesian analysis of the different order approximants, inorder to study their convergence properties as the orderof the approximation is increased. For an [ n/m ] Pad´eapproximant, the model parameters are the coefficients (cid:126)a = ( a , a , . . . , a n ) and (cid:126)b = ( b , b , . . . , b m ), with a totalof n + m + 1 coefficients to be determined. The posteriorprobability can be written as: P ( (cid:126)a,(cid:126)b | data) = 1 Z P (data | (cid:126)a,(cid:126)b ) P prior ( (cid:126)a,(cid:126)b ) , (3)where assuming Gaussian errors around the correctmodel parameters, the likelihood is given by: P (data | (cid:126)a,(cid:126)b ) = exp (cid:18) − χ (cid:19) ,χ = χ + χ µ B ,χ = T (cid:88) i =1 (cid:18) c i − ∂ i R nm ( µ B ; (cid:126)a,(cid:126)b ) ∂ ( µ B ) i (cid:12)(cid:12)(cid:12) µ B =0 (cid:19) σ c i ,χ µ B = L (cid:88) i =1 (cid:16) f i − R nm ( iµ IB,i ; (cid:126)a,(cid:126)b )) (cid:17) σ f i , (4)with T being the number of derivatives known at µ B = 0and L being the number of function values known for µ B < Z is a normalization constant. The Taylorcoefficients at µ B = 0 are clearly correlated, but theircorrelation matrix was not given in Ref. [2] so we ignorethe correlations. If the correlations between the Taylorcoefficients are known, including them in our method iscompletely straightforward. The data at different valuesof imaginary µ B come from different Monte Carlo runs,and are thus uncorrelated.The variables that the Bayesian analysis code uses forthe construction of the Pad´e approximants are not thecoefficients of the polynomial themselves. For the multi-point Pad´e approximants we use a number of interpolatedvalues at fixed node points in µ B /T . For the case ofPad´e-type approximants we use a smaller number of in-terpolated values at node points and a number of deriva-tives at µ B = 0. These are of course in a one-to-one cor-respondence with the polynomial coefficients, once therestriction b = 1 has been made in Eq. (2). Details ofthe implementation will be discussed in Section III B.An important part of our procedure is that we do notwork with the space of all Pad´e approximants of order[ n/m ], rather, the allowed approximants are restricted bythe prior, which always contains a factor that excludesspurious poles both in the interpolated and the extrapo-lated range.The prior also contains a further factor - the exact formof which for the two different Pad´e approximants will bediscussed in Section III B - which prevents extra oscilla-tions of the interpolants in the µ B < µ B /T range. We have checked that ourresults are not sensitive to the choice of the node points.This is also expected on mathematical grounds, since un-like polynomial interpolants, rational interpolants are notextremely sensitive to the choice of the node points usedfor the interpolation [30].Our numerical results will be based on the poste-rior distribution. For a fixed value of µ B /T , we studythe posterior distribution of the crossover temperature T c = R nm (ˆ µ B ) and chemical potential µ B = (cid:112) ˆ µ B T c . Thecenter point will in both cases be the median, while theslightly asymmetric error bars represent the central 68%of the posterior distribution of both quantities. We willcall these percentile based errors . In practice, the integra-tion over the prior distribution is carried out with simpleMonte Carlo algorithms. The statistics needed is suchthat the posterior distribution of the studied observablesdoes not change anymore, which we explicity checked tobe the case in our analysis. III. ANALYTIC CONTINUATION OF T c ( µ IB ) A. Convergence of Pad´e approximants in a chiraleffective model
Before applying the method described in Sec. II tothe actual QCD data, we study the analytic continua-tion within the chiral limit of the two flavor ( N f = 2)constituent quark-meson (CQM) model. We show thatin this model both the diagonal and the subdiagonal se-quences constructed from T c ( µ IB ) exhibit apparent con-vergence to the exact T c ( µ B ) curve. We also show thatthe Pad´e approximant knows nothing about the locationof the tricritical point (TCP), as this information is notencoded in T c ( µ B ) . Finally, we investigate the effect ofthe error on the analytical continuation.
1. Convergence in the absence of noise
In Ref. [31], using leading order large- N techniques re-sulting in an ideal gas approximation for the constituentquarks, the coefficients of the Landau-Ginzburg type ef-fective potential V eff = m Φ + λ eff Φ + . . . for the chiralorder parameter Φ were determined in the chiral limit tobe : m = m + (cid:18) λ
72 + g N c (cid:19) T + g N c π µ B , (5a) λ eff = λ − g N c π (cid:20) Ψ (cid:18)
12 + i ˆ µ B π (cid:19) + Ψ (cid:18) − i ˆ µ B π (cid:19) +2 + 2 ln 4 πTM (cid:21) . (5b)In the expressions above Ψ( x ) is the digamma function,ˆ µ B = µ B /T, N c is the number of colors, g = m q / Φ (with m q = m N / = f π /
2) is the Yukawacoupling between the pion and sigma mesons and theconstituent quarks, and m and λ are the renormal-ized mass and the self-coupling in the O ( N ) symmetricmesonic sector of the CQM model, which at the value M = 886 MeV of the renormalization scale take the val-ues m = − and λ = 400 . For µ B ≥ µ B − T plane, which is ob-tained from the condition m = 0. This line of secondorder points ends in a tricritical point with coordinatesdetermined by m = λ eff = 0 . For µ B > µ TCP B the chi-ral phase transition is of first order and m = 0 givesthe location of the first spinodal down to T = 0 . Themerit of the expressions in (5a) and (5b) is that the lineof second order phase transitions, which is actually anellipse in the µ B − T plane, can be determined analyti-cally together with the location of the TCP. This makesthe analytic continuation very simple, as we just have tochange ˆ µ B → − ˆ µ B in the expression T c (ˆ µ B ) = (cid:115) − m λ + 6 g N c (cid:0) µ B / (3 π ) (cid:1) , (6)obtained from (5a), to go from real to imaginary chemicalpotentials.We can sample T c at imaginary values of ˆ µ B and fitmultipoint Pad´e approximants to the sampled data (seeApp. A). Then, we can evaluate the Pad´e approximant atreal values of ˆ µ B and compare the value of analytic con-tinued T c with the exact values obtained from (6). This These expressions corresponds to Eqs. (13) and (14) of [31], justthat we used the relation ∂∂n (cid:16) Li n ( − e z ) + Li n ( − e − z (cid:17)(cid:12)(cid:12)(cid:12) n =0 = − γ − ln(2 π ) − (cid:2) Ψ (cid:0) (1 + iz/π ) / (cid:1) + Ψ (cid:0) (1 − iz/π ) / (cid:1)(cid:3) / , which canbe proven by comparing the high temperature expansion usedthere with the one given in [32]. T c [ M e V ] µ B [MeV] N=3, [1/1] 4, [1/2] 5, [2/2] 6, [2/3] 7, [3/3] 8, [3/4]exactTCP -10-5 0 5 10 0 50 100 150 200 250 300 350 400 ( T c app r o x / T c e x a c t - ) × µ B2 /T FIG. 1. Apparent convergence of the multipoint Pad´e approx-imants C N (solid lines) determined from T c values at µ IB inthe CQM model in comparison to the Taylor expansion of or-der N − µ B = 0 (dashed lines). In the main plot theparametric curves for the Pad´e approximants are obtained as( (cid:112) ˆ µ B C N (ˆ µ B ) , C N (ˆ µ B )) . The bivaluedness of the subdiagonalapproximants (and the curves of the odd-order Taylor expan-sion) reflects that µ B := (cid:112) ˆ µ B C N (ˆ µ B ) has a maximum. Theinset shows the percentage difference between the approxi-mated and exact values of T c , with the vertical dotted lineindicating the radius of converges of the Taylor expansion. comparison is presented in Fig. 1, where the inset showsthe percentage difference between C N (ˆ µ B ) and T c (ˆ µ B )using the same interpolating points as those shown inFig. 3 in the case of the QCD data. The main figureshows that the diagonal sequence converges from above,while the subdiagonal sequence converges from below tothe line m = 0, which for µ B < µ TCB B is the line ofcritical points and for µ B > µ TCB B is the first spinodal.Given that the sampling range is ˆ µ B ∈ ( − . , , theaccuracy of the [3 /
4] Pad´e approximant around the loca-tion of the TCP is remarkable, even though T c (ˆ µ B ) is arather simple function, as it represents an ellipse. Thisis even more so when one compares to the radius of con-vergence of the Taylor expansion around µ B = 0 whichis ˆ µ B ≈ .
2, as given by the pole in (6). We only referto the TCP because µ B (or ˆ µ B ) is rather large there; thePad´e approximant does not know about the existence ofthe TCP, as this is encoded in the quartic part of thetree-level potential and the second derivative of (cid:104) ¯ qq (cid:105) / Φ( q is the constituent quark field) with respect to Φ, whichjointly determine λ eff .
2. Effect of the error
Next, we investigate what happens when analytic con-tinuation is performed in the presence of noise. As areference point we start by generating T ic configurationswith a normal distribution characterized by mean calcu-lated from (6) and standard deviation corresponding to µ B [MeV] T c [ M e V ] σ T : 1.0%, σ c i = σ c i σ T : 0.1%, σ c i = σ c i σ T : 0.1%, σ c i = σ c i /10 σ T : 0.01%, σ c i = σ c i /10 N=4 µ B [MeV] T c [ M e V ] σ T : 1.0%, σ c i = σ c i σ T : 0.1%, σ c i = σ c i /10 σ T : 10 -3 %, σ c i = σ c i /10 σ T : 10 -4 %, σ c i = σ c i /10 N=5 µ B [MeV] T c [ M e V ] σ T : 1.0%, σ c i = σ c i σ T : 0.1%, no c i σ T : 0.1%, no c i , range: 2 ×σ T : 10 -3 %, σ c i = σ c i /10 σ T : 10 -4 %, σ c i = σ c i /10 σ T : 10 -5 %, σ c i = σ c i /10 µ B [MeV] T c [ M e V ] σ T : 1.0%, σ c i = σ c i σ T : 0.1%, σ c i = σ c i σ T : 10 -3 %, σ c i = σ c i /10 σ T : 10 -4 %, σ c i = σ c i /10 σ T : 10 -5 %, σ c i = σ c i /10 , range: 3 ×
200 400 600 800 1000N=7 µ B [MeV] T c [ M e V ] σ T : 1.0%, σ c i = σ c i σ T : 0.1%, no c i σ T : 10 -5 %, no c i σ T : 10 -5 %, no c i , range: 2 ×σ T : 10 -7 %, no c i σ T : 10 -7 %, no c i , range: 2 ×
200 400 600 800 1000N=8 µ B [MeV] T c [ M e V ] σ T : 1.0%, σ c i = σ c i σ T : 0.1%, σ c i = σ c i /10 σ T : 10 -3 %, σ c i = σ c i /10 σ T : 10 -5 %, σ c i = σ c i /10 σ T : 10 -6 %, σ c i = σ c i /10 , range: 3 ×σ T : 10 -7 %, no c i , range: 2 × FIG. 2. Result of a mock analysis showing the effect of theerror of the data and of the sampling range on the qualityof the analytic continuation obtained via multipoint Pad´e ap-proximants of various order and by including or omitting in-formation (the latter is denoted by ’no c i ’ in the key) on theerror of the Taylor coefficients in the evaluation of χ . Foradditional information see the main text. the relative error σ T ic /T ic = 1% and investigate to whatextent should we decrease the relative error in order toget close to the curves obtained in Fig. 1 in the absenceof noise. Note that σ T ic is by a factor of two larger thanthe average error of the QCD data at imaginary chemicalpotential.We determine the coefficients of the multipoint Pad´eapproximant for each generated configuration, evaluatethe approximant for positive ˆ µ B and, using the Bayesianmethod presented in Sec. II B, calculate χ including oromitting information on the Taylor coefficients, and thenstudy the posterior distribution of these values. Themethod is applied to the QCD data in the next subsec-tion, where it is referred to as Method 1. The controlpoints used to calculate χ have ˆ µ B corresponding to theQCD data at imaginary chemical potential and T c ob-tained from (6). We use a unique relative error of T c at allinterpolating and control points, whose value is indicatedin the key of Fig. 2 ( σ T ic used to generate T ic instancesis twice the indicated value). When Taylor coefficients c and c are also used in the calculation of χ , theirvalues c = − .
575 and c = 0 . σ c = 0 .
626 and σ c = 0 . . The sampling points in therange ˆ µ B ∈ ( − . ,
0] are those used previously to ob-tain Fig. 1. We also investigate the effect of changing thesampling range for fixed value of the error by increasingthe lower bound of the interval by the factor indicated inthe keys of Fig. 2. In the modified range the interpolatingpoints are equidistant from each other.Worth noticing in Fig. 2 is that in the presence of noisethe bands for T c ( µ B ) can deviate above some value of µ B from the curves of Fig. 1 by more than the estimatedstatistical error. This reflects the ill-posedness of theanalytic continuation problem. However, even with thelargest error used, the Pad´e sequence converges up to µ B ≈
600 MeV. For the mathematically curious, we alsoshow the effect of increasing the range of the interpola-tion points. As expected, convergence is accelerated bythe increase of the sampling range. This is of course notdirectly relevant for QCD, as the Roberge-Weiss transi-tion puts a limit on the available range for the interpo-lation points.
B. Analytic continuation of QCD data
The continuation of the critical line to real values of µ B using data at imaginary chemical potential available indifferent formulations of the QCD was in the recent yearsa topic of constant interest [33–36]. Pad´e approximantswere also used, but looking for the apparent convergenceof a sequence was seemingly not in the focus of theseinvestigations, with the exception of [33]. As an applica-tion of the method presented in Sec. II we shall addressthis issue and study the convergence of Pad´e series of theform [ p/p ] and [ p/p + 1] constructed:1. based only on interpolating points (multipointPad´e approximants), or2. using interpolating points and the expansion f ( t ) ≈ c + c t + c t around t = 0 , as explained inSec. III A (Pad´e-type approximant).In order to do so, we use the continuum extrapolatedvalues of T c recently determined on the lattice at µ B = 0and seven imaginary values of ˆ µ B = µ B /T , namelyˆ µ B ( j ) = i jπ , j = 0 , , , , , , . , , (7)given in Table II. of [3] and the Taylor coefficients κ and κ , appearing in the parametrization T c ( µ B ) = T c (0) (cid:2) − κ ( µ B /T c (0)) − κ ( µ B /T c (0)) (cid:3) , (8)also extrapolated to the continuum limit in [2].With the notation of Sec. II, the (assumed) function f ( t ), which corresponds to T c (ˆ µ B ), is known at sevenpoints τ j = ˆ µ B ( j ) <
0, corresponding to j (cid:54) = 0 in the listgiven in (7), and we also know c = T c (ˆ µ B ( j = 0)), aswell as c and c in terms of κ and κ . The values of κ , κ and T c (0) reported in [2] give through the explicitrelations given in App. B c = − .
878 and c = 0 . σ c = 0 .
626 and σ c = 0 . .
1. Numerical implementation of the Bayesian approach
In order to use the method presented in III A, we needto generate { T ic } instances at chosen interpolating points
155 160 165 170 175 180 185 -8 -7 -6 -5 -4 -3 -2 -1 0 T c [ M e V ] µ B2 /T W-Bp lattice dataN=3, [1/1] 5, [2/2] 6, [2/3] 8, [3/4] mP: Pt: σ Tc ∀ N N ≥ ≥ ≠ ≥ ∀ N ∀ N N=8 N ≥ ≥ ∀ N FIG. 3. Choice of the interpolating points, whose positionis indicated by the vertical dotted line, in comparison to theactual lattice data. The labels indicate the use of the inter-polating point in the multipoint Pad´e approximant (bottom)and the Pad´e-type approximant (top) of a given order, charac-terized by the number of independent parameters N . At thevalue of ˆ µ B corresponding to the lattice data points we showthe mean and the error of the T c computed from the mul-tipoint Pad´e approximants generated with importance sam-pling (for the sake of the presentation the abscissa is shifted).The band indicates the standard deviation of the normal dis-tribution used in Method 1 to generate { T ic } instances, i.e.it indicates the prior distribution, excluding the factor thatremoves the spurious poles. (also values of c and c in the case of the Pad´e-type ap-proximant) and then evaluate χ , defined in (4), using theactual lattice data as control points. The interpolatingpoints ˆ µ B,i used for the two types of Pad´e approximantsmentioned above are indicated in Fig. 3. The idea behindour choice was that each interpolating point of any of theused approximant fall in between two nearby lattice datapoints and be more or less equally distributed in the sam-pling range. The actual choice of the interpolating pointsis not important, however, in order to maximize the sam-pling range, one interpolating point is chosen close to thelattice data point with ˆ µ B ( j = 7) and, since we are inter-ested in analytic continuation through µ B = 0 , we alsochoose ˆ µ B ( j = 0) = 0 as an interpolating point.We use two methods to generate input for the Pad´eapproximants. In the first method ( Method 1 ) we simplygenerate T ic from normal distribution with mean obtainedby interpolating the mean of the lattice data and withthe standard deviation (SD) indicated in Fig. 3. In thiscase c and c , used in the Pad´e-type approximant, aregenerated from a normal distribution with mean and SDgiven by Eqs. (B2) and (B3), respectively. As a result, c and c are taken into account in the calculation of χ only when using the multipoint Pad´e approximant. Ac-cording to our prior, we only accept those configurationsfor which the corresponding Pad´e approximant is free ofspurious poles in the wide range ˆ µ B ∈ [ − π , π ]. Whenusing this method we calculate T c at some value of ˆ µ B as T c = R nm (ˆ µ B ) and the value of the real chemical po-tential as µ B = (cid:112) ˆ µ B T c and determine their percentilebased error using the weight e − χ / . The second method (
Method 2 ) for generating inputfor the Pad´e approximants is the importance samplingusing the Metropolis algorithm with “action” χ /
2. Theproposed value of T ic in the Markov chain is generatedusing a normal distribution for the noise with vanishingmean and SD of O (1) MeV. In the case of the Pad´e-typeapproximant we use normal distribution with standarddeviation σ c i , i = 1 , c i . Configurations for which the correspond-ing Pad´e approximant has spurious poles in the rangegiven above are excluded by assigning to them the value χ = ∞ . For the remaining configurations χ is calcu-lated using all the available lattice data according to theformulas in (4). The average and percentile based errorof T c = R nm (ˆ µ B ) and µ B for the Pad´e approximants werecalculated in the standard way with the configurationsprovided by the Metropolis algorithm.There are some peculiarities when doing importancesampling in this context. These are related to the spu-rious poles of the Pad´e approximants, which appear as“walls” of infinite action in the Metropolis update. Con-figurations with spurious poles are not guaranteed to beisolated points in the space of all configurations, rather,there can be regions in configurations space where allapproximants have a pole. One can easily stumble onan accepted configuration that is surrounded in most di-rections by configurations with a pole, thereby trappingthe algorithm. To avoid this problem it is a good ideato mark out a temperature range sampled by the algo-rithm during the random walk and assign infinite valuefor the action if a proposed T ic lies outside this range.I.e. even in the case of Method 2 a prior, like that inFig. 3 is used. An other reason to introduce this band isto exclude the Pad´e approximants from having featuresin the interpolated range not present in the data, even ifsuch an approximant has no pole and fits the data pointsacceptably. This is also a possibility, since Pad´e approx-imants are rather flexible. In practice a two times widerband than the one shown in Fig. 3 proved sufficient.Another observation is that in some cases it was veryhard to thermalize the system by updating the value of T c only at one interpolating point at a time. It proved moreuseful to propose in the Metropolis algorithm an updatedarray of T c values, as this procedure also substantiallyreduced the autocorrelation time.
2. Results for the analytic continuation
The first thing worth checking is the distribution of T c calculated from the Pad´e approximants at ˆ µ B valuescorresponding to the actual lattice data points. For themajority of the approximants and lattice data points, thedistribution is very close to a normal one with standard c (0) (1- κ µ /T c2 (0))=:T cT ( µ B ) T c [ M e V ] µ B [MeV] N=4, [1/2] 5, [2/2] 6, [2/3] 7, [3/3] 8, [3/4] -15-10-5 0 5 10 15 0 5 10 15 20 25 ( R m n / T c T - ) × µ B2 /T FIG. 4. Result of the analytic continuation via the diagonaland subdiagonal sequences of Pad´e-type approximants con-structed based on { T ic } instances generated using importancesampling. The inset shows the percentage difference with re-spect to the T c (0)(1 − κ µ B /T c (0)) curve plotted in the mainfigure, for which we used T c (0) = 158 .
01 MeV and κ = 0 . deviation compatible with that of the lattice data. Thelatter can be seen in Fig. 3 in the case of the multipointPad´e approximant, meaning that the selection of the T c instances based on χ works as expected. However, dif-ferent low order approximants seem to select, withingthe error, different ranges in the distribution of T c (and c i when the Pad´e-type approximant is used). This ismost visible in the case of the [2/2] approximant wherethe points posses a structure unseen in the lattice data.The multipoint Pad´e approximant [1/1] is the most con-strained by the likelihood, the error of T c being smallerthan the lattice one, while the [3/4] approximant is theleast constrained, matching closely the lattice error atall lattice points, and showing a wider range of the com-puted c and c coefficients. This loss of constraint isalso reflected by the χ histogram whose pick moves tohigher values when the number of parameters of the ap-proximant increases.Now we turn our attention to the extrapolation. Sinceboth Method 1 and
Method 2 used to generate input forboth the multipoint Pad´e and Pad´e-type approximantsresulted in very similar results for the analytically con-tinued T c ( µ B ) curve, we only present those obtained with Method 2 (importance sampling) in the case of the Pad´e-type approximant constructed using the first and sec-ond derivative of T c ( µ B ) at µ B = 0 . The fact that with
Method 1 the analytic continuation does not depend onthe approximant used, means that it makes no differencewhether we take into account the Taylor coefficients onlyin the approximant or only in the calculation of χ . Weremind that when importance sampling is used the Tay-lor coefficients are taken into account in the calculationof χ irrespective of the type of approximant, since oth-erwise the range in which c and c varies during therandom walk would not be constrained.Our main result on the analytic continuation is pre-sented in Fig. 4 in comparison with a simple parametriza-tion of the crossover line based on the Taylor coefficient κ . One sees that with the exception of the [2/2] type,the Pad´e approximants tend to give smaller T c with in-creasing µ B . Also, the behavior of the diagonal an sub-diagonal sequences follow different patterns, similar tothat observed in the model study in Fig. 1. Appar-ently, the Pad´e sequences converge, as, although [3/3]and [3/4] have overlapping error bars of similar size, thelatter moves towards the band laid out by the [2/3] ap-proximant. It remains to be seen if this pattern survivesthe possible addition of new lattice data points, whichwill further constrain the fit, and/or an increase in theprecision of the lattice data.
IV. CONCLUSIONS AND OUTLOOK
We presented a method for the numerical analytic con-tinuation of data available at imaginary chemical poten-tial that uses also the Taylor coefficients of an expansionaround µ B = 0 . Using lattice data that became availablerecently, we have investigated the continuation to real µ B of the crossover line with a sequence of Pad´e approx-imants, looking for apparent convergence as the numberof independent coefficients increases. Such an analysiswould have been less conclusive using the smaller dataset available at imaginary µ B in [36] and without takinginto account the lattice data for the Taylor coefficients.Our largest order Pad´e approximants is very close tothe simplest quadratic curve obtained with just the κ coefficient. This means that if the observed apparentconvergence is genuine, such a quadratic approximationmight be applicable in a rather large range of µ B . Wewould like to stress that, as discussed in the case of aneffective model in Sec. III A, our results on the analyticcontinuation tell nothing on the possible existence andlocation of the critical end point (CEP). It is also notpossible to clearly determine the value of µ B up to whichthe analytic continuation could be trusted.The Taylor and imaginary chemical potential methodsare usually considered to be competitors in the study offinite density QCD. This is somewhat unfortunate, as thetwo methods tend to provide complimentary information.With the Taylor method, lower order coefficients tend tobe more precise, while data at imaginary µ B tends torestrict higher order coefficients better, without giving avery precise value for the lower orders. For the case ofbaryon number fluctuations, this can clearly be seen bycomparing Fig. 3 of Ref. [27], where the signal for χ B and χ B is better, with Fig. 1 of Ref. [37], where χ B ismuch more precise. This means that joint analysis of suchdata might be a good idea also for the equation of state,where there are some indications - both from an explicitcalculation on coarser lattices [38–40] and phenomenolog-ical arguments [41, 42] - that the radius of convergence for temperatures close to the crossover is of the order µ B /T ≈
2, making a Taylor ansatz unusable beyond thatpoint. This makes it mandatory to try different ansatze,or resummations of the Taylor expansion, and one possi-ble choice could be the Pad´e approximation method usedhere.
ACKNOWLEDGMENTS
We would like to thank Sz. Bors´anyi, M. Giordano,S. Katz and Z. R´acz for illuminating discussions onthe subject and J. G¨unther for providing the raw lat-tice data of [36] in an early stage of the project. Thiswork was partially supported by the Hungarian NationalResearch, Development and Innovation Office - NKFIHgrants KKP126769 and PD 16 121064, as well as by theDFG (Emmy Noether Program EN 1064/2-1). A. P. issupported by the J´anos B´olyai Research Scholarship ofthe Hungarian Academy of Sciences and by the ´UNKP-20-5 New National Excellence Program of the Ministryof Innovation and Technology.
Appendix A: The multipoint Pad´e approximationmethod
Following Refs. [28] and [43], we briefly summarize theconstruction of the multipoint Pad´e approximant usedto analytically continue functions known only at a finitenumber of points of the complex plane. In our case thecontinuation is done along the real axis, from negative topositive values.When one knows the function at N points f i = f ( z i ) , i = 0 . . . N −
1, the rational function approximating f ( z ) is most conveniently given as a truncated continuedfraction C N ( z ) = A A ( z − z )1+ · · · A N − ( z − z N − )1 , (A1)where we used the notation x ≡ x . The task isto determine the N coefficients A i from the conditions C N ( z i ) = f i , i = 0 . . . N − . Note that only N − z i appear in (A1), z N − appears in the condition C N ( z N − ) = f N − . The coefficients can be obtained effi-ciently as A i = g i ( z i ) , i = 0 . . . N −
1, with the functions g i ( z ) defined by the recursion g p ( z ) = g p − ( z p − ) − g p − ( z )( z − z p − ) g p − ( z ) , ≤ p ≤ N − , (A2)with initial condition g ( z ) = f ( z ) , which means g ( z i ) = f i , when the function is know only in some discretepoints. Working out explicitly the condition A i = g i ( z i )for a few values of i , one sees that one needs to con-struct an upper triangular matrix t i,j using the recursion t i,j = ( t i − ,i − /t i − ,j − / ( z j − z i − ) , for j = 1 , . . . , N − i = 1 , . . . , j , starting from its first row t ,j = f j , j = 0 , . . . , N −
1. The diagonal elements are the coeffi-cients of C N : A i = t i,i . The relation of C N with the Pad´esequence is as follows: if N ≥ C N = [ p/p ]with p = ( N − /
2, while when N ≥ C N = [ p/p + 1] with p = − N/ . Writing C N ( z ) in the form C N ( z ) = N ( z ) /D ( z ) , thenumerator and denominator (at a given value of z ) canbe easily determined from the coefficients of the trun-cated continued fraction via the following three-term re-currence relation X n +1 = X n + ( z − z n ) A n +1 X n − , (A3)where for the numerator ( X = N ) one has X = 0 , X = A and for the denominator ( X = D ) one has X = X =1 , and the iteration goes from n = 0 up to and including n = N − . The coefficients a i ( b i ) of the numerator(denominator) can be easily obtained by calling the aboverecursion (A3) at a finite number of points z and solvinga system of linear equations. Appendix B: Relating c , with κ , In order to relate the coefficients c and c of the Tay-lor expansion T c ( t ) = c + c t + c t , with t = ˆ µ B , withthe coefficients κ and κ used by the HotQCD Collabo-ration, the expansion (8) has to be rewritten in terms ofˆ µ B : T c ( µ B ) T c (0) = 1 − κ T c ( µ B ) T c (0) ˆ µ B − κ T c ( µ B ) T c (0) ˆ µ B . (B1)Then, using that T c ( µ B ) /T c (0) ≈ − κ ˆ µ B and T c ( µ B ) /T c (0) = 1 + O (ˆ µ B ), we obtain c = − κ T c (0) and c = (cid:0) κ − κ (cid:1) T c (0) . (B2)Ignoring the covariance between κ and κ , which is notknown to us, the error associated to these Taylor coeffi-cients are σ c = (cid:2) T c (0) σ κ + κ σ T c (0) (cid:3) ,σ c = (cid:2) T c (0) (cid:0) σ κ + 16 κ σ κ (cid:1) + (cid:0) κ + 4 κ (cid:1) σ T c (0) (cid:3) . (B3)In c and c and their errors we use the data of theHotQCD Collaboration also for T c (0) . [1] C. Bonati, M. D’Elia, F. Negro, F. Sanfilippo andK. Zambello, “Curvature of the pseudocritical line inQCD: Taylor expansion matches analytic continuation,”Phys. Rev. D , no.5, 054510 (2018).[2] A. Bazavov et al. [HotQCD], “Chiral crossover in QCDat zero and non-zero chemical potentials,” Phys. Lett. B , 15-21 (2019).[3] S. Bors´anyi, Z. Fodor, J. N. G¨unther, R. Kara,S. D. Katz, P. Parotto, A. P´asztor, C. Ratti andK. K. Szab´o, “The QCD crossover at finite chemical po-tential from lattice simulations,” Phys. Rev. Lett. ,no.5, 052001 (2020).[4] A. N. Tikhonov and V. Y. Arsenin, Solutions of Ill-PosedProblems (V. H. Winston & Sons, Washington D.C.,1977).[5] I. Sabba Stefanescu, “On the stable analytic continua-tion with a condition of uniform boundedness,” J. Math.Phys. (1986), 2657.[6] C. M. Bender and S. A. Orszag, Advanced Mathemati-cal Methods for Scientists and Engineers I, AsymptoticMethods and Perturbation Theory (Springer, New York,1999).[7] G. A. Baker,
Essentials of Pad´e approximants (Acad.Press, New York, 1975).[8] L. Schlessinger, “Use of Analyticity in the Calculation ofNonrelativistic Scattering Amplitudes” Phys. Rev. ,1411, (1969).[9] J. Gelfgren, “Multipoint Pad´e approximants used forpiecewise rational interpolation and for interpolation tofunctions of Stieltjes’ type”, Research Report, Universityof Ume˚a, Dept. of Mathematics, no. 8, 1978. [10] J. Nuttall, “Convergence of Pad´e approximants of mero-morphic functions”, J. Math. Anal. Appl. , 101-117 (2013).[14] B. Beckerman, G. Labahn, A. C. Matos, “On rationalfunctions without Froissart doublets”, Numerische Math-ematik , 615-633 (2018).[15] A. Pilaftsis and D. Teresi, “Symmetry Improved CJT Ef-fective Action,” Nucl. Phys. B , no.2, 594-619 (2013).[16] G. Mark´o, U. Reinosa and Zs. Sz´ep, “Pad´e approximantsand analytic continuation of Euclidean Φ-derivable ap-proximations,” Phys. Rev. D , no.3, 036002 (2017).[17] R. A. Tripolt, P. Gubler, M. Ulybyshev and L. VonSmekal, “Numerical analytic continuation of Euclideandata,” Comput. Phys. Commun. , 129-142 (2019).[18] C. Brezinski, M. Redivo-Zaglia, “Pad´e-type rational andbarycentric interpolation”, Numerische Mathematik 125,89-113, (2013).[19] M. Asakawa, T. Hatsuda and Y. Nakahara, “Maxi-mum entropy analysis of the spectral functions in latticeQCD,” Prog. Part. Nucl. Phys. , 459-508 (2001).[20] A. Jakov´ac, P. Petreczky, K. Petrov and A. Velytsky,“Quarkonium correlators and spectral functions at zero and finite temperature,” Phys. Rev. D , 014506 (2007).[21] G. Aarts, C. Allton, J. Foley, S. Hands and S. Kim,“Spectral functions at small energies and the electricalconductivity in hot, quenched lattice QCD,” Phys. Rev.Lett. , 022002 (2007).[22] H. B. Meyer, “Transport Properties of the Quark-GluonPlasma: A Lattice QCD Perspective,” Eur. Phys. J. A , 86 (2011).[23] A. Rothkopf, “Improved Maximum Entropy Analysiswith an Extended Search Space,” J. Comput. Phys. ,106-114 (2013).[24] Y. Burnier and A. Rothkopf, “Bayesian Approach toSpectral Function Reconstruction for Euclidean Quan-tum Field Theories,” Phys. Rev. Lett. (2013),182003.[25] S. Borsanyi et al. , “Charmonium spectral functions from2+1 flavour lattice QCD,” JHEP , 132 (2014).[26] A. Rothkopf, “Bayesian techniques and applications toQCD,” PoS Confinement2018 , 026 (2018).[27] S. Borsanyi, Z. Fodor, J. N. Guenther, S. K. Katz,K. K. Szabo, A. Pasztor, I. Portillo and C. Ratti, “Higherorder fluctuations and correlations of conserved chargesfrom lattice QCD,” JHEP , 205 (2018)[28] G. A. Baker and J. L. Gammel, The Pad´e Approximantin Theoretical Physics , vol. 71 of Mathematics in Scienceand Engineering (Academic Press, New York, 1970).[29] K. P. Murphy,
Machine learning: a probabilistic perspec-tive (MIT Press, Cambridge, MA, USA, 2012).[30] L. N. Trefethen,
Approximation Theory and Approxima-tion Practice (Society for Industrial and Applied Mathe-matics, Philadelphia, 2012).[31] A. Jakov´ac, A. Patk´os, Zs. Sz´ep and P. Sz´epfalusy, “ T − µ phase diagram of the chiral quark model from a largeflavor number expansion,” Phys. Lett. B , 179 (2004).[32] A. Ayala, A. Bashir, J. J. Cobos-Martinez, S. Hernandez-Ortiz and A. Raya, “The effective QCD phase diagramand the critical end point,” Nucl. Phys. B (2015),77-86. [33] M. P. Lombardo, “Series representation: Pade’ approx-imants and critical behavior in QCD at nonzero T and µ ,” PoS LAT2005 , 168 (2006).[34] P. Cea, L. Cosmai, M. D’Elia, C. Manneschi and A. Papa,“Analytic continuation of the critical line: Suggestionsfor QCD,” Phys. Rev. D , 034501 (2009).[35] P. Cea, L. Cosmai, M. D’Elia, A. Papa and F. Sanfil-ippo, “The critical line of two-flavor QCD at finite isospinor baryon densities from imaginary chemical potentials,”Phys. Rev. D , 094512 (2012).[36] R. Bellwied, S. Bors´anyi, Z. Fodor, J. G¨unther,S. D. Katz, C. Ratti and K. K. Szab´o, “The QCD phasediagram from analytic continuation,” Phys. Lett. B ,559-564 (2015).[37] A. Bazavov et al. “Skewness, kurtosis, and the fifth andsixth order cumulants of net baryon-number distributionsfrom lattice QCD confront high-statistics STAR data,”Phys. Rev. D , no.7, 074502 (2020).[38] M. Giordano and A. P´asztor, “Reliable estimation of theradius of convergence in finite density QCD,” Phys. Rev.D , no.11, 114510 (2019).[39] M. Giordano, K. Kapas, S. D. Katz, D. Nogradi andA. Pasztor, “Radius of convergence in lattice QCD atfinite µ B with rooted staggered fermions,” Phys. Rev. D , no.7, 074511 (2020).[40] M. Giordano, K. Kapas, S. D. Katz, D. Nogradi andA. Pasztor, “Towards a reliable lower bound on the lo-cation of the critical endpoint,” [arXiv:2004.07066 [hep-lat]].[41] A. Connelly, G. Johnson, S. Mukherjee and V. Skokov,“Universality driven analytic structure of QCD crossover:radius of convergence and QCD critical point,”[arXiv:2004.05095 [hep-ph]].[42] S. Mukherjee and V. Skokov, “Universality driven ana-lytic structure of QCD crossover: radius of convergencein baryon chemical potential,” [arXiv:1909.04639 [hep-ph]].[43] H. J. Vidberg and J. W. Serene, “Solving the Eliash-berg equations by means of N-point Pad´e approximants,”J. Low. Temp. Phys.29