#### Controlling the sign problem in finite-density quantum field theory

Eur. Phys. J. C
Controlling the sign problem in finite-density quantum field theory
Nicolas Garron 0
Kurt Langfeld 0
0 Theoretical Physics Division, Department of Mathematical Sciences, University of Liverpool , Liverpool L69 3BX , UK
Quantum field theories at finite matter densities generically possess a partition function that is exponentially suppressed with the volume compared to that of the phase quenched analog. The smallness arises from an almost uniform distribution for the phase of the fermion determinant. Large cancellations upon integration is the origin of a poor signal to noise ratio. We study three alternatives for this integration: the Gaussian approximation, the “telegraphic” approximation, and a novel expansion in terms of theory-dependent moments and universal coefficients. We have tested the methods for QCD at finite densities of heavy quarks. We find that for two of the approximations the results are extremely close-if not identical-to the full answer in the strong sign-problem regime.
1 Introduction
The sign problem is known to be one of the most important
challenges of modern physics. In theoretical particle physics,
it prevents us from simulating finite-density QCD with
standard Monte-Carlo methods. Hence most of the QCD phase
diagram cannot be explored by first-principle techniques,
such as lattice QCD. Many reviews can be found; see for
example [
1–9
].
Dropping the phase factor of the quark determinant
exp{i φ} from the functional integral results in a theory, say
with partition function Z P Q , that is accessible by standard
importance sampling Monte-Carlo simulations. Very early
on, it became clear that Z P Q and the partition function of the
full theory Z are only comparable for the smallest values of
the chemical potential μ [
10
]. The deviation is quantified by
the so-called phase factor expectation value,
eiφ P Q = Z (μ)/Z P Q (μ) ∝ e− f V ,
(1)
where f is the free energy difference between the full and
the phase quenched theory and V is the volume (see e.g. [
10
]).
The knowledge of this phase factor would give access to the
partition function Z (μ) (we assume that Z P Q (μ) has been
obtained by standard methods). In this work, we study its
expectation value, eiφ P Q : it is a very small number,
generically very hard to measure due to the statistical noise, which
only decreases proportionally to the square root of the
number of Monte-Carlo configurations. Our approach is based
on the density-of-states method and in particular on the LLR
formulation [
11,12
], which is ideally suited to calculate
probability distributions of observables: it features an exponential
error suppression [12] which can result in an unprecedented
precision for the observable (see e.g. for an early
example [
13
]). It is based upon a non-Markovian random walk,
which immediately provides two main advantages: it bears
the potential to overcome the critical slowing down for
theories close to a first order phase transition [
9,14
], and it is
not restricted to theories with a positive probabilistic weight
for Monte-Carlo configurations. In fact, the method has been
successfully applied to the Z3 theory at finite densities [15]
and QCD at finite densities of heavy quarks [
16
]. In both
cases, the probability density ρ(φ) of the phase φ has been
obtained to very high precision. The phase factor expectation
value is then given by
eiφ P Q =
dφ ρ(φ) exp{i φ}
dφ ρ(φ)
(2)
Despite the high quality numerical result for ρ(φ), the
challenge remains to extract a very small signal from the above
Fourier transform. An approach, put forward in [
15,16
], is to
first represent the numerical data for ln ρ(φ) by a fit function
and then to calculate the Fourier transform of the fit function
(semi-)analytically. The method produces reliable results if
all the numerical data are well represented by the fit
function with a small number of fit parameters [
15,16
]. With the
advent of high precision data for ρ(φ), the main obstacle for
gaining access to quantum field theories at finite densities
is the above Fourier transform. The method used in [
15,16
]
hinges on the fact that a fit function which faithfully
represents the data could be found. This might not be generically
the case.
In this paper, we propose three alternatives to this direct
method. In Sect. 3 we present the first approach, called
Gaussian approximation. No fitting procedure is required;
instead the phase factor is computed directly from the data.
Within this framework, the integral in the numerator of (2)
is known analytically. The second approximation, presented
in Sect. 4 is what we call the “telegraphic” approximation.
This approach can be implemented either on the fit function
or directly on the data (although it might require new
simulations). The integral is replaced by a simple difference. In
Sect. 5, we introduce a third method, the “advanced moment
expansion”, which can be seen as a variant of a cumulant
expansion [
17–20
]. It is a systematic expansion in the
deviation from the uniform distribution and as such is expected to
work better in the strong sign-problem regime. We will
provide evidence that the universal coefficients decrease
exponentially with increasing order, providing a rapid
convergence if the moments are bounded. Although the convergence
is faster in the strong sign-problem regime, for the phase
factor expectation value we find an excellent agreement already
at the third order of the expansion, regardless of the strength
of the sign problem. In this case we still rely on a fitting
procedure for the density of states. However, the direct
computation of the Fourier transform (2) is not needed, only the
elementary moments are required. Before going through the
details of these methods, we present the framework and the
numerical details of our simulations in the next section. Our
conclusions are presented in Sect. 6.
2 Generalities and framework
2.1 Full theory and phase quenching
We consider a generic theory with a partition function
Z =
DUμ exp{β SYM[U ]} DetM [U ],
and with a complex “matter” determinant:
DetM [U ] = |DetM [U ]| exp{i φ[U ]}, φ ∈] − π, π ]. (4)
With the help of the density of states
ρ(s) =
DUμ exp{β SYM[U ]}|DetM [U ]|δ(s − φ[U ]),
(6)
(7)
(8)
(9)
(10)
(11)
(12)
the partition function can then be recovered by a
1dimensional Fourier transform:
Z =
ds ρ(s) exp{i s} =
ds ρ(s) cos(s).
We also introduce the so-called phase quenched counter part
by
The expectation values of an observable A in the full and in
the phase quenched theory are given as usual by
DUμ A exp{β SYM[U ]}DetM [U ],
DUμ A exp{β SYM[U ]}|DetM [U ]|,
(3)
(5)
=
ds ρ(s).
1
A = Z
1
A P Q = Z P Q
A =
Aeiφ P Q
eiφ P Q
implying the well-known relations
Z = Z P Q eiφ P Q .
In terms of the density, the phase factor expectation value is
given by (2).
2.2 Extensive density of states
For theories for which the imaginary part arises from a local
action an extensive phase x ∈] − ∞, ∞[ can be defined
as the sum of the local phases. This has been e.g. the case
for the finite density Z3 and for heavy dense QCD [
15,16
].
The expectation is that the extensive phase scales with the
volume. For fermionic theories with the phase φ[U ] arising
from the (non-local) quark determinant, the definition of an
extensive phase is not obvious. If the quark operator still
admits (complex) eigenvalues, a straightforward definition
would involve a sum of the phases of the eigenvalues. Another
possibility was as pointed out in [21]:
x [U ] = Im ln (Det M )
μ/T ∂(ln Det M )
= 0 Im ∂μ/ T
=
0
μ/T
Im tr M −1 ∂ M
∂μ/ T
μ=μ¯
μ=μ¯
d
d
μ
¯
T
μ
¯ .
T
The definition of an extensive phase factor has proven to be
important to achieve the precision needed for the Fourier
transform. If ρE (x ) denotes the corresponding probability
distribution, the phase factor expectation value in (2) is
obtained by
−a1 x 2 − a2 xV4 − a3 Vx 62 + · · · .
If we define a “scaling” variable by x = s/√V , the
deviation from a Gaussian distribution decreases with increasing
volume:
1 π
eiφ P Q = Z P Q n∈Z −π
2.3 Volume dependence of the density
We here consider the class of theories for which the phase of
the Gibbs factor is proportional to the chemical potential μ
and for which this is the only μ dependence. Scalar theories
do not fall into this class since the real part of the action
also acquires a μ dependence, but fermion theories in the
ab initio continuum formulation might fall into this class.
For these theories, let us study the dependence of ρE (s) on
the physical volume V . We make explicit the μ dependence
of the phase factor expectation value and point out that the
partition function is positive for all μ:
z(μ) =
ei μ φ P Q ≥ 0.
Note that we have z(0) = 1 and that from z(μ) ∈ R it follows
that
e−i μ φ P Q =
ei μ φ P Q
⇒ z(−μ) = z(μ).
Note that since z(μ) is obtained by a Fourier transform of ρ,
see (2), the density of states can be recovered from z(μ) by
the inverse Fourier transform (up to a normalisation constant
Z P Q ≥ 0)
ρE (s) = Z P Q
dμ z(μ) e−isμ.
2π
As argued in [
21
], z(μ) can be viewed as a partition function
with free energy density f (μ) (a necessary condition is that
z(μ) ≥ 0), leaving us with the volume dependence:
z(μ) = exp{− f (μ) V },
= exp{−[c1μ2 + c2μ4 + c3μ6 + · · · ] V },
1
eiφ P Q = Z P Q
+∞
dx ρE (x ) cos(x ).
−∞
The density of states ρ(s) can easily be recovered from the
extended density ρE (x ). To see this, we subtract from x a
multiple of 2π until s ∈ [−π, π [, s = x − 2π n, n ∈ Z, and
we split the integration domain in intervals of size 2π ,
where the coefficients ck are volume independent. Inserting
(20) into (18), we find, with an expansion in inverse powers
of V ,
ρE (s) = const. exp
s2 s4 s6
−a1 V − a2 V 3 − a3 V 5 + · · · ,
(13)
(14)
(15)
(16)
(17)
(18)
2.4 Numerical details
We use the data obtained in our previous work [
16
] but
have also generated new simulations for reasons that we
explain below. We summarise here the parameters used for
the numerical simulations and the methods to obtain the
density of states. The interested reader will find more details in
the aforementioned reference. The lattice parameters are
84 lattice,
β = 5.8,
κ = 0.12.
and we let the chemical potential μ vary between 1.0421 and
1.4321. We identified the “strong sign-problem region” as
being 1.1 < μ < 1.4. We for each value of μ, we split the
domain of the phase s ∈ [0, smax] in nint small interval of size
δs and on each interval k, we compute the LLR coefficients
ak . In practise we choose smax ∼ 36, δs = 0.896 and nint =
40, except for a few values of the chemical potential, for
which we need a better resolution. The corresponding values
are reported in Table 1.
We reconstruct the probability density function for
discrete values of the phase sk = kδs + δs /2, namely
ρE (sk ) = exp −
ai δs − ak δs /2 .
k−1
i=1
(22)
(23)
nnint
In [
16
] we performed a polynomial fit of ln(ρE ) and
computed (13) by a semi-analytic integration (we refer to this
method as “Exact”).
Although the fits are of very good quality and very stable,
for three values of the chemical potential, we have also ran
new simulations with δs = π/5. As shown below, these new
data allow us to compute ρ(s) directly from the data (without
relying on any fitting procedure) and will be very useful to
check the methods presented here. We have implemented this
technique for three different values of the chemical potential.
This is illustrated in Figs. 1, 2 and 3, where we see that the
different methods give compatible results.
Finally, we mention that we use around 1000
configurations and that the statistical errors are estimated with the
bootstrap method, using 500 samples. Naturally we have checked
that the errors are stable with respect to the number of
samples.
!
In fact, numerical results suggest that the probability
distribution is Gaussian to a good extent [
17,22–24
], which
would imply that only the cumulant φ2 c is non-vanishing.
It has been argued in [21] that higher cumulants are
suppressed by factors of the volume V and that, however, higher
order cumulants are important for the medium and high range
of chemical potentials. Throughout this paper, we define
the Gaussian approximation as the approximation of the
extended density of states by a normal distribution:
φ4 c − · · · .
(24)
ρE (s) ≈ const. exp{− s2}.
The phase factor expectation value (2) is then analytically
obtained:
(25)
(26)
from the standard expectation
eiφ P Q = exp
1
− 4
where the subscript E indicates that the expectation values
are defined with respect to the extended density ρE . We test
this approach for heavy-dense QCD with partition function
(7). We find the expectation value in (27) directly from the
data: we take the density obtained through (23) and compute
the expectation value s2 E using a trapezoidal
approximation. We obtain in this way an estimate for the phase
factor expectation value (26) without invoking any fitting
procedure. Our numerical findings are summarised in Fig. 4.
We find that the Gaussian approximation provides a
surprisingly good approximation over the whole range of chemical
potentials μ. Even in the strong sign-problem regime at
intermediate values μ, the cancellations are well emulated and
the approximate result only underestimates the true result by
roughly a factor 2.
4 The “telegraphic” approximation
4.1 Methodology
As can be seen in Fig. 3, ρ weakly depends on its arguments
in the strong sign-problem regime and for large volumes.
In this case, a Poisson re-summation of (15) should yield a
rapidly converging series:
dn e2πi ν n ρE (s + 2π n)
∞
dx ei ν x ρE (x ).
(29)
(30)
(31)
(32)
If the sum over ν is rapidly converging, we find
approximately
ρ(s)/c0 ≈ 1 + 2 eiφ P Q cos(s).
In the strong sign-problem regime, the amplitude of the
cosine is very small, and therefore we see that ρ(s) is almost
a constant. Equation (33) then offers the possibility to extract
the phase factor expectation value, i.e.,
1
eiφ P Q ≈ 4c0 [ρ(0) − ρ(π )].
Using (15), we therefore find
eiφ P Q ≈ π2 k∈∞Z(−d1x)kρρEE(x(k) π ) . (35)
−∞
We call this the telegraphic approximation. It emerges by
neglecting higher contributions cν of the Poisson sum. In
order to get a feeling for the resulting systematic error, we
adopt, for now only, the Gaussian approximation (25) and
find:
c2
c1 ≈
exp
1
− 4
3
This implies that the correction to ρ(s) in (33) is of order:
where we have used (26) and (32). At least in the strong
signproblem regime, for which eiφ P Q is very small, we expect
the telegraphic approximation to work very well.
We finally point out that the telegraphic approximation can
be improved in a systematic way. The order of the
approximation is defined by the number of harmonics entering the
density of states. E.g., in third order we have
(28)
ρ(s)/c0 ≈ 1 + 2 eiφ P Q cos(s) + c cos(2s)
+ d cos(3s),
(33)
(34)
(36)
(37)
with the unknowns eiφ P Q and c, d. We generate three
equations by evaluating ρ(s) at s = 0, π/3, π and solve the
linear set of equations for the unknowns. We are predominantly
interested in the phase factor:
1 1 1
eiφ P Q ≈ − 2 + 4 ρ(0)/c0 + 3 ρ(π/3)/c0
1
− 12 ρ(π )/c0,
which can easily be converted to a discrete sum over discrete
set of points of ρE (s) using (15).
4.2 Numerical implementation
Again, we use Heavy-Dense QCD to test this approximation.
Having in hands the density of state—either ρE obtained
from the fit or ρ from the date through (15)—it is
straightforward to implement numerically (35). If we take the results
from the fit, we find that this approximation provide results
extremely close to the “exact” ones: except for a few values of
μ in the weak sign-problem regime, the results (central value
and variance) are actually indistinguishable. For example, for
μ = 1.0821, we find
ln eiφ ePxQact = −1.992175 ± 2.910279 × 10−3,
ln eiφ aPpQprox = −1.992174 ± 2.910306 × 10−3.
(38)
(39)
We show our results for the various μ in Table 2 and Fig. 5.
It is essential that the numerical data are accurate enough
to produce significant results upon the cancellations in (35).
Clearly, the LLR approach in combination with the fit of the
density of states succeeds. In the following, we will test the
“naive” implementation for which we compute ρ(s) directly
from the data (without relying on any fitting procedure). For
this purpose, we use our new simulation data with δs = π/5.
In this case we have ρ(s) for s = π/10, 3π/10, . . ., but do
not have ρ(0) nor ρ(π ). Therefore we use a variant of (34):
Although in the strong sign-problem regime μ = 1.2921,
we could not extract a signal directly from the data
(without a representation of the density of states by a fit), for the
two other values of μ we find reasonable agreement. This
shows that the telegraphic approximation cannot master the
cancellations without further technical advances. We stress,
however, that the telegraphic approximation greatly assists
the LLR approach by reducing the Fouier integration of the
density of states to an alternating sum (see (35)).
5 The advanced moments approach
5.1 General formulation
The starting point is the expansion of the density of states:
The coefficients d j depend on the underlying theory, and
N0 ≥ 2 will define the order of the expansion. Our
conjecture is that the coefficients d j are suppressed by powers of
the volume with increasing j . For QCD, this conjecture is
supported by the strong coupling expansion and the hadron
resonance gas model [
21
]. There is also some numerical
evidence by the WHOT-QCD collaboration [
22–24
]. Last but
not least, this conjecture becomes true for the limited class
of theories considered in Sect. 2.2. Using (47) in (14), we can
express the phase factor expectation in terms of the
theorydependent coefficients d j :
1
eiφ P Q = Z P Q j=1 d j I2 j ,
where d0 has dropped out upon integration, and where
I2 j =
ds s2 j cos(s) =
(−1) j−l+1
2(2 j )! π 2l−1.
(2l − 1))!
The values I2k can be efficiently calculated by the recursion
I2k = −2(2k)π 2k−1 − (2k)(2k − 1)I2k−2,
with the initial condition I0 = 0. Our strategy to access
the coefficients d j in an actual numerical simulation is to
calculate combinations as the simple moments s2n . Using
the truncation (47) for a given N0, we find
s2n+2
with
2π 2i+2 j+1
Ai j = 2i + 2 j + 1 .
=
1
N0−1
Z P Q j=0
Anj d j
Keeping in mind that we have s2n+2 available from a
numerical simulation, the idea is to choose a set of n-values
and to consider (51) as a linear set of equations to obtain the
unknowns d j . Note that for n = −1, s2n+2 = 1 = 0
follows from the symmetry ρ(−s) = ρ(s) and does not
contain theory specific information. We hence choose n =
0, . . . , N0 − 1 and obtain
( A−1) jn s2n+2 .
Inserting this into (48), we obtain
d j
Z P Q
=
We now have at our fingertips the moment expansion of the
phase factor for a given order N0. We have not yet achieved
γN +1 n =
k(N )
n−1 −
α2k γkn / k(NN )
M2N +2
N
1
= k(N )
N n=0
= s2N +2
N
1
+ n=1 k(NN ) kn(N−)1 −
kn(N ) s2n+2 −
α2k
γkn s2n
k
n=1
α2k γkn
s2n ,
where we have changed the order of the double sum. We
therefore find the recursion
a systematic expansion, featuring increments of decreasing
size (when we increase the order N0). To this aim, we define
the first advanced moment M4 for N0 = 2 by
We then define recursively for N = 2, . . . , (N0 − 1):
N0−1
n=1
α2N +2 = k(NN ),
and finally we achieve the systematic expansion
We stress that the coefficients α2n+2 are universal, i.e., the
only dependence on the theory under investigations enters
via the moments Mk . Last but not least, we would like to
have an explicit representation of the advanced moments M
in terms of the simple expectations values sn . We define
M2k =
γki s2i .
k
i=1
By construction of the advance moments, we have the
normalisation γkk = 1. Although for high order N 1 the
intermediate coefficients γNi can become very large (we will show
this below), the field theories of interest, i.e., finite-density
quantum field theory in the strong sign-problem regime,
should give advanced moments within bounds. In this case,
the convergence is then left to the coefficients αn. Inserting
(62) into (59), we find after renaming of indices
(48)
(49)
(50)
(51)
(52)
(53)
(54)
(55)
N
N
k=2
k=max(n,2)
N
k=max(n,2)
γN +1 N +1 = 1 ,
where 1 ≤ n ≤ N and 2 ≤ N ≤ N0 − 1. The recursion
can be solved in closed form for i ∈ {2, . . . , N0} and j ∈
{1, . . . , N0}:
γii = 1, γ21 =
γi j =
k(1)
k0(1) , γi j = 0 for j > i,
1
k(ji−−11) − k(ji−−12) , i > j and i > 2.
k(i−1)
i−1
5.2 The first advanced moments
For illustration purposes, we will explicitly calculate the first
few advanced moments. The main task is to obtain the
coefficients ki(N ), which emerge from the solution of a linear set
of equations; see (55).
For the leading order N0 = 2, we find
⎛ π 3
( Ai j ) = 2 ⎜ 3
⎜⎝ π 5
5
π 5 ⎞
5
π 7 ⎠⎟⎟ ,
7
(I2 j ) =
0
−4π
.
k(1)
1
175
= α4 = − 2π 6 ,
k0(1)/ k1(1)
Hence, the first advanced moment, see (56), is given by
M4 =
s4
− 53 π 2 s2 .
At next to leading order, i.e., N0 = 3, we have
⎛ π 3
⎜ 3
( Ai j ) = 2 ⎜⎜⎜ π 5
⎜⎝⎜⎜ π57
7
π 5
5
π 7
7
π 9
9
π 7 ⎞
7 ⎟
11
⎛
π 9 ⎟⎟⎟ , (I2 j ) = ⎝
9 ⎟⎟
π 11 ⎠⎟
The solution of the corresponding linear system is given by
k(2)
0
k(2)
1
k(2)
2
= − 9845 2π 2π−6 33 ,
315 16π 2 − 231
= 4 π 8 ,
4851 2π 2 − 27
= − 8 π 10
= α6.
From (66), we then find for the coefficients γ
γ31 =
k(2)
0 − α4γ21
k(2)
2
= 251 π 4,
(66)
(67)
(68)
(69)
(70)
(71)
(72)
(73)
(74)
(75)
We have computed the moment coefficients up to order N0 =
5. We find for the coefficient matrix (k ≥ 2, i ≥ 1):
(76)
(77)
1
0 ⎟⎟⎟⎟⎟⎟⎟
0 ⎟⎟⎟⎟
1 ⎠
We finally perform a consistency check. For a truncation of
the density of states at order N0, all the moments up to M2N0
contribute to the phase factor expectation value at this order,
see (61). If we consider (47) as exact for the moment in the
sense that all simple moments s2n are calculated with this
density, then the phase factor expectation value is obtained
exactly by summing all contributions including the term
containing M2N0 . Since this result is already exact, all moments
M2k with k > N0 must vanish. For example, assume that the
density is given by
ρ(s) = d0 + d1s2,
(N0 = 2),
then e.g. M6 (and all higher moments need to vanish for all
choices for d0 and d1. This devises a consistency check. We
find for the present example
s6
10
n
20
30
s4
s2
= 15
Inserting these simple moments into M6, (76), we find that all
terms cancel and that M6 indeed vanishes for all choices of d0
and d1. If we consider, in a quantum field theory setting, the
expansion (47) as an expansion with respect to some inverse
power of the volume, the moments M2n are then suppressed
by these powers.
5.3 Convergence
For high orders N0, the coefficients γ in the definition (62)
of the advanced moments M2k can become very large. In
this section, we will assume that for functions ρ(s)
arising in a quantum field theory setting the moments remain
within bounds. This occurs due to cancellations between
simple moments s2i , as we will show below. In this case, the
expansion (61) of the phase factor expectation value in terms
of the advanced moments is dictated by behaviour of the
coefficients α2k for large k. These coefficients are universal:
they do not depend on the underlying theory, i.e., ρ(s). They
arise from the solution of the linear system (55), which reads
in a shorthand notation
k = A−1 I,
and it is this linear system that we are going to study in greater
detail. Since the matrix A in (52) is symmetric and positive,
we perform a Cholesky decomposition and solve for k:
A = L L T ,
L y = I,
L T k = y,
where L is a lower triangular matrix. Note that if the system
L y = b is solved at order N0 and if subsequently the order
N0 is increased, the first N0 components of the solution y are
unaffected by the increase due to the triangular form of L.
The same is true for the matrix L: increasing the order from
N0 to N0 + 1 does not affect the first N0 rows and columns.
We are interested in the N0 dependence of the last component
of k:
k(NN00) = yN0 /L N0,N0 .
The Cholesky decomposition gives
Lii =
Aii −
1
Li j = Li j ( Ai j −
i−1
k=1
j−1
Li2k ,
k=1
Lik L jk ),
i > j
(78)
(79)
(80)
(81)
We have solved this iteration analytically for values N0 up to
30. We find that for large n the data is well described. We find
In essence, the advanced moments approach from Sect. 5 is
an efficient numerical method to evaluate the Fourier
transFig. 6 The exponential rate K , see (82), as a function of the order n
of the expansion
that very quickly Lnn reaches an asymptotic regime which is
well described by
Ln+1,n+1 = K π 2 Lnn,
1
K = 4 ,
see Fig. 6. Asymptotically, we therefore find the exponential
increase:
L N0,N0 ∝
Unfortunately, we could not prove any of these asymptotic
behaviours analytically, but we have verified (85) by also
solving the linear system L T k = y for k. Our analytical
result for N0 = 2 to N0 = 32 is shown in Fig. 8. We find the
remarkable result that the expansion coefficients α2N0 are
exponentially decreasing with N0 suggesting a rapid
convergence of the advanced moment expansion as long as the
moments M2n are bounded.
5.4 Application to HDQCD (82) (83) (85)
1.004
of the LLR method, they are extracted with a very good
statistical precision.
We also observe that going from s2 to s8 , the relative
error increases very slowly. We turn now to the advanced
moments: since all the elementary moments are positive, the
relative signs in (88)–(94) imply that important cancellations
occur. At leading order (LO), we have
Fig. 8 The behaviour of the coefficient α2N0 = kNN00 of the expansion
in terms of advanced moments
form (14) for sufficiently smooth integrands ρ(s). In this
section, we test the method in the quantum field theory
context of QCD at finite densities of heavy quarks (HDQCD).
Our preliminary results have been reported in [
25
].
Here we are interested in the strong sign-problem region
(in which μ ∼ 1.3): in Fig. 3, we show that the density
is almost constant whereas for μ ∼ 1 and μ ∼ 1.4, the
density has variation of order 1 (see Figs. 1, 2). Hence, we
expect that the advanced moment expansion will have a better
convergence in the strong sign-problem regime.
From now on, we focus on the severe sign problem region,
μ = 1.2921. Once the density is known, we can compute the
elementary moments (again using our fit results and
semianalytic integration). They are reported in Table 3. By virtue
(86)
(87)
(88)
(89)
(90)
(91)
(92)
(93)
(94)
As expected, strong cancellations between the simple moments
occur making it mandatory to determine the simple moments
with high precision. The analysis has been carried out using
the bootstrap resampling method, and we point out that strong
correlations are at work to obtain the advanced moments at
the level of precision reported here. The numerical values
are also reported in Table 4. One should note that the overall
signs of the advanced moments oscillate; however, αi Mi is
a positive quantity, as can be seen in (61) or in the numerical
values.
The phase factor expectation value (61) is then given by
eiφ
When the order of the expansion increases, the statistical
error decreases and that the results converges quickly to the
“exact” answer
eiφ
= 2.189(324) × 10−6,
obtained by fitting the extensive density ρE and by carrying
out the Fourier transform using the fit, as in [
16
]. (In the latter
we quote 2.37(21) × 10−6, the small difference in the central
value comes from the fact that we use a different δs .) We
observe a rapid convergence here.
Since the phase factor is a small number, it is useful to
look at the logarithm of this quantity. We find
log eiφ P Q = −13.032 ± 0.152
(Full).
The advanced moment method yields
log eiφ P Q
= −13.445 ± 0.152 + O(α6 M6)
= −13.065 ± 0.152 + O(α8 M8)
= −13.033 ± 0.152 + O(α10 M10)
(LO),
(NLO),
(NNLO).
It is remarkable that not only the central value but also the
variance is very well approximated by our expansions. Indeed
for this value of μ, the full (relative) variance is already given
by the first order. Of course the quality of the approximation
depends on the variation of ρ (and therefore on the strength
of the sign problem).
We now vary the value of μ in the range 1 < μ < 1.4
and compare the results of the phase factor expectation value
obtained in [
16
] with the method proposed here. It is
interesting to note that even in the weak sign-problem region,
in which the density ρ fluctuates between 0 and 1, the NLO
and NNLO approximations already yield reasonable
approximations. This is illustrated in Figs. 9 and 10. Our numerical
results can be found in Table 5. We quote the “full answer” as
(95)
(96)
(97)
(98)
(99)
(100)
(101)
(102)
obtained in [
16
] and the relative difference with the method
presented here, for the first three orders. (Here we implement
the advanced moments method with the same δs as in [
16
].)
The NLO approximation works at the percent level over the
full available range, even in the weak sign-problem region.
6 Conclusions
There are two main possibilities in addressing finite-density
quantum field theory: (i) facing the large cancellations that
give rise to the smallness of the partition function or (ii) to
reformulate to an equivalent theory say by dualisation [
6
] or
by a complexification of the fields [
5
]. Method (ii) would
be preferred if the approach exists and if exactness can be
1.042
1.062
1.082
1.102
1.122
1.142
1.162
1.182
1.202
1.232
1.252
1.272
1.292
1.312
1.332
1.352
1.372
1.392
ln exp(i φ)
guaranteed. The appeal of method (i) is that it is universally
applicable if a way is found to control the cancellations.
A first success for direction (i) emerged with the advent of
Wang–Landau type techniques and, most notably, the LLR
method [
11
]: due to the feature of exponential error
suppression of the LLR approach [
12
], high precision data for the
density of states ρ (s) of finding a particular phase s over
many orders of magnitude has become available. The
partition function now emerges as Fourier transform of ρ (s). Due
to large cancellations, this Fourier transform is a challenge in
its own right. The recent success reported in [
15
] and in [
16
]
hinge on the ability to find a fit function for ln ρ (s) that well
represents hundreds of numerical data points with relatively
few fit parameters. This situation is unsatisfactory since the
quest for this fit function might not be always successful.
The present paper explores three methods to perform the
Fourier transform:
• The Gaussian approximation of the extensive density of
states ρE is most easily implemented, but hard to improve
in a systematic way. For the example of HDQCD, we
found this approximation yields the right order of
magnitude throughout and only misses the exact phase factor
by a factor of 2 when the sign problem is strongest.
• The telegraphic approximation yields the phase factor
through an alternating (discrete) sum of the extensive
density of states ρE . The relative systematic error is of the
order of the phase factor itself, which makes the
approximation excellent in the strong sign-problem regime.
• The advanced moment approach is a systematic
expansion of this Fourier transform with respect to the
deviations of ρ (s) from uniformity. The expansion therefore
works best in the strong sign-problem regime. The
expansion is independent of the quantum field theory setting
and can be applied to the Fourier transform of any
sufficiently smooth function ρ (s), s ∈ [−π, π ]. At the heart
of expansion are the so-called advanced moments. We
have thoroughly derived these moments and the
theoryindependent expansion coefficients α. We found evidence
that the expansion coefficients decrease exponentially
with increasing order, thus guaranteeing rapid
convergence if ρ (s) admits moments M that are bounded. We
have tested and validated the advanced moment
expansion in the context of HDQCD: we have confirmed that
the expansion converges very quickly. It works best in
the strong sign-problem region as expected, although at
third order the results agree with the “full” answer at the
sub-percent level even in the weak sign-problem regime.
Acknowledgements We are grateful to B. Lucini and A. Rago for
helpful discussions. NG and KL are supported by the Leverhulme Trust
(Grant RPG-2014-118) and, KL by STFC (Grant ST/L000350/1).
Open Access This article is distributed under the terms of the Creative
Commons Attribution 4.0 International License (http://creativecomm
ons.org/licenses/by/4.0/), which permits unrestricted use, distribution,
and reproduction in any medium, provided you give appropriate credit
to the original author(s) and the source, provide a link to the Creative
Commons license, and indicate if changes were made.
Funded by SCOAP3.
1. K. Fukushima , T. Hatsuda, The phase diagram of dense QCD. Rep. Prog. Phys . 74 , 014001 ( 2011 ). doi: 10 .1088/ 0034 -4885/74/ 1/014001. arXiv: 1005 . 4814
2. P. de Forcrand, Simulating QCD at finite density . In: PoS LAT2009 , p. 010 ( 2009 ). arXiv: 1005 . 0539
3. S. Gupta , QCD at finite density . In: PoS LATTICE2010 , p. 007 ( 2010 ). arXiv: 1101 . 0109
4. L. Levkova, QCD at nonzero temperature and density . In: PoS LATTICE2011 , p. 011 ( 2011 ). arXiv: 1201 . 1516
5. G. Aarts, Complex Langevin dynamics and other approaches at finite chemical potential . In: PoS LATTICE2012 , p. 017 ( 2012 ). arXiv: 1302 . 3028
6. C. Gattringer , New developments for dual methods in lattice field theory at non-zero density . In: PoS LATTICE2013 , p. 002 ( 2014 ). arXiv: 1401 . 7788
7. D. Sexty , New algorithms for finite density QCD . In: PoS LATTICE2014 , p. 016 ( 2014 ). arXiv: 1410 . 8813
8. S. Borsnyi , Fluctuations at finite temperature and density . In: PoS LATTICE2015 , p. 015 ( 2016 ). arXiv: 1511 . 06541
9. K. Langfeld , Density-of-states . In: Proceedings, 34th international symposium on lattice field theory (Lattice 2016 ), Southampton, 24 - 30 July 2016 ( 2016 ). arXiv: 1610 . 09856
10. M. Troyer , U.-J. Wiese, Computational complexity and fundamental limitations to fermionic quantum Monte Carlo simulations . Phys. Rev. Lett . 94 , 170201 ( 2005 ). doi:10.1103/PhysRevLett.94. 170201. arXiv:cond-mat/0408370
11. K. Langfeld , B. Lucini , A. Rago , The density of states in gauge theories . Phys. Rev. Lett . 109 , 111601 ( 2012 ). doi: 10 .1103/ PhysRevLett.109.111601. arXiv: 1204 . 3243
12. K. Langfeld , B. Lucini , R. Pellegrini , A. Rago , An efficient algorithm for numerical computations of continuous densities of states . Eur. Phys. J. C 76 ( 6 ), 306 ( 2016 ). doi: 10 .1140/epjc/ s10052-016-4142-5. arXiv: 1509 . 08391
13. K. Langfeld , J.M. Pawlowski , Two-color QCD with heavy quarks at finite densities . Phys. Rev. D 88 ( 7 ), 071502 ( 2013 ). doi: 10 .1103/ PhysRevD.88.071502. arXiv: 1307 . 0455
14. B. Lucini , W. Fall , K. Langfeld , Overcoming strong metastabilities with the LLR method . In: Proceedings of the 34th international symposium on lattice field theory (Lattice 2016 ), Southampton, 24 - 30 July 2016 ( 2016 ). arXiv: 1611 . 00019
15. K. Langfeld , B. Lucini , Density of states approach to dense quantum systems . Phys. Rev. D 90 ( 9 ), 094502 ( 2014 ). doi: 10 .1103/ PhysRevD.90.094502. arXiv: 1404 . 7187
16. N. Garron , K. Langfeld , Anatomy of the sign-problem in heavydense QCD . Eur. Phys. J. C 76 ( 10 ), 569 ( 2016 ). doi: 10 .1140/epjc/ s10052-016-4412-2. arXiv: 1605 . 02709
17. S. Ejiri, On the existence of the critical point in finite density lattice QCD . Phys. Rev. D 77 , 014508 ( 2008 ). doi: 10 .1140/epjc/ s10052-016-4412-2. arXiv: 0706 . 3549
18. Y. Nakagawa , S. Ejiri , S. Aoki , K. Kanaya , H. Ohno , H. Saito , T. Hatsuda , T. Umeda, Histogram method in finite density QCD with phase quenched simulations . In: PoS LATTICE2011 , p. 208 ( 2011 ). arXiv: 1111 . 2116
19. H. Saito , S. Aoki , K. Kanaya , H. Ohno , S. Ejiri , Y. Nakagawa , T. Hatsuda , T. Umeda, Finite density QCD phase transition in the heavy quark region . In: PoS LATTICE2011 , p. 214 ( 2011 ). arXiv: 1202 . 6113
20. H. Saito , S. Ejiri , S. Aoki , K. Kanaya , Y. Nakagawa , H. Ohno , K. Okuno , T. Umeda, Histograms in heavy-quark QCD at finite temperature and density . Phys. Rev. D 89 ( 3 ), 034507 ( 2014 ). doi: 10 . 1103/PhysRevD.89.034507. arXiv: 1309 . 2445
21. J. Greensite , J.C. Myers , K. Splittorff , The density in the density of states method . JHEP 10 , 192 ( 2013 ). doi: 10 .1007/ JHEP10( 2013 ) 192 . arXiv: 1308 . 6712
22. S. Ejiri , S. Aoki , T. Hatsuda , K. Kanaya , Y. Nakagawa , H. Ohno , H. Saito , T. Umeda, Numerical study of QCD phase diagram at high temperature and density by a histogram method . Central Eur. J. Phys . 10 , 1322 - 1325 ( 2012 ). doi:10.2478/s11534-012-0054-7 . arXiv: 1203 . 3793
23. S. Ejiri , Y. Nakagawa , S. Aoki , K. Kanaya , H. Saito , T. Hatsuda , H. Ohno , T. Umeda, Probability distribution functions in the finite density lattice QCD . In: PoS LATTICE2012 , p. 089 ( 2012 ). arXiv: 1212 . 0762
24. S. Ejiri, Phase structure of hot dense QCD by a histogram method . Eur. Phys. J. A 49 , 86 ( 2013 ). doi: 10 .1140/epja/i2013-13086-7. arXiv: 1306 . 0295
25. N. Garron , K. Langfeld , Tackling the sign problem with a moment expansion and application to Heavy dense QCD . In: Proceedings of the 34th international symposium on lattice field theory (Lattice 2016 ), Southampton, 24 - 30 July 2016 ( 2016 ). arXiv: 1611 . 01378