Analytic continuation of the critical line: suggestions for QCD
Abstract
We perform a numerical study of the systematic effects involved in the determination of the critical line at real baryonic chemical potential by analytic continuation from results obtained at imaginary chemical potentials. We present results obtained in theories free of the sign problem, such as twocolor QCD with finite baryonic density and threecolor QCD with finite isospin chemical potential, and comment on general features which could be relevant also to the continuation of the critical line in real QCD at finite baryonic density.
pacs:
11.15.Ha, 12.38.Gc, 12.38.AwI Introduction
The study of QCD at nonzero baryonic density by numerical simulations on a spacetime lattice is plagued by the wellknown sign problem: the fermionic determinant is complex and the Monte Carlo sampling becomes unfeasible.
One of the possibilities to circumvent this problem is to perform Monte Carlo numerical simulations for imaginary values of the baryonic chemical potential, where the fermionic determinant is real and the sign problem is absent, and to infer the behavior at real chemical potential by analytic continuation.
The idea of formulating a theory at imaginary was first suggested in Ref. (1), while the effectiveness of the method of analytic continuation was pushed forward in Ref. (2). Since then, the method has been extensively applied to QCD (3); (4); (5); (6); (7); (8); (9); (10) and tested in QCDlike theories free of the sign problem (11); (12); (13); (15); (14); (16) and in spin models (17); (18).
The stateoftheart is the following:

the method is well founded and works fine within the limitations posed by the presence of nonanalyticities and by the periodicity of the theory with imaginary chemical potential (19);

the analytic continuation of the (pseudo)critical line on the temperature – chemical potential plane is well justified, but a careful test in twocolor QCD (15) has cast some doubts on its reliability.
In particular, the numerical analysis in twocolor QCD of Ref. (15) has shown that, while there is no doubt that an analytic function exists which interpolates numerical data for the pseudocritical couplings for both imaginary and real across , determining this function by an interpolation of data at imaginary could be misleading. Indeed, in the case of polynomial interpolations, there is a clear indication in twocolor QCD that nonlinear terms in play a relevant role at real , but are less visible at imaginary , thus calling for an accurate knowledge of the critical line there and, consequently, for very precise numerical data.
The above described scenario could well be peculiar to twocolor QCD and strongly depends on the choice of parameters of Ref. (15). Therefore, in this work we perform a systematic study of the analytic continuation of the critical line

in twocolor QCD with a fermionic mass different than in Ref. (15),

in another signfree theory, SU(3) with a nonzero density of isospin.
The aim of this study is to single out some general features of the analytic continuation of the critical line and to understand if and to what extent they can apply also to the physically relevant case of QCD.
Ii Analytic continuation of the critical line in twocolor QCD
As in Ref. (13), to which we refer for further details about our numerical setup, we have performed Monte Carlo simulations of the SU(2) gauge theory with degenerate staggered fermions. The algorithm adopted has been the usual exact algorithm described in Ref. (21), properly modified for the inclusion of a finite chemical potential.
The general strategy is the following: first the critical line is sampled by determining, for a set of values, the critical couplings for which the susceptibilities of a bulk observable exhibit a peak. Since the location of the critical line may well depend on the chosen observable, we have used three probe observables: the chiral condensate, the Polyakov loop and the plaquette. The values are chosen to lie inside the socalled first RobergeWeiss (RW) sector, i.e. in the range delimited by and (see Ref. (13) for a detailed discussion on the phase diagram in the plane).
Then, the critical line is guessed by interpolating the values of for only. The validity of the interpolation is evaluated by comparing its analytic continuation to the region with the direct determinations of the critical coupling in this region.
The procedure above is repeated for each of the three probe observables considered.
ii.1 ,
The SU(2) gauge theory with degenerate staggered fermions of mass has been studied on a lattice in Ref. (15). We found that

there is a very weak dependence of the values of on the probe observable;

there is no better interpolation of data at than a polynomial of the form ;

the extrapolation to the region overshoots the direct determinations of in that region, the discrepancy becoming larger and larger for increasing ;

data for for both and can be nicely interpolated by an analytic function (a polynomial of third order in works nicely), this excludes the option that the critical line be not smooth and possibly nonanalytic in the thermodynamic limit.
1.01  0.51  1.13  
0.97  0.21  0.20  
0.30  0.66  0.79  
0.92  2.10  0.73  
0.37  0.20  0.24  
0.30  1.50  0.58  
0.03  1.01  0.14  
1.04  1.87  0.74  
0.30  0.73  0.68 
observable  /d.o.f.  

1.61  
1.34  
1.40 
The conclusion we drew in Ref. (15) was that terms of order and play a relevant role at , but are less visible at , thus calling for an accurate knowledge of the critical line in all the first RW sector.
Verifying if this scenario is peculiar to SU(2) with the parameters choice as in Ref. (15) is the main motivation of the present work. For this reason, in the following subsection we repeat the same analysis of Ref. (15) in SU(2) with a different numerical set up, whereas in the next section we turn to a completely different theory.
ii.2 ,
We have considered SU(2) with degenerate staggered fermions of mass on a lattice. The reason for setting a larger value for the quark mass is simple: at any fixed temperature the critical value of gets larger, so that the critical line could become less curved in the physical region . This should reduce the importance of higher orders in in the description of the critical line by a polynomial.
As in Ref. (15), the critical line is sampled by looking for peaks in the susceptibilities of Polyakov loop, chiral condensate and plaquette. In Fig. 1 we show as example the behavior in of the susceptibility of the three observables at different values of , together with the Lorentzian interpolation of the peak. The statistics for each data point is about 1020k trajectories of one Molecular Dynamics (MD) time length.
All determinations of are summarized in Table 1, from which we see a mild dependence on the probe observable for ; at , instead, determinations from the Polyakov loop are systematically below those from the chiral condensate and the plaquette, which instead agree within errors.
As in Ref. (15), the best interpolation of data at is a polynomial of the form . In Table 2 we report the values of the fit parameters and their uncertainties as obtained with the MINUIT minimization code. The extrapolation to compares very well with the direct determinations of in that region, as shown in Figs. 2. In these figures, and in the following of the same kind, the band around the best fit (dotted line) represents the 95% confidence level (CL) region.
We can draw some partial conclusions about the analytic continuation of the critical line in SU(2):

for our given accuracy, there seems to be no room for quartic and sextic terms in in the interpolation of data for imaginary ;

the extrapolation to works definitely better for larger masses, i.e. away from the chiral limit.
Both these features could be peculiar to the SU(2) theory. This motivates the need to repeat a similar analysis in a different theory, sharing with SU(2) the property of being free of the sign problem.
Iii Analytic continuation of the critical line in threecolor QCD at finite isospin chemical potential
The SU(3) gauge theory with a finite density of isospin charge (22) is a theory in which the chemical potential is for half of the fermionic species and for the other half. The partition function, which is even in and depends only on , can be written as follows:
(1) 
where the integration is over gauge link variables, is the pure gauge action and the fermion matrix (we adopt a standard staggered discretization). This leads to a real and positive measure, because of the property , and therefore to a theory free of the sign problem. This theory is obviously closer to real QCD than twocolor QCD, being yet unphysical, since it implies a zero baryonic density, while in nature a nonzero isospin density is always accompanied by a nonzero baryonic density; moreover the isospin charge is not a conserved number in the real world.
Nevertheless, for our purposes this theory is very convenient since it provides us with another theoretical laboratory for the method of analytic continuation.
Similarly to SU(2) with finite baryonic density, at imaginary values of the chemical potential (which we will call simply from now on) the theory exhibits RWlike transition lines (we refer to the following subsection for a discussion about the determination of their position). The first RW sector is given by the strip .
In our numerical analysis, we consider finite isospin SU(3) with degenerate staggered fermions of mass on a lattice. Differently from SU(2) with finite baryonic density, the critical line in the temperature – chemical potential plane is a line of (strong) first order transitions, over all the investigated range of values, . This is one of the reasons for working on a smaller volume than in the SU(2) case (tunneling between the different phases would have been sampled with much more difficulty on a larger volume) and clearly emerges from the distribution on the thermal equilibrium ensemble of the values of observables like the (real part of) the Polyakov loop, the chiral condensate, and the plaquette across the transition (see Fig. 3 as an example showing the typical twopeak structures). As a further evidence of the first order nature of the transition, we show in Fig. 4 the Monte Carlo run history of the (real part of) the Polyakov loop at and , which exhibits tunneling events between the two phases every few thousands trajectories, on average. Typical statistics have been around 10K trajectories of 1 MD unit for each run, growing up to 100K trajectories for 23 values around , for each , in order to correctly sample the critical behavior at the transition. The critical is determined as the point where the two peaks have equal height and, in all the cases considered, this point turned out to be the same for all the adopted observables. In particular this means that, at the transition, both the chiral condensate and the (real part of) the Polyakov loop “jump” at the same point (see Fig. 5, as an example).
If one looks at the “gap” between the two peaks at the transition (determined as the difference between the values of an observable at the maxima of the right and left distribution peaks), it turns out that the dependence of this gap on in the considered range is slightly observabledependent, being roughly constant for the Polyakov loop and increasing for the chiral condensate and the spatial plaquette, up to a turning point at , where the strength of the transition seems to reach a maximum. In particular, no critical endpoint for the first order line is detectable, within the explored range. We notice that in Ref. (23), where the same theory has been explored even if with a slightly larger quark mass (), a critical endpoint was reported, but for a real chemical isospin chemical potential () which is well outside the range explored by us.
When combining the first RWlike transition line with the numerical determinations for the critical ’s (see Sec. III.2), the resulting phase diagram in the () looks like in Fig. 7.
iii.1 RWlike transitions at finite isospin density
In order to understand the unphysical, RWlike phase transition taking place in the high temperature region of the QCD phase diagram in the presence of an imaginary isospin chemical potential, let us first remember the origin of the usual RW transition in the presence of an imaginary baryon chemical potential. That is linked to the dynamics of the Polyakov loop. Let us define , where is the baryon chemical potential, and let us call , with , the eigenvalues of the untraced Polyakov loop, which satisfy the condition . The oneloop effective potential of the untraced Polyakov loop in presence of the baryon chemical potential is then given by (19)
(2) 
The first term is the gluonic contribution, which has three degenerate minima, corresponding to SU() matrices belonging to the center of the group. In absence of fermions that means spontaneous breaking of center symmetry. For the fermion contribution is minimum for , and the center element corresponding to a real Polyakov loop is selected. For the minimum of the fermion contribution moves and as crosses the absolute minimum of the potential moves to a different center element (RW transition).
In the case of an isospin chemical potential instead, the chemical potential is for half of the fermionic flavors and for the other half. The effective potential therefore looks as follows:
(3) 
In this case the fermion contribution is left unchanged as , i.e. it is symmetric under complex conjugation of the Polyakov loop. However it easily verified that as crosses a critical value , the two degenerate minima corresponding to complex center elements becomes deeper than then one corresponding to a real center element, hence for the conjugation symmetry of the Polyakov loop is spontaneously broken and one of the two complex minima is selected.
Notice that the above effective potential is computed at one loop and for zero mass fermions. Higher order and finite quark mass corrections do not affect the location of the RW transition in the case of an imaginary baryon chemical potential, which is fixed at by symmetry. The same corrections may however influence the location of for isospin chemical potential, hence may be quark mass and temperature dependent. Indeed, as shown in Fig. 7, numerical simulations locate approximately as the point where the unphysical RWlike line meets the analytic continuation of the physical transition line.
A more detailed analysis of the QCD phase diagram in the presence of imaginary baryon and isospin chemical potentials goes beyond our present purposes and will be presented elsewhere (24).
iii.2 Results for the critical line at finite isospin
In Table 3 we summarize our determinations of the critical couplings for the finite isospin SU(3) theory on a 8 lattice with fermionic mass =0.1.
We observe from the very beginning that data for for both and can be globally fitted by an analytic function (a polynomial of third order in nicely works). The fit parameters of the global fit are given in the last line of Table 4, while Fig. 8 shows how the fit compares with the data.
The question is if there are interpolations of the critical couplings at only, that, when continued to , agree with the critical couplings directly determined in the latter region.
/d.o.f.  

6.33  0.  
2.24  0.  
0.88  0.  
0.97  0.  
0.94  0.  
1.07  0.  
1.10  0.2025 
We have tried several kind of interpolations of the critical couplings at . At first, we have considered interpolations with polynomials up to order (see Table 4 for a summary of the resulting fit parameters). We can see that data at are precise enough to be sensitive to terms beyond the order ; indeed, a good /d.o.f. is not achieved before including terms up to the order , in agreement with the outcome of the global fit discussed above. The extrapolation to for the polynomial of order is shown in Fig. 9 (right), in comparison with that obtained from the polynomial of order (Fig. 9 (left)). The extrapolation agrees with direct determinations of , within the 95% CL band, only if nonlinear terms in , at least up to , are taken into account.
/d.o.f.  

2.04  
0.50  
0.49  
0.44  
0.74  
0.55  
0.46  
0.54  
0.65 
Then, we have considered interpolations with ratios of polynomials of order up to (see Table 5 for a summary of the resulting fit parameters). In all but one case we got good fits to the data at , but only two extrapolations to compare well with numerical data in that region: the ratio of a 4th to 6th order polynomial and the ratio of a 6th to 4th order polynomial (see Fig. 10). It is interesting to observe that the two interpolations which “work” have in common the number of parameters.
/d.o.f.  

  
0.39  
0.29  
0.22  
0.39  
0.39  
0.39  
1.83  
3.72 
Both kinds of fits considered so far have evidenced that the role of
terms of order larger than cannot be neglected. Since the data
more sensitive to these terms are those farther from , while data
closer to should be enough to fix the coefficient of the
term in a polynomial interpolation, we followed a different strategy:
we performed a fit with a polynomial of the form
in the region
and varied in order to find the largest possible interval
in which the fit gives both a good /d.o.f. and a stable value
(within errors) for the linear coefficient ; this fixes the value of
the parameter with some accuracy (see Table 6)
At last, we have attempted a fit strategy never used so far in our studies on the present subject and, to our knowledge, nowhere in the literature on the analytic continuation from imaginary chemical potential. The idea is to write the interpolating function in physical units and to deduce from it the functional dependence of on , after establishing a suitable correspondence between physical and lattice units. The natural, dimensionless variables of our theory are , where is the critical temperature at zero chemical potential, and . The question that we want to answer is if fitting directly the dependence of on may lead to increased predictivity for analytic continuation. We shall name this kind of fits, made directly in terms of the dimensionless physical variables of the theory, as “physical” fits.
While is one of the dimensionless variables used in our
simulations, is not and must be deduced from the relation
,
where is the number of lattice sites in the temporal direction
and is the lattice spacing at a given
We have tried several different fitting functions and report two cases which work particularly well. The first is given by the following 3parameter function:
(4) 
leading to the following implicit relation between and :
(5)  
In Fig. 12(left) we compare the values of obtained using our determinations for with the interpolation according to the fit function (5) to data in the region : one can see that the extrapolation to the region behaves very well. In Figs. 12(right) we have rescaled the vertical axis in order to give the behavior of versus the isospin chemical potential. The values of the fit parameters are
(6) 
with /d.o.f.=1.33. We observe that is in nice agreement with the direct determination (see Table 3).
As an alternative function for the shape of the critical line, we have also tried the following
(7) 
which explicitly encodes the expected periodicity of the partition function for imaginary . The fit to data at imaginary is very good and its extrapolation to the real chemical potential side compares impressively well with data (see Fig. 13). The resulting fit parameters are
(8) 
with /d.o.f.=0.39. This function is a good candidate to parametrize the critical line for small values of .
In both cases, Eqs. (5) and 7, the physical fit has worked very well and with a reduced number of parameters with respect to our previous fits, leading to increased predictivity and consistency with data at real chemical potentials. In both cases one can easily check that the adopted functions are not appropriate for a continuation of the critical line down to the axis, but this is not the aim of our study since such extrapolation would be questionable anyway.
We have used other functional forms for , inspired by the GrossNeveu model (26) and by the random matrix theory (27), and found that they allow one to interpolate data at , but their extrapolation does not match data at .
Finally, in order to compare the efficiency of the different predictions of the critical line, we present in Table 7 the extrapolated value of at a common value of the real chemical potential () for all the interpolation methods adopted. We can see that the last two methods listed (constrained and physical fit) have a narrower error band, the physical fit having a better agreement with data.
Iv Discussion
The present study is part of a larger project which, based on the investigation of QCDlike theories which are free of the sign problem, is aimed at testing the validity of the method of analytic continuation and at improving its predictivity, in view of its application to real QCD.
In particular, we have presented results concerning the analytic continuation of the critical line in 2color QCD and in QCD with a finite density of isospin charge. We have detected some common features and developed some strategies, which could apply and be useful also for real QCD at finite baryon density. Let us briefly summarize them.

Nonlinear terms in the dependence of the pseudocritical coupling on in general cannot be neglected. A polynomial of order seems to be sufficient in all explored cases. The prediction for the pseudocritical couplings at real chemical potentials may be wrong if data at imaginary are fitted according to a linear dependence.

The coefficients of the linear and nonlinear terms in in a Taylor expansion of are all negative. That often implies subtle cancellations of nonlinear terms at imaginary chemical potentials () in the region available for analytic continuation (first RW sector). The detection of such terms, from simulations at only, may be difficult and requires an extremely high accuracy. As a matter of fact, the simple use of a sixth order polynomial to fit data at imaginary leads to poor predictivity, which is slightly improved if ratio of polynomials are used instead.

An increased predictivity is achieved if the following strategy is adopted: first, the linear term in is fixed from data at small values of only. Then a nonlinear fit (e.g. with a sixth order polynomial) is performed, keeping the term constrained.

We have proposed a new, alternative ansatz to parametrize the critical line directly in physical units in the plane (instead than in the plane) and given two explicit realizations. We have shown that this “physical” ansatz provides a very good description of the critical line, moreover with a reduced number of parameters, and leads to an increased predictivity, comparable to that achieved by the “constrained” fit.
We plan to apply our experience to the determination of the critical line for real QCD in the near future. The lesson from the present study is that the preferred extrapolation methods to be considered in that case are a nonlinear polynomial, at least of sixth order in , in which the term is constrained, or rather the new “physical” ansatz for the critical line that we have proposed. The two methods have shown comparable predictivity. The former requires dedicated simulations at small chemical potentials, in order to better fix the term, which could also be fixed by using alternative methods like reweighting or Taylor expansion. The latter is more appealing since it describes the critical line with a reduced number of parameters: this is very likely at the origin of its increased predictivity and is probably related to the fact that one is looking for a direct relation between two physical quantities ( and ) in this case, without going through bare lattice quantities.
V Acknowledgments
We acknowledge the use of the computer facilities at the INFN apeNEXT Computing Center in Rome and of the PC clusters of CNAFBologna and INFNBari. We thank Ph. de Forcrand and M.P. Lombardo for very useful discussions.
Footnotes
 A similar strategy has been followed also in Ref. (25).
 Strictly speaking the lattice spacing depends also on the bare quark mass, which in our runs slightly changes as we change since we fix . However in the following evaluation, which is only based on the perturbative 2loop function, we shall neglect such dependence.
References
 M.G. Alford, A. Kapustin, and F. Wilczek, Phys. Rev. D 59, 054502 (1999).
 M.P. Lombardo, Nucl. Phys. Proc. Suppl. 83, 375 (2000).
 Ph. de Forcrand and O. Philipsen, Nucl. Phys. B 642, 290 (2002); Nucl. Phys. B 673, 170 (2003).
 M. D’Elia and M.P. Lombardo, Phys. Rev. D 67, 014505 (2003); Phys. Rev. D 70, 074509 (2004).
 V. Azcoiti, G. Di Carlo, A. Galante and V. Laliena, Nucl. Phys. B 723, 77 (2005).
 H.S. Chen, X.Q. Luo, Phys. Rev. D 72, 034504 (2005).
 Ph. de Forcrand and O. Philipsen, JHEP 0701, 077 (2007); PoS LATTICE2007, 178 (2007).
 L.K. Wu, X.Q. Luo and H.S. Chen, Phys. Rev. D 76, 034505 (2007).
 M. D’Elia, F. Di Renzo, M.P. Lombardo, Phys. Rev. D 76, 114509 (2007).
 M. D’Elia and F. Sanfilippo, Phys. Rev. D 80, 014502 (2009).
 A. Hart, M. Laine, and O. Philipsen, Phys. Lett. B 505, 141 (2001).
 P. Giudice and A. Papa, Phys. Rev. D 69, 094509 (2004); Nucl. Phys. Proc. Suppl. 140, 529 (2005).
 P. Cea, L. Cosmai, M. D’Elia and A. Papa, JHEP 0702, 066 (2007); PoS LAT2006, 143 (2006).
 S. Conradi, M. D’Elia, Phys. Rev. D 76, 074501 (2007).
 P. Cea, L. Cosmai, M. D’Elia and A. Papa, Phys. Rev. D 77, 051501 (2008); PoS LATTICE 2007, 214 (2007).
 Y. Shinno and H. Yoneyama, arXiv:0903.0922 [heplat].
 S. Kim, P. de Forcrand, S. Kratochvila, and T. Takaishi, PoS LAT2005, 166 (2006).
 F. Karbstein, M. Thies, Phys. Rev. D 75, 025003 (2007).
 A. Roberge, N. Weiss, Nucl. Phys. B 275, 734 (1986).
 M.P. Lombardo, PoS LAT2005, 168 (2006).
 S. A. Gottlieb, W. Liu, D. Toussaint, R. L. Renken, and R. L. Sugar, Phys. Rev. D 35, 2531 (1987).
 D.T. Son and M.A. Stephanov, Phys. Rev. Lett. 86, 592 (2001); J.B. Kogut and D.K. Sinclair, Phys. Rev. D 66, 034505 (2002).
 Ph. de Forcrand, M.A Stephanov, U. Wenger, PoS LATTICE 2007, 237 (2007), arXiv:0711.0023 [heplat].
 M. D’Elia, C. Manneschi, F. Sanfilippo, in progress.
 P. de Forcrand and O. Philipsen, JHEP 0811, 012 (2008).
 K.G. Klimenko, Z. Phys. C 37, 457 (1988); S. Hands, A. Kocic and J.B. Kogut, Nucl. Phys. B 390, 355 (1993).
 M.A. Halasz, A.D. Jackson, R.E. Shrock, M.A. Stephanov and J.J.M. Verbaarschot, Phys. Rev. D 58, 096007 (1998).