arXiv:0807.1597 [hep-lat] Stochastic quantization at finite chemical potential

arXiv:0807.1597 [hep-lat]

Stochastic quantization at finite chemical potential

Gert Aarts and Ion-Olimpiu Stamatescu

Department of Physics, Swansea University, Swansea, United Kingdom
Institut für Theoretische Physik, Universität Heidelberg, Heidelberg, Germany
and FEST, Heidelberg, Germany
July 10, 2008

A nonperturbative lattice study of QCD at finite chemical potential is complicated due to the complex fermion determinant and the sign problem. Here we apply the method of stochastic quantization and complex Langevin dynamics to this problem. We present results for U(1) and SU(3) one link models and QCD at finite chemical potential using the hopping expansion. The phase of the determinant is studied in detail. Even in the region where the sign problem is severe, we find excellent agreement between the Langevin results and exact expressions, if available. We give a partial understanding of this in terms of classical flow diagrams and eigenvalues of the Fokker-Planck equation.

1 Introduction

A nonperturbative understanding of QCD at nonzero baryon density remains one of the outstanding problems in the theory of strong interactions. Besides the theoretical challenge, there is a clear phenomenological interest in pursuing these studies, due to the ongoing developments in heavy ion collision experiments, at RHIC, LHC and the planned FAIR facility at GSI.

The standard nonperturbative tool to study quarks and gluons, lattice QCD, cannot be applied in a straightforward manner, because the complexity of the fermion determinant prohibits the use of approaches based on importance sampling. This is commonly referred to as the sign problem. Since the start of the millenium a number of new methods has been devised. These include reweighting [1, 2, 3, 4], Taylor series expansion in [5, 6, 7, 8], imaginary chemical potential and analytical continuation [9, 10, 11, 12], and the use of the canonical ensemble [13, 14] and the density of states [15]. Except for the last two, these approaches are approximate by construction and aimed at describing the QCD phase diagram in the region of small chemical potential and temperatures around the crossover between the confined and the deconfined phase. In this paper we discuss an approach which is manifestly independent of those listed above: stochastic quantization and complex Langevin dynamics. How well this method will work is not known a priori. However, one of the findings of our study is that excellent agreement is found in the case of simple models, where comparison with results obtained differently is available. In particular we find that the range of applicability is not restricted to small chemical potential and, importantly, does not depend on the severity of the sign problem. The first results we present for lattice QCD at nonzero density are encouraging, although a systematic analysis has not yet been performed.

This paper is organized as follows. In the following section we briefly describe the idea behind stochastic quantization and the necessity to use complex Langevin dynamics in the case of nonzero chemical potential. In Secs. 3 and 4 we apply this technique to U(1) and SU(3) one link models. In both cases a comparison with exact results can be made. We study the phase of the determinant in detail. In the case of the U(1) model, we employ the possibility to analyse classical flow diagrams and the Fokker-Planck equation to gain further understanding of the results. In Sec. 5 we turn to QCD, using the full gauge dynamics but treating the fermion determinant in the hopping expansion. Our findings and outlook to the future are summarized in Sec. 6. The Appendix contains a brief discussion of the Fokker-Planck equation in Minkowski time for the one link U(1) model.

2 Stochastic quantization and complex Langevin dynamics

The main idea of stochastic quantization [16, 17] is that expectation values are obtained as equilibrium values of a stochastic process. To implement this, the system evolves in a fictitious time direction , subject to stochastic noise, i.e. it evolves according to Langevin dynamics. Consider for the moment a real scalar field in dimensions with a real euclidean action . The Langevin equation reads


where the noise satisfies


By equating expectation values, defined as


where is an arbitrary operator and the brackets on the left-hand side denote a noise average, it can be shown that the probability distribution satisfies the Fokker-Planck equation


In the case of a real action , the stationary solution of the Fokker-Planck equation, , will be reached in the large time limit , ensuring convergence of the Langevin dynamics to the correct equilibrium distribution. When the action is complex, as is the case in QCD at finite chemical potential, the situation is not so easy. It is still possible to consider Langevin dynamics based on Eq. (2.1) [18, 19, 20, 21]. However, due to the complex force on the right-hand side, fields will now be complex as well: . As a result, proofs of the convergence towards the (now complex) weight are problematic.

In the past, complex Langevin dynamics has been applied to effective three-dimensional spin models with complex actions, related to lattice QCD at finite in the limit of strong coupling and large fermion mass [22, 23, 24] (for applications to other models, see e.g. Ref. [25]). Our work has also partly been inspired by the recent application of stochastic quantization to solve nonequilibrium quantum field dynamics [26, 27, 28]. In that case the situation is even more severe. Nevertheless, numerical convergence towards a solution can be obtained under certain conditions, both for simple models and four-dimensional field theories. For an illustration we present some original results in the appendix.

Here we consider models with a partition function whose form is motivated by or derived from QCD at finite chemical potential. The QCD partition function reads


where is the bosonic action depending on the gauge links and is the complex fermion determinant, satisfying


Specifically, for Wilson fermions the fermion matrix has the schematic form


Here are lattice translations, , and is the hopping parameter. Chemical potential is introduced by multiplying the temporal links in the forward (backward) direction with () [29]. We use Eq. (2.7) as a guidance to construct the U(1) and SU(3) one link models considered next.

3 One link U(1) model

3.1 Complex Langevin dynamics

We consider a one link model with one degree of freedom, written as . The partition function is written suggestively as


where the “bosonic” part of the action reads


while the “determinant” is constructed by multiplying the forward (backward) link with (),


Due to the chemical potential, the determinant is complex and has the same property as the fermion determinant in QCD, i.e. . For an imaginary chemical potential , the determinant is real, as is the case in QCD.

Observables are defined as


In this model most expectation values can be evaluated analytically. We consider here the following observables:

  • Polyakov loop:


    where the partition function equals


    and are the modified Bessel functions of the first kind.

  • Conjugate Polyakov loop:


    At finite chemical potential, and are both real, but different.

  • Plaquette:


    Note that .

  • Density:


    At small chemical potential increases linearly with , while at large chemical potential exponentially fast.

We now aim to estimate these observables using numerical techniques. Due to the complexity of the determinant, they cannot be estimated using methods based on importance sampling. Instead, we attempt to obtain expectation values using stochastic quantization.

At nonzero chemical potential, the action is complex and it becomes necessary to consider complex Langevin dynamics. We write therefore , and consider the following complex Langevin equations


Here we have discretized Langevin time as , and the noise satisfies


The drift terms are given by


where the effective action reads


Explicitly, the drift terms are




Occasionally we will also consider this model by expanding in small , the hopping expansion, and take


This limit is motivated by the model of Heavy Dense Matter used in Ref. [36]. A direct application of our method to QCD in the hopping expansion is presented in Sec. 5.

In order to compute expectation values, also the observables have to be complexified. For example, after complexification , the plaquette reads


All operators we consider are now complex, with the real (imaginary) part being even (odd) under .

Figure 1: Real part of the Polyakov loop (left) and the conjugate Polyakov loop (right) as a function of for three values of at fixed . The lines are the analytical results, the symbols are obtained with complex Langevin dynamics.

The Langevin dynamics can be solved numerically. In Fig. 1 the real parts of the Polyakov loop and the conjugate Polyakov loop are shown as a function of for three values of at fixed . In Fig. 2 (left) the density is shown. The lines are the exact analytical results. The symbols are obtained with the stochastic quantization. We observe excellent agreement between the analytical and numerical results. For the results shown here and below, we have used Langevin stepsize and time steps. Errors are estimated with the jackknife procedure. The imaginary part of all observables shown here is consistent with zero within the error in the Langevin dynamics.333Analytically they are identically zero. This can be understood from the symmetries of the drift terms and the complexified operators, since the drift terms behave under as


while the imaginary parts are odd. Therefore, after averaging over the Langevin trajectory the expectation value is expected to reach zero within the error, which is what we observe. As an aside, we note that the symmetry of the drift terms under ,


relates positive and negative chemical potential.

Figure 2: Left: Real part of the density . Right: Real part of the plaquette versus . Results at positive (negative) have been obtained with complex (real) Langevin evolution.

At imaginary chemical potential the determinant is real, so that the complexification of the Langevin dynamics is not necessary. We demonstrate the smooth connection for results obtained at imaginary using real Langevin dynamics and results obtained at real using complex Langevin dynamics for the expectation value of the plaquette in Fig. 2 (right). Since the plaquette is even under , we show the result as a function of , so that the left side of the plot corresponds to imaginary chemical potential, while the right side corresponds to real chemical potential. On both sides excellent agreement with the analytical expression can be observed. We also note that the errors are comparable on both sides.

3.2 Phase of the determinant

At finite chemical potential the determinant is complex and can be written as


In order to assess the severity of the sign problem, we consider the phase of the determinant and study the behaviour of . An observable often used for this purpose [38, 39] is


where we used Eq. (2.6). At zero chemical potential the ratio equals one, while at large one finds in this model that


for nonzero . In expressing Eq. (3.23) as the expectation value obtained from the complex Langevin process, complex conjugation has to be performed after the complexification of the variables, as discussed above. In that case as a complex number is not the complex conjugate of . To avoid confusion we write in all relevant expressions. Notice that this implies that itself is also complex.

Figure 3: Left: Real part of . Right: Real part of .

In Fig. 3 (left) we show the real part of this observable as a function of . The imaginary part is again zero analytically and zero within the error in the Langevin process. The lines are obtained by numerical integration of the observable with the complex weight, while the symbols are again obtained from Langevin dynamics. We note again excellent agreement between the semi-analytical and the stochastic results. In particular, there seems to be no problem in accessing the region with larger where the average phase factor becomes very small. The numerical error is under control in the entire range. We find therefore that the sign problem does not appear to be a problem for this method in this model.444In QCD, the average phase factor is expected to go to zero exponentially fast in the thermodynamic limit.

Figure 4: Scatter plot of during the Langevin evolution for various values of at , . Note the different scale in the middle box.

In contrast to what could be inferred from the result above, expectation values of are not phase factors, due to the complexity of the action. This can be seen by considering


where the second and third equality follow from the cancelation of det in the definition of the expectation value and from being even in . We have also computed this observable using Langevin dynamics and the result is shown in Fig. 3 (right). For the Langevin parameters used here, we observe that the numerical estimate is consistent with 1, but with quite large errors when increases at small values of . We have found that increasing the Langevin time reduces the uncertainty. We conclude that at large chemical potential this ratio of determinants is the most sensitive and slowest converging observable we considered.

We observed above that the average phase factor becomes very small at large but that this does not manifest itself in all but one observable we consider. Insight into this feature can be gained by studying scatter plots of the phase factor during the Langevin process. In Fig. 4 we show the behaviour of during the Langevin evolution in the two-dimensional plane spanned by and . At zero chemical potential, and during the entire evolution. For increasing one finds more and more deviations from this, with an interesting structure at intermediate values of . Note that the resulting distribution is approximately invariant under reflection in , ensuring that the imaginary part of the expectation value vanishes within the error. Due to the wide distribution, the horizontal and vertical scales in the middle section of Fig. 4 are much larger than in the top and bottom part. However, the average phase factor remains well defined for all values of , as can be seen in Fig. 3. At large , the average phase factor becomes very small. However, the distribution is very narrow, see Fig. 4 (bottom). Therefore, although the average is close to zero, the error in the Langevin dynamics is very well under control.

3.3 Fixed points and classical flow

The excellent results obtained above can partly be motivated by the structure of the dynamics in the classical limit, i.e. in absence of the noise. As we demonstrate below, the classical flow and fixed point structure is easy to understand when and, most importantly, does not change qualitatively in the presence of nonzero chemical potential.

Classical fixed points are determined by the extrema of the classical (effective) action, i.e. by putting . We first consider the “bosonic” model and take . The drift terms are555For the bosonic model, there is of course no need to complexify the Langevin dynamics and one may take . This yields the same fixed points.


We see that there is one stable fixed point at and one unstable fixed point at . Moreover, the classical flow equation, , can be solved analytically, with the solution


where is the initial value at . We find therefore that the stable fixed point is reached for all , except when . On this line the solution reads


and the flow diverges to , except when starting precisely on the unstable fixed point . Note, however, that the noise in the direction will kick the dynamics of the unstable trajectories.

We now include the determinant, starting with the hopping expansion (3.18). Putting yields again one stable fixed point at and one unstable fixed point at , where


Note that in the strong coupling limit . We find therefore a simple modification of the bosonic model: in response to the chemical potential the two fixed points move in the vertical direction, but not in the direction.

We continue with the full determinant included. Consider the case with first, where real dynamics can be considered. Again we find the stable fixed point at and the unstable fixed point at . Provided that


there are no additional fixed points. In order to satisfy this condition, we take throughout. Using complex dynamics, while keeping , we find that the stable fixed point at remains, but that there are now three unstable fixed points at , given by and , with . Interestingly, the fixed-point structure is therefore different for real and complex flow.

Figure 5: Classical flow diagram in the plane for , , (top) and (bottom). The big dots indicate the fixed points at and . The small circles indicate a trajectory during the Langevin evolution, each dot separated from the previous one by 500 steps. Note the periodicity .
Figure 6: As in the previous figure, with (top) and (bottom).

Finally, we come to the full determinant at finite chemical potential. In this case the fixed points can only be determined numerically. We find a stable fixed point at and unstable fixed points at . The coordinates of these fixed points are determined by


where the upper (lower) sign applies to (). At there is only one solution, while at we find numerically that there are three (unstable) solutions for small chemical potential, two for intermediate and only one for large . Although the number of fixed points at depends on , we find that they are always unstable such that this has no effect on the dynamics, which is attracted to the stable fixed point at .

In Figs. 5 and 6 we show the classical flow diagrams in the plane. The direction of the arrows indicates , evaluated at . The lengths of the arrows are normalized for clarity. The fixed points are indicated with the larger black dots. In the bosonic model (), the analytical solution showed that the fixed point at is globally attractive, except when . At nonzero and , the fixed point at appears to be globally attractive as well, except again for . The small (blue) dots are part of a Langevin trajectory, each dot separated from the previous one by 500 steps. For vanishing , the dynamics takes place in the direction only; for increasing it spreads more and more in the direction. An interesting asymmetry around the classical fixed point in the direction can be observed. However, the dynamics remains well contained in the plane.

We conclude therefore that the complex Langevin dynamics does not change qualitatively in the presence of a chemical potential, small or large. We take this as a strong indication that the method is insensitive to the sign problem.

3.4 Fokker-Planck equation

The microscopic dynamics of the Langevin equation,


where is the (continuous) Langevin time, can be translated into the dynamics of a distribution , via the relation


From the Langevin equation, it follows that satisfies a Fokker-Planck equation,


where is the complex Fokker-Planck operator


The stationary solution of the Fokker-Planck equation is easily found by putting and is given by . In order to cast Eq. (3.34) as an eigenvalue problem, we write . The solution of the Fokker-Planck equation can then be written as


It is therefore interesting to study the properties of the Fokker-Planck equation and the nonzero eigenvalues .

In order to do so, we consider the model in the hopping expansion (3.18), with the action


Explicitly, the Fokker-Planck equation then reads


where primes/dots indicate / derivatives. Using periodicity, , we decompose


and we find




We note that this equation is completely real, such that all ’s are real. This is expected for the stationary solution, since from it follows that and therefore . The numerical solution of Eq. (3.40) is shown in Fig. 7 for a number of modes for (left) and 3 (right). The initial values for all . The number of modes is truncated, with . For large , exponentially fast. The zero mode is independent and equal to 1 by construction. We have verified that the other modes converge to the values determined by the stationary solution .

Figure 7: Solution of the complex Fokker-Planck equation: Langevin time dependence of the modes for various values of , with , , and (left) and (right).

The convergence properties can be understood from the eigenvalues of the Fokker-Planck operator. Writing gives the eigenvalue equation


Since all are real, also all eigenvalues are real. Although at first sight this may seem surprising, it is a consequence of the symmetry of the action and therefore it also holds in e.g. the full model.

The case corresponds to the stationary solution. Here we consider . First take . It follows from Eq. (3.42) that . As a result, the sequences for and split in two. Written in matrix form, they read


Approximating this matrix by a large but finite matrix, one can easily compute the eigenvalues numerically. We find that they are all positive and that the , sequences yield identical eigenvalues. In Fig. 8 (left) the four smallest nonzero eigenvalues are shown as a function of chemical potential. All eigenvalues are strictly positive and increase with . In Fig. 8 (right) the dependence on is indicated. At vanishing , the dependence cancels, since in that case . Also as a function of we observe that the eigenvalues are strictly positive.

Figure 8: Left: Four smallest nonzero eigenvalues of the complex Fokker-Planck equation as a function of with , . Right: Smallest nonzero eigenvalue as a function of for various values of at .

If the action and therefore the Langevin dynamics would be real, these results would be sufficient to explain the convergence of the observables towards the correct values as observed above, by employing Eq. (3.36) in the large limit. In the complex case we consider here, this is not immediately clear. Given the real Langevin equations,


one can also consider the real distribution , satisfying the Fokker-Planck equation


with the real Fokker-Planck operator


After complexification, expectation values should then satisfy


In contrast to the complex distribution , the distribution is real and has the interpretation as a probability distribution in the plane. The real and complex Fokker-Planck operators can be related, using


Here partial integration for finite is used as well as the analytic dependence of on . Eq. (3.48) suggests a relation between the spectrum of the complex and the real Fokker-Planck operator. However, we have not yet been able to show that a stationary solution of the real Fokker-Planck equation exists. We hope to come back to this issue in the future.

4 One link SU(3) model

4.1 Model

In this section we consider a one link model where the link is an element of SU(3). The partition function reads


with the bosonic part of the action666Note that for an SU(3) matrix, . Nevertheless, we write to allow for a straightforward complexification of the Langevin dynamics.


For the fermion matrix we take


with . We use the Pauli matrix rather than matrices to avoid factors of 2. The determinant has the product form


where the remaining determinants on the right-hand side are in colour space. In order to exponentiate the determinant, we use the identity, valid for SL(),


We find therefore that


with the quark and anti-quark contributions


Here we introduced the Polyakov loop and its “conjugate” ,


Note that . At large , the anti-quark contribution and no longer contributes. However, the term is crucial to preserve the symmetry (2.6) and is in particular relevant at imaginary and small real .

Observables are defined as


The observables we consider are the Polyakov loop , the conjugate Polyakov loop and the density . The latter is determined by


and reads


At zero chemical potential, the density vanishes while at large the density , the maximal numbers of (spinless) quarks on the site.

In this model, expectation values can be obtained directly by numerical integration, allowing for a comparison with the results from stochastic quantization presented below. Since we only consider observables that depend on the conjugacy class of , we only have to integrate over the conjugacy classes. These are parametrized by two angles . The reduced Haar measure on the conjugacy classes reads


where is a normalization constant. The matrix is parametrized as


such that




It is now straightforward to compute expectation values by numerical integration over and .

4.2 Complex Langevin dynamics

In contrast to in the U(1) model, the Langevin dynamics is now defined in terms of matrix multiplication. We denote and , where is again the Langevin time and consider the Langevin process,


Here () are the traceless, hermitian Gell-Mann matrices, normalized as . The noise satisfies


and the drift term reads


Differentiation is defined as


In particular, and .

The explicit expressions for the drift terms are




written in terms of


Let us first consider the case without chemical potential and take SU(3). Then it is easy to see that and therefore . Moreover, since the Gell-man matrices are traceless, . Therefore, if is an element of SU(3), it will remain in SU(3) during the Langevin process. The same results are found at finite imaginary chemical potential . At nonzero real on the other hand, we find that , although still holds. Therefore will be an element of SL() during the Langevin evolution. If is parametrized as , this implies that the gauge potentials are complex.

Figure 9: Real part of the Polyakov loop (left) and the conjugate Polyakov loop (right) as a function of for three values of at fixed .
Figure 10: Scatter plot of the Polyakov loop for , and .
Figure 9: Real part of the Polyakov loop (left) and the conjugate Polyakov loop (right) as a function of for three values of at fixed .

We solved the Langevin process (4.17) numerically. The matrix is computed by exponentiating the complex traceless matrix , employing Cardano’s method [37] for finding the eigenvalues. In Fig. 10 the real part of the Polyakov loop and the conjugate Polyakov loop are shown as a function of for three values of at fixed . The lines are the ‘exact’ results obtained by numerically integrating over the angles and , as discussed above. The symbols are obtained with complex Langevin dynamics, using the same Langevin stepsize and number of time steps as in the U(1) model ( and time steps). Errors are estimated with the jackknife procedure. Again, the imaginary part is zero analytically and consistent with zero within the error in the Langevin dynamics. Excellent agreement between the exact and the stochastic results can be seen.

Figure 11: Left: Real part of the density . Right: Deviation from SU(3) during complex Langevin evolution: as a function of Langevin step, for at .

Scatter plots of the Polyakov loop during the Langevin evolution are shown in Fig. 10 for four values of at and . Every point is separated from the previous one by 500 time steps. Note that the distribution is approximately symmetric under reflection , ensuring that within the error. We observe that the characteristic shape visible at becomes more and more fuzzy at larger , but the average remains well defined for all values of . The density is shown in Fig. 11 (left), with again good agreement between the exact and the stochastic results. We observe that saturation effects set in for smaller values of compared to the U(1) model, e.g.  already at .

During the complex Langevin evolution the matrix SL. In order to quantify how much it deviates from SU(3), we may follow the evolution of


It is easy to show that , with the equality in the case that SU(3).777Consider SL and . Using a polar decomposition, , with SU and a positive semidefinite hermitian matrix with , it is easy to show that , with the equal sign holding when SU. It provides therefore a good measure to quantify the deviation from SU(3). In Fig. 11 (right), we show this quantity during the Langevin evolution. We observe that the deviations from 1 are present but not too large. If is parametrized as , this implies the imaginary parts of the gauge potentials do not become unbounded.

4.3 Phase of the determinant

Figure 12: Left: Real part of . Right: Real part of .
Figure 13: Scatter plot of during the Langevin evolution, for , and .

As in the U(1) model, we study the phase of the determinant in the form


At zero chemical potential, the ratio is 1. Due to the SU(3) structure, however, the behaviour at large is qualitatively different. We find


such that


As a result the average phase goes to 1 at large and not towards 0 as in the U(1) model.888This difference can be traced back to Eq. (4.5). Therefore we expect the sign problem to become exponentially small in the saturation regime at large .

We have computed the average phase factor and the results are shown in Fig. 12 (left). The lines are again the ‘exact’ results. As is clear from this plot, the sign problem is quite mild for all values of , since the maximal deviation from 1 is less than 15%. In Fig. 12 (right) the ratio is shown. Here we observe a small but systematic deviation from 1, more pronounced at smaller and intermediate . However, we found that the deviation from 1 is reduced when continuing the Langevin evolution to larger and larger times. As in the U(1) model, this observable is the most sensitive and slowest converging quantity.

Scatter plots of the phase are presented in Fig. 13. At small chemical potential (top figure), there appears to be a similar structure as in the U(1) model, although not as pronounced. In the intermediate region (middle), the distribution is wider. At large (bottom), the distribution becomes narrow again, centered around .

5 QCD at finite chemical potential

5.1 Hopping expansion

In this section we leave the one link models behind and consider QCD at chemical potential. The SU(3) gauge field contribution to the euclidean lattice action is999We write rather than ; see footnote 6.


with . The plaquettes are defined as




The fermion matrix for Wilson fermions was already given in Eq. (2.7). The matrices satisfy and . We use periodic boundary conditions in space and antiperiodic boundary conditions in the euclidean time direction; the temperature and the number of time slices are related as (the lattice spacing ). The fermion matrix obeys


such that the determinant obeys Eq. (2.6).

We consider the heavy quark limit, where all spatial hopping terms are ignored and only the temporal links in the fermion determinant are preserved. We write therefore


where we defined and the (conjugate) Polyakov loops are


In the final line of Eq. (5.5) the determinant refers to colour space only. The sign appears because of the antiperiodic boundary conditions.

This approximation is motivated by the Heavy Dense Model considered e.g. in Refs. [30, 31, 32, 33, 34, 35, 36], in which the limit


was taken. However, here also the backward propagating quark, with the inverse Polyakov loop, is kept in order to preserve the relation (2.6).

Using Eq. (4.5), the determinant can now be written as