# Thermodynamics of two flavor QCD from imaginary chemical potentials

###### Abstract

We study QCD thermodynamics in presence of two independent imaginary chemical potentials coupled to two degenerate flavors of staggered quarks. Analytic continuation is used to determine non-linear susceptibilities, to test the Hadron Resonance Gas (HRG) model below the zero density critical temperature, , and to determine the average phase factor of the fermion determinant. Deviations from HRG predictions, of the order of a few percent, are clearly visible for temperatures . The determination of non-linear susceptibilities, using different interpolating functions for analytic continuation, gives consistent results and in agreement with Taylor expansion computations, apart from some systematic effects at or right above . Results for the average phase factor are compared with the predictions of Chiral Perturbation Theory; below we are able to distinguish the contribution of different hadron states, which is positive (i.e. tends to mitigate the sign problem) in the case of baryons.

###### pacs:

11.15.Ha, 12.38.Aw, 12.38.Mh## I Introduction

The study of QCD at finite temperature and baryon density has increasing phenomenological interest related to the physics of heavy ion experiments and compact astrophysical objects. The main open questions regard the location and nature of phase transitions in the QCD phase diagram, as well as the properties of strongly interacting matter around the transitions. A reliable answer to these questions requires at treatment of QCD at a non-perturbative level: unfortunately lattice QCD simulations, which are the only available tool for a non-perturbative study of the theory based on first principles, are not possible at finite baryon chemical potential, because of the well known sign problem: the QCD fermion determinant becomes complex and the probability interpretation of the QCD Euclidean action, necessary for standard importance sampling Monte-Carlo, is lost.

A number of strategies have been developed to partially circumvent that problem, like reweighting techniques glasgow ; fodor ; density , the use of an imaginary chemical potential either for analytic continuation muim ; immu_dl ; azcoiti ; chen ; giudice ; cea ; sqgp ; conradi ; cea2 or for reconstructing the canonical partition function rw ; cano1 ; cano2 , Taylor expansion techniques taylor1 ; taylor2 ; gagu1 ; gagu2 ; gagu3 and non-relativistic expansions hmass1 ; hmass2 ; hmass3 .

The aim of the present work is that of exploiting the method of analytic continuation from an imaginary chemical potential to study the properties of hadronic matter around the deconfinement transition in QCD with two light flavors (). As an improvement with respect to previous studies based on analytic continuation, we introduce to independent chemical potentials, and , coupled to the two different quark flavors. That is equivalent to the introduction of two independent chemical potentials, and , coupled respectively to the baryon and to the isospin charges and .

Our strategy will be to determine the dependence of the free energy on the two chemical potentials, apart from constant terms, by measuring its first derivatives with respect to and (quark number densities) for imaginary values of the two variables, and by then fitting them by suitable functions, to be continued within proper analyticity domains.

One of our aims is the study of generalized susceptibilities with respect to different conserved charges of the model (baryonic, isospin). These quantities are of significant phenomenological interest and have been determined till now mostly by the Taylor expansion method. We shall compare our results with those obtained by previous studies and comment on the efficiency and systematic effects of analytic continuation. In the confined region, i.e. below the critical temperature , we shall be able to perform a high precision test of the Hadron Resonance Gas (HRG) model, leading to the uncover of violations close to . Finally, the knowledge of the dependence of the free energy on the two independent chemical potentials will allow us a study of the average phase factor, which gives a direct measurement of the severeness of the sign problem.

Our study is made for QCD with two flavors of unimproved staggered quarks and is based on a standard RHMC algorithm. The choice of parameters is taken from Ref. gagu2 . The paper is organized as follows: In Section II we describe the model that we have investigated as well as the relevant physical observables; we also discuss the symmetries of the model, which are important for the choice of the free energy interpolating functions to be used for analytic continuation. In Section III we report the technical details of our numerical simulations. In Section IV we present results obtained below and compare them to the predictions of the HRG model. In Section V we report results obtained above . In Section VI and VII we discuss results obtained respectively for generalized susceptibilities and for the analytic continuation of the average phase factor. Finally, in Section VII, we draw our conclusions.

## Ii QCD with two independent chemical potentials and analytic continuation.

QCD with two continuum degenerate flavors is described, in the (rooted) staggered fermion discretization of the theory, by the following partition function

(1) |

where is the discretized pure gauge action (standard Wilson plaquette action in our case) and is the staggered fermion matrix describing 4 continuum flavors. Periodic (antiperiodic) boundary conditions are assumed for gauge (fermion) fields along the Euclidean time direction.

The introduction of two independent chemical potentials, and , coupled to the number operators of each quark family leads to the following expression for the grand canonical partition function:

(2) |

where the fermion matrix in the standard staggered formulation at finite chemical potential reads:

(3) | |||||

Here and refer to lattice sites, is a unit vector on the lattice, are staggered phases; and are respectively the chemical potential and the quark mass in lattice units.

The two chemical potentials can be rewritten in terms of a quark number chemical potential (or equivalently a baryon chemical potential ) and of an isospin chemical potential .

While the original theory is invariant under both charge conjugation and isospin rotations, the theory in presence of finite chemical potentials obviously is not. However the original invariance is reflected in the fact that the free energy must be an even function of and separately, or equivalently it must be invariant under the two following transformations and , which are easily verified to be symmetries of the partition function in Eq. (2). That places strong constraints on its possible functional dependence.

In presence of a finite chemical potential becomes complex and . Therefore, apart from the case (), the integrand in Eq. (2) is complex and cannot be interpreted as a probability distribution over gauge fields, so that standard importance sampling techniques cannot be applied (sign problem).

Positivity is recovered if the chemical potentials and are taken as purely imaginary: in this case numerical simulations are feasible and results can be used to fit the functional dependence of relevant observables.

Due to the above mentioned symmetries of the free energy, analytic continuation is actually a continuation from negative to positive values of and . Of course it is expected to be applicable as long as no phase transitions are met along the continuation path.

It is convenient for the following discussion to introduce the variables

and

where is the number of lattice sites in the temporal direction.

It can be easily shown that the introduction of an imaginary chemical potential is equivalent to a twist in the temporal boundary conditions for fermions by an angle . Hence both determinants appearing in Eq. (2) are periodic functions, respectively of and , with period , so that the free energy itself is a periodic function of these variables.

In terms of and that means again periodicity with period in both variables, plus invariance under . However, following the argument given by Roberge and Weiss in Ref. rw , it is possible to prove that a transformation , where is the number of colors and is an integer, can be cancelled by a change of variables in the functional integration in which all temporal links at a given time slice get multiplied by a center element (center transformation). Hence the free energy is expected to be a periodic function of with period instead of ( in our case). An analogous change of variables does not work for translations in , which rotate the link variables appearing in each determinant in a different way, therefore the period in is really .

For temperatures below the zero density critical temperature, , no phase transitions are expected, as in the case, in the whole plane. Therefore, due to the discussed periodicity and required symmetries, the most natural parametrization of the free energy is in terms of a trigonometric series as follows:

(4) |

with and both integers; moreover and must have the same parity because of the invariance under . Further constraints on the number of terms appearing in Eq. (4) may be predicted by particular effective models of strong interactions below , like for instance the HRG model to be discussed in Section IV. In such regime, information valid for analytic continuation can be gathered in the whole plane.

For we expect instead phase transitions in the plane, corresponding either to the continuation of the physical deconfinement transition or to the generalization of Roberge-Weiss (RW) transitions. Therefore a limited region around is available for the purpose of analytic continuation to real chemical potentials, and we shall write an expression for the free energy valid in that region which respects the predicted symmetries under and separately. In particular the free energy will be expressed as a polynomial like

(5) |

with non negative integers, or as a ratio of polynomials of the same kind

(6) |

The latter is an example of Chisholm approximant, i.e. the generalization to the case of two independent variables of usual Padè approximants, which have revealed to be better suited for analytic continuation in some cases mpl05 ; cea ; shinno .

Some of the quantities we are interested in are generalized susceptibilities with respect to the different chemical potentials, which for are defined as follows

(7) |

where is the pressure. Analogous susceptibilities are defined in terms of and

(8) |

The free energy symmetries discussed above imply precise constraints on the susceptibilities computed at zero chemical potentials. In particular we have only if and are both even, while if is even and .

Such quantities encode all relevant information about fluctuations of conserved charges, which are generally considered to be sensitive probes for the properties of the thermal medium produced in heavy ion collisions. They have been computed mostly in the Taylor expansion approach taylor1 ; taylor2 ; gagu1 ; gagu2 ; gagu3 , where they are expressed as average values at of operators which are more and more complex and computationally demanding as the order grows, since they require more and more matrix inversions. It is therefore sensible to explore the consistency and the efficiency of different strategies. In the analytic continuation approach we determine numerically the functional dependence, for imaginary values of the chemical potentials, of the first derivatives of the free energy. In terms of adimensional quantities, which are most conveniently determined on the lattice, they are given by

(9) |

where and are respectively the quark number and isospin charge operators, with

(10) | |||||

for . In terms of the susceptibilities defined in Eq. (7) and . Such first derivatives, which are purely imaginary for imaginary chemical potentials, can be measured quite efficiently (only one matrix inversion is needed for the noisy estimation of the trace) and, apart from constant terms, encode all information about the dependence of the free energy on . Information gathered at imaginary values of can then be analytically continued to real values of , in particular higher order derivatives at can be extracted.

In comparison to the Taylor expansion approach, the great advantage related to the much simpler observables can be compensated by the need for multiple simulations at different values of the chemical potentials. Moreover, this procedure involves some systematic dependence on the function chosen to interpolate data at imaginary ’s, which should be eventually checked by comparing results obtained with different functions. We shall compare trigometric expansions with polynomials below , polynomials with ratio of polynomials above .

## Iii Parameter details and numerical setup

Since we want to compare our results for the generalized susceptibilities with those obtained by the Taylor expansion approach, we have chosen for this study a subset of the parameters used in Ref. gagu2 , which is reported in Table 1. That corresponds to five different temperatures with a standard staggered lattice discretization on lattices and a fixed value (on the corresponding lattices) for the pion mass, MeV (actually and ). The critical temperature reported in Ref. gagu2 is MeV.

0.9 | 0.02778 | 5.26 | 95 | 2300 | |
---|---|---|---|---|---|

0.951 | 0.02631 | 5.275 | 95 | 2460 | |

1 | 0.025 | 5.2875 | 95 | 3500 | |

1.048 | 0.0238 | 5.30 | 24 | 3120 | |

1.25 | 0.02 | 5.35 | 77 | 2270 |

In particular, we have made simulations on a lattice using a RHMC algorithm. Our spatial size corresponds to about 6.6 inverse pion masses, hence finite size effects are not expected to be important.

For we have made simulations on a grid of about 100 different pairs , in the range : because of the above described periodicity, this surely contains all possible information available at imaginary chemical potentials (actually in a redundant way, which however is a benefit for checking the reliability of our statistical analysis). Since susceptibilities are calculated at null values of and , more points were taken in a restricted region around the origin, in order to perform fits of low-degrees polynomials in and around the origin easily. Morover, we have decided to perform a more accurate study of HRG model along the axis , therefore we have chosen further points there.

For we have performed a preliminary study aimed at finding the position of transition lines, with the purpose of delimiting the region at imaginary chemical potentials available for analytic continuation. Further information about this region are given in Section V.

For each we have produced about 2-3K thermalized trajectories of 1 Molecular Dynamics time length each. More details about the amount of pairs explored and average numbers of generated configurations are given in Table 1.

Quark densities have been measured by using noisy estimators. It is possible to minimize the total error of these observables (sum of statistical and noise fluctuations) at fixed simulation time by choosing an appropriate number of random vectors used for each noisy estimation. Assuming that noise and statistical fluctuations are independent of each other, the optimal number of random vectors to be used for each configurations is given by

(11) |

where is the variance of the observable (quark density) over different configurations, is the variance of the different estimates of the observable over a fixed configuration, is the time needed to generate a new configuration and is the time needed to perform one noisy estimate of the observable. We have measured those quantities in preliminary runs and we have found that, with our numerical setup, this number is around 30 for all explored parameter sets. Notice that Eq. (11) does not take into account the autocorrelation among configurations and thus overestimates ; we have however directly checked, by comparing different choices of , that the efficiency is almost stable for . We have always chosen in our production runs.

Simulations have been done on two PC farms in Genoa and in Bari. The complete collection of our data is not reported here, but is at disposal for interested readers.

## Iv Results at : precision test on the Hadron Resonance Gas model

The thermal medium below the critical temperature is generally believed to be well described as a gas of free hadron resonances (HRG model). This model provides a good description of thermal conditions at freeze-out redlich ; redlich2 ; andronic and has received theoretical support from lattice QCD simulations kareta . Deviations from the model have been recently detected close to in a lattice study based on the Taylor expansion method taylor2 .

In the HRG model the free energy is expressed as the sum of free particle energies. In particular, the free energy for species of spin , mass , baryon number and isospin , is given by

(12) | |||||

where , the upper (lower) sign applies to mesons (baryons) and

(13) |

The expression in Eq. (12) is an approximation in the case of unstable particles, for which an integration over a Breit-Wigner distribution in the particle mass would be more appropriate. The Bessel function is exponentially suppressed for large values of the argument, , hence for we can keep just the first term in the expansion, corresponding to the Boltzmann approximation in which quantum statistics effects are neglected. Summing up over all known particles and resonances and grouping together all charge conjugation and isospin partners we get

(14) | |||||

where and

Such prediction is easily continued to imaginary chemical potentials, where hyperbolic functions get transformed into trigonometric functions, in particular we have

(15) | |||||

(16) | |||||

(17) | |||||

where . The average quark densities are always purely imaginary for imaginary chemical potentials, for that reason we shall simply write and in the following, meaning implicitely that their imaginary part is taken.

Predictions from the HRG model to be tested in lattice QCD simulations can be classified as follows:

1) The free energy has a particularly simple form since, on the basis of known hadron resonances, only , , , are different from zero in previous equations. That means a further strong restriction on the expected form of the free energy at low temperatures: a necessary condition for the HRG model to be valid is that only the few lowest terms of the Fourier expansion in Eq. (4) give contribution;

2) Also the numerical values of the coefficients can be predicted from the known experimental resonance mass spectrum.

Latter prediction is easily affected by lattice artifacts and by the unphysical quark masses used in simulations, which change the actual hadron spectrum on the lattice. The former, instead, is expected to be more robust and less sensitive to discretization details. The method of analytic continuation is particularly well suited for lattice QCD tests of the HRG model, since it gathers information, below , from the whole range of possible imaginary chemical potentials, so that the number of terms actually contributing to the Fourier expansion in Eq. (4) can be checked with great precision: this idea has been followed in earlier studies limited to the axis immu_dl ; cano1 , in which the presence, within errors, of a single Fourier contribution, corresponding to , has been verified. In this respect the aim of our work is to extend such studies by increasing precision and by exploring also the region.

### iv.1

We start by discussing results obtained at . Let us first look at the axis: is zero in this case, while in general can be Fourier expanded as:

(18) |

and the HRG model predicts contribution only from the lowest harmonic, . Indeed a simple sine term, corresponding to , is perfectly compatible with our data, as showed in Fig. 1 and reported in Table 4. A second term with is therefore not necessary, at least within the precision of our data, even if a two sine fit leads to a smaller with a within three standard deviations (see again Table 4). As shown in Table 5, completely equivalent results are obtained if, instead of fitting our data, we compute the coefficients by explicit Fourier transform,

(19) |

where the integration is performed numerically by linear interpolation of consecutive data points.

Next we consider data for and obtained in the whole range of and explored, which are shown in Figs. 2 and 3, and try to fit them according to the expressions in Eqs. (16) and (17), considering more and more parameters till an acceptable value for the test is obtained. Fit results are reported in Table 6: a reasonable value of is obtained if a term with quantum numbers and is allowed for, besides those corresponding to usual meson () and baryons ( or ). Such term does not correspond to any known or even possible exotic hadron jaffe , but it is easily recognized as the first term, , neglected in Eq. (12) in the Boltzmann approximation in the case of pions: this is actually the first correction taking into account quantum statistics effects for pions, i.e. the fact that they are bosons, and corresponds to a two-pion exchange. With a pion mass as that used in our simulations, MeV, such term would mimic a coefficient , in very good agreement with the value obtained in our fit. Notice that terms with are negligible in our discretization setup, but would not be so, already at this temperature, in the case of physical pion masses. As for the data at , allowing for a term with leads to a lower value of , but is not strictly necessary, at least within the precision of our data.

Our conclusion is therefore that at numerical data do not contradict, within errors, the prediction coming from the HRG model and regarding the number of terms actually contributing to the free energy, apart from marginal evidence for a term which however is not strictly needed to fit data. Other deviations can be ascribed to the crudeness of the Boltzmann approximation for pions and are indeed well accounted for by the first neglected term.

Of course if one looks at the numerical value of the coefficients,
checking the agreement with experimental data is less trivial:
taking into account all
non-strange (since we are considering ) hadron resonances reported
in the Particle Data Book PDG ,
we would expect, for instance,
^{1}^{1}1More precisely we considered
all mesons of widely accepted existence, marked with a dot in
the meson summary table.,
which is roughly twice the value we have obtained ().
A more careful comparison is made using the unphysical pion and masses
realized in our lattice simulations ( MeV and
MeV gagu2 ): the coefficient becomes
including all resonances, taking into account just
pions and particles, and including
just pions (the errors here take roughly into account the
uncertainties given for the lattice estimate of the masses in Ref. gagu2 ),
i.e. much closer to our numerical result or even perfectly compatible
in the last case. We notice that,
since already masses are beyond the UV scale of our lattice
( MeV), it is perfectly reasonable that the
contribution from higher
resonances is not properly take into account.
That also clearly shows that a comparison of the numerical values
of the fitted coefficients with the HRG model prediction is
unavoidably affected by the systematics of the lattice discretization.

### iv.2

Once again we first look at results obtained for at , which are shown in Fig. 4. In this case two Fourier terms, corresponding to and , are necessary to fit our data. The second term is small, giving a contribution of the order of the total signal, but our data are precise enough to detect it; indeed a of order 2 is obtained if a single sine fit is tried (see Table 4).

In this case the presence of the term cannot be simply ascribed to
a violation of the Boltzmann approximation: assuming a mass
of order 1 GeV for the lightest baryon, the first neglected
term should lead to a signal a factor
smaller than what we get; moreover it should be negative,
as appropriate for a two-fermion exchange term.
The presence in the thermal medium of baryon-baryon bound states,
like deuterons, is a viable hypothesis: however assuming a mass
difference GeV between those states and
the lowest baryon states, one would expect a suppression
factor of the order
at this temperature, i.e. much
smaller than what we have obtained ^{2}^{2}2Notice however that also
for this states lattice artifacts due to the low UV cutoff,
MeV, could be important..
A simpler explanation is
that at this temperature corrections to the HRG model, induced by
non-trivial interactions close to the phase transition,
start to be important.

That is confirmed by analyzing the complete set of data for and as a function of and : fit results are reported in Table 6. Also in this case a term with is needed, but its value comes out to be about twice than expected from the first term neglected in the Boltzmann approximation for pions. In order to get a reasonable value for it is necessary to introduce also terms corresponding to (in agreement with results at ) and terms with and isospin up to . We interpret this again as a violation of the HRG model.

Regarding the numerical values obtained for the fitted coefficients, we obtain for instance , to be compared with if only pions are taken into account, including pions and mesons, including all meson resonances. The same considerations made for and regarding this comparison also apply here.

### iv.3

Finally let us briefly discuss results obtained at . Since at this temperature we stay in the confined phase as we switch an imaginary chemical potential, however small, it is still sensible to test predictions from the HRG model. However it is sufficient to look at results obtained for at (Fig. 5) to realize that violations to the model are important: in this case inclusion of the first three harmonics () is necessary to obtain a reasonable value for (see Table 4). This fact is confirmed by fits to the complete set of data for and which are reported in Table 6: the value decreases as more and more terms in the expansion in Eq. (4) are added.

In conclusion, within the current precision of our data, corrections to the HRG model are clearly detectable starting from .

Our best fits reported in Tab. 6, which are marked in the field by a *, provide us with a valid parametrization of the free energy (apart from a constant term). We shall make use of these parametrizations in the following Sections to derive generalized susceptibilities at and to study the analytic continuation of the average phase of the fermionic determinant. Systematic effects are expected at , where the of our best fit is somewhat bigger than 1.

In order to check for systematic effects related to the choice of the interpolating function we have also performed polynomial fits in a limited range of chemical potentials: our results are reported in Table 7. Fits chosen for analytic continuation are again marked by a * in the field.

## V Results at

The range of imaginary chemical potentials available for analytic continuation is limited, above , either by the presence of unphysical phase transitions related to center group dynamics (RW transitions) or by transitions corresponding to the analytic continuation of the deconfinement surface present at real chemical potentials. A full account of the high temperature phase structure in presence of two different imaginary chemical potentials will be given elsewhere delmansan ; in the present context we are just interested in the location of such transitions for the two temperatures explored, i.e. and . To that aim we have performed preliminary simulations on a small lattice to get a rough idea of the phase structure at these temperatures and thus delimit a safe region for analytic continuation, where to perform simulations on the larger lattice.

As for , in Fig. 6 we show the behaviour of the modulus of the Polyakov loop and of the chiral condensate as a function of at and as a function of at . It is clear that along both axes a transition is met where the system gets back into a phase with confinement and chiral symmetry breaking: at those points the system is crossing the analytic continuation of the pseudo-critical surface, present also at real chemical potentials. On the same symmetry grounds as for the deduction of general properties of the free energy in Section II, one expects that for small chemical potentials such pseudo-critical surface must be of the form

As clear from Fig. 6 the transition happens at approximately equal points along both axes (), i.e. . Also the observables (Polyakov loop and chiral condensate) seem to be, within a good approximation, universal functions of , where , at least not too far from the origin . Only imaginary chemical potentials strictly within the deconfined region can be considered for analytic continuation: Fig. 6 suggests us to take , with .

The phase structure is less trivial at . In Fig. 7 we plot the behaviour of the modulus of the Polyakov loop and of the chiral condensate as a function of in three cases: fixed , fixed and fixed . We observe again an approximate universal dependence on for relatively small values of this variable. Along the axis a transition is met, at , which clearly belongs to the pseudocritical confinement/deconfinement surface. Along the axis instead the system always stays in the deconfined phase and the Roberge-Weiss transition is met at where the system enters a different sector, as also apparent from the behaviour of the Polyakov loop phase shown in Fig. 8. What happens along the axis is less clear: presumably there one meets a pseudo-critical point close to the junction between the Roberge-Weiss transition and the pseudo-critical deconfinement surface. In this context we are only interested in delimiting a region safe for analytic continuation: from Fig. 7 it is clear that points with are surely contained in that region and this has been our conservative choice.

In this temperature regime we have tried to fit our results for and as a function of according to polynomials derived from the general expansion for the free energy given in Eq. (5) and truncated to a given order, or according to expressions derived from a parametrization of the free energy given in terms of ratios of polynomials as in Eq. (6).

At a fourth order polynomial provides a good fit, while coefficients are largely indetermined if a sixth order polynomial is used: not enough information can be extracted from the limited region available for analytic continuation. A marginally good fit is obtained with the ratio of two second order polynomials, but a fourth order polynomial at the numerator seems preferable.

At a sixth order polynomial or the ratio between fourth and second order polynomial are instead the best interpolating functions.

## Vi Generalized susceptibilities

Best fits to our data provide us with a parametrization for the dependence of the free energy on the chemical potentials, from which generalized susceptibilities can be extracted. Results can be considered reliable as long as different interpolations provide consistent results.

In Table 2 we report results obtained for , , and (definined in Eq. (8)) from free energy best fits marked by a * in the tables. Results obtained for , and are reported also in Figs. 11, 12 and 13 respectively, where they are compared with analogous results obtained using the Taylor expansion method in Ref. gagu2 .

The following general features can be observed. Different extrapolations provide always consistent results for and . A good agreement with Taylor expansion results can be observed as well, apart from the case.

For we observe a discrepancy between different interpolations only for and ; the agreement with Taylor expansion is less good around .

For different extrapolations disagree or are at most marginally compatible in the whole range of temperatures: with the current precision of our data, we cannot get reliable results for sixth or higher order susceptibilities.

In general, the comparison among different interpolation methods and with Taylor expansion results is good, apart from the region around . This in not unexpected: right above the region of imaginary chemical potentials usable for analytic continuation is small and restricted by the continuation of the pseudo-critical line, so that poor information is available. Moreover, right at we could not get best fits to the free energy dependence with a less than 1.5, therefore we do not have a completely satisfactory parametrization of the free energy for this temperature and systematic effects related to analytic continuation may be more important.

We have reported in Table 1 the total number of Dirac matrix multiplications needed in our numerical simulations at each temperature. We infer, from a rough estimate, that the effort for measurement purposes in our case (which is more or less half of the total) is approximately two orders of magnitude larger than what needed (again for measurement purposes) in Ref. gagu2 . The increased effort leads to corresponding smaller errors (about one order of magnitude) only for the lowest susceptibilities ( and ), while for higher order susceptibilities the Taylor expansion method seems to be more efficient. One has to consider, however, that our numerical simulations were not designed to be optimized for the computation of susceptibilities, and that in our case we obtain a complete parametrization of the free energy dependence in terms of and , which is usable for different purposes.

In Table 3 we report also results obtained for the susceptibilities with respect to quark and isospin chemical potentials and defined in Eq.(8). In Fig. 14 we show in particular the values of and for all temperatures, as obtained from polynomial fits: notice that is always larger than below , meaning that isospin charge fluctuations can be excited more easily (mainly in the form of pions) than baryon charge fluctuations below , while in the deconfined region the two susceptibilities become almost equal, as appropriate for a system made up mostly of quark-like degrees of freedom.

0.9 | 0.2925(20) | -0.0535(17) | 1.287(24) | 9.5(3) | |
---|---|---|---|---|---|

0.289(3) | -0.0588(24) | 1.17(7) | 5.6 1.2 | ||

gagu2 | 0.311(19) | -0.057(15) | 1.495(75) | 11.2 7.0 | |

0.951 | 0.439(4) | -0.058(3) | 2.32(8) | 22(2) | |

0.434(4) | -0.062(3) | 2.16(8) | 14(2) | ||

gagu2 | 0.423(21) | -0.080(17) | 3.16(26) | -29 11 | |

1 | 0.759(7) | -0.039(5) | 5.09(13) | 61(3) | |

0.734(7) | -0.060(5) | 4.27(13) | 31(2) | ||

gagu2 | 0.946(20) | -0.0331(72) | 6.51(20) | -5.3 10.7 | |

1.048 | 1.557(6) | -0.032(5) | 3.4(3) | - | |

1.557(7) | -0.033(6) | 3.3(4) | 1 24 | ||

gagu2 | 1.55(16) | -0.0385(98) | 4.33(23) | -69 16 | |

1.25 | 1.8470(12) | -0.0130(9) | 1.960(20) | 0.64(23) | |

1.8473(11) | -0.0121(7) | 1.968(16) | 2.78(25) | ||

gagu2 | 1.84(12) | -0.0138(85) | 2.181(31) | 5.5 1.7 |

0.9 | 0.478(6) | 0.692(4) | 4.92(25) | 4.05(7) | 1.94(3) | |
---|---|---|---|---|---|---|

0.461(6) | 0.696(9) | 4.15(17) | 4.1(4) | 1.76(14) | ||

0.951 | 0.762(9) | 0.993(10) | 8.2(3) | 8.3(5) | 3.45(11) | |

0.744(8) | 0.992(10) | 6.95(23) | 7.5(4) | 3.36(15) | ||

1 | 1.440(14) | 1.597(16) | 19.8(5) | 16.8(8) | 7.47(21) | |

1.348(13) | 1.589(18) | 13.9(3) | 15.6(8) | 6.46(23) | ||

1.048 | 3.052(17) | 3.178(15) | 7.5 | 12.8 | 5.5(4) | |

3.045(21) | 3.176(15) | 2(5) | 11(4) | 6.6(9) | ||

1.25 | 3.668(3) | 3.720(3) | 4.59(13) | 4.75(12) | 3.67(3) | |

3.671(3) | 3.7188(24) | 4.72(11) | 4.68(9) | 3.681(17) |

## Vii Phase of the fermionic determinant

As we have recalled in Section II, the complex phase of the fermion determinant, , hinders numerical simulations in presence of a real baryon chemical potential . The problem is however milder in case the fluctuations of the phase around zero, over the gauge configurations which are typical of the statistical ensemble, are small: in that case efficient numerical methods, like reweighting, can be used. A typical measure of the severeness of the sign problem is therefore given by the average of the phase factor (or some power of it), computed for convenience over the ensemble at finite isospin density. In particular in our case we can define:

(20) | |||||

As clear from Eq. (20), a way to determine is to take the average of the ratio of two determinants over the ensemble at real isospin chemical potential: that is feasible but computationally demanding, especially at large volumes. Studying the analytic continuation of at imaginary values of ,

(21) |

is an alternative: it has been shown splitt1 ; splitt2 that, in the full QCD case, the average phase factor is analytic around , and an efficient numerical method for the evaluation of the ratio of partition functions appearing in Eq. (21) has been proposed in Ref. conradi .

In the present context we adopt a much faster and cheaper approach: having measured and fitted first derivatives with respect to both chemical potentials, we have a complete knowledge, apart from constant terms, of the dependence of the free energy on and , so that computing the ratio in Eq. (21) is straightforward. Let us consider for instance the low temperature case, where we have used the HRG parametrization in Eq. (15) that we rewrite:

(22) | |||||

then

(23) |

Non-zero coefficients , apart from the constant which does not enter in the computation of the average phase factor, have been obtained by fitting our numerical data. The expression can then be easily continued to real chemical potentials obtaining:

(24) | |||||

The same procedure applies to other functional forms used in our fits: the comparison of different extrapolations to real chemical potentials based on different fitting functions, when available, gives a measure of the systematic effects involved in analytic continuation. Notice that in the case of the HRG parametrization we can distinguish the different contributions to the average phase factor, hence to the sign problem, coming from different particle species: this feature will be useful in our analysis.

In Figs. 15 and 16 we report, as a function of , results obtained respectively at and using HRG inspired and polynomial interpolations. Where visible, the two lines reported for each extrapolation delimit the 90% confidence level region and give an estimate of our uncertainties: a good agreement between HRG inspired and polynomial extrapolations can be appreciated.

It is interesting to make a direct comparison of our results with predictions coming from chiral perturbation theory (PT). The average phase factor has been computed to one loop order of PT in Ref. splitt3 . According to the results reported in Section VI of Ref. splitt3 , our spatial lattice size is big enough () to justify taking the thermodynamical limit at fixed T of the one loop PT result, which coincides with the prediction of a HRG model including only pions:

(25) |

with

(26) |

This prediction (assuming in our case MeV and MeV) is reported in Figs. 15 and 16 as a solid line. It is apparent that the agreement of PT with the analytic continuation of our data is not satisfactory. In particular analytic continuation provides a higher value for , meaning a milder sign problem. To better understand the origin of this discrepancy, we have tried to compute the average phase factor from our HRG model best fit, but neglecting all contributions to the free energy with , which cannot be taken into account by PT, i.e. taking only contributions from and in Eq. (24). Results are shown in Figs. 15 and 16: in this case the agreement with PT is almost perfect for , and acceptable for . This is expected since, as we have discussed in Section IV.1, the coefficients and obtained by our fits are compatible within errors, at , with those predicted if only pions are taken into account: of course that may be an accident and the contribution of higher meson resonances should be better understood.

Anyway, an outcome of our analysis, which is in agreement with HRG model expectations, is that contributions to the average phase factor coming from physical states with are significant and tend in general to make the sign problem less severe.

In Fig. 17 we report the analytic continuation of the average phase factor obtained at all temperatures from a polynomial fit: of course results reported in the figure must be intended to be valid for chemical potentials bounded, below , by the deconfinement critical line present at real chemical potentials. As expected, at fixed chemical potential the sign problem is much milder for . This can be put again in connection with the fact that states with , which are more easily created above , tend to mitigate the sign problem.

## Viii Conclusions

In this paper we have studied QCD thermodynamics, exploiting analytic continuation from two imaginary chemical potentials coupled to baryon and isospin charges. Simulations have been performed at five temperatures around the critical value MeV, using a lattice with a standard staggered action and a fixed pion mass MeV.

We have computed free energy first derivatives with respect to the chemical potentials (quark number densities) and interpolated them by suitable functions, in order to perform analytic continuation. In particular we have tested HRG predictions below , reconstructed generalized susceptibilities at zero chemical potentials and determined the analytic continuation of the average phase factor.

We have checked that HRG model predictions are in very good agreement with our numerical results for . Small but clearly detectable deviations start to be visible at , in agreement with similar results reported in Ref. taylor2 . They appear, in a HRG inspired parametrization of the free energy, as contributions from unphysical states with higher values of baryon or isospin charges, which are of the order of a few percent at and above 10 % at .

Regarding the computation of generalized susceptibilities, analytic continuation gives consistent results which are in agreement with those obtained by the Taylor expansion method, apart from temperatures in correspondence or right above , where the range of imaginary chemical potentials available for analytic continuation is small and larger systematic effects are expected. Poor information has been obtained for susceptibilities beyond sixth order.

We have obtained consistent determinations, by analytic continuation with different interpolating functions, of the average phase factor. In particular below , in the case of HRG inspired interpolations, we have been able to distinguish the contribution to the average phase factor coming from the different hadron states: results from analytic continuation are consistent with PT results, below , if one takes into account only meson contributions. Baryons give contributions to the average phase factor which in general tend to make the sign problem less severe. The sign problem is much milder for , and this can be put again in connection with the fact that states with , which are more easily created above , tend to mitigate the sign problem.

Our results should be refined and could be improved in several respects. Simulation closer to the continuum limit and possibly closer to the physical quark mass spectrum would clarify the comparison with HRG predicitions, as well as that with PT for the average phase factor. An improvement in the determination of generalized susceptibilities could be obtained by combining analytic continuation with other techniques: for instance fixing lowest order terms in a polynomial expansion by the Taylor expansion method or by reweighting could lead to enhanced predictivity for analytic continuation. We shall continue our investigation along those lines in the future.

## Acknowledgments

We thank F. Becattini, F. Karsch and K. Splittorf for very useful discussions, as well as R. Gavai and S. Gupta for very interesting discussions and for providing us with their numerical results for generalized non-linear susceptibilities. Numerical simulations have been performed on two PC farms in Genoa and in Bari provided by INFN.

0.1536(14) | - | - | 28/21 |

0.1514(15) | 0.0046(14) | - | 17/20 |

0.2413(25) | - | - | 42/21 |

0.2383(26) | 0.0102(22) | - | 21/20 |

0.3865(4) | - | - | 248/21 |

0.395(4) | 0.048(3) | - | 50/20 |

0.389(4) | 0.048(3) | 0.018(4) | 23/19 |

0.9 | 0.1521(11) | 0.0052(11) | - |
---|---|---|---|

0.951 | 0.2387(19) | 0.0101(17) | - |

1 | 0.392(3) | 0.0503(27) | 0.018(3) |

0.2284(11) | - | 0.0110(6) | 0.0202(3) | - | - | - | - | 284/187 |

0.2157(18) | 0.0050(6) | 0.0115(6) | 0.0198(3) | - | - | - | - | 206/186 |

0.2156(18) | 0.0051(6) | 0.0111(6) | 0.0197(3) | - | - | 0.00043(13) | - | 196/185 * |

0.2862(13) | - | 0.0199(7) | 0.0305(4) | - | - | - | - | 640/187 |

0.258(2) | 0.0114(6) | 0.0212(7) | 0.0292(4) | - | - | - | - | 281/186 |

0.257(2) | 0.0117(6) | 0.0203(8) | 0.0290(4) | - | - | 0.00084(18) | - | 259/185 |

0.256(2) | 0.0114(6) | 0.0210(8) | 0.0264(6) | 0.0017(3) | - | 0.00088(18) | - | 230/184 |

0.257(2) | 0.0106(7) | 0.0212(8) | 0.0265(6) | 0.0009(4) | 0.0006(2) | 0.00090(18) | - | 222/183 * |

0.3775(15) | - | 0.0412(7) | 0.0456(4) | - | - | - | - | 1798/187 |

0.3322(21) | 0.0219(7) | 0.0391(7) | 0.0465(4) | - | - | - | - | 808/186 |

0.3269(21) | 0.0225(7) | 0.0363(7) | 0.0464(4) | - | - | 0.00372(24) | - | 562/185 |

0.3184(22) | 0.0246(7) | 0.0349(7) | 0.0396(6) | 0.0053(3) | - | 0.00436(24) | - | 330/184 |

0.3208(22) | 0.0218(8) | 0.0342(7) | 0.0391(6) | 0.0038(4) | 0.0019(3) | 0.00445(24) | - | 288/183 |

0.3214(22) | 0.0220(8) | 0.0344(8) | 0.0393(6) | 0.0042(4) | 0.0015(3) | 0.0031(6) | 0.010(4) | 281/182 * |

0.34 | 0.479(3) | 0.1892(22) | - | - | - | - | - | - | - | 6014/108 |

0.34 | 0.659(5) | 0.412(4) | -2.66(7) | -1.22(4) | -2.40(4) | - | - | - | - | 237/105 |

0.34 | 0.696(9) | 0.461(6) | -4.1(4) | -1.76(14) | -4.15(17) | 23(6) | 8(3) | 12(3) | 30(3) | 113/101 * |

0.34 | 0.636(3) | 0.301(3) | - | - | - | - | - | - | - | 8483/108 |

0.34 | 0.897(5) | 0.651(5) | -3.69(6) | -2.00(5) | -3.76(5) | - | - | - | - | 388/105 |

0.34 | 0.992(10) | 0.744(8) | -7.5(4) | -3.36(15) | -6.95(23) | 62(7) | 28(4) | 26(4) | 54(4) | 125/101 * |

0.34 | 0.838(5) | 0.394(4) | - | - | - | - | - | - | - | 12955/108 |

0.34 | 1.340(9) | 1.099(9) | -5.75(8) | -3.33(8) | -6.67(7) | - | - | - | - | 972/105 |

0.34 | 1.589(18) | 1.348(13) | -15.6(8) | -6.46(23) | -13.9(3) | 157(13) | 65(6) | 47(5) | 121(5) | 239/101 * |

0.12 | 3.029(8) | -2.941(9) | - | - | - | - | - | - | - | 228/36 |

0.12 | 3.178(15) | 3.052(16) | -12.8 1.9 | -5.5(4) | -7.5 1.9 | - | - | - | - | 39/33 * |

0.12 | 3.178(24) | 3.05(3) | -10(6) | -7.6 1.5 | -4(6) | 639(1005) | 125(148) | 224(135) | -991(1037) | 34/29 |

0.3 | 3.2810(11) | 3.2438(12) | - | - | - | - | - | - | - | 170616/111 |

0.3 | 3.7156(18) | 3.6677(2) | -4.67(4) | -3.555(9) | -4.74(4) | - | - | - | - | 142/111 |

0.3 | 3.720(3) | 3.668(3) | -4.75(12) | -3.67(3) | -4.59(13) | 0(3) | 1.7(5) | 1.5(5) | 6(3) | 123/107 * |