Apparent convergence of Pad\'e approximants for the crossover line in finite density QCD

We propose a novel Bayesian method to analytically continue observables to real baryochemical potential $\mu_B$ in finite density QCD. Taylor coefficients at $\mu_B=0$ and data at imaginary chemical potential $\mu_B^I$ are treated on equal footing. We consider two different constructions for the Pad\'e approximants, the classical multipoint Pad\'e approximation and a mixed approximation that is a slight generalization of a recent idea in Pad\'e approximation theory. Approximants with spurious poles are excluded from the analysis. As an application, we perform a joint analysis of the available continuum extrapolated lattice data for both pseudocritical temperature $T_c$ at $\mu_B^I$ from the Wuppertal-Budapest Collaboration and Taylor coefficients $\kappa_2$ and $\kappa_4$ from the HotQCD Collaboration. An apparent convergence of $[p/p]$ and $[p/p+1]$ sequences of rational functions is observed with increasing $p.$ We present our extrapolation up to $\mu_B\approx 600$ MeV.


I. INTRODUCTION
Despite considerable effort invested so far, the phase diagram of QCD in the temperature(T )-baryon chemical potential(µ B ) plane still awaits determination from first principles. At the moment, the only solid information available is the curvature of the crossover temperature [1][2][3], together with some upper bound on the absolute value of the next Taylor coefficient of order O(µ 4 B ) [2,3]. These results come either from the evaluation of Taylor coefficients with lattice simulations performed at µ B = 0 or via simulations performed at imaginary µ B , where the sign problem is absent, with the Taylor coefficients obtained from a subsequent fit.
Whether the input data are the Taylor coefficients or the values of a function at several values of the imaginary chemical potential, fact is that the numerical analytic continuation needed to extrapolate the crossover to real µ B is a mathematically ill-posed problem [4,5]. This means that although the analytic continuation of a function sampled inside some domain D is uniquely determined by the approximant used, the extention of a function differing on D by no matter how small an amount can lead to arbitrarily different values at points outside D. That is to say: analytic continuation is unique, but is not a continuous function of the data. For such illposed problems, the only way to achieve convergence in the results is to use some kind of regularization. This makes sure that the noise in the data is not overemphasized by the analytic continuation. As the noise is reduced, the regularizing term is made weaker. This leads to a kind of double limit when the regularization and * apasztor@bodri.elte.hu † szepzs@achilles.elte.hu ‡ gmarko@physik.uni-bielefeld.de the noise are taken to zero together. The simplest kind of regularization for analytic continuation is the use of some ansatz, which is assumed to describe the physics both in the range where data is available, and in the range where one tries to extrapolate. The conservative view is to use for analytic continuation few-parameter approximants, which all fit the data well, and perform the continuation only in a range where they do not deviate much from each other, assessing the systematic error of the continuation from this deviation. Here we pursue a more adventurous approach, by considering a sequence of approximants of increasing functional complexity, and trying to observe whether they converge or not.
In the absence of physically motivated ansatz, a good guess is to to study the [p/p] (diagonal) and [p/p + 1] (subdiagonal) Padé sequences, as these are only slightly more complicated to work with than polynomials, but have far superior convergence properties. Ordinary Padé approximants (i.e. rational functions constructed using approximation-through-order conditions to match the Taylor expansion of a function at a given point) are known to converge uniformly on the entire cut plane for functions of Stieltjes type [6] (which have a cut on the negative real axis). For this class of functions the subdiagonal sequence of multipoint (or N −point) Padé approximants [7], also know as the Schlessinger point method in the context of scattering theory [8], is also convergent (see [9] and references therein). For a meromorphic function, on the other hand, Padé approximants are known to converge in measure [10,11], i.e. almost everywhere on the complex plane, in stark contrast to polynomial approximations, which stop converging at the first pole of such a function.
While the convergence properties of Padé approximants in exact arithmetic are often very good, even in cases where the mathematical reason for the convergence is not fully understood yet, these approximations tend to be very fragile in the presence of noise. This often manifests itself in spurious poles, whose residue goes to zero as the noise level is decreased, as well as spurious zero-pole pairs (called Froissart doublets [12]). The distance between the zero and the pole goes to zero as the noise decreases, eventually leading to the annihilation disappearance of the doublet. There is a large body of mathematical literature devoted to the removal of these spurious poles. Procedures which do so typically involve some further regularization, like in Ref. [13], where this is based on singular value decomposition, or monitoring the existence of Froissart doublets for later removal, like in Ref. [14]. In cases where the noise level on the data cannot be arbitrarily decreased 1 , the exclusion of spurious poles is mandatory if one wants to go to higher order approximants in the analysis.
When dealing with numerical analytic continuation, we need to select the approximant from a class of possible functions (i.e. a model) and a method to take into account the data (i.e. a fitting method). For the former we use two types of rational approximants, the classical multipoint Padé approximants recently used for analytic continuation in Refs. [15][16][17] and a slight generalization of the the Padé-type approximant introduced and studied recently in [18]. The parameters of the multipoint Padé approximant are determined solely in terms of the interpolating points and information on Taylor coefficients, if it exists, can be taken into account in the second, data fitting step. In contrast, the Padé-type approximant allows for a joint use of interpolating points and Taylor coefficients in determining the parameters of the approximant. Although the focus in [18] was on the diagonal sequence [p/p] of Padé-type approximants, the method can be easily generalized to construct the subdiagonal sequence as well. For the data fitting step we use a Bayesian analysis. The likelihood function ensures that approximants are close to both the data on the Taylor coefficients at µ B = 0 and the data at purely imaginary µ B , while a Bayesian prior makes sure that spurious poles are excluded from the extrapolation. Considering two different types of Padé approximants is a non-trivial consistency check, mainly because the exact form of the prior distribution will be different for the two cases, as the number of interpolation point where the function values will be restricted is different.
We note that while Bayesian methods -especially variations of the Maximum Entropy method with different entropy functionals -for the analytical continuation to real time are quite commonly used in lattice QCD [19][20][21][22][23][24][25][26], as far as we are aware, such methods have not been applied 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 1 E.g. if the noise is only coming from machine precision in floating point arithmetic, the Froissart doublets can often be removed by simply using multiple precision arithmetic for the "naive" algorithm.
derivatives of the pressure around µ B = 0 from data at imaginary µ B . One must note however, that the mathematical problem in that paper is that of numerical differentiation, which is distinct from the analytic continuation problem discussed here. This paper, therefore, is the first attempt of using this class of mathematical techniques to the analytic continuation problem in finite density QCD. The paper is organized as follows. In Section II we introduce the mathematic tools used for our analysis. First we treat the novel Padé-type approximants in the absence of noise. Since the traditional multipoint Padé approximants are quite well known, they are relegated to Appendix A. Next, we discuss the Bayesian analysis in the presence of noise in a general manner that includes both the multipoint Padé and the mixed Padé approximant case. In Section III we turn to physical applications. We first demonstrate the effectiveness of Padé approximants in a chiral effective model. Finally, we perform a joint analysis of the continuum extrapolated lattice data on the Taylor coefficients at µ B = 0 and the crossover line at imaginary µ B . Appendix B summarizes the formulas relating our notational conventions on the Taylor coefficients to those found elsewhere in the literature.

II. NUMERICAL METHOD FOR ANALYTIC CONTINUATION
A. Padé-type rational approximants in the absence of noise Using the notation of [18] the mathematical formulation of analytic continuation is as follows. Assuming the existence of a continuous real function f : R → R, we would like to know its value for t > 0 given that: around t = 0 are known.
A widely used method to tackle the problem is to fit a rational fraction 2 to the set of function values and/or available derivatives. When only derivatives at t = 0 are used, one obtains the ordinary Padé approximant in which the coefficients of both the denominator and the numerator are fully determined by the Taylor coefficients by imposing the approximation-through-order conditions 3 R n m (t) = f (t) + O(t m+n+1 ). This condition implies that the Taylor expansion of the Padé approximant around t = 0 agrees with the Taylor expansion of the function up to and including the order of the highest Taylor coefficient known.
When only the set of function values at the interpolating points are used, one obtains the so-called multipoint Padé approximant, which is particularly useful in numerics in its continuous fraction formulation because the coefficients can be determined easily from recursion relations. For odd number of points N = 2k + 1, k ≥ 0, one obtains the approximant [k/k], while for even number of points N = 2k, k ≥ 1, one obtains the approximant [k+1/k]. More details can be found in Appendix A, where the construction of the multipoint approximant C N is summarized.
As mentioned in Ref. [28] (see p. 16), an obvious modification of the multipoint Padé approximation can be given if any number of successive derivatives exists at the points where the value of the function is known. Recently such a modification, called Padé-type rational approximant with n = m = k was constructed in [18], with k being the degree up to (and including) which the expansion of the approximant matches the Taylor expansion of the function. The denominator of this [k/k] approximant is fixed by function values at arbitrarily chosen interpolating points and the coefficients of the numerator are obtained by imposing the approximation-through-order conditions.
It is easy to generalize the construction used in [18] to obtain Padé-type approximants for which l = k. With k + 1 coefficients of the Taylor expansion (including the value of the function at zero) one can construct many Padé-type approximants of this type, one just has to satisfy the relation n + m = k + l. In this case k + 1 coefficients of the numerator N n (t) are determined from the approximation-through-order conditions, meaning that strictly speaking R n m satisfies by construction n ≥ k, and the remaining n + m − k = l coefficients are fixed by function values at l number of interpolating points via In what follows we shall use Padé-type approximants of the form R p p and R p p+1 , with p ≥ 1, satisfying 2p = k + l and 2p + 1 = k + l, respectively. To construct for example R 3 4 (t) using c 0 , c 1 , and c 2 , one equates (1) with (2) and after cross multiplication one matches the coefficients of t 0 , t 1 , and t 2 . This gives a 0 = c 0 , a 1 = c 1 + b 1 c 0 , and a 2 = c 2 + b 1 c 1 + b 2 c 0 , which are common for all 3 One equates the expansion of f (t) with R n m (t), cross multiply and then equates the coefficients of t on both sides of the equation. approximants with n ≥ 2. Using these expressions for a 0 , a 1 , and a 2 in N 3 (t), one sees that the five conditions f i D 4 (τ i ) = N 3 (τ i ), i = 1, . . . , 5 represents a system of linear equations for the five unknown a 3 , b 1 , b 2 , b 3 , b 4 , which can be easily solved numerically with some standard linear algebra algorithm. As for the R 1 2 (t) approximant, this is constructed very similarly to the original Padé approximant, only the condition on the third derivative (unknown in our case) is replaced by R 1 2 (τ 1 ) = f 1 , where τ 1 is an interpolating point.

B. Bayesian approach
The Bayesian approach [29] that considers the data sample fixed and the model parameters as random variables gives a perspective on the curve fitting problem which is particularly suited for a meta-analysis of data with noise.
We do not include Padé approximants of different order in one large meta-analysis, rather we perform a separate Bayesian analysis of the different order approximants, in order to study their convergence properties as the order of the approximation is increased. For an [n/m] Padé approximant, the model parameters are the coefficients a = (a 0 , a 1 , . . . , a n ) and b = (b 1 , b 2 , . . . , b m ), with a total of n + m + 1 coefficients to be determined. The posterior probability can be written as: where assuming Gaussian errors around the correct model parameters, the likelihood is given by: with T being the number of derivatives known at µ B = 0 and L being the number of function values known for µ 2 B < 0. Z is a normalization constant. The Taylor coefficients at µ B = 0 are clearly correlated, but their correlation matrix was not given in Ref. [2] so we ignore the correlations. If the correlations between the Taylor coefficients are known, including them in our method is completely straightforward. The data at different values of imaginary µ B come from different Monte Carlo runs, and are thus uncorrelated.
The variables that the Bayesian analysis code uses for the construction of the Padé approximants are not the coefficients of the polynomial themselves. For the multipoint Padé approximants we use a number of interpolated values at fixed node points in µ 2 B /T 2 . For the case of Padé-type approximants we use a smaller number of interpolated values at node points and a number of derivatives at µ B = 0. These are of course in a one-to-one correspondence with the polynomial coefficients, once the restriction b 0 = 1 has been made in Eq. (2). Details of the implementation will be discussed in Section III B.
An important part of our procedure is that we do not work with the space of all Padé approximants of order [n/m], rather, the allowed approximants are restricted by the prior, which always contains a factor that excludes spurious poles both in the interpolated and the extrapolated range.
The prior also contains a further factor -the exact form of which for the two different Padé approximants will be discussed in Section III B -which prevents extra oscillations of the interpolants in the µ 2 B < 0 range, which are not warranted by the data. This is enforced by using a prior distribution of the interpolated values at the node points at fixed µ 2 B /T 2 range. We have checked that our results are not sensitive to the choice of the node points. This is also expected on mathematical grounds, since unlike polynomial interpolants, rational interpolants are not extremely sensitive to the choice of the node points used for the interpolation [30].
Our numerical results will be based on the posterior distribution. For a fixed value of µ B /T , we study the posterior distribution of the crossover temperature T c = R n m (μ 2 B ) and chemical potential µ B = μ 2 B T c . The center point will in both cases be the median, while the slightly asymmetric error bars represent the central 68% of the posterior distribution of both quantities. We will call these percentile based errors. In practice, the integration over the prior distribution is carried out with simple Monte Carlo algorithms. The statistics needed is such that the posterior distribution of the studied observables does not change anymore, which we explicity checked to be the case in our analysis.

III. ANALYTIC CONTINUATION OF Tc(µ I B )
A. Convergence of Padé approximants in a chiral effective model Before applying the method described in Sec. II to the actual QCD data, we study the analytic continuation within the chiral limit of the two flavor (N f = 2) constituent quark-meson (CQM) model. We show that in this model both the diagonal and the subdiagonal sequences constructed from T c (µ I B ) exhibit apparent convergence to the exact T c (µ B ) curve. We also show that the Padé approximant knows nothing about the location of the tricritical point (TCP), as this information is not encoded in T c (µ B ). Finally, we investigate the effect of the error on the analytical continuation.

Convergence in the absence of noise
In Ref. [31], using leading order large-N techniques resulting in an ideal gas approximation for the constituent quarks, the coefficients of the Landau-Ginzburg type ef- . . for the chiral order parameter Φ were determined in the chiral limit to be 4 : In the expressions above Ψ(x) is the digamma function, obtained from (5a), to go from real to imaginary chemical potentials.
We can sample T c at imaginary values ofμ B and fit multipoint Padé approximants to the sampled data (see App. A). Then, we can evaluate the Padé approximant at real values ofμ B and compare the value of analytic continued T c with the exact values obtained from (6). This 4 These expressions corresponds to Eqs. (13) and (14) of [31], just that we used the relation ∂ ∂n Lin(−e z ) + Lin(−e −z n=0 = −γ − ln(2π) − Ψ (1 + iz/π)/2 + Ψ (1 − iz/π)/2 /2, which can be proven by comparing the high temperature expansion used there with the one given in [32].   comparison is presented in Fig. 1, where the inset shows the percentage difference between C N (μ 2 B ) and T c (μ 2 B ) using the same interpolating points as those shown in Fig. 3 in the case of the QCD data. The main figure shows that the diagonal sequence converges from above, while the subdiagonal sequence converges from below to the line m 2 eff = 0, which for µ B < µ TCB B is the line of critical points and for µ B > µ TCB B is the first spinodal. Given that the sampling range isμ 2 B ∈ (−7.35, 0], the accuracy of the [3/4] Padé approximant around the location of the TCP is remarkable, even though T c (μ 2 B ) is a rather simple function, as it represents an ellipse. This is even more so when one compares to the radius of convergence of the Taylor expansion around µ B = 0 which isμ 2 B ≈ 44.2, as given by the pole in (6). We only refer to the TCP because µ B (orμ 2 B ) is rather large there; the Padé approximant does not know about the existence of the TCP, as this is encoded in the quartic part of the tree-level potential and the second derivative of qq /Φ (q is the constituent quark field) with respect to Φ, which jointly determine λ eff .

Effect of the error
Next, we investigate what happens when analytic continuation is performed in the presence of noise. As a reference point we start by generating T i c configurations with a normal distribution characterized by mean calculated from (6) and standard deviation corresponding to the relative error σ T i c /T i c = 1% and investigate to what extent should we decrease the relative error in order to get close to the curves obtained in Fig. 1 in the absence of noise. Note that σ T i c is by a factor of two larger than the average error of the QCD data at imaginary chemical potential.
We determine the coefficients of the multipoint Padé approximant for each generated configuration, evaluate the approximant for positiveμ 2 B and, using the Bayesian method presented in Sec. II B, calculate χ 2 including or omitting information on the Taylor coefficients, and then study the posterior distribution of these values. The method is applied to the QCD data in the next subsection, where it is referred to as Method 1. The control points used to calculate χ 2 haveμ 2 B corresponding to the QCD data at imaginary chemical potential and T c obtained from (6). We use a unique relative error of T c at all interpolating and control points, whose value is indicated in the key of Fig. 2 (σ T i c used to generate T i c instances is twice the indicated value). When Taylor coefficients c 1 and c 2 are also used in the calculation of χ 2 , their values c 1 = −1.575 and c 2 = 0.0267 is determined from the Taylor expansion of (6), as for the reference value of their error, indicated in the key of Fig. 2, we use the error of the QCD data obtained from Ref. [2], namely σ 0 c1 = 0.626 and σ 0 c2 = 0.627. The sampling points in the rangeμ 2 B ∈ (−7.35, 0] are those used previously to obtain Fig. 1. We also investigate the effect of changing the sampling range for fixed value of the error by increasing the lower bound of the interval by the factor indicated in the keys of Fig. 2. In the modified range the interpolating points are equidistant from each other. Worth noticing in Fig. 2 is that in the presence of noise the bands for T c (µ B ) can deviate above some value of µ B from the curves of Fig. 1 by more than the estimated statistical error. This reflects the ill-posedness of the analytic continuation problem. However, even with the largest error used, the Padé sequence converges up to µ B ≈ 600 MeV. For the mathematically curious, we also show the effect of increasing the range of the interpolation points. As expected, convergence is accelerated by the increase of the sampling range. This is of course not directly relevant for QCD, as the Roberge-Weiss transition puts a limit on the available range for the interpolation 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 in different formulations of the QCD was in the recent years a topic of constant interest [33][34][35][36]. Padé approximants were also used, but looking for the apparent convergence of a sequence was seemingly not in the focus of these investigations, with the exception of [33]. As an application of the method presented in Sec. II we shall address this issue and study the convergence of Padé series of the form [p/p] and [p/p + 1] constructed: given in Table II. of [3] and the Taylor coefficients κ 2 and κ 4 , appearing in the parametrization also extrapolated to the continuum limit in [2].
With the notation of Sec. II, the (assumed) function f (t), which corresponds to T c (μ 2 B ), is known at seven points τ j =μ 2 B (j) < 0, corresponding to j = 0 in the list given in (7), and we also know c 0 = T c (μ B (j = 0)), as well as c 1 and c 2 in terms of κ 2 and κ 4 . The values of κ 2 , κ 4 and T c (0) reported in [2] give through the explicit relations given in App. B c 1 = −1.878 and c 2 = 0.0451 with errors σ c1 = 0.626 and σ c2 = 0.627.

Numerical implementation of the Bayesian approach
In order to use the method presented in III A, we need to generate {T i c } instances at chosen interpolating points  (also values of c 1 and c 2 in the case of the Padé-type approximant) and then evaluate χ 2 , defined in (4), using the actual lattice data as control points. The interpolating pointsμ 2 B,i used for the two types of Padé approximants mentioned above are indicated in Fig. 3. The idea behind our choice was that each interpolating point of any of the used approximant fall in between two nearby lattice data points and be more or less equally distributed in the sampling range. The actual choice of the interpolating points is not important, however, in order to maximize the sampling range, one interpolating point is chosen close to the lattice data point withμ 2 B (j = 7) and, since we are interested in analytic continuation through µ B = 0, we also chooseμ 2 B (j = 0) = 0 as an interpolating point. We use two methods to generate input for the Padé approximants. In the first method (Method 1 ) we simply generate T i c from normal distribution with mean obtained by interpolating the mean of the lattice data and with the standard deviation (SD) indicated in Fig. 3. In this case c 1 and c 2 , used in the Padé-type approximant, are generated from a normal distribution with mean and SD given by Eqs. (B2) and (B3), respectively. As a result, c 1 and c 2 are taken into account in the calculation of χ 2 only when using the multipoint Padé approximant. According to our prior, we only accept those configurations for which the corresponding Padé approximant is free of spurious poles in the wide rangeμ 2 B ∈ [−π 2 , 60π 2 ]. When using this method we calculate T c at some value ofμ 2 B as T c = R n m (μ 2 B ) and the value of the real chemical potential as µ B = μ 2 B T c and determine their percentile based error using the weight e −χ 2 /2 .
The second method (Method 2 ) for generating input for the Padé approximants is the importance sampling using the Metropolis algorithm with "action" χ 2 /2. The proposed value of T i c in the Markov chain is generated using a normal distribution for the noise with vanishing mean and SD of O(1) MeV. In the case of the Padé-type approximant we use normal distribution with standard deviation σ ci , i = 1, 2 to generate the noise for the Taylor coefficients c i . Configurations for which the corresponding Padé approximant has spurious poles in the range given above are excluded by assigning to them the value χ 2 = ∞. For the remaining configurations χ 2 is calculated using all the available lattice data according to the formulas in (4). The average and percentile based error of T c = R n m (μ 2 B ) and µ B for the Padé approximants were calculated in the standard way with the configurations provided by the Metropolis algorithm.
There are some peculiarities when doing importance sampling in this context. These are related to the spurious poles of the Padé approximants, which appear as "walls" of infinite action in the Metropolis update. Configurations with spurious poles are not guaranteed to be isolated points in the space of all configurations, rather, there can be regions in configurations space where all approximants have a pole. One can easily stumble on an accepted configuration that is surrounded in most directions by configurations with a pole, thereby trapping the algorithm. To avoid this problem it is a good idea to mark out a temperature range sampled by the algorithm during the random walk and assign infinite value for the action if a proposed T i c lies outside this range. I.e. even in the case of Method 2 a prior, like that in Fig. 3 is used. An other reason to introduce this band is to exclude the Padé approximants from having features in the interpolated range not present in the data, even if such an approximant has no pole and fits the data points acceptably. This is also a possibility, since Padé approximants are rather flexible. In practice a two times wider band than the one shown in Fig. 3

proved sufficient.
Another observation is that in some cases it was very hard to thermalize the system by updating the value of T c only at one interpolating point at a time. It proved more useful to propose in the Metropolis algorithm an updated array of T c values, as this procedure also substantially reduced the autocorrelation time.

Results for the analytic continuation
The first thing worth checking is the distribution of T c calculated from the Padé approximants atμ 2 B values corresponding to the actual lattice data points. For the majority of the approximants and lattice data points, the distribution is very close to a normal one with standard deviation compatible with that of the lattice data. The latter can be seen in Fig. 3 in the case of the multipoint Padé approximant, meaning that the selection of the T c instances based on χ 2 works as expected. However, different low order approximants seem to select, withing the error, different ranges in the distribution of T c (and c i when the Padé-type approximant is used). This is most visible in the case of the [2/2] approximant where the points posses a structure unseen in the lattice data. The multipoint Padé approximant [1/1] is the most constrained by the likelihood, the error of T c being smaller than the lattice one, while the [3/4] approximant is the least constrained, matching closely the lattice error at all lattice points, and showing a wider range of the computed c 1 and c 2 coefficients. This loss of constraint is also reflected by the χ 2 histogram whose pick moves to higher values when the number of parameters of the approximant increases. Now we turn our attention to the extrapolation. Since both Method 1 and Method 2 used to generate input for both the multipoint Padé and Padé-type approximants resulted in very similar results for the analytically continued T c (µ B ) curve, we only present those obtained with Method 2 (importance sampling) in the case of the Padétype approximant constructed using the first and second derivative of T c (µ B ) at µ B = 0. The fact that with Method 1 the analytic continuation does not depend on the approximant used, means that it makes no difference whether we take into account the Taylor coefficients only in the approximant or only in the calculation of χ 2 . We remind that when importance sampling is used the Taylor coefficients are taken into account in the calculation of χ 2 irrespective of the type of approximant, since otherwise the range in which c 1 and c 2 varies during the random walk would not be constrained.
Our main result on the analytic continuation is presented in Fig. 4 in comparison with a simple parametrization of the crossover line based on the Taylor coefficient κ 2 . One sees that with the exception of the [2/2] type, the Padé approximants tend to give smaller T c with increasing µ B . Also, the behavior of the diagonal an subdiagonal sequences follow different patterns, similar to that observed in the model study in Fig. 1. Apparently, the Padé sequences converge, as, although [3/3] and [3/4] have overlapping error bars of similar size, the latter moves towards the band laid out by the [2/3] approximant. It remains to be seen if this pattern survives the possible addition of new lattice data points, which will further constrain the fit, and/or an increase in the precision of the lattice data.

IV. CONCLUSIONS AND OUTLOOK
We presented a method for the numerical analytic continuation of data available at imaginary chemical potential that uses also the Taylor coefficients of an expansion around µ B = 0. Using lattice data that became available recently, we have investigated the continuation to real µ B of the crossover line with a sequence of Padé approximants, looking for apparent convergence as the number of independent coefficients increases. Such an analysis would have been less conclusive using the smaller data set available at imaginary µ B in [36] and without taking into account the lattice data for the Taylor coefficients.
Our largest order Padé approximants is very close to the simplest quadratic curve obtained with just the κ 2 coefficient. This means that if the observed apparent convergence is genuine, such a quadratic approximation might be applicable in a rather large range of µ B . We would like to stress that, as discussed in the case of an effective model in Sec. III A, our results on the analytic continuation tell nothing on the possible existence and location of the critical end point (CEP). It is also not possible to clearly determine the value of µ B up to which the analytic continuation could be trusted.
The Taylor and imaginary chemical potential methods are usually considered to be competitors in the study of finite density QCD. This is somewhat unfortunate, as the two methods tend to provide complimentary information. With the Taylor method, lower order coefficients tend to be more precise, while data at imaginary µ B tends to restrict higher order coefficients better, without giving a very precise value for the lower orders. For the case of baryon number fluctuations, this can clearly be seen by comparing Fig. 3 of Ref. [27], where the signal for χ B 6 and χ B 8 is better, with Fig. 1 of Ref. [37], where χ B 4 is much more precise. This means that joint analysis of such data might be a good idea also for the equation of state, where there are some indications -both from an explicit calculation on coarser lattices [38][39][40] and phenomenological 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 that point. This makes it mandatory to try different ansatze, or resummations of the Taylor expansion, and one possible choice could be the Padé approximation method used here.

ACKNOWLEDGMENTS
We would like to thank Sz. Borsányi, M. Giordano, S. Katz and Z. Rácz for illuminating discussions on the subject and J. Günther for providing the raw lattice data of [36] in an early stage of the project. This work was partially supported by the Hungarian National Research, Development and Innovation Office -NKFIH grants KKP126769 and PD 16 121064, as well as by the DFG (Emmy Noether Program EN 1064/2-1). A. P. is supported by the János Bólyai Research Scholarship of the Hungarian Academy of Sciences and by theÚNKP-20-5 New National Excellence Program of the Ministry of Innovation and Technology.
Appendix A: The multipoint Padé approximation method Following Refs. [28] and [43], we briefly summarize the construction of the multipoint Padé approximant used to analytically continue functions known only at a finite number of points of the complex plane. In our case the continuation is done along the real axis, from negative to positive 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 continued fraction where we used the notation 1 1+ x ≡ 1 1+x . The task is to determine the N coefficients A i from the conditions C N (z i ) = f i , i = 0 . . . N − 1. Note that only N − 1 values of z i appear in (A1), z N −1 appears in the condition C N (z N −1 ) = f N −1 . The coefficients can be obtained efficiently as A i = g i (z i ), i = 0 . . . N − 1, with the functions g i (z) defined by the recursion with initial condition g 0 (z) = f (z), which means g 0 (z i ) = f i , when the function is know only in some discrete points. Working out explicitly the condition A i = g i (z i ) for a few values of i, one sees that one needs to construct an upper triangular matrix t i,j using the recursion t i,j = (t i−1,i−1 /t i−1,j −1)/(z j −z i−1 ), for j = 1, . . . , N −1 and i = 1, . . . , j, starting from its first row t 0,j = f j , j = 0, . . . , N − 1. The diagonal elements are the coefficients of C N : A i = t i,i . The relation of C N with the Padé sequence is as follows: if N ≥ 1 is odd, then C N = [p/p] with p = (N − 1)/2, while when N ≥ 1 is even, then C N = [p/p + 1] with p = −1 + N/2.
Writing C N (z) in the form C N (z) = N (z)/D(z), the numerator and denominator (at a given value of z) can be easily determined from the coefficients of the truncated continued fraction via the following three-term recurrence relation where for the numerator (X = N ) one has X 1 = 0, X 0 = A 0 and for the denominator (X = D) one has X 1 = X 0 = 1, and the iteration goes from n = 0 up to and including n = N − 2. The coefficients a i (b i ) of the numerator (denominator) can be easily obtained by calling the above recursion (A3) at a finite number of points z and solving a system of linear equations.