Applying Complex Langevin Simulations to Lattice QCD at Finite Density

TL;DR

This paper investigates the use of complex Langevin equations for simulating lattice QCD at finite density, assessing conditions for accurate results and exploring how parameters affect the simulation's fidelity.

Contribution

It demonstrates that complex Langevin simulations can produce correct results at certain parameters and explores how proximity to the SU(3) manifold influences accuracy.

Findings

Good agreement with expected results at small and large chemical potentials.

Accuracy improves with weaker coupling and smaller quark mass.

Proximity to the SU(3) manifold correlates with correct physics outcomes.

Abstract

We study the use of the complex-Langevin equation (CLE) to simulate lattice QCD at a finite chemical potential ($\mu$) for quark-number, which has a complex fermion determinant that prevents the use of standard simulation methods based on importance sampling. Recent enhancements to the CLE specific to lattice QCD inhibit runaway solutions which had foiled earlier attempts to use it for such simulations. However, it is not guaranteed to produce correct results. Our goal is to determine under what conditions the CLE yields correct values for the observables of interest. Zero temperature simulations indicate that for moderate couplings, good agreement with expected results is obtained for small $\mu$ and for $\mu$ large enough to reach saturation, and that this agreement improves as we go to weaker coupling. For intermediate $\mu$ values these simulations do not produce the correct physics. We compare our results with those of the phase-quenched approximation. Since there are indications that correct results might be obtained if the CLE trajectories remain close to the $SU(3)$ manifold, we study how the distance from this manifold depends on the quark mass and on the coupling. We find that this distance decreases with decreasing quark mass and as the coupling decreases, i.e. as the simulations approach the continuum limit.

Tables

Table 1. Table 1: Plaquettes from simulations at β = 5.6 𝛽 5.6 \beta=5.6 , μ = 0 𝜇 0 \mu=0 for various types of simulation and choice of parameters.

	lattice	$\beta$	$\mu$	$dt$	$d{t}_{adaptive}$	cools	start	plaquette
RHMC	${12}^4$	5.6	0.0				ordered	0.43552(5)
RHMC	${16}^4$	5.6	0.0				ordered	0.43556(2)
RLE	${16}^4$	5.6	0.0	0.01	0.00103	0	ordered	0.43566(3)
CLE	${12}^4$	5.6	0.0	0.01	0.000309	5	ordered	0.43667(9)
CLE	${12}^4$	5.6	0.0	0.01	0.000299	15	disordered	0.43672(12)
CLE	${12}^4$	5.6	0.0	0.01	0.000294	100	ordered	0.43686(10)
CLE	${16}^4$	5.6	0.0	0.01	0.000183	5	ordered	0.43681(6)
CLE	${16}^4$	5.6	0.0	0.005	0.000094	5	$dt=0.01$	0.43682(7)

Table 2. Table 2: Chiral condensate from simulations at β = 5.6 𝛽 5.6 \beta=5.6 , μ = 0 𝜇 0 \mu=0 for various types of simulation and choice of parameters.

	lattice	$\beta$	$\mu$	$dt$	$d{t}_{adaptive}$	cools	start	$⟨\overline{\psi }\psi ⟩$
RHMC	${12}^4$	5.6	0.0				ordered	0.2133(11)
RHMC	${16}^4$	5.6	0.0				ordered	0.2158(3)
RLE	${16}^4$	5.6	0.0	0.01	0.00103	0	ordered	0.2160(5)
CLE	${12}^4$	5.6	0.0	0.01	0.000309	5	ordered	0.1993(15)
CLE	${12}^4$	5.6	0.0	0.01	0.000299	15	disordered	0.1996(29)
CLE	${12}^4$	5.6	0.0	0.01	0.000294	100	ordered	0.1985(17)
CLE	${16}^4$	5.6	0.0	0.01	0.000183	5	ordered	0.1968(8)
CLE	${16}^4$	5.6	0.0	0.005	0.000094	5	$dt=0.01$	0.1997(9)

Table 3. Table 3: Unitarity norms F ​ ( U ) 𝐹 𝑈 F(U) from CLE simulations at β = 5.6 𝛽 5.6 \beta=5.6 , μ = 0 𝜇 0 \mu=0 for various choices of parameters.

	lattice	$\beta$	$\mu$	$dt$	$d{t}_{adaptive}$	cools	start	$F\left(U\right)$
CLE	${12}^4$	5.6	0.0	0.01	0.000309	5	ordered	0.1436(28)
CLE	${12}^4$	5.6	0.0	0.01	0.000299	15	disordered	0.1479(41)
CLE	${12}^4$	5.6	0.0	0.01	0.000294	100	ordered	0.1496(31)
CLE	${16}^4$	5.6	0.0	0.01	0.000183	5	ordered	0.1533(15)
CLE	${16}^4$	5.6	0.0	0.005	0.000094	5	$dt=0.01$	0.1487(20)

Table 4. Table 4: Comparison between local observables from CLE simulations at β = 5.7 𝛽 5.7 \beta=5.7 , m = 0.025 𝑚 0.025 m=0.025 on 16 4 superscript 16 4 16^{4} and 16 3 × 36 superscript 16 3 36 16^{3}\times 36 lattices for μ = 0.2 , 0.3 , 0.8 𝜇 0.2 0.3 0.8 \mu=0.2,0.3,0.8 .

plaquette
$\mu$	${16}^4$	${16}^3\times 36$
0.2	0.42365(4)	0.42364(3)
0.3	0.42368(4)	0.42365(3)
0.8	0.43686(6)	0.43706(3)

Equations

$$-i\left(\frac{d}{dt}{U}_l\right)U_l^{-1}=-i\frac{\delta }{\delta U_l}S(U)+{\eta }_l$$

$$⟨{\eta }_l^a(t){\eta }_{l^{\prime }}^b(t^{\prime })⟩={\delta }^{ab}{\delta }_{l{l}^{\prime }}\delta (t-t^{\prime })$$

$$S(U,\mu )=\beta \sum _{\square }\left\{1-\frac{1}{6}Tr[UUUU+(UUUU{)}^{-1}]\right\}-\frac{N_f}{4}\mathrm{Tr}\{ln[M(U,\mu )]\}$$

$$U^{(n+1/2)}$$

$$X_0$$

$$U^{(n+1)}$$

$$X_1$$

$$⟨{\eta }_l^{a(m)}{\eta }_{l^{\prime }}^{b(n)}⟩=\left(1-\frac{3}{2}dt\right){\delta }^{ab}{\delta }_{l{l}^{\prime }}{\delta }^{mn}$$

$$\mathrm{Tr}\{ln[M(U,\mu )]\}\to {\xi }^{†}\{ln[M(U,\mu )]\}\xi ,$$

$$⟨{\xi }^{\ast i(m)}(x){\xi }^{j(n)}(y)⟩={\delta }^{ij}{\delta }_{xy}{\delta }^{mn}$$

$$f_{\max }=\stackrel{\mathrm{MAX}}{l,i,j}\left|f_{ij}(l)\right|$$

$$d{t}_{\mathrm{adaptive}}=\frac{dt}{f_{\max }}$$

$$F(U)=\frac{1}{4V}\sum _{x,\mu }Tr\left[U^{†}U+(U^{†}U{)}^{-1}-2\right]\ge 0$$

Full text

I Introduction

QCD at finite baryon/quark-number density describes nuclear matter, the constituent of neutron stars and the interiors of heavy nuclei. We are interested in calculating the phase diagram for nuclear matter. Knowing the properties of nuclear matter can yield an equation-of-state and a better description of neutron stars. In addition it can yield information on the interaction of nuclear matter with particles passing through it. Nuclear matter at high temperatures undergoes a transition to a quark-gluon plasma. Such transitions can be observed in relativistic heavy-ion collisions. While for low densities this transition is expected to be a crossover, it is believed that for high densities it becomes a first-order transition. The critical end-point is the second-order transition (critical point) where this change occurs.

Finite density QCD has a sign problem which prevents the direct application of standard lattice QCD simulation methods, which rely on importance sampling. When finite density is implemented by use of a quark-number chemical potential, $\mu$, the sign problem manifests itself by making the fermion determinant complex, with a real part of indefinite sign. Simulations using the complex Langevin equation (CLE) Parisi:1984cs ; Klauder:1983nn ; Klauder:1983zm ; Klauder:1983sp can accommodate such complex actions. However, the CLE can only be shown to yield correct values for observables if the space over which the fields evolve is compact, the drift (force) term is holomorphic in the fields, and the solutions are ergodic Aarts:2009uq ; Aarts:2011ax ; Nagata:2015uga ; Nishimura:2015pba ; Nagata:2016vkn ; Aarts:2017vrv ; Seiler:2017wvd ; Aarts:2017hqp ; Nagata:2018net .

For QCD, implementation of the CLE requires extending the gauge-field manifold from $SU(3)$ to $SL(3,C)$. Keeping the action holomorphic in the fields except on a space of measure zero requires making the action a function of the gauge fields and their inverses only. Then the drift term is meromorphic in the gauge fields, having poles at the zeros of the fermion determinant. Early attempts at simulations were frustrated by runaway behaviour which could not be controlled by adaptively decreasing the updating interval. Recently it was discovered that, for small enough couplings, this behaviour could be tamed by ‘gauge-cooling’, gauge transforming the fields after each update to keep them as near as possible to the $SU(3)$ manifold Seiler:2012wz . When this is done, the gauge fields appear to evolve over a compact manifold, at least for weak enough coupling. It remains to be determined for what if any range of quark masses, lattice spacings, chemical potentials and temperatures, and for what choices of lattice actions, the CLE produces correct values for chosen observables, despite the presence of poles in the drift term and/or in the operators defining these observables.

Extensive studies have been performed of heavy dense lattice QCD at finite $\mu$ using the CLE. Aarts:2008rr ; Aarts:2013uxa ; Aarts:2014bwa ; Aarts:2016qrv ; Langelage:2014vpa ; Rindlisbacher:2015pea . Some CLE simulations of lattice QCD have been performed at finite $\mu$ with lighter quark masses, both on small lattices and at finite temperatures Sexty:2013ica ; Aarts:2014bwa ; Fodor:2015doa ; Nagata:2016mmh ; Tsutsui:2018jva ; Scherzer:2018udt . Preliminary work directed towards zero temperatures has been reported Ito:2018jpo . In heavy dense lattice QCD at finite chemical potential, the zeros of the fermion determinant (poles in the drift term) only affect the CLE results very close to the transition, provided one excludes the regions where the real part of the fermion determinant is negative, when necessary. Alternatively, good results for heavy quarks can be obtained by restricting the length of CLE trajectories to keep them close to the $SU(3)$ manifold. For lighter quarks, CLE simulations are found to produce good results at high temperatures (above the finite temperature phase transition), but the situation at lower temperatures is less clear.

We simulate zero temperature lattice QCD with 2 flavours (tastes) of staggered quarks at $\mu$ values from zero to saturation, using the CLE. Here we are interested in the phase transition from hadronic to nuclear matter which should occur at $\mu\sim m_{N}/3$. Random matrix theories (RMT) related to QCD at finite $\mu$ suggest that when CLE simulations fail, they produce the results of the phase-quenched model (the theory where the fermion determinant is replaced by its magnitude), which has a transition to a superfluid state with a pion-like condensate at $\mu\approx m_{\pi}/2$. This has been observed by Mollgard and Splittorff Mollgaard:2013qra who simulate the Osborn RMT Osborn:2004rf ; Bloch:2012bh and find that the CLE fails for small masses, approaching phase-quenched results for small-enough masses. They suggest a solution in a subsequent paper Mollgaard:2014mga . Bloch et al. Bloch:2017sex using the Stephanov RMT Stephanov:1996ki (which has a non-trivial phase structure) find that the CLE generates phase-quenched results. (Note that other random matrix CLE simulations seem more optimistic Nagata:2016alq .) For this reason we also perform RHMC simulations of the phase-quenched theory at the same values of $\beta=6/g^{2}$ and quark mass $m$ over the same range of $\mu$ values and on the same lattice size as the full theory, for comparison. Hence it is important that we choose $m$ such that $m_{N}/3$ is significantly larger than $m_{\pi}/2$. We perform our simulations at $\beta=5.6$, $m=0.025$ (in lattice units) on a $12^{4}$ lattice and at $\beta=5.7$, $m=0.025$ on a $16^{4}$ lattice. Preliminary results were reported at Lattice 2015–2018. See Sinclair:2018rbk and its references to our earlier talks.

At $\beta=5.6$, the CLE measurement of the plaquette for $\mu=0$ exhibits a systematic error of $\approx 0.31$ %, while at saturation it shows a systematic error of $\approx 1.43$ %. These should be compared with the increase in value of the plaquette over this range which is $\approx 9.2$ %. At $\beta=5.7$ the CLE plaquette measurement has a systematic error of $\approx 0.16$ %, while at saturation the systematic error is $\approx 0.3$ %. The increase in the known value of the plaquette over this range is $\approx 6.6$ %. We note that for both $\beta$s, the plaquette values at saturation show excellent agreement with those obtained from CLE simulations of $SU(3)$ lattice gauge theory in the absence of quarks, as expected.

At $\beta=5.6$, the CLE measurement of the chiral condensate at $\mu=0$ lies $\approx 6.9$ % below the correct value, while at $\beta=5.7$ the chiral condensate predicted by the CLE is $\approx 1.22$ % lower than the known value. For $\mu<m_{\pi}/2$ there is a similar improvement in the CLE-predicted chiral condensate between $\beta=5.6$ and $\beta=5.7$. At $\beta=5.7$ and $0.5\leq\mu\leq 0.9$ the CLE produces values of the chiral condensate, quark-number density and plaquette in good agreement with the phase-quenched theory. However, the Wilson Line/ Polyakov Loop indicates that the fermion determinant is still complex. Although the CLE results appear to be approaching the correct physics for very small and large $\mu$ as the coupling decreases towards the continuum limit, they still fail to produce the expected physics in the transition region. The transition to nuclear matter appears to start at $\mu$ even less than $m_{\pi}/2$ instead of $\mu\approx m_{N}/3$.

Since there are indications that the CLE might produce correct results when its trajectories remain close to the $SU(3)$ manifold, we perform a systematic study of how the average distance to this manifold (measured using the ‘unitarity norm’) depends on the quark mass and the coupling. We find that this norm decreases with decreasing quark mass and with decreasing coupling, i.e. as we approach the continuum limit.

In section 2 we present the formulation of the CLE, we use. Section 3 describes our simulations at $\beta=5.6$ and $\beta=5.7$ and presents results. In section 4 we present our simulations to determine how the unitarity norm depends on quark mass $m$ and lattice coupling $g$. Section 5 gives a summary, discussion and conclusions.

II Complex Langevin for finite density Lattice QCD

If $S(U)$ is the gauge action after integrating out the quark fields, the Langevin equation for the evolution of the gauge fields $U$ in Langevin time $t$ is:

[TABLE]

where $l$ labels the links of the lattice, and $\eta_{l}=\sum_{a}\eta^{a}_{l}\lambda^{a}$. Here $\lambda_{a}$ are the Gell-Mann matrices for $SU(3)$. $\eta^{a}_{l}(t)$ are Gaussian-distributed random numbers normalized so that:

[TABLE]

We note in passing that the discretized Langevin Equation is the limiting case of the Hybrid Molecular Dynamics method where each trajectory has only a single update.

The complex-Langevin equation has the same form except that the $U$s are now in $SL(3,C)$. $S$, now $S(U,\mu)$ is

[TABLE]

where $M(U,\mu)$ is the unimproved staggered Dirac operator with quark-number chemical potential $\mu$, for a single staggered fermion field (corresponding to 4 continuum flavours). Note: backward links are represented by $U^{-1}$ not $U^{\dagger}$. Note also that we have chosen to keep the noise term $\eta$ real.

To simulate the time evolution of the gauge fields we use the partial second-order formalism of Fukugita, Oyanagi and Ukawa. Ukawa:1985hr ; Fukugita:1986tg ; Fukugita:1988qs For an update of the fields by a ‘time’ increment $dt$, this gives:

[TABLE]

where $\gamma=\frac{1}{2}+\frac{1}{4}dt$ and the Gaussian noise $\eta$ is normalized such that:

[TABLE]

To proceed, we replace the spacetime trace with a stochastic estimator $\xi$

[TABLE]

$\xi$ is a vector over space-time and colour of Gaussian random numbers, normalized such that:

[TABLE]

which means, in particular, that the $\xi$s in $X_{0}$ and $X_{1}$ are independent, unlike the $\eta$s. After performing $\frac{\delta}{\delta U}$ of $ln(M)$ it is useful to use the cyclic property of the trace to rearrange the terms proportional to $U$ and $U^{-1}$ prior to introducing the stochastic estimators, so that this operator is antihermitian when $\mu=0$ and $U$ is unitary. That way, in this special case, the complex Langevin equation becomes the real Langevin equation.

We apply adaptive updating. If $f_{ij}(l)$ are the components of the drift term, we define

[TABLE]

where $l$ runs over the links of the lattice. $i=1,2,3$, $j=1,2,3$ are the colour indices. Then, if $f_{\rm max}>1$, we replace the input updating increment $dt$ by the adaptive increment

[TABLE]

for the current update. Because the Dirac operator is often ill-conditioned, we use 64-bit floating point precision throughout.

After each update, we adaptively gauge fix to the gauge which minimizes the unitarity norm:

[TABLE]

which equals [math], if and only if $U$ is unitary. $V$ is the space-time volume of the lattice. $F(U)$ is a measure of the lattice averaged distance of the gauge fields from the $SU(3)$ manifold. We implement gauge cooling after each updating following the method described in Seiler:2012wz , equations 8–10. Here $\epsilon=dt$ the input increment for our CLE updating, and we choose $\alpha=1/4$. We make this updating adaptive by the following ansatz. If $G(n)=\sum_{a}\lambda_{a}G_{a}(n)$, $G(n)_{norm}=\sqrt{{\rm Tr}(G(n)^{2})}$, then if $\alpha G_{norm}(n)>1$, we replace $G(n)$ by $G(n)/[\alpha G_{norm}(n)]$ for this gauge-cooling step.

III CLE simulations of lattice QCD at finite $\mu$ and zero temperature

We simulate two-flavour lattice QCD at $\beta=5.6$, $m=0.025$ on a $12^{4}$ lattice, and at $\beta=5.7$, $m=0.025$ on a $16^{4}$ lattice, at $0\leq\mu\leq 1.5$. Note $\mu=1.5$ is well into the saturation domain, where all fermion levels are filled and the fermion number density, normalized to 1 staggered fermion field (4 flavours/tastes)is 3. The input $dt$ is chosen to be $0.01$. Our runs for individual $\mu$ values vary between $10^{6}$ and $3\times 10^{6}$ updates in length. Discarding the first fifth of each run for equilibration, this makes the length of each run in Langevin time units vary between 80 and over 1000, with the shortest runs being in the saturation regime. For most of the runs, we choose 5 gauge-cooling steps after each CLE update.

To test our choice of $dt$ and the number of gauge-cooling steps per update, as well as to observe the finite-lattice-size effects we ran a number of test runs at $\beta=5.6$, $\mu=0$ where we could compare our CLE simulations with those performed using the exact RHMC algorithm and with the corresponding real Langevin equation (RLE). The results of these simulations are summarized in tables 1,2,3

There is good agreement between the RLE and the exact RHMC observables when we used an input $dt=0.01$. We note that changing the number of gauge-cooling steps has little effect on the chiral condensates, plaquettes and unitarity norms for the CLE simulations. Neither does changing the starting conditions for the runs. This was also found to be true at $\mu=0.5$ where we ran CLE simulations from an ordered start with 5 gauge-cooling steps per update and from a disordered start with 10 gauge-cooling steps per update. The chiral condensate shows a small but significant finite volume effect in going from a $12^{4}$ to a $16^{4}$ lattice. Reducing $dt$ from $0.01$ to $0.005$ on the $16^{4}$ lattice increases the CLE measurement of the chiral condensate. Considering only the leading dependence on $dt$ and hence on $dt_{adaptive}$ which is linear, we predict $\langle\bar{\psi}\psi\rangle\approx 0.2025$ in the limit $dt\rightarrow 0$. We note that the difference between this and the true (RHMC) value is $\approx 0.013$, compared with the difference between the $dt=0.01$ value and the true value on the $12^{4}$ lattice is $\approx 0.014$. Since we know that one needs a smaller $dt$ on a larger lattice we believe that $dt=0.01$ is adequate for the $12^{4}$ lattice, at least at $\mu=0$. For $\beta=5.7$ on a $16^{4}$ lattice at $\mu=0$, the RHMC simulations predict that the chiral condensate, the observable most sensitive to simulation methods and parameters, $\langle\bar{\psi}\psi\rangle=0.1752(2)$, the RLE simulations with $dt=0.01$ give $\langle\bar{\psi}\psi\rangle=0.1756(3)$, while the CLE simulations with $dt=0.01$ and 5 gauge-cools/update yield $\langle\bar{\psi}\psi\rangle=0.1731(10)$. While the CLE probably still has some systematic error, we consider it to be small enough to be acceptable. Hence we choose $dt=0.01$ and 5 gauge-cools/update for all CLE simulations except those mentioned in the tables. Here we implicitly assume that adaptive rescaling of $dt$ is sufficient to tame any wild fluctuations which are produced by going to non-zero $\mu$.

For comparison, we perform RHMC simulations of the phase-quenched theory with the same parameters and lattice sizes. At each $\beta$ and $\mu$ we run 10,000 length-1 trajectories (except at $\mu=0$, where we ran 20,000 trajectories). The observables for the phase-quenched theory should remain constant at their $\mu=0$ values up to $\mu\approx m_{\pi}/2$ in the limit that the explicit symmetry-breaking term’s coefficient $\lambda$ vanishes, up to finite temperature and volume corrections, whereas in the full theory, these observables should remain constant until $\mu\approx m_{N}/3$. (See Kogut:2002zg for definition of $\lambda$). Note that QCD with an isospin chemical potential $\mu_{I}$ is identical to the phase-quenched theory with chemical potential $\mu=\mu_{I}/2$. This is because the 2-flavour phase-quenched theory is a theory with 1 quark which is a colour triplet and 1 conjugate-quark which is a colour anti-triplet and has the opposite parity from the regular quark. These can form a pion-like state of 1 quark and 1 conjugate quark, which can form a colourless quark-number breaking condensate. If we write $\Psi$ as a ‘flavour’ doublet with components the normal quark and the anti-conjugate-quark, then the symmetry breaking term in the Lagrangian is $i\lambda\bar{\Psi}\gamma_{5}\tau_{2}\Psi$.) For $\beta=5.6$, $m=0.025$, $m_{\pi}/2\approx 0.21$ and $m_{N}/3\approx 0.33$ Bitar:1990cb while for $\beta=5.7$, $m=0.025$, $m_{\pi}/2\approx 0.194$ and $m_{N}/3\approx 0.28$ Brown:1991qw ; Schaffer:1992rq , so these 2 transitions should be distinguishable. At saturation, the quark-number density should be 3 (one quark of each colour at each site), the chiral condensate should vanish, and the fermions should decouple from the gauge fields, so gauge observables such as the plaquette should have their pure gauge values.

Figure 1 shows the average plaquettes as functions of $\mu$ for $\beta=5.6$, $m=0.025$ on a $12^{4}$ lattice and $\beta=5.7$, $m=0.025$ on a $16^{4}$ lattice, from our CLE simulations. The corresponding results for the phase-quenched theory are shown for 2 different values of the symmetry breaking parameter $\lambda$ indicating that there is relatively little $\lambda$ dependence. For $\mu=0$ at $\beta=5.6$, while the difference between the CLE value and the exact value is not large, it is significant. At $\beta=5.7$ the relative difference is almost a factor of 2 smaller. At saturation, for both $\beta$s, the plaquette for the phase-quenched theory is identical to that of the pure gauge theory within statistical errors. For both $\beta$s the plaquette value for the full theory at saturation predicted by the CLE, agrees within statistical errors with that predicted by CLE simulations of the pure gauge theory. The values predicted by the CLE for pure gauge theory differ significantly from the correct value (see section 4). However, this difference is much smaller for $\beta=5.7$ than for $\beta=5.6$.

The chiral condensates for $\beta=5.6$ and $\beta=5.7$ are shown as functions of $\mu$ in figure 2 for both our CLE simulations of QCD at finite $\mu$ and RHMC simulations of the phase-quenched theory. At $\mu=0$ the CLE value of this condensate at $\beta=5.6$ is $\approx$ 6.6% too low, while at $\beta=5.7$ it is $\approx$ 1.2% too low, a considerable improvement. As $\mu$ is increased from zero, the condensates for both $\beta$s start to fall for $\mu<m_{\pi}/2$ rather than for $\mu\approx m_{N}/3$. Therefore, even though $\beta=5.7$ shows an improvement over $\beta=5.6$, it is still a worse approximation to the physics than is the phase-quenched theory. It is still unclear if going to a much weaker coupling will produce the correct results or those of the phase-quenched theory. At both $\beta$s the condensate reaches its saturation value of zero at large $\mu$.

Figure 3 shows the quark-number densities for our $\beta=5.6$ and $\beta=5.7$ simulations. On this scale, the CLE results for the full theory and RHMC results for the phase-quenched theory look almost identical, except on their approach to saturation. Figure 4, shows an expanded version of the low $\mu$ region which indicates that the apparent agreement is only because the number densities are so close to zero for small $\mu$. The results for the 2 theories are closer for $\beta=5.7$ than for $\beta=5.6$. This would be expected if the CLE produces phase-quenched results in the weak-coupling limit, but also could indicate that it produces correct results in that limit. Note that it is difficult, if not impossible, to determine the positions of the transitions from these quark-number density graphs.

Figure 5 shows the average unitarity norms for the CLE simulations at $\beta=5.6$ and $\beta=5.7$ as functions of $\mu$. In both cases this norm decreases from its value at $\mu=0$ as $\mu$ is increased reaching a minimum around $0.5$ or $0.6$ before increasing to a maximum at saturation, which is the value for CLE simulations of the pure gauge theory at the same $\beta$ value. The values of the unitarity norms for $\beta=5.7$ lie below the values of those at $\beta=5.6$ for the same $\mu$. It has been suggested, on the basis of simulations of lattice QCD in the heavy-dense limit that there is some value of the unitarity norm around $0.1$ below which the CLE will produce correct results Aarts:2016qrv 111We have investigated the suggestion from this paper that one should terminate the simulation at a given $\beta$ and $\mu$ before the fluctuations in the observables increase substantially. For our simulations, while this improves the results at $\mu=0$, it does not prevent the precocious onset of the transition. In addition, such restricted simulations are short and their length decreases with increasing $\mu$, essentially vanishing above the transition, so it is doubtful that the system has time to equilibrate.. From our simulations at $m=0.025$ there does appear to be such a value for $\mu=0$, between $0.13$ and $0.15$. However, if there is such a value for $\mu>0$, it decreases with increasing $\mu$, lying below the value of this norm for $\beta=5.7$, at least through the transition region.

At $\beta=5.7$, the value of the unitarity norm drops by roughly a factor of 2 from $\mu=0$ where it is $\approx 0.131$ to $\mu=0.3$ at the upper end of the transition region, where it is $\approx 0.067$. From there it falls by an order of magnitude, reaching a minimum between $\mu=0.5$ where it is $\approx 0.0058$ and $\mu=0.6$ where it is $\approx 0.0054$. It is tempting to suggest that the CLE will produce correct results for $\mu$s around this rather broad minimum. Moreover, for $0.5\leq\mu\leq 0.9$ the plaquette, the chiral condensate, and the quark-number density are in good agreement with those of the phase-quenched approximation. (Above $\mu=0.9$ we are in the regime controlled by saturation a lattice artifact.) Either, in the large $\mu$ region, the full and phase-quenched theories give the same physics or this is a sign that the CLE breakdown produces phase-quenched results as suggested by random matrix theory. We shall have more to say about this shortly.

To test whether the assumption that the CLE is valid for $\beta=5.7$, $m=0.025$ for $\mu$ near the unitarity-norm minimum is reasonable, even though it fails for small $\mu$, we examine the distribution of values of the chiral condensate for $0\leq\mu\leq 0.9$. (For $\mu>0.9$ the system is influenced by saturation, a lattice artifact.) The chiral condensate is chosen since it has poles at the same places as those of the drift term. Hence if these poles are approached too closely, which invalidates the CLE, it should show large (non-Gaussian) excursions in its distribution. Figure 6 presents the histogram of values of $\langle\bar{\psi}\psi\rangle$ at $\mu=0$, which is typical of the distributions for small $\mu$. This histogram has long tails, with a few outliers, which indicates that the poles in the drift term as well as those in $\bar{\psi}\psi$ are being approached. The fact that good results are obtained for $\mu=0$ indicates that although the trajectory of the system approaches the poles this does not necessarily cause a breakdown of the CLE. However, as $\mu$ is increased, and the tails slowly become more prominent, the CLE does start to deviate from the correct physics. In figure 7 we show histograms at $\mu=0.2$, close to $m_{\pi}/2$ and hence at the beginning of the transition region if the physics were that of the phase-quenched theory, and at $\mu=0.3$ close to $m_{N}/3$ and thus to the transition expected for QCD at finite $\mu$. At $\mu=0.2$ the tails appear to be close to maximal. Beyond this, they decrease, and are noticeably smaller at $\mu=0.3$. Figure 8 shows the histograms for $\mu=0.6$ at the minimum of the unitarity norm, and $\mu=0.8$ in the large $\mu$ regime. By $\mu=0.6$ the non-gaussian tails have almost vanished and remain insignificant over the range $0.5\leq\mu\leq 0.9$ as can be seen in the histogram at $\mu=0.8$. In this high $\mu$ regime we have also looked at the histograms of quark-number density distributions (since the quark-number operator also has poles at the zeros of the fermion determinant), which are also well-behaved over this range of $\mu$s.

Figure 9 shows the real parts of the Wilson Lines (Polyakov Loops) and the Inverse Wilson Lines as functions of $\mu$ for our $\beta=5.6$, and $\beta=5.7$ from our CLE simulations. The imaginary parts of these quantities are small, consistent with zero. We include the Wilson Lines from the corresponding phase-quenched RHMC simulations with $\lambda=0.001$ for comparison. What we notice is that once $\mu$ is above the transition, if not before, the Wilson lines start to diverge from their inverses. This means that the fermion determinant is complex. The fact that the inverse Wilson line lies above the Wilson line up to some $\mu$ between $0.8$ and $0.9$ indicates that the phase of the Wilson line and that of the fermion determinant are positively correlated. Such a correlation was anticipated in deForcrand:1999cy , and observed in simulations of heavy-dense lattice QCD Seiler:2012wz ; Aarts:2008rr ; Langelage:2014vpa ; Rindlisbacher:2015pea . It is interesting to note that the crossover point where the phase of the determinant vanishes lies near to the quark-number density of 1.5 where the fermi states are half-filled. This behaviour was predicted and observed in heavy-dense lattice QCD Rindlisbacher:2015pea .

The Wilson line and the inverse-Wilson line increase with increasing $\mu$ reaching a maximum just before the effects of saturation start to be felt. At $\beta=5.7$ the fact that local observables measured in CLE simulations agree with phase-quenched results, despite the fact that the fermion determinant is complex, is possible because the phase of the fermion determinant can be determined by a few low lying eigenvalues of the Dirac operator, whereas other quantities are determined by the distribution of eigenvalues. This was discussed most clearly in Bloch:2017sex . We also note that the CLE Wilson and inverse-Wilson lines do not show any significant effect from the pseudo-transition at $\mu=m_{\pi}/2$, unlike the Wilson line for the phase-quenched theory, which might indicate that the CLE will not eventually yield phase-quenched results in the transition region.

Because the presence of the chemical potential introduces an asymmetry between space and time which tends to amplify the effects of temperature, the question arises as to whether the $12^{4}$ lattice for $\beta=5.6$ and the $16^{4}$ lattice for $\beta=5.7$, while being good approximations to zero temperature at $\mu=0$, remain good approximations to zero temperature as $\mu$ is increased. The increase in value of the Wilson Line as $\mu$ is increased, discussed in the previous paragraph, emphasizes the possibility that some of the $\mu$ dependence of observables other than the Wilson line could actually be finite temperature effects. To check this we performed test CLE runs on a $16^{3}\times 36$ lattice with $\beta=5.7$ and $m=0.025$ at $\mu=0.2$ and $\mu=0.3$, which bracket the transition region, and at $\mu=0.8$, a large value of $\mu$, but not so large as to be influenced by saturation. At all 3 of these $\mu$ values the Wilson lines are consistent with zero. The values of the the plaquette, chiral condensate, and quark-number density and unitarity norm measured in these simulations are compared with those obtained for the corresponding $16^{4}$ lattice simulations in table 4.

For $\mu=0.2,0.3$ there is good agreement between these local observables for $N_{t}=16$ and $N_{t}=36$. For $\mu=0.08$ there are small but noticeable differences in the plaquettes and chiral condensates between the 2 lattice sizes, while the quark-number densities and unitarity norms are in good agreement. We conclude that the finite temperature effects from using $N_{t}=16$ are acceptably small over the whole range of $\mu$ values.

IV Dependence of the unitarity norm on quark mass and coupling

It has been observed that the CLE is better behaved if its trajectories remain near to the $SU(3)$ manifold, i.e. if the unitarity norm remains small. It is therefore of interest to determine how the unitarity norm depends on the parameters of the theory.

First we consider how the unitarity norm depends on the quark mass $m$. In the previous section, we have observed that, for fixed $m$ and $\beta$, the local maximum for small $\mu$ occurs at $\mu=0$. The global maximum occurs at saturation. Since this maximum is the value for the pure gauge theory at this $\beta$, it does not depend on $m$. Therefore we shall determine the $m$ dependence of the unitarity norm at $\mu=0$, which is thus relevant for small $\mu$.

Since the unitarity norm at $\mu=0$ appears to decrease with increasing $\beta$ and we shall be interested in $\beta\geq 5.6$, we study the $m$ dependence with $\beta$ fixed at $\beta=5.6$. We perform CLE simulations with $0.01\leq m\leq\infty$. For $m=0.01$ we use a $16^{4}$ lattice. For $m=0.025$ we use both $16^{4}$ and $12^{4}$ lattices, while for $m=0.05$, $m=0.1$, $m=0.25$, $m=0.5$ and $m=\infty$ we use $12^{4}$ lattices. Except for $m=0.025$ and $m=\infty$ on $12^{4}$ lattices, each run uses $3\times 10^{6}$ updates.

In figure 10 we plot the unitarity norm versus $1/m$ at $\beta=5.6$. We note that it decreases by almost an order of magnitude as $m$ is decreased from infinity to $0.01$.

For a given $\beta$, the unitarity norm has its maximum for $\mu$ at saturation, for which it is mass independent. Hence we choose to perform a CLE at saturation for each $\beta$, that is we simulate the pure gauge theory for that $\beta$, knowing that this will give the upper bound to the unitarity norm for that $\beta$. Moreover, since there are no quarks, we do not have to worry that we should really change $m$ when we change $\beta$ so as to keep on a line of constant physics. Since there are no quarks, these simulations are fast, and for given $\beta$ one can use a much smaller lattice than for the same $\beta$ with light quarks. Without quarks, the drift term is holomorphic in the gauge fields, so the CLE should produce results correct to order $dt^{2}$, provided the fields evolve on a compact domain.

We perform CLE simulations for a selection of $\beta$s in the range $5.6\leq\beta\leq 7.0$. $\beta=5.6$ and $\beta=5.7$ simulations were performed on $12^{4}$ lattices, $\beta=5.8$, $\beta=5.9$, $\beta=6.0$ and $\beta=6.2$ simulations were performed on $16^{4}$ lattices, $\beta=6.5$ simulations were performed on a $24^{4}$ lattice and $\beta=7.0$ simulations were performed on a $32^{4}$ lattice. First we find that at $\beta=5.7$, if we start on the $SU(3)$ manifold, the CLE trajectory stays on the $SU(3)$ manifold for at least $10^{7}$ updates. However, if we start slightly off this manifold, the system evolves away from the $SU(3)$ manifold to a region where the unitarity norm fluctuates around a stable, non-zero value. This value appears to be independent of how far the starting point is from the $SU(3)$ manifold. We assume similar behaviour for $\beta>5.7$ and start the simulations at larger $\beta$s away from the $SU(3)$ manifold. At each $\beta$, we have at least 1 run with a non-$SU(3)$ start of $5\times 10^{6}$ updates or more. The measured unitarity norms are plotted as a function of $\beta$ in figure 11. These decrease as $\beta$ increases, that is as the coupling $g$ decreases. In fact over the range of $\beta$s considered, this norm decreases by more than an order of magnitude.

In figure 12 we plot the average plaquette from our CLE simulations and compare it to the ‘exact’ value from a Monte-Carlo simulation. Except for $\beta=5.6$ these values are close and get closer as $\beta$ is increased. This difference for the lower $\beta$s is too large to be due to statistical errors or the inexact nature of Langevin simulations. (This is surprising since the drift term is holomorphic in the fields, and the region spanned by these CLE simulations appears to be bounded.) Hence it is a systematic of the CLE, which must be due to the distribution of the plaquette values not falling off fast enough at the boundaries or non-ergodicity which could indicate that there are other regions of the $SL(3,C)^{4V}$ space which are not accessible to these simulations. For $\beta=5.6$ where this systematic error is unacceptably large, we have observed one large excursions in some of our runs. At all $\beta$s the distributions of plaquette values do have ‘tails’ or ‘skirts’ as do those of the unitarity norms.

Because of the size of the systematic errors for pure gauge theory at $\beta=5.6$, we have checked to see if any of this can be due to inadequate gauge-cooling by running simulations with 5, 7 and 100 cooling steps with $dt=0.01$. The plaquette values from all 3 are in good agreement, as are the unitarity norms. We also performed a simulation with 100 cooling steps per update but with $dt=0.005$, and find good agreement with the plaquettes and unitarity norms from our $dt=0.01$ simulations. For $\beta=5.7$ we have performed simulations with 5, 6, 7, 8 and 10 gauge-cooling steps, with different starting configurations. For the case with 5 gauge-cooling steps per update, where we used an ordered start, the system stayed in $SU(3)$ (RLE) and the plaquette was very close to that from the (exact) Monte Carlo value. For the other 4 choices, the plaquettes were in agreement and only slightly above the exact value. The agreement between the CLE and exact simulations improves with increasing $\beta$.

V Summary, discussion and conclusions

We have performed CLE lattice simulations of lattice QCD at zero temperature and finite $\mu$ at $\beta=5.6$, $m=0.025$ and at $\beta=5.7$, $m=0.025$. Neither $\beta$ shows the expected physics in the transition region. Whereas one expects that the transition from hadronic to nuclear matter should occur at $\mu\approx m_{N}/3$, these simulations show transitions for $\mu<m_{\pi}/2$. For the weaker coupling ($\beta=5.7$) for $\mu$ close to zero, the CLE produces values of the observables which are close to the correct results, and considerably better than those for $\beta=5.6$. For $\mu$ large enough to produce saturation, where all available fermion states are filled, both $\beta$s indicate that the quarks have decoupled leaving us with a pure $SU(3)$ gauge theory. In both cases the plaquette observable agrees with that from a CLE simulation of the pure gauge theory. This value is much closer to the exact result for the weaker coupling. For $\beta=5.7$ at very small $\mu$ and for a significant range of $\mu$s above $m_{N}/3$, there is good agreement between these CLE simulations and the exact (RHMC) simulations of the phase-quenched theory.

For $\mu>m_{N}/3$ we have not seen any sign of new exotic phases, such as a colour-superconducting phase. Nor is there any indication of a transition to quark-matter below saturation. There is also no indication of a difference between the full theory and its phase-quenched approximation for $\mu\geq 0.5$ while not so high that saturation is a dominant influence, except that the Wilson Line indicates that the fermion determinant is complex for the full theory. The apparently gaussian distribution of chiral condensate measurements (and of quark-number density measurements) suggests that the CLE should be reliable in this domain. A difference between the full theory and the phase-quenched theory might be expected because of the existence of a pion-like superfluid phase in the latter. If the full theory does have local observables identical to the phase-quenched theory in this region, and this is not an artifact of the CLE, it might be possible to check this by reweighting from the phase-quenched theory to the full theory in the high $\mu$ domain.

Because there is some indication that the CLE is more likely to produce correct results if the trajectories stay close to the $SU(3)$ manifold, we have performed CLE simulations over a range of quark masses at a fixed coupling, and over a range of couplings at infinite quark mass (pure $SU(3)$ gauge theory). What we find is that the average distance of the trajectories from $SU(3)$ as determined by the unitarity norm decreases as the coupling and quark mass are decreased, that is as we approach the continuum limit. This gives us some hope that the CLE might produce the correct physics in this limit. However, random matrix theories suggest that it might produce phase-quenched results. Since the simulations described in section 3 do not rule out either possibility, further simulations at weaker couplings, which require larger lattices, are needed. Since our conclusion that the CLE fails to observe the transition from hadronic to nuclear matter at the expected $\mu$ value is based on measurements of the chiral condensate and quark-number density, both of which are expectation values of operators with poles at the same place as those in the drift term, another possibility needs to be considered. That is the possibility that the simulation is correct and produces the correct values for the expectation values of non-singular operators but fails for such singular operators.

Other methods have been suggested to try and obtain correct results from CLE simulations. One is to add terms to the action which cause the CLE to avoid the poles, with coefficients, which when taken to zero, yield the original action Nagata:2018mkb . The question then is can these coefficients be taken small enough to allow them to be continued to zero, without them losing their effectiveness in avoiding the poles. Preliminary results on very small lattice look promising, but it remains to be seen if this method will work on larger lattices.

A second method is based on the observation that one needs to keep the unitarity norm small for the CLE to work. This is achieved by changing the dynamics of the CLE by adding a force in the direction of decreasing unitarity norm to the drift term, with a coefficient which can be made arbitrarily small Attanasio:2018rtq . This additional force should be irrelevant (in the re-normalization group sense) so that it will vanish as the lattice spacing goes to zero. It should also vanish when the gauge fields lie on the $SU(3)$ manifold. Such a force will not be holomorphic in the gauge fields, so that adding it to the drift term could completely destroy any relationship between the CLE and the physics contained in the original functional integral, so careful testing is needed.

Another possible reason for the failure of the CLE has been discussed by Block and Schenk Bloch:2017jzi . Their claim is that part of the problem with the usual application of the CLE to lattice QCD at finite $\mu$ is the use of stochastic estimators for the traces of the inverses of the poorly conditioned Dirac operator. They promote the use of newer methods to replace these stochastic estimators, which produce better results.

The use of other actions should be investigated. For example, the phase structure of lattice QCD at zero temperature should be studied as a function of $\mu$ using the CLE with Wilson fermions, since Wilson fermions preserve the continuum order of the zeros of the fermion determinant, which is equal to the number of flavours. For staggered fermions this order is decreased by a factor of 4 by ‘taste’ breaking, and is only recovered in the continuum limit.

Application of the CLE to finite temperature QCD at finite $\mu$, in particular with regard to the effect of finite $\mu$ on the transition of hadronic/nuclear matter to a quark-gluon plasma, should be studied, continuing the pioneering work of Fodor et alFodor:2015doa . It has been observed that the CLE should work better in this case, since finite temperatures move the zeros of the fermion determinant away from the CLE trajectories. In the case of 2-flavour staggered quarks, our studies indicate that this is only likely to work (at least for $m=0.025$) if $\beta=5.6$ lies in the low temperature phase. This requires $N_{t}\gtrsim 12$.

Acknowledgments

DKS is supported in part by US Department of Energy contract DE-AC02-06CH11357. Some of the computing for this research was performed on the Blues and Bebop clusters at the Laboratory Computing Resource Center at Argonne National Laboratory. Part of the computing was performed on the Edison and Cori CRAY computers at NERSC through an ERCAP allocation. Some of the computing for this research was provided at the Stampede and Stampede 2 clusters at TACC, the Bridges cluster at PSC and the Comet cluster at SDSC, through an allocation provided through XSEDE. DKS would like to thank the CSSM and CoEPP at the University of Adelaide for their hospitality while some of this research was being performed.