Lattice QCD equation of state at finite chemical potential from an alternative expansion scheme

Taylor expansion of the equation of state of QCD suffers from shortcomings at chemical potentials $\mu_B \geq (2-2.5)T$. First, one faces difficulties inherent in performing such an expansion with a limited number of coefficients; second, higher order coefficients determined from lattice calculations suffer from a poor signal-to-noise ratio. In this work, we present a novel scheme for extrapolating the equation of state of QCD to finite, real chemical potential that can extend its reach further than previous methods. We present continuum extrapolated lattice results for the new expansion coefficients and show the thermodynamic observables up to $\mu_B/T\le3.5$.


I. INTRODUCTION
The phase diagram of Quantum Chromodynamics (QCD) is an open field of investigation, which is currently at the center of intense efforts from the theoretical and experimental communities alike. At vanishing baryon density, the transition between confined and deconfined matter is known to be an analytic crossover [1]. This knowledge was gained through lattice simulations, which represent a systematically improvable method to solve equilibrium QCD. Although at finite baryon density lattice QCD faces a sign problem, numerous results have been published for moderate chemical potentials in recent years [2,3]. New techniques that allow for direct simulations at finite chemical potential in the presence of a sign problem are the subject of intense investigation. Promising results are available from Lefschetz thimbles [4][5][6], the Complex Langevin equation [7][8][9] or reweighting-based methods [10]. However, these novel approaches cannot be applied to large scale QCD simulations with physical quark masses yet.
There are many indirect lattice methods to study QCD at finite density. These are based on the known analytic feature of the QCD free energy at zero baryo-chemical potential (µ B ). The conceptually simplest method is the Taylor expansion, where the leading µ B -derivatives of the relevant observables are calculated [11][12][13][14]. These derivatives were often calculated for chiral observables, and the µ B dependence of the transition temperature could be extracted [15][16][17][18].
An additional twist to the Taylor method is the observation that these coefficients can be efficiently calcu- * Corresponding author: parotto@uni-wuppertal.de lated by simulating at imaginary values of the chemical potential(s), besides µ B = 0 [19,20]. The use of imaginary chemical potentials is motivated by the analytic crossover at µ B = 0, from which the smooth behavior of the thermodynamic observables as a function of µ 2 B follows. This approach has been popular for more than a decade [17,[21][22][23][24][25][26]. This possibility opened an avenue towards finite density physics, which is often referred to as analytical continuation. This name suggests that, besides the computation of the Taylor expansion coefficients, other extrapolation schemes can be established. The success of this method was most visible in the study of the QCD transition line, where continuum extrapolated results are available for the leading µ B dependence [17,[27][28][29], and recently also for the next-to-leading coefficient [30]. On the other hand, the Taylor series is just one of the possible expansion or continuation schemes. For example, the Padé summation was also considered in the context of QCD thermodynamics [31][32][33][34][35] .
The knowledge of the features of the QCD phase diagram from lattice simulations is currently limited to small chemical potentials and data are mostly available only in the transition region. We have to mention, though, that at higher temperatures resummed perturbation theory has provided a quantitative description of the chemical potential dependence of several observables [36][37][38]. Dedicated lattice studies have bridged the gap between the transition region and perturbative temperatures and found perfect agreement [39,40].
On the experimental side, several heavy-ion collision programs are in place, with the explicit intent of mapping out the phase structure of strongly interacting matter. An important tool for the theoretical interpretation of the results from these experiments are hydrodynamic simulations, which describe the evolution of the system created in the collisions. These simulations need the equation of  [20], at µB/T = 1, 2, 3, as a function of the temperature. Different colors correspond to the order to which the expansion is carried out.
state of QCD as an input, in the whole range of temperatures and densities covered in the experiments. Recently, a Bayesian analysis based on a systematic comparison between heavy-ion data and theoretical predictions showed that the posterior distribution over possible equations of states is compatible with the one calculated on the lattice [41]. For this reason, the equation of state at finite density is a crucial ingredient to understand the nature of strongly interacting matter, and to support the experimental program.
The equation of state at vanishing chemical potential has been known now for several years over a broad range of temperatures [42][43][44]. The first continuum extrapolated extension to finite µ B using the Taylor method in Ref. [13] was followed by several works with the intent of extending these results to higher µ B by adding more terms in the Taylor series [14,45]. Currently, even the eighth µ B -derivative of the QCD pressure is available with modest precision from lattice simulations [20,46]. Very recently, similar results were found by solving a QCD-assisted effective theory with functional methods [47].
In this work, we propose a new extrapolation scheme to finite density QCD. We intend to remedy some shortcomings of the Taylor-based equation of state, e.g. the extrapolation through a crossover boundary. We discuss this issue and suggest a solution in Section II. The formalism for our method is worked out in Section III. We then compute the new observables from lattice simulations in Section IV and perform their continuum estimate. Using this lattice input, we construct the finite density thermodynamic functions in Section V. We conclude and discuss further aspects in Section VI.

II. MOTIVATION
The knowledge of the equation of state from lattice simulations commonly consists of the established µ B = 0 result [43,44] and the Taylor expansion coefficients of the pressure around µ B = 0 where χ B j are the j-th derivatives of the normalized pressure: Besides diagonal coefficients, one can also define offdiagonal correlators between different conserved charges in QCD. Correlators between baryon number and strangeness, which we will need in our procedure, are defined as follows Such correlators have phenomenological relevance [48] and they can also be used to extrapolate the equation of state of QCD in the full, four-dimensional phase diagram at finite T, µ B , µ S , µ Q [49,50]. We will use thê µ i = µ i /T shorthand notation in this manuscript. Currently, results for the expansion coefficients are available up to order O(µ 6 B ) [20,46]. The region of validity of the resulting expansion is usually determined by the range in chemical potential within which an apparent convergence is achieved with the available coefficients. This is stated to be µ B /T < ∼ 2 − 2.5 [14,45]. High order derivatives of the pressure are notoriously difficult to calculate, as they suffer from a low signal-tonoise ratio. This is because their direct determination involves large cancellations of different terms containing derivatives of the Dirac operator [39]. Moreover, studies of chiral models revealed that the structure of the temperature dependence of such observables becomes more and more complex when higher orders are considered [51]. It was pointed out in Ref. [20] that the linear µ 2 Bdependence of the crossover temperature may explain the basic structure. This may explain why including one more term in a truncated Taylor series will not always improve the convergence. On the contrary, pathological behaviornamely non-monotonicity in the T -or µ B -dependenceappears in the extrapolated thermodynamic quantities at chemical potentials beyond µ B /T < ∼ 2 − 2.5. This is due to the fact that, for large enough µ B /T , the observables at finite chemical potential are dictated by the µ B = 0 temperature dependence of the last coefficient included in the expansion. Hence, the structures appearing around the QCD transition temperature in higher order coefficients are "translated" into the finite-µ B behavior of e.g., the entropy, baryon density, etc.
Another inherent problem with the Taylor expansion is the fact that it is carried out at constant temperature. This means that the values of the coefficients at µ B = 0 and a certain temperature T , determine the equation of state at the same T at finite µ B , while the pseudo-critical temperature T pc might have varied considerably. While a sufficiently large number of expansion coefficients would lead to smooth extrapolated functions, even though the Taylor coefficients themselves present a complex structure around the transition temperature, the problem here is rather practical. A scheme that could work with fewer coefficients would be much preferable from the numerical cost point of view.
In Fig. 1 we show the baryon density n B (T ) obtained from a Taylor expansion with the coefficients in Ref. [20], atμ B = 1, 2, 3. The extrapolation is shown including an increasing number of coefficients, to show the effect of higher-order ones. The leading-order (LO) and higher truncations refer to ∼μ B ∂n B (T )/∂μ B , or While atμ B = 1 apparent convergence is achieved at the NLO level, for higher chemical potential this is not the case. Especially atμ B = 3, the inclusion of all the coefficients in Ref. [20] causes unphysical non-monotonic behavior. Ultimately, a pathological behaviour could also come from a finite radius of convergence. Incidentally, recent estimates on coarse lattices [52,53], but also universality arguments [54] place the convergence in the same ball-park in µ B .
In this work, we present an alternative summation scheme which can better cope with the fact that the QCD transition temperature presents a µ B -dependence.
We start from the observation that we made while working with imaginary values of the chemical potentials in an earlier work on analytical continuation. In the upper panel of Fig. 2 we show temperature scans of the quantity n B (T )/μ B = χ B 1 (T,μ B )/μ B for several fixed imaginary µ B /T ratios. The 0/0 limit at µ B = 0 can be easily resolved and equals χ B 2 (T ). The T -dependence of the normalized baryon density at finite chemical potential appears to be simply shifted/rescaled towards higher temperatures from the µ B = 0 results for χ B 2 . This behavior is more apparent in the vicinity of the transition, where the slope of these T rescaled using κ=0.0205 curves is larger. At very large, as well as at very low temperatures a simple shift cannot describe the physics, since the curves become extremely flat. A simple rescaling of temperatures can be described as: where the actual temperature difference can be expressed through a µ B -dependent rescaling factor that we write for simplicity as In the lower panel of Fig  A similar behavior is observed for other quantities too. We show in Fig. 3 the first and second order fluctuations of strangeness at imaginary baryon chemical potentials. In analogy with Eq. (4) one has: where T is defined analogously to Eq. (5), albeit with different κ parameters. In order to motivate our alternative summation scheme, let us first consider a crude approximation that we will later refine. We take Eq. (5) at face value and use it together with Eq. (4) to obtain a well defined χ B 1 (T,μ B ) function. We need a χ B 2 (T, 0) function as well, which we borrow from a deliberately simple fit f (T ) = a+b arctan(c(T −d)) to our data on a 48 3 ×12 lattice. In principle, we could not only calculate χ B 1 (T,μ B ) at arbitraryμ B but, blindly believing Eq. (5), one could calculate the higher µ B -derivatives as well. While this will not describe Nature precisely, it can serve as a test for the Taylor expansion method. To this end, we took severalμ B -derivatives of our χ B 1 (T,μ B ) function and calculated its truncated Taylor series for the lowest four orders. Here LO means just plotting χ B 2 (T, 0). We compared our mock curve (labelled as "full") against its Taylor expansion at three real values of the chemical potential (see Fig. 4).
Forμ B = 1, 2, the summation up to LO and NLO is sufficient to perfectly reproduce the function. However, as the chemical potential is increased, the Taylor expansion carried as far as the NNNLO does not reproduce the original function. On the one hand, convergence is achieved more slowly; on the other hand, spurious effects appear, which generally manifest themselves in a non-monotonicity of the function. These spurious effects in truncated Taylor series were pointed out also in Refs. [55][56][57].
The picture emerging from this simple analysis is rather suggestive, especially when compared to the results shown in Fig. 1 (right panel) obtained from actual lattice data. We also note that this simple analysis does not suffer from the additional complications of signal extraction for higher order expansion coefficients, which in turn play a relevant role in the real-data analysis.

III. FORMALISM
At vanishing chemical potential, it is possible to express the normalized baryon density as a Taylor expansion: As we saw in Fig. 2, the behavior of at finite chemical potential clearly resembles that of χ B 2 (T,μ B ), although shifted/rescaled in temperature. As long as χ B 1 /μ B is a monotonic function of T , the finite density physics can be encoded into the T (T,μ B ) function. A straightforward, but systematic generalization of Eq. (5) reads: In the above equation, we introduced the new parameters κ BB 2 (T ) and κ BB 4 (T ), which describe the shift/rescaling of the temperature of χ B 1 /μ B at finite µ B . Analogous parameters will be introduced below to for the case of and κ BS 4 ) and of χ S 2 (κ SS 2 and κ SS 4 ) at finite µ B . In a way, this formalism replaces the fixed temperature µ B expansion by a fixed-observable temperature expansion.
Having now two expressions, Eqs. (7) and (4), for the same quantity we require their equality at each order in theμ B expansion at µ B = 0, having: which in turn yields: A similar treatment can be applied to the other observables. For the second order fluctuations including baryon number and strangeness, one can consider: 120 χ BS 51 (T, 0)+· · · (11) and: Similarly as before, one can show that: and: κ SS 4 (T ) = 1 24χ S 2 (T ) Before we embark into the discussion of the lattice analysis, let us have an impression on the discussed quantities. Eq. (10) gives a way to directly determine κ BB 2 (T ) and κ BB 4 (T ) using only µ B = 0 data. This approach might be subject to numerical problems, especially in the case of κ BB 4 (T ), which is obtained as the difference of two competing terms. Notice, too, that Eq. (10) contains temperature-derivatives of the χ(T ) coefficients, which may pose a numerical challenge, unless the coefficients are known at sufficient statistics and resolution in T .
On lattices where high statistics data taking is feasible, we can investigate Eq. (10), at least for κ BB 2 . In the top panel of Fig. 5 we compare the numerator and (rescaled) denominator of Eq. (10), while their ratio κ BB 2 (T ) is shown in the bottom panel. In the entire transition region the ratio is consistent with a constant, because the peak in χ B 4 (T ) is replicated in the temperature dependence of χ B 2 (T ). As opposed to the Taylor coefficients, κ BB 2 (T ) shows a very mild temperature dependence. Finally we remark that very similar equations have been already used in Ref. [14] to calculate "lines of constant physics" to O(µ 4 B ) order. In this reference the pressure, energy density and entropy were calculated using the Taylor method, and in a further step lines were drawn on the µ B − T phase diagram, where these quantities are constant in some normalization. The obtained κ 2 coefficients are closely related to ours.
Contrary to Ref. [14], we use Eq. (8) as the definition of a truncation scheme rather than to investigate a Taylor expanded result. In a way, we work in the opposite direction: we will use lattice data at zero and imaginary µ B to obtain the coefficients in Eq. (8), which then can be used to either calculate the Taylor coefficients or, even better, to extrapolate the equation of state at finite µ B with no reference to the Taylor coefficients themselves.
Also, we used imaginary chemical potentials not only to calculate the coefficients, but also to study the single observable first, on which the analysis is based. We base our description of the entire chemical potential-dependence of the QCD free energy function on χ B 1 (T,μ B )/μ B . It is essential to rely on one observable only, in order to guarantee thermodynamic consistency: entropy, pressure and energy density will obey the known thermodynamic relations (see Section V) only if they come from the same truncation scheme.

A. Lattice details
In this work, we use the lattice action and the parameters described in Ref. [39]. The action benefits from treelevel Symanzik improvement in the gauge sector and four levels of stout smearing for the staggered flavors. The up and down quarks are degenerate. The resulting light pair of quarks, as well as the strange and charm quarks assume their respective physical mass.
We performed simulations at µ B = 0 on 32 3 × 8, 40 3 × 10, 48 3 × 12 and 64 3 × 16 lattices in a temperature range of 130 − 300 MeV, and up to 500 MeV on larger volumes. These simulations were complemented at imaginary values of the chemical potential in the temperature range 135 − 245 MeV for the lattice resolutions N τ = 8, . . . , 12. The µ B = 0 data were simulated at µ S = 0, some of these ensembles were already used in Ref. [20].
In addition, we performed a high-statistics run on the cheap and coarse 24 3 ×8 lattice, mainly to produce Fig. 5. These data did not enter the continuum extrapolation. κ BB 2 (T ), shown in Fig. 5. To extract κ BB 4 (T ) using the same strategy, a precise result on χ 6 B (T ) would be necessary.
Instead, we utilize imaginary chemical potential simulations as follows. We perform simulations at imaginary values of the baryon chemical potential: withμ Q =μ S = 0. We simulate temperatures in the range T = 135 − 245 MeV. For each temperature T and chemical potential µ B we determine the temperature T for which Eq. (7) holds, hence defining a function T (T,μ B ). Rearranging the terms in Eq (8), and in a similar way for the other observables, one can write: with the proxy quantity The proxy Π(T,μ B ) can be determined using lattice simulations. The relatively precisely known function χ B 2 (T, 0) is first interpolated in temperature. Afterwards, we only need to measure χ B 1 (T, µ B )/μ B on an ensemble with imaginary µ B . Equating this to χ B 2 (T , 0) gives us T (T, µ B ), while we have to take care of the propagation of the statistical errors.
Having determined Π(T,μ B ) for several imaginary chemical potentials and several lattice spacings for each given temperature, one can perform a polynomial fit inμ 2 B and obtain the expansion coefficients. This can be done separately for each lattice. However, we prefer to combine theμ B and continuum fits in one twodimensional fitting procedure. This combined fit is repeated for every temperature, in steps of 5 MeV.
In order to estimate the systematic uncertainties associated to our results, we perform a number of analyses at each temperature. There are several ambiguous points that need to be considered. Most obviously, one could choose to include the κ ij 6 (T ) term in the fit or not, or consider the fit of 1/Π instead. Also the range in imaginary µ B is arbitrary: we consider Im µ B ≤ 2.0 or Im µ B ≤ 2.4. When different lattice spacings are fitted together in a continuum extrapolation, one selects the bare parameters such that the ensembles correspond to the same physical temperature. This choice is, however, ambiguous too, as the scale setting may be based on various observables. In our case, we consider f π or w 0 to this purpose. As we mentioned before, χ B 2 (T, 0) is subject to an interpolation, performed through basis splines. The same is true for χ B 1 (T,μ B ) at finite chemical potentials. Since the location of the node points is also arbitrary, we include three versions at µ B = 0 and two at imaginary µ B , each with perfect fit quality. Finally, in the continuum extrapolation the coarsest lattice, 32 3 ×8, has either been used or dropped. The listed options can be considered in arbitrary combinations. In total we carry out all 144 fits to perform a continuum extrapolation of κ ij 2 (T ) and κ ij 4 (T ). After dropping the fits with a Q-value below a percent, we use uniform weights to produce histogram out of these (somewhat less than) 144 results for each temperature. The width of the histogram defines the systematic error. In the plots we show combined errors, where we assume that statistical and systematic errors add up in quadrature. This systematic error estimation procedure has been used and described in more detail in several of our works, most recently in Ref. [20].
In the top panel of Fig. 6 we show the results of the temperature-by-temperature fit procedure for the parameters κ BB 2 (T ) and κ BB 4 (T ), alongside the corresponding expectations from the Hadron Resonance Gas (HRG) model. We find that, within errors, κ BB 2 (T ) has hardly any dependence on the temperature, while κ BB 4 (T ) is everywhere consistent with zero at our current level of precision. Nonetheless, a clear separation of almost one order of magnitude appears between these two coefficients. We also note that good agreement with the HRG results is found up to at least T = 160 MeV.
In the central and bottom panels of Fig. 6 we show our results for the parameters κ BS 2 , κ BS 4 and κ SS 2 , κ SS 4 respectively, together with their HRG determinations. While for χ BS 11 barely any temperature dependence is observed, as in the case of χ B 2 , for χ S 2 a much stronger T -dependence is clearly visible. As in the case of κ BB 4 , κ BS 4 is consistent with zero throughout the whole temperature range we consider. On the other hand, κ SS 4 rises above zero for temperatures T > ∼ 220 MeV. Note that the error bars in these plots are highly correlated. The correlation is mostly systematic. The apparent 'waves' are often statistically not as significant as it seems at a first glance. For comparison we refer to the direct result on our coarsest lattice in Fig. 5, where these 'waves' are absent. Before moving on to calculate thermodynamic observables at finite real chemical potential, we construct a smoother version of our final results for κ ij 2 and κ ij 4 , in order to limit the influence of numerical effects on the final observables, but also to facilitate the temperature derivatives that enter the entropy formula. Thanks to the very mild dependence of the coefficients on the temperature, we perform a polynomial fit of order 5 for the κ ij 2 , and of order 2 for the even less T -dependent κ ij 4 . The very good fit qualities did not motivate higher order polynomials, using fourth or sixth order for κ ij 2 hardly changes the result. In order to stabilize the low-temperature behavior, we include in the fit two points from the HRG model for T = 120, 130 MeV, to which we associate an uncertainty of 5% for κ ij 2 and of 300% for κ ij 4 . The choice of these particular values for the uncertainties is uniquely guided by the necessity of placing a constraint on the low-T behavior, while avoiding to drive the fit too strongly. For this reason, these arbitrary errors are chosen to be smaller, but comparable to the lattice ones. We note that the fits we perform take fully into account the correlations between results at different temperatures, systematic as well as statistical. Thus we encode all errors into the (correlated) errors of the coefficients of a polynomial.
In Fig. 7 we show the results of the fits in darker color,  with the fitted data in lighter shades. The HRG points included in the fit are shown as well.
C. Continuum result on χ B 2 (T ) and its temperature derivative In order to determine the value of thermodynamic quantities at finite real chemical potential, a continuum result for the basic quantity χ B 2 (T, 0) is required in addition to κ BB 2 (T ) and κ BB 4 (T ). For some observables like the entropy, the temperature derivative of χ B 2 (T ) has to be calculated as well. We have already published this quantity in Ref. [39], but not its temperature derivatives. To make this work self-contained, we update this determination with the updated statistics using lattices up to 64 3 × 16.
We have divided the temperature range into two parts: the transition region and the higher, nearperturbative temperatures. For the lower temperature part, we interpolate using basis splines using the T ∈ [130 MeV, 300 MeV] data points. For the quantity T dχ B 2 (T )/dT we extended the data range with the HRG curve, so that the numerical derivative has a level arm at the lowest temperatures. The basis splines are cubic splines in the given temperature range. We fit the splines in temperature and in 1/N 2 τ in one step. The fitted function in this range is thus: where b i (t j ) = δ ij and t j are the knots for the spline with j = 1 . . . n. n and the t j set are the only arbitrary inputs. We combine the result from four sets of values for them, so that we can estimate the systematics. Only the lattices with N τ = 10, 12 and 16 enter the continuum extrapolation. For χ B 2 (T, 0; N τ ) in Eq. (18), the tree-level improvement was already applied [39]. The lattice resolutions and the spline knots are selected such that the fit qualities (Q value) are distributed as expected. We also took into account the scale setting ambiguity by combining the results with w 0 − and f π −based scale settings. The systematics of the T -derivative T dχ B 2 (T, 0)/dT and of χ B 2 (T, 0) were determined separately. In the high temperature regime, the smooth monotonic behavior is not well described with cubic splines. Instead we performed a high order polynomial fit in the inverse temperature 1/T . The known analytical form of the convergence to the Stefan-Boltzmann limit is, of course, not this polynomial. We do not wish to enforce the perturbative behavior at the intermediate temperatures that we describe. Thus, the constant in the 1/T description is not exactly the Stefan-Boltzmann limit, and for this reason, we cannot regard this as a basis for an extrapolation. However, this simple approach allows the interpolation and the calculation of the derivative. We used lattice data in the range T = 180 − 450 MeV.
High-and low-temperature regions overlap and the two fitting methods give consistent results between T = 200 − 280 MeV for both quantities. Thus, we simply concatenate the resulting functions at T = 260 MeV. We show the final version of the T dχ B 2 (T, 0)/dT and χ B 2 (T, 0) functions in Fig. 8.

V. THERMODYNAMICS AT REAL CHEMICAL POTENTIAL
Once n B is determined, we have everything we need to extract the other thermodynamic quantities. The integration constant for the pressure is obviously the pressure itself at µ B = 0.
We note here that, on the lattice, we always deal with dimensionless thermodynamic quantities, which correspond to the physical ones divided by suitable powers of the temperature. E.g., the dimensionful baryon density is n B = T 3n B (we will hereafter use the hat to indicate dimensionless quantities).
From the baryon densityn B (μ B , T ), the pressure is obtained through simple integration: The entropy density is defined as: For dimensionless quantities: where in the last step we converted the derivative at constant µ B into a derivative at constantμ B .
The T -derivative of the pressure involves the chain rule is calculated at µ B = 0 as already described.
The dimensionless expression for the energy densitŷ = /T 4 follows as: Finally, we present our results for the finite real chemical potential extrapolation of several thermodynamic quantities. The various panels of Fig. 9 show the baryon density, pressure, entropy, energy density, strangeness density and χ S 2 forμ B = 0−3.5. Alongside our results, we show predictions from the HRG model for T < 150 MeV, which we find in very good agreement with our extrapolation for all observables, at all values of the chemical potential.
We also note that in all cases, the observables do not suffer from pathological behavior. The uncertainties are under control for our range of chemical potentials, which highly improves on the results currently achievable via Taylor expansion.
We devote the two panels of Fig. 10 to the comparison of our results for the baryon density (left) and energy density (right) to the simplified case where κ BB 4 is neglected. We can appreciate how the inclusion of the next-to-leading-order parameter came at the cost of an increased uncertainty at larger chemical potential. This does not come unexpected, as we saw from our results that κ BB 4 was compatible with zero at all temperatures. In the case of the energy density, which is dominated by the µ B = 0 contribution, hardly any effect is visible.

VI. CONCLUSIONS
In this work, we proposed an alternative summation scheme for the equation of state of QCD at finite real chemical potential, designed to overcome the shortcomings which are characteristic of the Taylor expansion approach. Combining simulations at both zero and imaginary chemical potentials, we determined the LO and NLO parameters describing the chemical potential dependence of the baryon density, which we could then extrapolate to large real chemical potentials.
By combining this new element, and previously pub-lished results for the EoS at vanishing density, we could reconstruct all thermodynamic variables at chemical potentials as large as µ B /T = 3.5 with rather limited uncertainty. Systematic as well as statistical errors were considered in the analysis. These results, although still limited in precision at the level of the parameters κ ij 2 and κ ij 4 , suggest that the avenue we pursue in this work is rather promising for the description of QCD thermodynamics at finite chemical potential. Moreover, our procedure is systematically improvable with sufficient computing power, and might prove to be a better strategy than existing "canonical" approaches.
In this work we limited ourselves to the case where the strange and electric chemical potentials are set to zero. We reserve for future work the exploration of the phenomenologically relevant case of strangeness neutrality and fixed electric charge to baryon ratio.