Stochastic quantization at finite chemical potential
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
- 2 Stochastic quantization and complex Langevin dynamics
- 3 One link U(1) model
- 4 One link SU(3) model
- 5 QCD at finite chemical potential
- 6 Summary and outlook
- A Fokker-Planck equation in Minkowski time
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 . 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. ). 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 () . 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:
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.
Note that .
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
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 .
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.
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
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.
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.
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.
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 .
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.
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
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
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
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.
We solved the Langevin process (4.17) numerically. The matrix is computed by exponentiating the complex traceless matrix , employing Cardano’s method  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.
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
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
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.
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