Anatomy of the sign-problem in heavy-dense QCD

Abstract

QCD at finite densities of heavy quarks is investigated using the density-of-states method. The phase factor expectation value of the quark determinant is calculated to unprecedented precision as a function of the chemical potential. Results are validated using those from a reweighting approach where the latter can produce a significant signal-to-noise ratio. We confirm the particle–hole symmetry at low temperatures, find a strong sign problem at intermediate values of the chemical potential, and an inverse Silver Blaze feature for chemical potentials close to the onset value: here, the phase-quenched theory underestimates the density of the full theory.

1 Introduction

Monte Carlo simulations of quantum chromo dynamics (QCD) at finite baryon densities would provide direct insights into cold, but dense matter as it occurs in compact stars. They would also trigger the evolution of effective theories. To date, there are numerous proposals for such theories and models. Those rise from exact solvable models that mimic certain aspects of QCD (see the Gross–Neveu model [1, 2]) or are motivated by certain limits of QCD: The limit of many colours has led to the proposal of the “quarkyonic phase” [3, 4]. Reducing the gluon sector to the essence of the centre elements has revealed that “centre-dressed quarks” obey Bose statistics and can undergo Bose–Einstein condensation in the dense, but still confined phase (see “Fermi–Einstein condensation” in [5]). Since heavy-ion collision experiments probe matter at high temperatures, but—at best—at moderate densities, the essential input for understanding cold-dense baryonic matter has to come from first principles computer simulations. Standard Monte Carlo simulations, however, fail since the Gibbs factor is complex at non-vanishing chemical potentials and, thus, lacks the interpretation of a probabilistic weight for lattice configurations. This problem does not exclusively relate to dense QCD, but is generic for dense matter quantum field theories. It has become known as the notorious “sign-problem” over the last three decades.

The recent years, however, have seen significant progress in the numerical studies of complex action systems, both with Monte Carlo methods and techniques that do not rely on importance sampling. Among the most promising methods are the complexification of the fields in a Langevin based approach [6, 7], worm or flux algorithms [8, 9] to simulate the dual theory if it happens that this theory is real [10–13] and the use of techniques that explicitly exploit the cancellations of classes of fields [14].

Among the alternatives to conventional Monte Carlo sampling, the so-called density-of-states simulations (for early results for the gauge and spin systems see [15, 16]): this approach performs Monte Carlo updates according to the number of states for a given (complex) action and employs the pioneering techniques introduced by Wang and Landau [17] to refine the density-of-states during simulation. Once this quantity has been determined, the partition function and derived expectation values of observables can be computed semi-analytically by integrating the density-of-states with the appropriate (potentially complex) Boltzmann weight. More recently, a Wang–Landau type method originally introduced for continuous systems has been put forward in [18–20]. This method features an exponential error suppression and allows one to calculate the density-of-states over many orders of magnitude [21]. At least for the $Z_3$ spin model at finite densities, the achieved precision of the density-of-states has been high enough to solve the strong sign problem by direct integration [22].

Heavy-dense QCD (HDQCD) emerges in the limit in which the quark mass and chemical potential are simultaneously large [23, 24]. This theory has a non-trivial phase diagram in the plane of temperature and chemical potential, which qualitatively agrees with the one expected for real QCD: e.g., at vanishing chemical potential, there is a thermal deconfinement transition as the temperature is increased with the transition being first order for very heavy quarks and a crossover for slightly lighter but still heavy quarks [25]. The gluonic part of HDQCD is given by the SU(3) Yang–Mills theory, and a dualisation that could leave us with a real theory at presence of a chemical potentials is not known. So far, this rules out any flux or worm-type algorithms and makes it a significant testing ground for the density-of-states techniques. We point out that HDQCD has been simulated with complex Langevin method providing results for bench-marking our findings [25, 26]. We also refer the reader to [27] for a recent study of HDQCD using reweighting and a mean-field approximation.

In this paper, we study HDQCD with the density-of-states approach detailed in [22]. The theory is real in the limit of vanishing and of large chemical potentials and for chemical potential equalling the heavy quark mass. Although the phase-quenched approximation sketches a qualitatively correct picture for this reason, we do find a strong sign problem for chemical potentials close to the mass threshold.

2 Heavy-dense QCD and the generalised density-of-states approach

2.1 HDQCD: definitions and features

The partition function of QCD with the quarks field integrated out is a functional integral over SU(3) unitary matrices only:

$$\begin{aligned} Z(\mu ) = \int \mathcal{D }U_\mu \exp \{ \beta S_\mathrm {YM}[U] \} \hbox {Det} M(\mu ) , \end{aligned}$$

(1)

where we use the Wilson formulation of the Yang–Mills action:

$$\begin{aligned} S_\mathrm {YM}[U]= & {} \frac{1}{3} \sum _{x,\mu >\nu } \hbox {Re} \hbox {tr}\Bigl [ U_\mu (x) U_\nu (x+\mu ) \nonumber \\&\times U^{\dagger } _\mu (x+\nu ) U^{\dagger } _\nu (x) \Bigr ] . \end{aligned}$$

(2)

The so-called quark determinant possesses the property

$$\begin{aligned} (\hbox {Det} M(\mu ))^*= \hbox {Det} M(- \mu ), \quad (\mu \in {\mathbb {R}}), \end{aligned}$$

(3)

which implies that QCD at vanishing chemical potential, i.e., $\mu =0$, is a real theory. For large quark mass m and simultaneously large chemical potential $\mu$, the quark determinant factorises into [23–27]

$$\begin{aligned} \hbox {Det} M(\mu )= & {} \prod _{\mathbf {x}} {\det } ^2\left( 1 + h \mathrm {e}^{\mu / T} P(\mathbf {x})\right) \nonumber \\&\times {\det } ^2 \left( 1 + h \mathrm {e}^{ - \mu / T} P^\dagger (\mathbf {x})\right) , \end{aligned}$$

(4)

where $T = 1/ N_t a$ is the temperature with a the lattice spacing and $N_t$ the number of lattice points in temporal direction. The parameter h is related to the quark hopping parameter $\kappa$ and $P(\mathbf {x})$ is the Polyakov line starting at position $\mathbf {x}$ and winding around the torus in temporal direction:

$$\begin{aligned} h = (2 \kappa )^{N_t},\quad P(\mathbf {x}) = \prod _{t=1}^{N_t} U_4 (\mathbf {x},t). \end{aligned}$$

(5)

The determinants at the right hand side of (3) extend over colour indices only. Introducing the heavy quark mass m by

$$\begin{aligned} m a =- \ln (2 \kappa ), \end{aligned}$$

(6)

we find that $h = \exp \{ - m/T \}$, yielding for (4)

$$\begin{aligned} \hbox {Det} M(\mu )= & {} \prod _{\mathbf {x}} {\det }^2 ( 1 + \mathrm {e}^{(\mu -m) / T} P(\mathbf {x})) \nonumber \\&\times {\det }^2 ( 1 + \mathrm {e}^{ - (\mu + m) / T} P^\dagger (\mathbf {x})). \end{aligned}$$

(7)

Inspection of the latter equation easily confirms that

$$\begin{aligned} \hbox {Det} M(\mu =0) \in {\mathbb {R}}. \end{aligned}$$

(8)

For non-vanishing $\mu$, we will indeed find that the determinant is complex (albeit the imaginary part can be very small; see below). However, we are going to show that the partition function is nevertheless real, i.e., the imaginary part of Z vanishes upon the integration over gauge configurations. This can be most easily seen by adopting the Polyakov gauge where

$$\begin{aligned} U_4 (t \not =1,\mathbf {x}) = 1, \quad P(\mathbf {x}) = U_4 (t=1,\mathbf {x}). \end{aligned}$$

The partition function takes the form

$$\begin{aligned} Z = \int \mathcal{D}U_\mu \mathrm {e}^{\beta S_\mathrm {YM}} f (U_4(1,\mathbf {x}), U^\dagger _4(1,\mathbf {x})), \end{aligned}$$

where f is a real and analytic function. Given that the Haar measure and the action are real, we find upon the substitution $U_4(1,\mathbf {x}) \rightarrow U^\dagger _4(1,\mathbf {x})$ that

$$\begin{aligned} Z (\mu ) = Z^*(\mu ). \end{aligned}$$

(9)

For positive chemical potentials and for low temperatures, i.e.,

$$\begin{aligned} \mu \ge 0, \quad \frac{m}{T} \gg 1, \end{aligned}$$

(10)

we can neglect quark excitations from the Dirac sea. Formally, the second determinant in (7) equals unity to a very good approximation, and we find

$$\begin{aligned} \hbox {Det} M (\mu )\approx & {} \prod _{\mathbf {x}} {\det }^2 (1 + \mathrm {e}^{(\mu -m) / T} P(\mathbf {x})). \end{aligned}$$

(11)

For any unitary matrix $P \in SU(3)$, we find that

$$\begin{aligned} \det (1 + c P) = 1 + c \hbox {tr}P + c^2 \hbox {tr}P^\dagger + c^3 . \end{aligned}$$

(12)

This implies that the quark determinant is also real for $\mu = m$ (i.e., $c=1$) (see also [27]):

$$\begin{aligned} \hbox {Det} M (\mu =m) \in {\mathbb {R}}. \end{aligned}$$

(13)

Let us now study the case of large chemical potentials, i.e., $\mu \gg m$. Starting from (11), we obtain

$$\begin{aligned} \hbox {Det} M (\mu )= & {} \mathrm {e}^{2N_c V (\mu -m) / T} \nonumber \\&\times \prod _{\mathbf {x}} {\det }^2 ( 1 + \mathrm {e}^{- (\mu -m) / T} P^\dagger (\mathbf {x})), \end{aligned}$$

(14)

where $N_c=3$ is the number of colours, $V=\sum _{\mathbf {x}}$ is the spatial volume and where we have used the fact that P is a unitary matrix, i.e., $P P^\dagger =1$, $\det P =1$. It is convenient to introduce the scaled chemical potential relative to the mass threshold:

$$\begin{aligned} t = \frac{ \mu - m }{T} . \end{aligned}$$

(15)

Using (11) in the functional integral (1), the partition function only depends on t and obeys the relation

$$\begin{aligned} Z (t) \approx \mathrm {e}^{2N_c V t } Z (-t) \quad (m\gg T) , \end{aligned}$$

(16)

where we have used the fact that Z is real (see (9)). As usual, we define the quark density by

$$\begin{aligned} \sigma (t) = \frac{T}{V} \frac{ \partial \ln Z(\mu ) }{\partial \mu } = \frac{1}{V} \frac{ \partial \ln Z(t) }{\partial t } . \end{aligned}$$

(17)

Using (16), we find the duality

$$\begin{aligned} \sigma (t) \approx 2N_c - \sigma (-t) \quad (m\gg T) . \end{aligned}$$

(18)

For negative t, the chemical potential is below the mass threshold and the density $\sigma (t)$ rapidly approaches zero with decreasing t. This implies with the help of (18) that, for large t, the density rapidly approaches the saturation density:

$$\begin{aligned} \sigma (t) \mathop {\rightarrow }\limits ^{t \rightarrow \infty } 2N_c . \end{aligned}$$

(19)

As a side-remark, we point out that in this regime, i.e., $\mu \gg m$, the quark determinant becomes a (real) constant (see (14)),

$$\begin{aligned} \hbox {Det} M (\mu ) \approx \mathrm {e}^{2N_c V (\mu -m) / T} , \end{aligned}$$

and the partition function at large $\mu$ is given by that of pure SU(3) Yang–Mills theory up to a multiplicative constant.

2.2 Reweighting simulations

If the imaginary part of the quark determinant is small, i.e., for $\mu \approx 0$ or $\mu \approx m$ or $\mu \gg m$, the standard reweighting procedure can produce reliable results. Using a polar decomposition of the determinant, the partition function (1) can be rewritten as

$$\begin{aligned} Z(\mu ) = \int \mathcal{D }U_\mu \mathrm {e}^{ \beta S_\mathrm {YM}[U] } \vert \hbox {Det} M(\mu ) \vert \exp \{ i \phi [U] \} . \end{aligned}$$

(20)

We here introduce the partition function of the phase-quenched theory by

$$\begin{aligned} Z_\mathrm {PQ}(\mu ) = \int \mathcal{D }U_\mu \mathrm {e}^{ \beta S_\mathrm {YM}[U] } \vert \hbox {Det} M(\mu )\vert . \end{aligned}$$

(21)

Sometimes, the phase-quenched theory is referred to as QCD with an iso-spin chemical potential. Indeed, rewriting e.g.

$$\begin{aligned} \vert \hbox {Det} M(\mu ) \vert ^2= & {} \hbox {Det} M(\mu ) \hbox {Det}^*M(\mu ) \nonumber \\= & {} \hbox {Det} M(\mu ) \hbox {Det} M(-\mu ) , \nonumber \end{aligned}$$

the phase-quenched theories can be interpreted as (in this case) a 2-flavour quark theory with a chemical potential coupling to the flavours with opposite sign.

The Monte Carlo simulation based upon reweighting generates a Markov chain of configurations $\{U_\mu \}$ of the phase-quenched theory (21). The expectation value of any observable A is then obtained by

$$\begin{aligned} \langle A \rangle = \frac{ \langle A \exp \{ i \phi [U] \} \rangle _\mathrm {PQ}}{\langle \exp \{ i \phi [U] \} \rangle _\mathrm {PQ} } . \end{aligned}$$

(22)

For a successful reweighting approach, it is essential that the phase factor expectation value, i.e.,

$$\begin{aligned} \langle \exp \{ i \phi [U] \} \rangle _{PQ} = \frac{ Z(\mu ) }{ Z_\mathrm {PQ}(\mu ) } , \end{aligned}$$

(23)

is of significant size. This would ensure a good signal-to-noise ratio. However, it has been known for a long time (see e.g. [28]) that the full and phase-quenched theories have a difference in their free energy density, say $\Delta f$. Using the triangle inequality, one also finds that

$$\begin{aligned} Z_\mathrm {PQ}(\mu ) \le Z(\mu ) . \end{aligned}$$

Hence, the ratio of their partition function in (23) is exponentially suppressed with the volume V:

$$\begin{aligned} \langle \exp \{ i \phi [U] \} \rangle _{PQ} = \exp \{ - \Delta f V \} , \quad \Delta f \ge 0 . \end{aligned}$$

Consequently, reweighting simulations are restricted to the parameter space for which the quark determinant is almost real, i.e.,

$$\begin{aligned} \Delta f (\mu ) \approx \mathcal{O}( 1/V ) . \end{aligned}$$

2.3 Density-of-states method

The density-of-states method belongs to the class of Wang–Landau type simulations [17]. It has been argued in [21] that the LLR version [18] possesses an exponential error suppression that allows one to estimate a strongly suppressed phase factor expectation value (23) with good relative precision. This has been demonstrated for the first time for the $Z_3$ spin model at finite densities [22].

Central to all Wang–Landau methods is the density-of-states, which is defined in the present case of HDQCD by

$$\begin{aligned} \rho (s) = \int \mathcal{D}U_\mu \delta ( s - \phi [U]) \mathrm {e}^{\beta S_\mathrm {YM}[U] } \vert \hbox {Det} M\vert . \end{aligned}$$

(24)

Using this definition, the phase factor expectation value (23) can be obtained by Fourier transform

$$\begin{aligned} \langle \mathrm {e}^{i \phi } \rangle = \frac{ \int \mathrm {d}s \rho (s) \exp \{ i s\} }{ \int \mathrm {d}s \rho (s)}. \end{aligned}$$

(25)

Since the final answer is potentially a very small number, the density-of-states method needs to overcome two issues here: (i) $\rho (s)$ must be calculated to high precision over the whole range of s. This is where standard histogram methods fail: they do not produce enough statistics in certain regions of s (overlap problem). (ii) The smallness of $\langle \exp \{i \phi \} \rangle$ arises from cancellations implying that the numerical integration must be carried out with extreme care. The LLR algorithm generically overcomes the issue (i), and we refer to the literature for details (most notably see [29] for a thorough discussion of the theoretical framework). To resolve issue (ii), we will adopt the approach that proved successful in the case of the $Z_3$ spin model [22], and we will present details in the result section.

Fig. 1

Left The phase factor expectation value $\langle \mathrm {e}^{i\phi } \rangle$ as a function of the chemical potential $\mu$ (simulation parameters in (29)); black symbols the reweighting approach; red symbols the LLR approach as a preview. Right detail of the graph

We finally point out that the quark density $\sigma (\mu )$ can be calculated once good results for the phase factor expectation value are available. This arises from the observation that (23) leads to

$$\begin{aligned} \sigma (\mu )= & {} \frac{T}{V} \frac{\mathrm{d}}{\mathrm{d}\mu } \ln \langle \mathrm {e}^{i \phi } \rangle (\mu ) + \sigma _\mathrm {PQ}(\mu ), \end{aligned}$$

(26)

where we have introduced the phase-quenched quark density by

$$\begin{aligned} \sigma _\mathrm {PQ}(\mu )= & {} \frac{T}{V} \frac{ \mathrm{d} \ln Z_\mathrm {PQ} (\mu ) }{ \mathrm{d} \mu }. \end{aligned}$$

(27)

3 Results from reweighting

Throughout this paper, we use discretised space-time employing a $N^4$ cubic lattice and the Wilson action (2). We work in the Polyakov gauge, i.e., all links are updated except

$$\begin{aligned} U_4(\mathbf {x},t\not =1) = 1 . \end{aligned}$$

This implies that the remaining time-like links are identified with the Polyakov line:

$$\begin{aligned} U_4(\mathbf {x},t=1) = P(\mathbf {x}) . \end{aligned}$$

Using the gauge invariance of the quark determinant, it is apparent that $\hbox {Det} M$ does only depend on $\hbox {tr}P^n(\mathbf {x})$ [30–32]. We use the local hybrid-Monte Carlo (LHMC) simulation algorithm (with respect to the angles of the algebra) for the update of configurations according to the phase-quenched partition function (21). We have validated and fine-tuned the algorithm by comparing some of the results with those obtained by the standard Cabibbo–Marinari method. The LHMC update shows shorter auto-correlation times (e.g. for the topological charge). The simulation parameters are

$$\begin{aligned} N=8,\,\beta = 5.8,\,\kappa = 0.12,\,N_\mathrm {conf} = 12{,}000 \end{aligned}$$

(28)

where $N_\mathrm {conf}$ is the number of the independent configurations for the Monte Carlo estimators. Errors are obtained by a bootstrap analysis. Our findings from the reweighting approach are shown in Fig. 1. The chemical potentials are chosen symmetrically around the mass threshold, which is (using $\kappa = 0.12$, into (6))

$$\begin{aligned} a m \approx 1.427 . \end{aligned}$$

Our numerical findings are in line with the theoretical predictions in Sect. 2.1: the phase factor expectation value approaches 1 for small and large values of $\mu$ and for $\mu$ close to the mass threshold. Because of the particle–hole duality (18), we can confine ourselves to discussing only the case $\mu \le m$. It is remarkable that on a quantitative level the reweighting approach produces reliable results for $\mu$ as large as 1. Note, however, that, for the intermediate values, i.e.,

$$\begin{aligned} 1.15 \mathop {_\sim }\limits ^{<} \mu \mathop {_\sim }\limits ^{<} 1.4 , \end{aligned}$$

we do encounter a sign problem with the signal being much smaller than the noise.

Let us discuss the implications for the quark density $\sigma (\mu )$. We start with a discussion of the phase-quenched density. Since the only $\mu$ dependence is in the quark operator, we find

$$\begin{aligned} \sigma _\mathrm {PQ}(\mu )= & {} \frac{1}{Z_\mathrm {PQ}(\mu )} \int \mathcal{D }U_\mu \mathrm {e}^{ \beta S_\mathrm {YM}[U] } \vert \hbox {Det} M(\mu )\vert \nonumber \\&\times \frac{\partial }{\partial \mu } \ln \vert \hbox {Det} M(\mu ) \vert , \end{aligned}$$

(29)

where for HDQCD $\hbox {Det} M$ is given in (4). As detailed in Sect. 2.1, HDQCD is real for vanishing chemical potential, for $\mu =m$ and for large $\mu$ implying that $\sigma = \sigma _\mathrm {PQ}$ for these limiting cases. This signals that the phase-quenched density shows the correct behaviour for small $\mu$, the correct onset at $\mu =m$ and the correct asymptotic value given by saturation. It is therefore expected that $\sigma _\mathrm {PQ}(\mu )$ qualitatively reflects the $\mu$ dependence of the full density $\sigma$. This is indeed verified by our direct evaluation of (29) shown in Fig. 2. Although phase quenching produces qualitatively correct results, we cannot conclude that the sign problem is weak (see below).

Fig. 2

Quark density $\sigma _\mathrm {PQ}(\mu )$ of the phase-quenched theory as a function of the chemical potential $\mu$

Regardless of the quantitative details, we can draw some interesting conclusions for the density using the identity (26). For small chemical potentials, e.g., $\mu \le 1.1$, the phase factor expectation value is decreasing. Consequently, the correction

$$\begin{aligned} \frac{T}{V} \frac{\mathrm{d}}{\mathrm{d}\mu } \langle \mathrm {e}^{i \phi } \rangle \end{aligned}$$

is negative implying that the phase-quenched result overestimates the true result $\sigma$. This is usually referred to as “Silver Blaze problem”. With a smoothness assumption of $\langle \mathrm {e}^{i \phi } \rangle$, we expect that its derivative with respect to $\mu$ vanishes at $\mu _{c1}$ with $1.15< \mu _{c1} < 1.4$. For this chemical potential, we find agreement:

$$\begin{aligned} \sigma _\mathrm {PQ} (\mu _{c1}) = \sigma (\mu _{c1}) . \end{aligned}$$

For $\mu _{c1}< \mu < m$, the derivative of $\langle \mathrm {e}^{i \phi } \rangle$ is positive. We here find an inverted Silver Blaze behaviour: close to the mass threshold, the phase-quenched theory underestimates the value of the density. We stress, however, that a study involving several volumes and temperatures would be needed to decide whether this effect has a role to play for phenomenology. This is left to future work.

4 LLR results

4.1 Foundations of the LLR simulation

Our aim is to calculate an approximation of the density-of-states $\rho (s)$ for the imaginary part s of the quark determinant. We divide the domain of support for $\rho$ into intervals $[s_k, s_k+ \mathrm{d}s]$. Under physically motivated assumptions, $\rho (s)$ is a smooth function such that a Taylor expansion over these intervals yields a valid approximation. Central to the LLR approach are the Taylor coefficients (also called LLR coefficients)

$$\begin{aligned} a_k := \frac{ \mathrm{d} \ln \rho }{\mathrm{d}s } \vert _{s=s_k+ \mathrm{d}s/2}, \end{aligned}$$

(30)

which will be the target of our numerical simulations below. With these coefficients at our fingertips, we use a piece-wise linear approximation for $\ln \rho$ and derive the approximation:

$$\begin{aligned} \rho (s)= & {} \rho _0 \left( \prod _{k=1}^{N-1} \mathrm {e}^{a_k \mathrm{d}s} \right) \exp \left\{ a_N \left( s - s_N \right) \right\} , \end{aligned}$$

(31)

where, for a given s, the upper boundary N is chosen such that

$$\begin{aligned} s _N \le s \le s_N + \mathrm{d}s , \quad s_k = s_0 + k \mathrm{d}s . \end{aligned}$$

The goal of the LLR method is to calculate the coefficients from a stochastic non-linear equation. A key ingredient of this equation is the restricted and reweighted expectation values [18] with a being an external variable (not to be confused with the lattice spacing):

(32)

$$\begin{aligned} N_k= & {} \int \mathcal{D} U_\mu \vert \hbox {Det}M \vert \mathrm {e}^{\beta S_{YM} } \theta _{[s_k,\delta _s]}(\phi [U]) \nonumber \\&\times \mathrm {e}^{-a \phi [U] } , \end{aligned}$$

(33)

where we have introduced the modified Heaviside function,

$$\begin{aligned} \theta _{[s_k,\delta _s]} (\phi ) = \left\{ \begin{array}{ll} 1 &{}\quad \hbox {for} \, s_k \le \phi \le s_k + \delta _s\\ 0 &{}\quad \hbox {otherwise . } \end{array} \right. \end{aligned}$$

For the particular choice

$$\begin{aligned} W[\phi ] = \phi - s_k - \frac{\delta _s}{2} =: \Delta \phi \end{aligned}$$

we showed that

$$\begin{aligned} \left\langle \left\langle \Delta \phi \right\rangle \right\rangle _{k} (a) = 0 \quad \hbox {for} \quad a = a_k . \end{aligned}$$

(34)

The latter equation is a non-linear equation to obtain a. For instance, this can be done by using the fixed point iteration:

$$\begin{aligned} a_k^{(n+1)} = a_k^{(n)} + \frac{12}{\delta _s^2} \left\langle \left\langle \Delta \phi \right\rangle \right\rangle _{k} \left( a^{(n)}_k\right) . \end{aligned}$$

Note that the expectation value $\left\langle \left\langle \Delta \phi \right\rangle \right\rangle _{k}$ is not known exactly. An estimate, however, can be obtained by standard Monte Carlo simulations. The issue here is that the statistical error interferes with convergence of the fixed point iteration. The mathematical framework to obtain a solution was developed by Robbins and Monro. They showed that the under relaxed iteration

$$\begin{aligned} a_k^{(n+1)}= & {} a_k^{(n)} + \alpha _n \frac{12}{\delta _s^2} \left\langle \left\langle \Delta \phi \right\rangle \right\rangle _{k} \left( a^{(n)}_k \right) \end{aligned}$$

(35)

$$\begin{aligned} \sum _n \alpha _n\rightarrow & {} \infty , \quad \sum _n \alpha ^2_n = \hbox {finite} , \end{aligned}$$

(36)

converges to the correct answer. Moreover, if the iteration is truncated at $N=N_\mathrm {cut}$ and independently repeated many times, the final values $a_k^{(N_\mathrm {cut})}$ are normal distributed with the true value $a_k$ as mean. This paves the way to a bootstrap analysis to obtain an error estimate for our estimate for $a_k$. A common choice is ($0 < \gamma \le 1$)

$$\begin{aligned} \alpha _n = \left\{ \begin{array}{ll@{\quad }l} 1 &{} \quad \hbox {for} \quad &{} 0 \le n \le n_t ,\\ 1/(n-n_t)^\gamma &{} \quad \hbox {for} \quad &{} n > n_t , \end{array} \right. \end{aligned}$$

(37)

where the iterations with $n \le n_t$ are considered as thermalisation steps, and for which the limiting case $\gamma =1$ is the optimal choice for error suppression.

Table 1 Simulation parameters for one particular value s

Once the Taylor coefficients are obtained for the range s of interest, the generalised density-of-states $\rho (s)$ can be calculated in the usual way:

$$\begin{aligned} \ln \rho (s)= & {} -\sum _{k=1}^{n-1} a_i \delta _s- a_n \delta _s/2 \end{aligned}$$

(38)

$$\begin{aligned}&n \,\hbox {such that:}\, s_n \le s < s_{n+1} . \end{aligned}$$

(39)

Our final target is phase factor expectation value, which can be obtained by means of two LLR integrals (details of the numerical method will be presented in Sect. 4.4 below):

$$\begin{aligned} \langle \mathrm {e}^{i\phi } \rangle =\frac{\displaystyle \int _0^{s_\mathrm{max}} \rho (s) \cos (s) \mathrm{d}s }{\displaystyle \int _0^{s_\mathrm{max}} \rho (s) \mathrm{d}s } \end{aligned}$$

(40)

Since $\rho (s)$ is rapidly decreasing, we will find that it is not difficult to find a reliable cutoff $s_\mathrm{max}$.

4.2 Thermalisation

We find that the thermalisation is most demanding for small interval sizes $\delta _s$ and for chemical potentials near the onset value. In order to provide insight into the thermalisation history, we present here some results for the simulation parameters listed in Table 1.

The thermalisation history for 40 independent random starts is shown in Fig. 3. Between each iterations, we performed 40 sweeps at a fixed parameter $a_k^{(n)}$ in order to let the system equilibrate.

We see a decrease of the width of the error band with increasing iteration number n, which is due to the Robbins Monro underrelaxation. In the production runs for the results below, we have chosen $n_t=200$ and a maximum of 1, 000 iterations. We then make use of the Robbins Monro feature that the final values for $a_k$ are normal distributed with the correct mean. For the statistical analysis, we repeated each iteration 40 times and use the copies for $a_k$ for the bootstrap analysis.

Fig. 3

Thermalisation history (simulations parameters are in Table 1)

Fig. 4

Strong sign-problem regime: the LLR coefficient a(s) as a function of s for $\mu =1.3321$ (left panel). Right detail of the graph

For a consistency check and to analyse the effect of the Robbins Monro parameter $\gamma$, we calculated the average $a_k$ for different values of $\gamma$. We find

$$\begin{aligned} \begin{array}{l|lllll} \gamma &{} 0.6 &{} 0.7 &{} 0.8 &{} 0.9 &{} 1.0 \\ - a_k &{} 3.287 &{} 3.334 &{} 3.256 &{} 3.288 &{} 3.300 \\ \mathrm{err} [10^{-2}] &{} 5.397 &{} 3.495 &{} 2.356 &{} 1.656 &{} 1.082 \end{array} \end{aligned}$$

We did not observe any ergodicity issues and found that the limiting case $\gamma =1$ is most effective for error reduction as expected.

4.3 Probability distribution of the imaginary part

According to Fig. 1, we will distinguish three parameter regimes, depending on the choice of the chemical potential $\mu$:

• Low density regime for $\mu \mathop {_\sim }\limits ^{<} 1.1$: this regime might be accessible by a Taylor expansion with respect to $\mu$ and simulations using reweighting.

• Regime with a strong sign problem for $1.1 \mathop {_\sim }\limits ^{<} \mu \mathop {_\sim }\limits ^{<} 1.4$: this regime is beyond the scope of standard Monte Carlo methods and will be specifically targeted with the LLR method below.

• Dense regime for $1.4 \mathop {_\sim }\limits ^{<} \mu \le m\approx 1.427$: the system possesses a significant quark density, which reaches half of the saturation density for $\mu =m$.

Because of the duality (16), we do not need to explicitly explore the regime $\mu > m$. We stress that the above regime boundaries have been chosen in an ad hoc way. We are not aware of any physical phenomenon that would define these boundaries in a rigorous way. The different regimes above, however, have quite distinct features as we will reveal in this section by exploring the density-of-states.

To this aim, we have calculated the LLR coefficients $a_k$ (30) over a range of imaginary parts s for given chemical potentials. The simulation parameters again have been

$$\begin{aligned} 8^4 \quad \beta =5.8 \quad \kappa =0.12 . \end{aligned}$$

Note that the LLR method becomes exact in the limit of vanishing interval size $\delta _s$.

In practice, we check that our result for $a_k$ does not dependent on $\delta _s$. We illustrate this fact for $\mu = 1.3321$, which belongs to the interesting regime of a strong sign problem. Our findings are shown in Fig. 4. We find that the coefficients are quite insensitive to size of $\delta _s$. This also holds for the other regimes. Note that a smaller $\delta _s$ requires more intervals to cover the same (integration) domain for s. We found that $\delta _s= 0.896$ is a good compromise between accuracy and computational effort, and it is this value which we have used in most simulations.

Fig. 5

Low density regime: the LLR coefficient a(s) as a function of s for several values of the chemical potential $\mu$ between 1.0421 and 1.1421. The error bars are smaller than the symbols

Fig. 6

Strong sign-problem regime (i): a(s) for several values of the chemical potential $\mu$ between 1.2021 and 1.2921. Note that the y-scale differs from the previous plot. We observe that a(s) is an increasing function of $\mu$ for any value of $s>0$

Fig. 7

Strong sign-problem regime (ii): a(s) for several values of the chemical potential $\mu$ between 1.3121 and 1.3721. In this range, we observe that a(s) is a decreasing function of $\mu$

Figures 5, 6, 7 and 8, left panel, show the LLR coefficient as a function of s for various values of the chemical potential. We stress that in these figures, the error bars are present but smaller than the symbols. Error bars are obtained from 40 independent sets of a that are subjected to 500 bootstrap samples. Figure 5 shows the low density regime. We find a slight modulation of a(s) with s, which did not occur for $\mu =1.3321$ (see Fig. 4). In Figs. 6 and 7, we summarise our findings for a(s) for a range of chemical potentials that mostly belong to the strong sign-problem regime. We observe a quite distinct behaviour: the curvature of the curves increases with increasing chemical potential. For the largest values of $\mu$ shown in Fig. 7, we enter the dense phase. Our largest values of $\mu$ are shown in Fig. 8, left panel. Here, we observe that a(s) starts to show an oscillatory behaviour. Needless to say that we have checked that these oscillations are independent of the choice of $\delta _s$ and statistically significant. This is illustrated in Fig. 8, left panel, where we show the coefficient a(s) for the chemical potential $\mu =1.4321$, which is slightly above the mass threshold of $m=1.42711$.

Fig. 8

Left Dense regime (i): a(s) for several values of the chemical potential $\mu$ close to the mass threshold m. Right Dense regime (ii): a(s) as a function of s near “half-filling” (slightly above m) for three different values of $\delta _s$

Table 2 Fit results for $\mu =1.3321$ and $\delta s = 0.29867$. We show the fit coefficients for different truncations $A_i$, the corresponding $\chi ^2$ per degree-of-freedom and the result of the integration. Missing results imply that the corresponding coefficient is fixed to zero. In the last rows, the results are obtained by numerical integration either with or without folding

4.4 The LLR integration

Once the coefficient a(s) has been extracted, we are in a position to calculate the phase factor expectation value $\langle \mathrm{e} ^{i \phi }\rangle$ for a given value of $\mu$ by means of (40). The straightforward method would be to make use of the piece-wise linear interpolation (39) and to control the systematic errors in the Riemann sense by making $\delta _s$ smaller. It was already noted in [22] for the case of the $Z_3$ theory at finite densities that this method does not muster enough precision at an affordable size $\delta _s$ to obtain a good signal-to-noise ratio. Instead of seeking convergence in the Riemann sense, we expand $\ln \rho$ in terms of basis functions $f_n(s)$:

$$\begin{aligned} \ln \rho (s) = \sum _{n}^{N_\mathrm {max}} c_n f_n(s) . \end{aligned}$$

(41)

The approximation now occurs by the truncation of the above sum at $N_\mathrm {max}$. Here, we follow the strategy of compressed sensing (see e.g. [33]) and choose the basis in such a way that a minimal number of coefficients $c_n$ represents the data at given accuracy and $\chi ^2$ per degree-of-freedom (dof) of the fit. It is quite remarkable that a basis with simple powers of s, i.e.,

$$\begin{aligned} f_n(s) = s ^{n}. \end{aligned}$$

(42)

already produces very good results, at least for the $Z_3$ theory [22]. Equation (42) is also our choice here for HDQCD. Note that coefficients $c_n$, with n are incompatible with the theory’s reflection symmetry $\ln \rho (-s) = \ln \rho (s)$ and are therefore set to zero.

In summary, our approach is:

• Using the numerical estimates $a_k$, we build the function $P(s)=\ln (\rho (s))$ according to

  $$\begin{aligned} P(s)= & {} -\sum _{k=1}^{n-1} a_i \delta _s- a_n \delta _s/2 , \end{aligned}$$

  (43)
  $$\begin{aligned} s= & {} s_n + \delta _s/2 = n \delta _s+ \delta _s/2 , \end{aligned}$$

  (44)

  where in the last equation, we choose $s_0=0$ as a starting point.

• We fit the result to an even-power polynomial

  $$\begin{aligned} P(s) = \sum _{i=0}^{\mathrm{deg}/2} {c_{2i}} s^{2i} . \end{aligned}$$

  (45)

• From the fit result, we reconstruct the density

  $$\begin{aligned} \rho (s) = \exp (P(s)) . \end{aligned}$$

  (46)

• Finally, we semi-analytically compute the LLR integral

  $$\begin{aligned} \langle \mathrm {e}^{i\phi } \rangle =\frac{\displaystyle \int _0^{s_\mathrm{max}} \rho (s) \cos (s) \mathrm{d}s }{\displaystyle \int _0^{s_\mathrm{max}} \rho (s) \mathrm{d}s }. \end{aligned}$$

  (47)

We have performed various checks in order to ensure that our procedure is stable. First, we have tried different truncations: we denote by $A_i$ a fit to a polynomial of degree i in which all the coefficient $c_{2i}$ are free parameters. We also performed some fits with $c_0$ fixed to 0, we call them $\tilde{A}_i$. Some details of our fit procedure for the finest $\delta s$ can be found in Table 2 for the specific value of $\mu =1.3321$. By comparing $\tilde{A}_2$ with $A_2$ and $A_6$ with $\tilde{A}_6$, we see that constant term $c_0$ has very little effect on the other fit parameters. All in all, we observe that the fit procedure is robust, however, our data are clearly best fitted by a degree-6 polynomial. Adding higher degrees gives compatible results with larger errors (see $A_8$). We also present the fit results for $\delta _s=0.29867$ in Fig. 9.

Fig. 9

Left Natural logarithm of the density-of-states as a function of s for $\mu =1.3321$, in the strong sign-problem regime. We show both the data points with their error bars (blue cross) and the fit results (red solid line) for 120 intervals in s, corresponding to $\delta _s\sim 0.299$. Right same as left, but on a linear scale

Fig. 10

Left Result for the phase factor expectation value as a function of the cut $s_\mathrm{max}$. Results are shown for $\mu =1.3321$ and $\delta s = 0.29867$. Right $\langle \mathrm {e}^{i\phi } \rangle$ for $\mu =1.3321$ as a function of $\delta _s^2$. The results of the extrapolation is $1.185 (18) \times 10^{-5}$, statistical error only

Since we are looking for a very small signal emerging after large cancellations, even the trivial identity

$$\begin{aligned} \int _0^{s_\mathrm{max}} \rightarrow \frac{1}{2} \int _{-s_\mathrm{max}}^{s_\mathrm{max}} \quad \hbox {(folding)} \end{aligned}$$

(48)

might perform differently upon its numerical implementation. In order to check the robustness of our results, we implemented both integrals. In Table 2, the first integral (from 0 to $s_\mathrm{max}$) is denoted by (i) and the second (from $-s_\mathrm{max}$ to $s_\mathrm{max}$) is marked by (ii). We see that the difference is smaller than the statistical error.

We have also checked that the results do not depend on the cutoff $s_\mathrm{max}$, which is expected since $\rho (s)$ is rapidly decreasing. This is illustrated in Fig. 10, left panel, where we have changed the value of $s_\mathrm{max}$ before performing the fit of $\ln (\rho )$, in other words we have varied the value of n in the functional form Eq. 43. We have also checked that the integral itself does not depend on $s_\mathrm{max}$.

Finally we investigate the $\delta _s$ dependence. We have already seen that the LLR coefficients exhibit very little dependence, but it remains to be checked that the same holds for the LLR integrals leading to the phase factor expectation value. In fact, we expect the artefacts to be dominated by order $\delta s^2$ terms [29]. Using $\mu =1.3321$ (from the severe sign-problem region), we carried out simulations with several different values of $\delta s$, reconstructed the LLR coefficients and finally performed the LLR integrals to obtain values of $\langle \mathrm {e}^{i \phi } \rangle$ for this set of $\delta _s$. We then performed a linear extrapolation in $\delta s^2$. Our findings are summarised in Fig. 10, right panel: we indeed find a very small $\delta s^2$ dependence. In fact, the final results for $\langle \mathrm {e}^{i\phi } \rangle$ are more or less independent of $\delta _s$ within statistical error bars. Our numerical findings for $\langle \mathrm {e}^{i\phi } \rangle$ for different truncations can be found in Table 3.

Table 3 Result for the phase factor expectation value for $\mu =1.3321$ as a function of $\delta _s$ for various fit Ansätze. The integral has been computed from $-s_\mathrm{max}$ to $s_\mathrm{max}$ (with folding)

4.5 The phase factor expectation value

We have repeated the analysis outlined in the previous subsection for several values of the chemical potential in the low density region, in the strong sign-problem and in the dense regimes (see Sect. 4.3 for a more formal definition of these regimes). The numerical results are given in the appendix. Each regime has its own challenges:

In the low density regime, the LLR coefficients a(s) are rapidly increasing with s. This implies a rather narrow density-of-states $\rho (s)$, which might approximate a Dirac function $\delta$ if $\mu$ approaches zero. Here, a careful fine-tuning of $\delta _s$ and of the upper integration limit $s_\mathrm {max}$ would be in order. Since this regime is easily accessible by the reweighting approach, we did not further pursue an optimal choice of parameters, but used a generic choice of parameters for a validation of the method only.

In the strong sign-problem regime, our method works best: the results are very robust against the parameter choice. The LLR coefficients show a monotonic behaviour as a function of s, and the choice of even powers of s for the base functions $f_n(s)$ in (42) is converging rapidly: a few non-vanishing coefficients represents hundreds of data points with a $\chi ^2 /$dof well below one.

The dense regime is obtained if the chemical potential takes values close to the heavy quark mass, i.e., its onset value. The sign problem in this regime is mild, and good results are obtained by the reweighting approach. The coefficients a(s) show oscillations around a significant (negative) mean value. Upon reconstructing the density-of-states (see (31)), we find still find a monotonic decreasing $\rho (s)$ (by virtue of the mean value of a), but clearly a significant number of base functions $f_n(s)$ is needed to grasp the oscillatory behaviour, and the method loses its appeal. Insights into the cause of the oscillations would help to develop a new set of base functions $f_n(s)$ that, again with few coefficients, would grasp the essence of the numerical data. For the present paper, we do present LLR results for this regime as well, but observe that the representation of the data with the base functions $f_n(s)=s^{2n}$ failed. Further work in this direction is needed, which we will be presented elsewhere.

Fig. 11

Natural logarithm of $\langle \mathrm {e}^{i\phi } \rangle$ for different values of $\mu$, only the statistical errors are shown. The colour code is as follows: the plain blue points (between $\mu =1.0621$ and $\mu = 1.3721$) have a $\chi ^2$ per degree-of-freedom of order one, the light blue points between 10 and 50, and the white points larger than 50

Finally, we point out our rationale for the approximation of the numerical data for $\ln \rho (s)$ in terms of $f_n(s)$: if few base functions can approximate the data well ($\chi ^2 /\mathrm{dof} < 1$), the bootstrap analysis for the final value of the phase factor expectation values yields small statistical errors, and if the final result is insensitive to the interval size $\delta _s$, we are confident that the LLR approach solves the sign problem in this regime. We have presented evidence for HDQCD in the cases for which the reweighting method can still produce statistical significant results. We also note that if the base function fit fails in the sense that it produces a $\chi ^2 /\mathrm{dof} \ge 100$, it does not necessarily fail to produce a result for the phase factor expectation value close to the true value: it might that fit fails at a large scale in a region of the integration parameter s that is irrelevant to the final result of the integration. We indeed have observed this for dense regime: although the fit fails according to the obtained $\chi ^2 /\mathrm{dof}$, the final results is close to the value known from the reweighting method.

Table 4 Fit result of the phase factor expectation value for the low values of $\mu$ and different truncations

Table 5 Fit result of the phase factor expectation value for the middle-low values of $\mu$ and different truncations

Table 6 Fit result of the phase factor expectation value for the middle-high values of $\mu$ and different truncations

We finally present our main numerical finding. We are interested in $\ln \langle \exp \{ i \phi \} \rangle$ since it is this quantity that enters in e.g. the calculation of the quark density (see (17)):

$$\begin{aligned} \sigma (\mu ) = \frac{T}{V} \frac{ \partial }{\partial \mu } \ln \langle \mathrm{e}^{i\phi } \rangle + \frac{T}{V} \frac{ \partial }{\partial \mu } \ln Z_{PQ}(\mu ) \end{aligned}$$

(49)

Our result for $\ln \langle \exp \{ i \phi \} \rangle$ as a function of the chemical potential $\mu$ is shown in Fig. 11. Further details, such as the quality of the fits, are given in Tables 4, 5, 6 and 7 in the appendix. We have also added these LLR results to the Fig. 1 of Sect. 3 to validate the LLR method against the reweighting data and to demonstrate the quality of the LLR data in the strong sign-problem regime.

Table 7 Fit result of the phase factor expectation value for the high values of $\mu$ and different truncations. We do not give the $\chi ^2$-values larger than 5000

5 Conclusions

We have thoroughly studied QCD with a chemical potential for heavy quarks using the density-of-states approach (LLR version [18, 22]). This approach allows for a determination of the probability distribution of the imaginary part of the quark determinant featuring exponential error suppression. The partition function appears as Fourier transform of this probability distribution. We have bench-marked the LLR results against results from the standard reweighting procedure (in the regime where the latter produces a viable signal-to-noise ratio) and find excellent agreement. We stress, however, that our approach yields an error that is typically smaller by five orders of magnitude.

Due to an (approximate) particle–hole duality at low temperatures, the phase factor expectation value $\langle \exp \{ i \phi \} \rangle (\mu )$ is symmetric around the onset chemical potential $\mu = m$ for which $\langle \exp \{ i \phi \} \rangle = 1$. This suggests an inverted Silver Blaze behaviour: close to the mass threshold, the phase-quenched quark density underestimates the result of the full theory.

Depending on the chemical potential, we found three different regimes which exhibit a different qualitative behaviour of the density-of-states $\rho (s)$:

• (i) In the low density regime, where the theory is almost real, the domain of support of $\rho (s)$ is limited to small values of s as expected.

• (ii) For intermediate values of $\mu$, we find a strong sign problem with $\langle \exp \{ i \phi \} \rangle (\mu )$ reaching values as low as $10^{-6}$ for a small lattice size of $8^4$ (see (28) for the simulation parameters).

• (iii) For chemical potentials close to the onset value, the theory is almost real again. By contrast to the low density regime, however, the density-of-states for the imaginary part, i.e., $\rho (s)$, has a large domain of support, and the corresponding LLR coefficients a(s) show an oscillatory behaviour. It is exceedingly difficult to control the errors of the Fourier transform that is needed to access the phase factor expectation value. Further studies to explore the nature of the oscillations of a(s) is left to future work. We point out, however, the regime close to onset is accessible by reweighting.

In summary, we find that the LLR approach to the probability distribution of the imaginary part of the quark determinant is a viable tool for the whole range of chemical potentials (with a possible exemption near the onset transition). At least for the moderate lattice size explored in this paper, the approach does solve a strong sign problem.

Appendix A: Numerical details

The tables below present details of the fit of the base function expansion depending on the truncation (see Sect. 4.4 for details). In boldface is the fit used for the final results presented in Fig. 11.