Complex Langevin and other approaches to the sign problem in quantum many-body physics

TL;DR

This paper reviews complex Langevin and other methods addressing the sign problem in quantum many-body physics, discussing theoretical foundations, various approaches, and recent applications to quantum field theory.

Contribution

It provides a comprehensive overview of complex stochastic quantization and compares it with alternative methods for solving the sign problem in quantum many-body systems.

Findings

Complex Langevin can be applied to quantum field theories with sign problems.

Various methods like Lefschetz thimbles and density-of-states improve convergence.

The paper highlights open challenges and practical solutions in the field.

Abstract

We review the theory and applications of complex stochastic quantization to the quantum many-body problem. Along the way, we present a brief overview of a number of ideas that either ameliorate or in some cases altogether solve the sign problem, including the classic reweighting method, alternative Hubbard-Stratonovich transformations, dual variables (for bosons and fermions), Majorana fermions, density-of-states methods, imaginary asymmetry approaches, and Lefschetz thimbles. We discuss some aspects of the mathematical underpinnings of conventional stochastic quantization, provide a few pedagogical examples, and summarize open challenges and practical solutions for the complex case. Finally, we review the recent applications of complex Langevin to quantum field theory in relativistic and nonrelativistic quantum matter, with an emphasis on the nonrelativistic case.

Figures

Tables

Table 1. Table 1 : Properties and expressions for the Langevin method with real-valued fields versus complex-valued fields.

	Real Langevin	Complex Langevin
Fields	$\varphi (\overrightarrow{x},\tau )$	${\varphi }_R(\overrightarrow{x},\tau )+i{\varphi }_I(\overrightarrow{x},\tau )$
\hdashline Action	$S\left[\varphi \right]$	$S\left[{\varphi }_R+i{\varphi }_I\right]$
\hdashline Field update	$\Delta \varphi =K\left[\varphi \right]\Delta t+\eta \left(t\right)$
\hdashline Drift	$K\left[\varphi \right]=-\frac{\delta S\left[\varphi \right]}{\delta \varphi (\overrightarrow{x},\tau )}$
\hdashline Noise	$⟨\eta ⟩$ = 0	$⟨{\eta }_R⟩=⟨{\eta }_I⟩=0$
	$⟨\eta \left(t\right)\eta \left(t^{\prime }\right)⟩=2\Delta \delta \left(t-t^{\prime }\right)$

Equations

$$\mathcal{Z}=\sum _{C}W(C).$$

$$\mathcal{Z}=\text{Tr}\left[{\mathrm{e}}^{-\beta (\hat{H}-\mu \hat{N})}\right],$$

$$\mathcal{Z}=\int \mathcal{D}\varphi {\mathrm{e}}^{-S[\varphi ]},$$

$$S[\varphi ]=\int \mathrm{d}\tau {\mathrm{d}}^{d}x\left\{{\varphi }^{\ast }\left({\partial }_{\tau }+\mathcal{H}\right)\varphi +V[\varphi ]\right\},$$

$${\mathrm{e}}^{-S_{\mathrm{int}}[\varphi ]}=\int \mathcal{D}\sigma {\mathrm{e}}^{-S_{\mathrm{aux}}[\varphi ,\sigma ]-S_0[\sigma ]},$$

$$S_{\mathrm{int}}[\varphi ]=\int \mathrm{d}\tau {\mathrm{d}}^{d}xV[\varphi ],$$

$$\mathcal{Z}=\int \mathcal{D}\sigma {\mathrm{e}}^{-S_{\mathrm{HS}}[\sigma ]},$$

$${\mathrm{e}}^{-S_{\mathrm{HS}}[\sigma ]}=\frac{{\mathrm{e}}^{-S_0[\sigma ]}}{\det M[\sigma ]}.$$

$$\det M[\sigma ]=\det M_1[\sigma ]\det M_2[\sigma ]\cdots \det M_{N_f}[\sigma ].$$

$${\mathrm{e}}^{-S[\varphi ]}=\det M[\varphi ]{\mathrm{e}}^{-S_g[\varphi ]},$$

$$\det M[\varphi ]=\det M_1[\varphi ]\det M_2[\varphi ]\cdots \det M_{N_f}[\varphi ].$$

$$⟨\hat{N}⟩=\frac{\partial \ln \mathcal{Z}}{\partial (\beta \mu )}=\frac{1}{\mathcal{Z}}\text{Tr}\left[\hat{N}{\mathrm{e}}^{-\beta (\hat{H}-\mu \hat{N})}\right],$$

$$⟨\hat{N}⟩=\int \mathcal{D}\varphi P[\varphi ]\left[-\frac{\partial S[\varphi ]}{\partial (\beta \mu )}\right],$$

$$⟨\hat{N}⟩=\int \mathcal{D}\varphi P[\varphi ]\mathrm{Tr}\left[M^{-1}\frac{\partial M}{\partial (\beta \mu )}\right],$$

$$P[\varphi ]=\frac{{\mathrm{e}}^{-S[\varphi ]}}{\mathcal{Z}}.$$

$$⟨\mathcal{O}⟩=\int \mathcal{D}\varphi P[\varphi ]\mathcal{O}[\varphi ]\simeq \frac{1}{{\mathcal{N}}_{\varphi }}\sum _{\{\varphi \}}\mathcal{O}[\varphi ],$$

$$⟨\mathcal{O}⟩=\frac{1}{\mathcal{Z}}\int \mathcal{D}\varphi |P[\varphi ]|{\mathrm{e}}^{i\theta [\varphi ]}\mathcal{O}[\varphi ]=\frac{\int \mathcal{D}\varphi |P[\varphi ]|{\mathrm{e}}^{i\theta [\varphi ]}\mathcal{O}[\varphi ]}{\int \mathcal{D}\varphi |P[\varphi ]|}\frac{\int \mathcal{D}\varphi |P[\varphi ]|}{\int \mathcal{D}\varphi |P[\varphi ]|{\mathrm{e}}^{i\theta [\varphi ]}}=\frac{⟨⟨\mathcal{O}[\varphi ]{\mathrm{e}}^{i\theta [\varphi ]}⟩⟩}{⟨⟨{\mathrm{e}}^{i\theta [\varphi ]}⟩⟩},$$

$$⟨⟨{\mathrm{e}}^{i\theta [\varphi ]}⟩⟩=\frac{\int \mathcal{D}\varphi |P[\varphi ]|{\mathrm{e}}^{i\theta [\varphi ]}}{\int \mathcal{D}\varphi |P[\varphi ]|}=\frac{\mathcal{Z}}{{\mathcal{Z}}_{pq}}={\mathrm{e}}^{-\beta (\Omega -{\Omega }_{pq})},$$

$$⟨⟨{\mathrm{e}}^{i\theta [\varphi ]}⟩⟩={\mathrm{e}}^{-\beta V(\omega -{\omega }_{pq})},$$

$$\frac{\Delta }{⟨⟨{\mathrm{e}}^{i\theta [\varphi ]}⟩⟩}\sim \frac{{\mathrm{e}}^{\beta V(\omega -{\omega }_{pq})}}{\sqrt{{\mathcal{N}}_s}}.$$

$$\hat{H}=-t\sum _{s=\uparrow ,\downarrow }\sum _{⟨\mathbf{i}\mathbf{j}⟩}{\hat{c}}_{s,\mathbf{i}}^{†}{\hat{c}}_{s,\mathbf{j}}+U\sum _{\mathbf{i}}({\hat{n}}_{\uparrow ,\mathbf{i}}-1/2)({\hat{n}}_{\downarrow ,\mathbf{i}}-1/2),$$

$${\mathrm{e}}^{-\tau \hat{H}}\simeq {\mathrm{e}}^{-\tau \hat{T}}{\mathrm{e}}^{-\tau \hat{V}},$$

$${\mathrm{e}}^{-\tau U({\hat{n}}_{\uparrow }-1/2)({\hat{n}}_{\downarrow }-1/2)}=\frac{{\mathrm{e}}^{-\tau U/4}}{2}\sum _{\varphi =\pm 1}{\mathrm{e}}^{-{\lambda }_d\varphi ({\hat{n}}_{\uparrow }-{\hat{n}}_{\downarrow })},$$

$${\mathrm{e}}^{-\tau U({\hat{n}}_{\uparrow }-1/2)({\hat{n}}_{\downarrow }-1/2)}=\frac{{\mathrm{e}}^{-\tau U/4}}{\sqrt{2\pi }}{\int }_{-\infty }^{\infty }\mathrm{d}\varphi {\mathrm{e}}^{-{\varphi }^2/2-{\lambda }_c\varphi ({\hat{n}}_{\uparrow }-{\hat{n}}_{\downarrow })},$$

$$P[\varphi ]=\det M_{\uparrow }[\varphi ]\det M_{\downarrow }[\varphi ],$$

$${\mathrm{e}}^{-\tau U({\hat{n}}_{\uparrow }{\hat{n}}_{\downarrow }-{\hat{n}}_{\uparrow }/2-{\hat{n}}_{\downarrow }/2)}={\int }_{-\infty }^{\infty }\mathrm{d}\varphi p[\varphi ]{\mathrm{e}}^{\varphi ({\hat{n}}_{\uparrow }+{\hat{n}}_{\downarrow })},$$

$${\int }_{-\infty }^{\infty }\mathrm{d}\varphi p[\varphi ]$$

$${\int }_{-\infty }^{\infty }\mathrm{d}\varphi p[\varphi ]{\mathrm{e}}^{\varphi }$$

$${\int }_{-\infty }^{\infty }\mathrm{d}\varphi p[\varphi ]({\mathrm{e}}^{\varphi }{)}^2$$

$$\mathcal{Z}=\int \mathcal{D}\varphi {\mathrm{e}}^{-S[\varphi ]}=\int \mathrm{d}s\rho (s){\mathrm{e}}^{-is}=2\int \mathrm{d}s\rho (s)\cos (s),$$

$$\rho (s)=\int \mathcal{D}\varphi \delta (S_I[\varphi ]-s){\mathrm{e}}^{-S_R[\varphi ]}.$$

Full text

††journal:  

1. 1 Introduction

2. 1.1 The challenge of many-body quantum mechanics: memory and statistics

3. 1.2 Path integrals and the sign problem

4. 2 Approaches to the sign problem: from reweighting to complex Langevin

5. 2.1 Reweighting

6. 2.2 Alternative Hubbard-Stratonovich transformations

7. 2.3 Density-of-states methods

8. 2.4 Dual variables for bosons

9. 2.5 Dual variables for fermions and fermion bags

10. 2.6 Majorana fermions

11. 2.7 Imaginary asymmetry

12. 2.8 Lefschetz thimbles

13. 2.9 Path optimization method

14. 3 The Langevin method for real and complex variables

15. 3.1 Complex Langevin: origins and modern re-emergence

16. 3.2 Stochastic quantization: path integrals and the Langevin process

17. 3.3 A practical guide to real Langevin

18. 3.3.1 Toy problem I: a pedagogical example of real Langevin

19. 3.4 A practical guide to complex Langevin

20. 3.4.1 Toy problem II: a pedagogical example of complex Langevin

21. 3.5 Formal aspects and justification

22. 3.5.1 Practical monitoring of the boundary terms

23. 3.6 Challenges and some solutions

24. 3.6.1 Mathematical shortcomings

25. 3.6.2 Instabilities

26. 3.6.3 Convergence to incorrect values

27. 3.6.4 Practical solutions

28. 3.6.5 Schwinger-Dyson representation

29. 4 Applications in relativistic physics: from toy models to finite density QCD

30. 4.1 QCD-inspired toy models

31. 4.1.1 XY model

32. 4.1.2 Polyakov and SU(3) spin models

33. 4.1.3 Random matrix theory (RMT)

34. 4.1.4 Thirring model

35. 4.1.5 Heavy-dense QCD (HDQCD) and hopping parameter expansion

36. 4.1.6 Low-dimensional QCD at non-zero chemical potential

37. 4.2 $3+1$ dimensional QCD at finite chemical potential

38. 4.3 Other applications in relativistic field-theories

39. 4.3.1 Real-time dynamics

40. 4.3.2 Theories with topological terms

41. 4.3.3 Matrix models for superstring theory

42. 5 Applications in non-relativistic matter: ultracold atomic gases

43. 5.1 Non-relativistic bosons

44. 5.2 Non-relativistic fermions in one dimension

45. 5.2.1 1D fermions in the ground state

46. 5.2.2 1D fermions at finite temperature

47. 5.3 Non-relativistic fermions in three dimensions: the polarized unitary Fermi gas

48. 6 Summary and Outlook

1 Introduction

1.1 The challenge of many-body quantum mechanics: memory and statistics

In the early days of quantum mechanics it was quickly discovered that the Schrödinger equation could be solved analytically for hydrogen and hydrogen-like atoms in a straightforward manner [1]. However, each new particle added to the problem came at a dauntingly steep price, leaving the vast majority of the periodic table of the elements unattainable due to the complexity of the equations and the accompanying high cost of computation. Indeed, the presence of more than two interacting particles yields equations that are analytically intractable. The quantum few-body problem thus appeared to be very difficult, and the chances of solving the quantum many-body problem seemed dire. At the heart of the problem is the fact that, while the influence of the massive atomic nucleus on the much lighter electrons can be approximated (as a static external field à la Born-Oppenheimer), addressing the Coulomb interaction among the electrons is far more challenging. Dirac famously remarked in 1929 that, while the underlying physical laws were then completely known, “the difficulty is only that the exact application of these laws leads to equations much too complicated to be soluble”[2]. This difficulty could hardly be overemphasized then, and remains a challenge to this day. In facing that challenge, a wide variety of algorithms was – and continues to be – developed by specialists around the world to fit the paradigms of their specific area of physics or chemistry.

The most common first-principles approaches to the quantum many-body problem can be roughly divided into two sets: memory intensive and statistics intensive. The former include methods such as exact diagonalization (see e.g. [3, 4]) and coupled cluster (see e.g. [5, 6]), while the latter include a set of stochastic techniques generally known as quantum Monte Carlo (QMC) methods. Within that QMC set, this review focuses on a large class of approaches for which the many-body problem is expressed in the language of second-quantization or quantum field theory, such that expectation values of operators are written as a path integral over continuous fields living on a spacetime lattice. That formulation is in fact very general – it is natural in relativistic quantum field theory as well as nuclear and condensed matter physics, either in the form of low-energy effective field theories (see e.g. [7, 8, 9]), or as a reformulation of traditional Hamiltonians like the Hubbard model (see e.g. [10, 11]). Regardless of the application, the computational cost of path-integral QMC methods scales at face value (see below) polynomially with particle number and basis size (i.e. the size of the spacetime lattice), which makes them exceptionally well-suited for the many-body problem.

An essential component of QMC techniques is that they rely on a stochastic process governed by the Metropolis accept-reject algorithm [12], which itself requires a well-defined probability measure to guarantee convergence to the correct result. Simply put, the algorithm requires that the partition function $\mathcal{Z}$ be written as a sum of positive weights $W(C)$ (which play the role of the probability mentioned above) over some set of configurations $C$:

[TABLE]

The Metropolis algorithm, by construction, provides samples of the configurations $C$ distributed according to $W(C)$. Under many circumstances, however, a serious issue arises for this kind of algorithm, which has hindered computation in a wide range of situations: the infamous sign problem. In those cases, $W(C)$ does not have a well-defined sign or even becomes complex (as explained in further detail below). Unfortunately, by far most systems of interest suffer from such a problem: high-$T_{c}$ superconductors (due to strong repulsive interaction away from half filling, see e.g. [13]), nuclear structure (strong repulsive core, finite spin-isospin polarization, see e.g. [14, 15]), and quantum chromodynamics (finite quark density see e.g. [16, 17, 18]), to name only a few.

Over the last few decades, many ideas have been proposed to overcome the sign problem in quantum many-body physics and field theory. This review covers some of them briefly and focuses on the so-called complex Langevin (CL) approach, as applied to the calculation of equilibrium properties of quantum many-body systems in relativistic and nonrelativistic physics, with an emphasis on the latter. The next section sketches out the basic path-integral formalism involved in Metropolis-based and stochastic quantization approaches, with the goal of showing where and how the sign problem arises.

1.2 Path integrals and the sign problem

The central quantity of the field theoretical approach to the quantum many-body problem is the partition function, which in the grand canonical ensemble is given by

[TABLE]

where $\hat{H}$ is the Hamiltonian of the system, $\hat{N}$ the particle number operator, $\beta$ the inverse temperature, $\mu$ the chemical potential, and the trace is over all multiparticle states (i.e. Fock space). Note that $\mathcal{Z}$ is shown here for a single particle species, but is straightforwardly generalized to multiple chemical potentials, etc. As written, $\mathcal{Z}$ contains the thermodynamic information of the system: by differentiating with respect to $\beta$ and $\mu$ one obtains expectation values of the Hamiltonian and the particle number operators. More detailed information, such as momentum distributions and other correlation functions, can be obtained by adding sources to the Hamiltonian as is common in quantum field theory. The direct evaluation of Eq. (2) is impossible for interacting systems, as it requires a priori knowledge of the full energy spectrum.

The path-integral approach to the many-body problem provides an alternative route. Either using an operator-based approach or coherent states, one arrives at an expression for $\mathcal{Z}$ which is written generically as

[TABLE]

where $\phi$ is a field living in $(d+1)$-dimensional spacetime and represents any degrees of freedom in the system. One thus replaces the problem of evaluating Eq. (2) with that of calculating the above path integral. In practice, boundary conditions in the spatial directions can be chosen in a variety of ways, but those in the time direction, which is compact and runs in the range $\tau\in[0,\beta)$, are set by the quantum statistics of the problem: bosonic fields will obey periodic boundary conditions and fermionic fields anti-periodic.

From this point on, our discussions will focus on nonrelativistic systems unless otherwise specified. In purely bosonic theories the action $S[\phi]$ will typically take a simple local form such as

[TABLE]

where $\mathcal{H}$ represents the noninteracting Hamiltonian (including external trapping potentials) and $V[\phi]$ represents the interactions (i.e. terms cubic and beyond in $\phi$).

For real bosonic variables $\phi$, the action $S[\phi]$ is also real and therefore ${\rm e}^{-S[\phi]}$ can be used as a probability measure in a stochastic process. For complex $\phi$, however, the fact that $\partial_{\tau}$ is an antisymmetric operator results in a complex $S[\phi]$, which is the source of a sign problem in this formulation (see below) and has a counterpart in relativistic bosons at finite chemical potential (see Secs. 2.4 and 4). The problem of the antisymmetry of $\partial_{\tau}$ can be circumvented for an even number of species with attractive interactions, which however render bosons unstable (but not fermions, see below).

To make contact with the fermionic case discussed below (see also Fig. 1), it is instructive to rewrite the interaction using a Hubbard-Stratonovich (HS) transformation, although it is not strictly needed in the bosonic case. Consider for instance the case of a complex field $\phi$. Schematically, one introduces an auxiliary field $\sigma$ such that

[TABLE]

where

[TABLE]

$S_{\rm aux}[\phi,\sigma]$ is a quadratic functional of both $\phi$ and $\sigma$, and $S_{0}[\sigma]$ is a pure-$\sigma$ term; both $S_{\rm aux}$ and $S_{0}$ depend on the specific choice of HS transformation. Since the action is now quadratic in $\phi$, the corresponding path integral can be carried out, which results in a $\sigma$-dependent determinant, i.e. the full partition function can now be expressed as

[TABLE]

where

[TABLE]

Usually it is possible to factor the determinant into a determinant for each particle species (flavors), such that if $N_{f}$ flavors are present then

[TABLE]

Naturally, calculations with bosons are not carried out using the action of Eq. (8). The formulation based on $S[\phi]$ of Eq. (4) is considerably easier to work with as there are no determinants involved. However, the appearance of the boson determinant in Eq. (8) shows that, if $\det M$ is real and an even number of species is present, then the sign problem can be avoided if $S_{0}[\sigma]$ is real. Unfortunately, that situation is not relevant for bosons as it corresponds to attractive interactions, which make a many-boson system unstable. On the other hand, the determinant-based representation of Eq. (8) has a direct counterpart in the fermion case, which we discuss next.

In theories with fermions, the action will require the much more complicated (non-linear, non-local) form based on determinants because the fermionic analogue of Eq. (4) is written in terms of anticommuting objects, i.e. Grassmann numbers, which are not amenable to numerical computation. We therefore assume that fermionic degrees of freedom (i.e. said Grassmann variables) have been integrated out. Taking such a step requires a HS transformation of some kind to decouple the interaction, i.e. one introduces auxiliary fields to obtain a quadratic action in the fermion fields, which are then integrated and result in a fermion determinant. Assuming such steps have already been taken, we have the schematic form,

[TABLE]

where $M$ encodes the dynamics of the fermions (quarks, electrons, atoms) in the external field $\phi$, and $S_{g}[\phi]$ is the “pure HS” part of the action (often called “pure gauge” part in QED and QCD); for the latter, the form of $S_{g}$ will depend on the kind of HS transformation utilized (see Sec. 2.2). Parameters like the fermion mass and chemical potential appear in $M$. In particular, in many cases it is possible to choose a HS transformation that decouples $N_{f}$ species (i.e. flavors) of fermions such that, as in the bosonic case described above,

[TABLE]

The above path integral formulation, being a rewriting of $\mathcal{Z}$, inherits the usual mechanisms to access expectation values of operators, namely differentiating $\mathcal{Z}$ with respect to a chosen parameter. For instance, the average particle number is given by

[TABLE]

such that in the bosonic case of Eq. (4),

[TABLE]

while in the fermionic case of Eq. (10),

[TABLE]

where in either case

[TABLE]

In evaluating Eqs. (13) or (14), the natural course of action is to sample field configurations $\phi$ according to the probability $P[\phi]$ and evaluate the quantities of interest that appear between square brackets. It is for that reason that the identification of $P[\phi]$ as a probability measure is a central aspect of conventional, Metropolis-based approaches to the evaluation of expectation values in quantum systems with many degrees of freedom. More specifically, in those cases where the sign (or phase) of $P[\phi]$ does not depend on $\phi$, one samples $\phi$ according to $P[\phi]$ using the Metropolis algorithm (combined with a suitable field updating procedure, e.g. Wolff, worm, or hybrid Monte Carlo algorithms) to obtain a set of $\mathcal{N}_{\phi}$ decorrelated samples $\{\phi\}$, which in turn are used to estimate expectation values as

[TABLE]

for a given operator $\mathcal{O}$.

As mentioned above, by far for most systems of interest in physics face a sign problem, as the sign (or more generally complex phase) of $P[\phi]$ varies with $\phi$. Then, $P[\phi]$ simply cannot be interpreted as a probability and the Metropolis algorithm is not applicable.

For nonrelativistic fermionic systems, the sign problem typically happens at finite polarization (i.e. chemical potential asymmetry) or when interactions contain a repulsive component. The problem therefore affects essentially all of condensed matter, nuclear physics, and quantum chemistry. There are notable exceptions, such as a large class of systems in one spatial dimension and the Hubbard model at half filling, for which the sign problem can be eliminated completely. For relativistic fermions, such as quarks at finite chemical potential, the sign problem has obstructed the investigation of the phase diagram of QCD.

The case of bosons is markedly different from that of fermions. Here, the nonrelativistic case presents a sign problem even in the absence of interactions or chemical potentials: it is the asymmetry of the single time derivative, see Eq. (4), that creates the problem, as we will explain in further detail in Sec. 3. This is to be contrasted with the relativistic case, which develops a sign problem when a chemical potential is turned on (see however Sec. 2.4).

The remainder of this review is organized as follows. Sec. 2 reviews a broad (but by no means complete) set of approaches to the sign problem. Sec. 3 introduces the formal aspects of stochastic quantization and the complex Langevin method in more detail, including pedagogical examples as well as a brief discussion of the challenges and shortcomings in the mathematical underpinnings. Secs. 4 and 5 review the recent and emerging applications of CL in relativistic and nonrelativistic physics, with an emphasis on the latter. Finally, Sec. 6 concludes the review with a summary and outlook.

2 Approaches to the sign problem: from reweighting to complex Langevin

There have been multiple approaches suggested to solve or ameliorate the sign problem. Some of these methods aim at solving the problem directly, typically by rewriting the partition function in new and clever ways that remove the sign problem entirely. Other approaches involve rewriting the original problem so it can be solved stochastically but with controlled sign fluctuations. Below we present a selection of those methods in a logical sequence that starts with the simplest idea, namely reweighting, and concludes with complex plane methods. Along the way, we present an elementary discussion of each method and comment on their advantages and shortcomings, which often result in valuable insights on the nature of the sign problem.

In Fig. 2 we propose a visual organization of the various approaches to the sign problem. Although we do not follow the proposed taxonomy in this review in a linear fashion, we do find it helpful to organize the information in this manner. On the left column of that figure we list “new variables” methods as those that attempt to tackle the sign problem by switching from the conventional path-integral formulation to a new set of variables. We begin by reviewing a classic work that looks for sign problem free HS transformations in Sec. 2.2. That work shows that a choice of HS transformation may not solve the problem but may help in addressing it (and certainly choosing the wrong one can be a recipe for trouble). Dual-variable and Majorana-fermion representations succeed in completely solving the sign problem in many cases, as shown in Secs. 2.4 and 2.6.

The middle column of Fig. 2 lists the set of what we call “statistical” approaches, which attempt to tackle the sign problem in a head-on manner. The simplest of those methods is by far the oldest and most commonly applied one across all areas of physics: reweighting, which we describe in Sec. 2.1. More recent statistical approaches, commonly referred to as “density-of-states” methods, proceed by probing the shape of the probability distribution at fixed action phase, and then integrating over that phase at the end; we describe those ideas in Sec. 2.3

Finally, the right column lists “complex plane” methods, which also come in different flavors. Sec. 2.7 reviews the use of imaginary asymmetries in the parameters of a given theory (e.g. chemical potential, mass imbalance) to carry out calculations without a sign problem, necessarily followed by some kind of analytic continuation to return to the real physical values. Complex Langevin methods, the focus of this review, start in Sec. 3.1, while in Sec. 2.8 and Sec. 2.9 we mention ideas based on contour deformations (Lefschetz thimbles and the path optimization method, respectively). All of those methods rely on complexifying the integration variables; based on that idea, there exist constructive approaches (mentioned in Sec. 3.1) that aim to define a real action in such a complex space.

2.1 Reweighting

The simplest (and likely oldest [19]) idea to overcome the sign problem is that of reweighting, which amounts to sampling $\phi$ using the magnitude of $P[\phi]$ as a probability measure. In such an approach, one rewrites the expectation value of $\mathcal{O}$ as

[TABLE]

where ${\rm e}^{i\theta[\phi]}$ is the phase of $P[\phi]$, and the double angle bracket denotes an expectation value taken with respect to $|P[\phi]|$.

Reweighting thus provides a way forward for systems that have a sign problem: simply compute the numerator and denominator of Eq. (17)) and then take their ratio. In practice, however, both the numerator and the denominator of Eq. (17)) vanish exponentially as the physical extent of the spacetime lattice is increased. The phase average is

[TABLE]

where ${\mathcal{Z}}_{pq}$ is the partition function of the “phase-quenched” theory111In the phase-quenched theory, the probability distribution is replaced by its absolute value, such that it is a positive-definite measure., and $\Omega_{pq}$ and $\Omega$ are the corresponding grand thermodynamic potentials of the phase-quenched and original theory. Both $\mathcal{Z}$ and ${{\mathcal{Z}}_{pq}}$ are real quantities, but since ${\mathcal{Z}}_{pq}$ is a sum over nonnegative real numbers, while ${\mathcal{Z}}$ accounts for the phase, we necessarily have ${\mathcal{Z}}\leq{{\mathcal{Z}}_{pq}}$. More importantly, thermodynamic potentials are extensive quantities in the spatial volume $V$ of the system, such that we may write the above in terms of intensive potentials $\omega$ and $\omega_{pq}$ as

[TABLE]

which exposes the exponential nature of the sign problem in the thermodynamic ($V\to\infty$) and ground-state ($\beta\to\infty$) limits. This can be seen more clearly by examining the statistical uncertainty $\Delta$, which in a Monte Carlo calculation with ${\mathcal{N}}_{\mathrm{s}}$ samples decreases as ${\mathcal{N}}_{\mathrm{s}}^{-1/2}$.In the case of re-weighting, the relative statistical uncertainty on the average phase is overpowered by the exponential behavior coming from Eq. (19):

[TABLE]

This last equation shows the difficulty in approaching the sign problem with a simple technique such as re-weighting: an exponentially large number of samples is needed in order to determine the average phase with any reasonable accuracy as the volume of spacetime is increased. Viewed through the lens of this simple idea, the sign problem may be regarded as the reappearance of an exponential type of computational wall, which affects non-stochastic methods (see Introduction) in the guise of memory requirements and statistical methods in the form of a signal-to-noise problem.

2.2 Alternative Hubbard-Stratonovich transformations

The partition function $\mathcal{Z}$ of Eq. (2) is naturally a sum of positive quantities ${\rm e}^{-\beta(E-\mu N)}$. The path-integral representation of $\mathcal{Z}$, while exact, introduces a large number of degrees of freedom to represent the same quantity. It therefore seems natural to expect that such a formulation would require massive cancellations (i.e. the sign problem) to yield correct physical answers. On the other hand, there are many ways to choose a HS representation, which may in turn yield different kinds of cancellations (i.e. more or less dramatic, by some measure). Even in the absence of a sign problem, different kinds of continuous or discrete HS transformations display varying behavior (see e.g. [20, 21]). It therefore makes sense to ask whether efficient representations exist, i.e. HS transformations which can substantially reduce the difference $\omega-\omega_{pq}$ in Eq. (19) or even eliminate it completely.

As an example, consider the fermionic Hubbard model given by

[TABLE]

where $t$ is the nearest-neighbor hopping, $U>0$ is the repulsive coupling, $\hat{c}^{(\dagger)}_{s,\bf i}$ is the annihilation (creation) operator for particles of spin $s=\uparrow,\downarrow$ at location ${\bf i}$, and $\hat{n}_{s,{\bf i}}$ is the corresponding density operator.

The work of Ref. [22] showed that a general HS transformation for the above model resulting in positive weights (i.e. a real action $S[\phi]$) does not exist. While attractive interactions – such as in the negative-$U$ Hubbard model – feature no sign problem, repulsive interactions and in general any finite polarization (i.e. non-zero chemical potential asymmetry) do yield sign oscillations. More specifically, the problem arises because the determinant in Eq. (10) becomes a product of two determinants which are real or can be made real by choosing a proper HS transformation, but which will generally have different signs. Reference [22] showed that it is possible to isolate the origin of the signs in such a way that the determinants are real and identical, i.e. one ends up with a square of a determinant, and the signs are not eliminated but can be predicted. This remarkable property is illustrated below.

We begin by implementing a Trotter-Suzuki factorization of the Boltzmann weight with imaginary time step $\tau$, such as

[TABLE]

where $\hat{T}$ contains the hopping terms and $\hat{V}$ the on-site interaction, as they appear in Eq. (21). It is to address the latter that an HS transformation is used. The two most common HS representations, used in calculations of the repulsive Hubbard model, proceed by writing (omitting the spatial indices)

[TABLE]

[TABLE]

where $\phi$ is the auxiliary field, $\lambda_{d}$ is set by $\cosh(\lambda_{d})={\rm e}^{\tau U/2}$, and $\lambda_{c}=\sqrt{\tau U}$. Both of these “density-channel” transformations successfully decouple the two spin species $\uparrow$ and $\downarrow$, and the resulting determinants are real, but they are generally different from each other, such that

[TABLE]

will generally vary in sign with $\phi$. Here we omit the pure-$\phi$ part for the continuous case, which is real and positive anyway. It should be pointed out that there are more general ways than the above factorized form that result in a sign problem free situation; for an exploration of more general conditions based on time-reversal invariance, see Refs. [23, 24] and further discussion in Sec. 2.6.

A more general transformation that aims to preserve the up-down symmetry of the Hubbard Hamiltonian, and therefore provide the square of a real determinant in $P[\phi]$, can be written as

[TABLE]

where we want $p[\phi]$ to be real and positive and, evaluating both sides at the eigenvalues of the density operators (i.e. setting $\hat{n}_{s}\to 0,1$), we see that

[TABLE]

Unfortunately, the last two equations can only be satisfied simultaneously if ${\rm e}^{\tau U/2}\leq 1$, i.e. if $U\leq 0$, which is not the case we are interested in here as there is then no sign problem.

The above shows that, at least within the rather general form proposed, it is not possible to generate the square of a determinant and avoid the sign problem at the same time for the repulsive Hubbard model. On the other hand, if $p[\phi]$ is allowed to vary in sign, then there is no constraint on $U$ and we obtain the square of a real determinant. In that case, the sign problem comes not from the fermion determinant but from $p[\phi]$, which means that it is completely predictable as soon as $\phi$ is known, without computing determinants. Such predictability of the sign or phase of the determinant has not been exploited in the literature beyond the work of Ref. [22], but it could be of interest in the context of the density-of-states methods discussed in the next section. In those methods, the knowledge of the precise form of the imaginary part of the action as a functional of the field is essential and has been used with some success to characterize a class of relativistic field theories.

Generally speaking, the fact that there exists a family of HS transformations representing the same partition function, especially if they are non-trivially related to each other, provides in effect a variety of calculations that can be used as checks against each other. As explored in Ref. [25], the density-channel decompositions mentioned above can be replaced by their “pairing channel” counterparts (of which a new family exists, with discrete and continuous members, as above), which display different sign properties. Other kinds of useful HS transformations have been discussed in Refs. [26, 27]

Crucially, the availability of different HS transformations with different sign behavior shows that the sign problem is not an intrinsic property of a given Hamiltonian, but rather depends on the decoupling scheme. Therefore, the search for a link between the physics of a given system and the sign problem should be taken with caution, as such a link may be entirely an artifact of the formulation of the problem. An interesting example in that regard is the elimination of the sign problem by way of a fermionic reformulation of a bosonic problem in the case of a frustrated Kondo model coupled to fermions [28], followed by a HS transformation on the resulting fermionic interaction (see also [29, 30]).

It is worth pointing out, however, that there is a link between the sign problem and phase transitions. Indeed, with the path integral formulation at hand, one can reasonably argue that the sign problem can be expected to be severe close to a critical point. One way to visualize that concept is in terms of the Lee-Yang zeros of the partition function $\mathcal{Z}$, written as a path integral (i.e. Eq. (7). When sampled over the relevant configurations of $\phi$, the integrand ${\rm e}^{-S[\phi]}$ must reflect the existence of an accumulation point of roots of $\mathcal{Z}$ when approaching the phase transition. By itself, that property would not pose a problem. However, the natural scale of the integrand is the exponential of an extensive quantity; therefore, ${\rm e}^{-S[\phi]}$ must oscillate dramatically in order to generate the large collection of zeros (in the thermodynamic limit), and the corresponding high sensitivity to the parameter values, around the transition point. For cases that do not have a sign problem, the integrand must necessarily tend to zero when approaching a phase transition (again, in the thermodynamic limit), which is often reflected in the appearance of zero modes in fermion matrices.

2.3 Density-of-states methods

The density-of-states (DoS) approaches are a class of methods that attempt to tackle the sign problem in a head-on manner, as opposed to rewriting the partition function in terms of new variables or straightforward reweighting (although it may be argued that DoS methods are actually a kind of reweighting). The original idea of sampling the density of states as an alternative to Metropolis-based methods is due to Wang and Landau [31] and has been applied to a wide variety of systems including gauge theories [32, 33], but its generalization to systems with a sign problem was explored later on in Refs. [34, 35, 36, 37] (see also Refs. [38, 39]). The result of those explorations is now known in the literature as the logarithmic linear regression (LLR) algorithm or the functional fit approach (FFA), both of which are very closely related but differ on specific details. We will restrict ourselves here to those approaches (which have also been reviewed recently in Ref. [16]) but it is worth pointing out that DoS methods have also been applied to finite density QCD in different forms which involve histograms of the phase of the fermion determinant (see e.g. [40, 41, 42]).

The idea common to all DoS approaches is that, in the presence of a sign problem where the action can be decomposed into real and imaginary parts $S=S_{R}+iS_{I}$, the partition function can be written as

[TABLE]

where we used the fact that the partition function is real and took the real part in the last step, and

[TABLE]

The determination of $\rho(s)$ is then carried out by combining two ingredients: first, propose a functional form that can account for its variation over vast orders of magnitudes; second, carry out restricted calculations at constant or approximately constant imaginary action $S_{I}[\phi]$ in order to determine the coefficients in the proposed functional form for $\rho(s)$. A key aspect of the method is that, using these elements, it can deliver exponential accuracy in the calculation of $\rho(s)$.

In the approach of Refs. [35, 43], the parametrization of $\rho(s)$ is done in a piecewise-linear fashion:

[TABLE]

for $s\in I_{n}$, $I_{n}=[s_{n},s_{n+1}]$, where the partitioning and sizes of the intervals $I_{n}$, i.e. the set of numbers $\{s_{j}\}$, can be chosen at will to reflect the desired precision in describing $\rho(s)$. By requiring continuity of $\rho(s)$ and a normalization condition $\rho(0)=1$, the constants $A_{n}$ can be determined as a function of $k_{n}$:

[TABLE]

where $\Delta_{j}=s_{j+1}-s_{j}$ is the size of the $j$-th interval. In order to determine the constants $k_{n}$, the FFA uses restricted expectation values defined by

[TABLE]

where the restricted partition function is

[TABLE]

with $\theta_{n}(x)=1$ for $x\in I_{n}$ and [math] otherwise.

With the above piecewise-linear parametrization of $\rho(s)$, the restricted partition function and expectation values can be computed in closed form. It turns out that

[TABLE]

where $D_{n-1}=\sum_{j=0}^{n-1}\Delta_{j}$ and

[TABLE]

These equations encode a crucial aspect of the method: once the intervals $I_{n}$ are chosen (such that the $\Delta_{n}$ are fixed constants), each of the functions $Y_{n}(\lambda)$ is entirely determined by the single parameter $k_{n}$ and must follow the shape dictated by $h(x)$. Thus, the function $Y_{n}(\lambda)$ is a kind of response function in which the source parameter $\lambda$ is coupled to the imaginary part of the action $S_{I}$ to constrain the value of $k_{n}$ for each $n$. If the one-parameter fit to $h(x)$ is unsatisfactory, that signals a poor choice of the discretization $\{s_{j}\}$, such that a more refined mesh is likely needed. An example of the typical shape of $Y_{n}(\lambda)$ is shown in the left panel of Fig. 3 for several values of $n$ for the SU(3) spin model of Ref. [44].

Once the $k_{n}$ are known, one may reconstruct $\rho(s)$ and calculate the full partition function as a Fourier transformation via Eq. (30)). This last step will be sensitive to the large oscillations due to the $\cos(s)$ factor. As an example, the density of states obtained in Ref. [44] for the SU(3) spin model is shown in the right panel of Fig. 3, for several values of the nearest-neighbor coupling $\tau$. The variations of $\rho(s)$ over many orders of magnitude are evident from that figure.

While we have focused here on the FFA, the LLR method [33] accomplishes exponential error suppression by calculating the slopes $k_{n}$ of the distribution using a fixed-point iteration method. The latter is based on the work of Robbins and Monro [45], which guarantees that the values obtained in subsequent iterations are Gaussian-distributed around the exact answer. The above-mentioned exponential error suppression amounts to a constant relative error in the determination of the density of states over the full domain of the phase, which is crucial in order to carry out the Fourier integral in Eq. (30)).

The FFA and the LLR methods have been used to analyze several models on the lattice at finite density such as the $Z_{3}$ spin model [43], the SU(3) gauge theory with static color sources [44, 46], and two-color QCD with heavy quarks [37]. One of the most interesting advantages of DoS methods is that they are extremely parallelizable on modern computers. These methods require a large set of calculations, e.g. several independent calculations as a function of $\lambda$ and $n$ to determine $k_{n}$, but each of those is an independent run and can therefore be done in a perfectly scalable fashion. As long as the number of $\lambda$ and $n$ points does not grow exponentially with the spacetime volume (there is not indication thus far in the literature that that is the case, due largely to the smoothness and lack of sharp features in $\ln\rho(s)$), the computational cost will scale better than that of reweighting methods.

On the other hand, DoS approaches present a paradigm that is quite different from conventional MC methods when it comes to calculating different observables. For observables $\mathcal{O}$ that do not depend explicitly (and only) on the action $S[\phi]$, it will not be enough to know $\rho(s)$. In such a case, a source term $j\mathcal{O}$ would have to be included in the action and a family of densities $\rho(s,j)$ would need to be calculated at least for small $j$. Once $\rho(s,j)$ is thus obtained, numerical differentiation of $\ln\mathcal{Z}[j]$ yields the desired $\langle O\rangle$ in the limit $j\to 0$. Such an approach may seem slightly cumbersome or costly, but it amounts to multiple applications of the same idea which moreover retains the full parallelizability property mentioned above.

2.4 Dual variables for bosons

Dualization is another approach involving rewriting the partition function $\mathcal{Z}$ so as to eliminate or ameliorate the sign problem. The dual variables are a new set of variables, typically discrete, that may yield a representation of $\mathcal{Z}$ entirely in terms of positive quantities. While the concept of duality is in itself an old one, the use of dual variables in quantum Monte Carlo calculations first appeared in the 1980’s: for instance, Ref. [47] showed that that the strong-coupling limit of QCD could be represented as a system of dimers. Later on, Ref. [48] used dual variables to analyze the bosonic Hubbard model. The concept was later applied to fermions as well, which we review in the next section. For pure bosonic theories at finite density, it was shown in Ref. [49] that dual variables are not just an alternative representation: they successfully solve the sign problem for both relativistic as well as non-relativistic systems. Moreover, the dual variables can be efficiently sampled using the worm algorithm [50, 51]. Since the early 2010’s, a few groups have pursued the study of several quantum field theories (from simple models to effective theories of QCD) at finite temperature and density (see e.g. [52, 53, 54, 55, 56, 57, 58, 59]).

Following the notation and steps of Ref. [16], we show how to introduce dual variables first in relativistic bosons and then in the non-relativistic case. The lattice action for a relativistic complex-valuedued field $\phi_{x}$ is

[TABLE]

where $x$ is a spacetime lattice point and $\hat{\nu}$ denotes a unit vector in the $\nu$-th direction ($\nu=4$ being the imaginary-time direction). At finite $\mu$, the quantity in square brackets ceases to be real, which exposes the sign problem. Exponentiating the action, as one would normally do to compute the partition function, we note that

[TABLE]

At each spacetime-Lorentz point $x,\nu$, the offending term becomes a factor that can be rewritten by expanding each of the exponentials in a Taylor series as

[TABLE]

where $\mathcal{N}_{n,\bar{n}}=\prod_{x,\nu}{1}/({n_{x,\nu}!\bar{n}_{x,\nu}!})$ and the sum $\sum_{\{n,\bar{n}\}}$ denotes a sum over all configurations of the Taylor indices $n_{x,\nu}\geq 0$ and $\bar{n}_{x,\nu}\geq 0$. Using the above and the polar form $\phi_{x}=r_{x}{\rm e}^{i\theta_{x}}$, we obtain for the partition function

[TABLE]

where

[TABLE]

which is a non-negative, local function of the index configuration $\{n,\bar{n}\}$ (note also that the exponent of $r_{x}$ is strictly positive), and

[TABLE]

which results in Kronecker delta functions for each $x$ imposing constraints on the configurations. In principle, the job is done at this point: we have shown that there is a discrete-field representation of the partition function as a sum over positive quantities. It is useful, however, to take a few more steps towards simplifying the calculation, specifically towards implementing the constraints imposed by the function $T[n,\bar{n},x]$. To that end, one parameterizes the sum and difference of $n$ and $\bar{n}$ via two new ‘dual’ variables $k,\ell$ defined via

[TABLE]

which take values over all integers and all non-negative integers, respectively. We finally obtain

[TABLE]

where

[TABLE]

[TABLE]

The constraint $\nabla_{v}k_{x,\nu}=0$ enforces that the field $k_{x,\nu}$ be solenoidal, such that the flux of $k_{x,\nu}$ is conserved.

Pure gauge theories can also be written in terms of dual variables and, as above, results in a new representation for the partition function in terms of purely positive terms. This approach yielded the first real and positive dualization of abelian gauge theories with a so-called $\theta$ term [60], which is a term in the action coupled to a topological charge (see Refs. [61, 62] for updates on that work). Because that term is necessarily complex, it results in a sign problem in conventional path integral representations. The current challenge for this approach is the extension to non-abelian gauge theories and the inclusion of fermions (see however next section).

The case of non-relativistic bosons can also be addressed with dual variables with some modifications with respect to the relativistic case. We show some of steps of that derivation here as a pedagogical example; they closely follow the relativistic case. As we will show in more detail in Sec. 5, the problem appears in non-relativistic bosons not because of the chemical potential but because of the asymmetry in the time derivative: there are only particles and no antiparticles. (The physical source of the problem is thus the same as in the relativistic case: the breaking of time-reversal invariance.) Starting with the lattice action for the complex-valued field $\phi_{x}$ in $3+1$ dimensions (although this example can be easily generalized to $d+1$ dimensions), we have

[TABLE]

Then,

[TABLE]

where we now expand the exponentials of the derivative terms in a power series, such that

[TABLE]

where $\mathcal{N}_{n,m,\bar{m}}=\prod_{x}{1}/{n_{x}}!\prod_{x,k}{1}/(2^{m_{x,k}+\bar{m}_{x,k}}{m_{x,k}!\bar{m}_{x,k}!})$, and $n_{x}$ is a site variable, whereas $m_{x,k}$ and $\bar{m}_{x,k}$ are link variables in the spatial directions. As in the previous example, one may now write the fields in terms of their polar representation $\phi_{x}=r_{x}{\rm e}^{i\theta_{x}}$ to obtain constraints for the integer fields ${n,m,\bar{m}}$ and eventually arrive at a sum of positive definite terms for $\mathcal{Z}$. Explicitly,

[TABLE]

where

[TABLE]

which is a non-negative, local function of the index configuration $\{n,m,\bar{m}\}$ (also, as before, the exponent of $r_{x}$ is strictly positive), and

[TABLE]

which results in Kronecker delta functions that impose constraints on the configurations.

It is then clear that the dual-variable formulation avoids the sign problem for non-relativistic bosons, as first noted in Ref. [49]. It is worth noting, however, that other cases such as coupling to angular momentum, are not obviously solvable with this technique.

While this method completely solves the sign problem in the cases shown above (and some others, e.g. bosons with non-abelian spin-orbit coupling), the calculation of specific observables acquires a new degree of complexity due to the dramatic change of variables. This is merely an algebraic inconvenience but, in practice, the change from the original fields $\phi$ to the discrete fields $n,m,\bar{m}$ implies that any operator expression in the $\phi$ language needs to be re-derived (e.g. by inserting sources in the original action or using the parameters of the theory).

In this section we have focused on exact alternative representations of scalar bosons. However, such representations have also been found for gauge fields [63, 64] which are valid also away from the strong coupling limit.

2.5 Dual variables for fermions and fermion bags

One of the first uses of dual variables for fermions in Monte Carlo calculations was unrelated to the sign problem: it was found in Ref. [47] that QCD in the strong-coupling limit could be represented as a system of dimers, which inspired multiple studies [65, 66, 67, 68, 69, 70]. However, dual variables by themselves do not necessarily avoid the fermion sign problem. As first described in Ref. [71], in what is now known as the ‘meron cluster’ approach, one must sum analytically over configurations in a given cluster (where the type of configuration cluster must be cleverly identified) and then stochastically over clusters. When the clusters are properly chosen, they contain configurations that may vary in sign but such that the overall contribution of a cluster is of constant sign across clusters.

The fermion bag approach of Refs. [72, 73] (extended to continuous time in Ref. [74]; see also Refs. [75, 76, 77, 78]) extends the meron cluster approach to a larger class of theories. Rather than following those derivations (based on Grassmann numbers), here we connect with them from a different perspective. We begin by re-writing the partition function as

[TABLE]

where $\beta=\tau N_{\tau}$. Using a Suzuki-Trotter decomposition, we may approximate

[TABLE]

where $\hat{T}$ is the kinetic energy and $\hat{V}$ the interaction. As an example, we specialize to the case where $\hat{V}=-U\sum_{x}\hat{n}_{\uparrow}(x)\hat{n}_{\downarrow}(x)$, and then we have

[TABLE]

where $B={\rm e}^{\tau U}-1$. We may use the above at each point in imaginary time by inserting this expression into Eq. (55)) to obtain

[TABLE]

where $\hat{\mathcal{T}}_{t}\equiv{\rm e}^{-\tau(\hat{T}-\mu\hat{N})}$ and we note that the product over the time slices actually factorizes across flavors as

[TABLE]

Below we will show the following trace-determinant identity

[TABLE]

where $W_{s}[\{m\}]$ is the free fermion matrix for spin $s$ in which the rows and columns ${x_{0},t_{0}}$ for which $m_{x_{0},t_{0}}=1$ are dropped (see below for details on the form of $W_{s}[\{m\}]$). Using the above in the definition of $\mathcal{Z}$ yields

[TABLE]

where we finally have completely re-written the full partition function as a sum over configurations of the monomer field $m_{x,t}$.

For unpolarized non-relativistic systems, $W_{s}[\{m\}]$ is real and takes on the same value for $s=\uparrow,\downarrow$, such that there is no sign problem, as long as $B\geq 0$, i.e. for attractive interactions. Thus far, the same conditions apply for the auxiliary field formulation of the problem: repulsive interactions ($B<0$) or polarization ($W_{\uparrow}[\{m\}]\neq W_{\downarrow}[\{m\}]$) would lead to a sign problem. This formulation, however, lends itself to an interpretation of the sum in terms of clusters known as fermion bags, which are disjoint regions ${\mathcal{B}}_{i}$ of the discrete field $m_{x,t}$ within which $m_{x,t}=0$ (see Fig. 4). As the corresponding interaction vertices (i.e. insertions of $U$) are thus absent in $W_{s}[\{m\}]$, fermions are free to move about inside the bag. Using that property, Ref. [72] argued that, while the contributions to a given fermion bag configuration may vary in sign, the overall contribution of each bag to the full partition function is actually positive. Thus, if one is able to add up the terms within each bag, then it is possible to use Metropolis-based importance sampling to sum over all possible fermion bag configurations.

The right-hand side of Eq. (59)) can also be calculated using Wick’s theorem, if interpreted as the expectation value of a time-dependent operator in a noninteracting system. In that weak-coupling interpretation, it is possible to show that

[TABLE]

where $M_{\uparrow}$ is the noninteracting spacetime fermion matrix and $G_{\uparrow}[\{m\}]$ is a propagator matrix whose size depends on the monomer configuration $m_{x,t}$ and which contains noninteracting propagators connecting the monomer sites where $m_{x,t}=1$. As explained in Ref. [73], the equality between Eqs. (59) and (61) represents a duality relation between strong coupling and weak coupling; in the former case, a large number of monomers appear and Eq. (59)) is easier to calculate than Eq. (61)), which becomes easier at weak coupling.

In combination with the hopping expansion (which amounts to expanding the exponential of the kinetic energy rather than the potential energy), one arrives at other useful sign problem-free representations of fermionic partition functions. One such example is well-known and is the case of non-relativistic fermions in 1D with two-body interactions [79]. Another more recently discovered case is that of non-relativistic fermions in 1D with four-body interactions, shown and used in Refs. [80, 81, 82], where also a fermion-bag type idea is used to sum over configuration clusters. Interestingly, baryons at strong coupling can also be described by bags where three quarks propagate coherently as a single free fermion (i.e. a baryon) inside bags, while the complementary domain displays quark- and di-quark-type excitations [83]. Other relevant examples can be found in Refs. [66, 84, 85].

Proof of trace-determinant identity

For completeness, we outline the proof of the trace-determinant identity Eq. (59)) used above. We have not seen this way to approach the proof anywhere else, and since we find it particularly clear, we include it here. First we quote the auxiliary identity

[TABLE]

which is a reformulation of a well-known identity in the operator formulation of quantum Monte Carlo (see e.g. [86]). Here, the left-hand side trace is over Fock space, the entries of the Fermi matrix on the right-hand side are themselves matrices (i.e. the above is shown in block form), and $\hat{A}_{k}=\sum_{i,j}[A_{k}]_{ij}\hat{c}^{\dagger}_{i}\hat{c}_{j}$ and $\hat{B}_{k}=\sum_{i,j}[B_{k}]_{ij}\hat{c}^{\dagger}_{i}\hat{c}_{j}$ are (generally non-commuting) one-body operators. For our purposes, $k$ represents a particular time slice, $\hat{B}_{k}=j_{x,k}\hat{n}_{x,k}$, such that $[{\rm e}^{B_{k}}]_{x,x^{\prime}}=\delta_{x,x^{\prime}}{\rm e}^{j_{x,k}}$, and $\hat{A}_{k}$ encodes the kinetic energy, i.e. it is actually a $k$-independent operator. Our focus is on the $B_{k}$ factors.

The crucial step in proving Eq. (59)) is in differentiating both sides of Eq. (69)) with respect to $j_{x,k}$ in as many points $\{x,k\}$ as needed (namely the points where $m_{x,k}=1$; recall each point appears only once in the matrix) to match the desired insertions of $\hat{n}(x)$ at the desired time slices $k$. To carry out the differentiations on the right-hand side, we first use the Laplace cofactor expansion of the determinant and set the corresponding $j_{x,k}$ to zero once the derivative is taken. Once all the desired derivatives are applied, what remains is the corresponding cofactor determinant. The latter is the determinant of the matrix on the right-hand side of Eq. (69)) where the rows and columns of a given ${\rm e}^{B_{k}}$ containing the differentiated points $\{x,k\}$ are simply dropped, and the sources $j$ in any remaining terms are set to zero; that prescription defines the square matrix $W_{s}[\{m\}]$, whose number of rows (and columns) is reduced from the original matrix by the number of non-zero monomers. Note that there will, typically, be more than one spatial point $x$ affected within a given temporal block ${\rm e}^{B_{k}}$; similarly, any temporal block not affected by the derivatives will turn into an identity matrix once the sources are set to zero. Note also that any overall signs are unimportant because they cancel against the corresponding expression for the other fermion species.

2.6 Majorana fermions

In Ref. [87] a fermion representation was introduced that does not display a sign problem for a broad class of systems. Those developments were precipitated in part by the work of Huffman and Chadrasekharan of Ref. [78] and were further investigated by several authors in different ways (see in particular Refs. [88, 89]). The main result of that line of research, which we will explain in this section, is that there is a new class of systems that do not have a sign problem, and that that class goes beyond the well-known time-reversal-symmetric situation.

To understand the main principle behind this new class of systems, note that for a typical non-relativistic, single-species Fermi system, the discretization of the time direction into $N_{\tau}$ slices followed by a Hubbard-Stratonovich transformation yield a partition function of the form

[TABLE]

where

[TABLE]

which satisfies $\det M[\sigma]=\det\left(1+U[\sigma]\right)$, where

[TABLE]

and the ${A_{k}}$ factors encode the kinetic or potential energy contributions (the latter in the form given by the choice of Hubbard-Stratonovich transformation). It is then clear that the question of which systems display a sign problem amounts to asking under what conditions expressions of the form $\det\left(1+U[\sigma]\right)$ have a constant sign. That question was “crowd sourced” by L. Wang on the MathOverflow website and then analyzed in detail by multiple authors, leading to the work of Ref. [88] (see also Ref. [24]), whose discussion we parallel next.

Under specific conditions on the $A_{k}$ matrices, reviewed below, the product $U[\sigma]$ lies in the split-orthogonal group $O(n,n)$, which is defined as the set of real matrices $R$ such that $R^{T}\eta R=\eta$, where the metric is

[TABLE]

where $1$’s and $-1$’s appear $n$ times each. It is easy to see that $|\det R|=1$ and it is also possible to show that, writing $R$ in terms of $n\times n$ blocks

[TABLE]

that $|\det R_{11}|\geq 1$ and $|\det R_{22}|\geq 1$, which defines four different sectors in $O(n,n)$ labeled by the signs of these two determinants, typically denoted $O^{\pm\pm}(n,n)$. Of those, only $O^{++}(n,n)$ contains the identity and forms a subgroup. Crucially, it can be shown that, if $U[\sigma]\in O^{++}(n,n)$ then $\det\left(1+U[\sigma]\right)\geq 0$, and if $U[\sigma]\in O^{--}(n,n)$ then $\det\left(1+U[\sigma]\right)\leq 0$. If $U[\sigma]$ is in not in those sectors, then the determinant vanishes.

Finally, we see that, if the generating matrices $A_{k}$ are in the algebra of $O(n,n)$, i.e. if $\eta A_{k}\eta=-A_{k}^{T}$, then their exponentials will be group elements and the product of such exponentials will be in $O(n,n)$, in particular it will be in $O^{++}(n,n)$. Furthermore, parametrizing one such algebra generator as

[TABLE]

the condition of it being in the algebra of $O(n,n)$ implies $X=B^{T}$, $C^{T}=-C$, and $D^{T}=-D$, such that the diagonal of $A$ must be entirely zero. Below we will consider bipartite systems and $C$ and $D$ will contain matrix elements corresponding to same-lattice indices, whereas $B$ and $B^{T}$ will connect different sublattices.

As an example of a model whose auxiliary-field representation satisfies the above constraints, Ref. [88] considers the spinless $t$-$V$ model on a bipartite lattice:

[TABLE]

where $\hat{c}^{\dagger}_{\bf i}$ and $\hat{c}_{\bf i}$ are the creation and annihilation operators and $\hat{n}_{{\bf i}}$ is the number density operator at site $\bf i$. The angle brackets denote a sum over nearest neighbors, which are assumed to belong to different sublattices. The two key properties of this model are: a) the kinetic matrix $K$ only connects terms across the two sublattices and is zero on the diagonal; b) the interaction term can be decoupled using a Hubbard-Stratonovich transformation which again only connect different sublattices (and to that end the constant $\Gamma$ plays a crucial role). Specifically, Ref. [88] proposed the following auxiliary-field representation:

[TABLE]

where $\lambda=\text{acosh}(1+V/(2\Gamma))$, which is a real number for $V,\Gamma\geq 0$, i.e. repulsive interactions. This is easily checked by using the fact that, when $\bf i\neq\bf j$, all positive even powers of the combination $(\hat{c}^{\dagger}_{\bf i}\hat{c}_{\bf j}+\hat{c}^{\dagger}_{\bf j}\hat{c}_{\bf i})$ take the same operator value, namely

[TABLE]

and all odd powers of $(\hat{c}^{\dagger}_{\bf i}\hat{c}_{\bf j}+\hat{c}^{\dagger}_{\bf j}\hat{c}_{\bf i})$ also take on one operator value, which then vanishes upon summing over $\sigma=\pm 1$.

As usual, there are two different kinds of matrices $A_{k}$: one for the kinetic energy operator $A^{(K)}$, and one for the potential energy $A^{(V)}$ (which includes the Hubbard-Stratonovich field). Since $K$ is a real symmetric matrix and only connects different sublattices, $A^{(K)}$ features $C=D=0$ and $B=B^{T}$ [c.f. Eq. (80)]. On the other hand, the potential energy factor resulting from Eq. (82) is also real and symmetric and is designed to connect different sublattices only, such that also in this case $C=D=0$ and $B=B^{T}$. Thus, with the above choice one is within the purview of the theorem of Ref. [88] and there is no sign problem for the spinless $t$-$V$ model on a bipartite lattice. As anticipated in a previous section, the choice of Hubbard-Stratonovich transformation (of which there exist an infinite number for any given system) can determine the appearance or not of a sign problem, and that is the case here.

It is possible and useful to recast the above discussion in terms of Majorana variables. Reference [87] showed that, by writing the fermion operators in the original Hamiltonian as

[TABLE]

where $\hat{\gamma}_{i,{\bf j}}$ are Majorana fermions, it is possible to avoid the sign problem in certain classes of spinless fermion models on bipartite lattices. Moreover, using this type representation, it is possible to generalize the conclusions obtained using the split-orthogonal group regarding the types of systems that do not display a sign problem. That generalization was carried out in Ref. [89]. There, it was shown that a system displays no sign problem if it admits a Majorana decomposition in which the usual kinetic and potential energy factors (after the Hubbard-Stratonovich transformation) take the bilinear form

[TABLE]

where the vector $\hat{\gamma}$ contains the operators $\hat{\gamma}_{i,{\bf j}}$ above (in some order), and

[TABLE]

where $C=-C^{T}$ is complex antisymmetric and $B$ is Hermitian positive (or negative) semidefinite. Based on the above result, Ref. [89] showed that not only models like the above $t-V$ case can be made sign problem free, but also cases in which coupling to a pairing channel is present can be shown to have no sign problem using Majorana fermions. A precursor to that result appeared early on in Ref. [90].

In a second theorem, Ref. [89] also showed that there is another class of systems which, though partially overlapping with the above, represents a new class of sign problem-free systems as a result of Majorana-Kramers positivity. For the latter to hold, operators $S$ and $P$ must exist such that $V$ satisfies

[TABLE]

where $S$ is a real antisymmetric matrix satisfying $S^{2}=-I$ and $S^{T}=-S$, $P$ is a symmetric or antisymmetric Hermitian matrix satisfying $P^{2}=I$ and $PS=-SP$. The first of the above equations ensures that $V$ is time-reversal symmetric, which by itself is not a sufficient condition to avoid the sign problem. The second equation enforces a Kramers degeneracy that ensures that there is no sign problem.

A characterization of the classes of sign problems that can be addressed with this technique can be found in Refs. [89, 91] (see also [24] for a recent review). Notably, exponentiating the generator $A$ of Eq. (86) yields elements of $O(n,n)$, which is a real group, while the exponentials of $V$ yield elements of $O(m,\mathbb{C})$. Reference [92] presented a more general approach to such a group-theoretic characterization (so far the most general, to the best of our knowledge), based on Lie semigroups. Various applications were explored in Refs. [93, 94, 95, 96, 97].

Rather than a new method or an algorithm, one may regard the Majorana representation as a way to discover systems that do not have a sign problem (and for which it is not otherwise obvious that this is the case in conventional fermion formulations). Once that property is established, conventional algorithms can be used to carry out the calculation.

2.7 Imaginary asymmetry

In non-relativistic physics, fermions at finite polarization or mass imbalance present a sign problem because although the determinant in Eq. (11) factorizes and the factors are real, they will not typically be equal for all values of the HS field. A trick to overcome that problem can be borrowed from condensed matter theory [98] which is to make the asymmetry imaginary. For instance, if there are two spin flavors with corresponding chemical potentials $\mu_{\uparrow,\downarrow}=\mu\pm\delta\mu$, then taking $\delta\mu\to i\delta\mu_{R}$, where $\delta\mu_{R}$ is real, such that $\delta\mu$ is purely imaginary, makes the fermion determinants complex conjugates of one another and their product is real and positive, i.e.,

[TABLE]

This trick allows efficient sampling of the modified path integral via conventional Monte Carlo approaches. Below, we shall refer to imaginary asymmetry methods in general by the abbreviation “iHMC”.

Naturally, a caveat of this approach is that one must return to the real asymmetry axis (see left panel of Fig. 5). This usually involves a fit of an ansatz to the numerical data which then needs to be analytically continued. At this point, some degree of uncontrolled approximation enters the analysis as the ansatz for the fit is not unique. Nevertheless, the fact that investigations on the imaginary asymmetry axis are done in an entirely nonperturbative and controlled way is certainly an attractive feature.

In relativistic physics, the equivalent of the above is the introduction of a finite chemical potential that breaks the particle-antiparticle symmetry, thus inducing a finite difference between the densities of particles and antiparticles. The associated breaking of the charge-conjugation symmetry creates a sign problem which can be cured by rendering the chemical potential entirely imaginary, not only the difference of the chemical potentials as in non-relativistic field theories. This idea was originally put forward in Ref. [99] in the late 1990s and was very successfully employed in lattice QCD shortly thereafter [100, 101, 102]. Since then, this approach has proven to be very valuable for studies of thermodynamics as well as the phase structure of QCD, see, e.g., Refs. [103, 104, 105, 106, 107, 108, 109] for more recent results. For a discussion of the analytic continuation and suitable functional parametrizations of the data in the case of QCD, we refer the reader to Refs. [110, 111, 108] and to Ref. [112] for a detailed analysis of this issue with the aid of an exactly solvable field-theoretical model.

In Fig. 6 (left panel), for illustration purposes, we show results for the dimensionless baryon density $(1/T^{3}){(\rm d}p/{\rm d}\mu_{\rm B})$, rescaled by a factor of $T/\mu_{\rm B}$, for (2+1)-flavor QCD, see Ref. [108]. Here, $\mu_{B}$ is the baryon chemical potential. To obtain the results for $(\mu_{B}/T)^{2}>0$ (physical case) in Fig. 6 (left panel), various functional forms have first been fitted to the points on the negative side of the horizontal axis and then analytically continued to the positive side. The invariance of QCD under a sign flip of the baryon chemical potential has been used in constructing the fit functions. Whereas the predictions for $(\mu_{B}/T)^{2}>0$ are impressively independent of the used functional form at small chemical potential, the uncertainty grows large when the chemical potential is increased, illustrating the limitations of this approach. In any case, we note that it is not possible to study the zero-temperature limit of QCD at finite baryon chemical potential with this approach. Indeed, because of the Roberge-Weiss symmetry [113], only the regime $|\mu_{\rm B}|/T<\pi$ is accessible.

Inspired by these studies, the same principle has been followed in non-relativistic physics in recent years. More specifically, it has been applied to chemical potential and mass asymmetries [114, 115, 116, 117, 118]. An example is shown in the center and right panel of Fig. 6 where the magnetization $m$ (in units of the noninteracting density $n_{0}$) of one-dimensional non-relativistic spin-$1/2$ fermions at several values of the chemical potential $\mu$ (given in units of the inverse temperature $\beta=1/T$) is depicted. In this case, calculations avoiding the sign problem are possible at imaginary chemical potential asymmetry $h_{I}$ (center panel). To obtain the results for real asymmetry $h$ (right panel), an ansatz for the magnetization in terms of a function is first fitted to the data for imaginary asymmetry (center panel) and then analytically continued to real asymmetry $h$.

Here, $h=\mu_{\uparrow}-\mu_{\downarrow}$ is the difference in the chemical potentials of the spin-up and spin-down component. Similarly to the case of finite baryon chemical potential in QCD, it is not possible to compute the magnetization as a function of chemical potential asymmetry at zero temperature with this approach. In fact, it is only possible to reach values in the regime $|h|/T<\pi$ because of the $2\pi$-periodicity of non-relativistic fermions at finite temperature [115]. As already stated above, the basic idea of “taking a detour in the complex plane of the parameter space” is not limited to chemical potentials. In fact, the energy equation of state of 1D non-relativistic two-component fermions coming with different masses has been successfully studied in Ref. [118] using an imaginary mass difference. Other interesting applications, such as imaginary angular velocity coupled to angular momentum, remain unexplored to date but could in principle also be studied in this way.

It should be pointed out that, for small asymmetries, one can avoid the problem of analytic continuation completely by performing a Taylor expansion of the path integral around vanishing asymmetry, where calculations can be carried out without a sign problem. This approach has been successfully employed in many finite-temperature lattice QCD studies at finite baryon chemical potential to compute the equation of state and extract the phase structure, at least at small baryon chemical potential. See Refs. [119, 120, 121, 122] for ground-breaking studies of 2- and (2+1)-flavor QCD with this approach and Refs. [123, 124] for recent state-of-the-art results. For example, we can exploit the relation between the baryon density $n_{\rm B}$ and the QCD partition function $\mathcal{Z}$:

[TABLE]

where $V$ is the spatial volume and the temperature is introduced to render the expression dimensionless. At finite baryon chemical potential $\mu_{\rm B}$, the partition function is invariant under $\mu_{\rm B}\to-\mu_{\rm B}$ and therefore the density can be written as a series in odd powers of $({\mu_{B}}/T)^{2}$. The coefficients of this series can then be computed rigorously with stochastic methods by realizing that they are directly related to derivatives of $\ln{\mathcal{Z}}$ with respect to $\mu_{\rm B}$ evaluated at $\mu_{\rm B}=0$. Moreover, the pressure equation of state can eventually be obtained by integrating the density with respect to $\mu_{\rm B}$ since $p=(T/V)\ln{\mathcal{Z}}$.

Essentially, this amounts to the computation of static response functions. This approach can indeed be very efficient provided that these functions can be calculated in a statistically controlled manner, see also Refs. [121, 125] for a discussion of the reliability of this approach. Presently, the effort has been pushed to sixth order in the baryon chemical potential [123]. Of course, a similar approach can also been applied to non-relativistic theories to study equations of state as well as the phase structure and appears to be a worthwhile endeavor, also to cross-check results obtained by using imaginary asymmetries.

2.8 Lefschetz thimbles

In the presence of a phase problem, i.e. when the action $S[\phi]$ has real and imaginary parts: $S[\phi]=S_{R}[\phi]+iS_{I}[\phi]$, it may be possible to deform the integration contour away from the real line $\phi\in(-\infty,\infty)$ and into the complex $\phi$ plane (see right panel of Fig. 5) such that $S_{I}[\phi]$ is constant, or approximately so. If our problem concerns a simple one-site model, such that the partition function is a one-dimensional integral

[TABLE]

then observables take the form

[TABLE]

The goal of the Lefschetz thimbles method is to achieve a deformation of the integration region into a new region $\mathcal{C}$ such that we may alculate the above as

[TABLE]

[TABLE]

where $\phi$ is now regarded as a complex variable and we have used the fact that $S_{\rm I}[\phi]$ is constant along $\mathcal{C}$, i.e. we have cancelled it in both Eqs. (93) and (94).

Such a contour deformation changes neither the theory nor the observables, as the integrand in the path integral of any theory of interest will be an analytic function of $\phi$. If such a contour can be determined either a priori or dynamically during a calculation, the sign problem could potentially be solved or at least tamed. As a Monte Carlo method, the idea can be traced back to the work of Ref. [126] where the so-called shifted contour auxiliary-field Monte Carlo method was put forward for electronic systems (see also Ref. [127, 128]).

To extend the case of the simple one-dimensional integral discussed above to general QFTs, one complexifies the field variable in accordance with complex Langevin (which will be explored in the remainder of this review). In such a complexified configuration space, the Lefschetz thimbles approach aims to find the stationary points for which $\delta S_{\rm I}[\phi]/\delta\phi=0$, as those points feature reduced phase oscillations for $\exp\left(-{\rm i}S_{\rm I}[\phi]\right)$. Such an approach is of course the generalization of the saddle-point (or critical point) method of evaluating complex integrals and the higher-dimensional version is often referred to as Picard-Lefschetz theory. The corresponding deformed, high-dimensional integration contours of steepest descent are called Lefschetz thimbles.

In practice, the locations of such (stable) points of steepest descent are found by evolving the field along a fictitious time $t$, which is similar in spirit to the fictitious Langevin time (although unrelated to the imaginary time $\tau$). This propagation proceeds according to the holomorphic gradient flow equation

[TABLE]

or, more explicitly,

[TABLE]

which, as we shall see below, is remarkably similar to the CL equations. In fact, the above expression corresponds, up to a sign in the imaginary part, to the drift term in the Langevin equations Eq. (135), given by $d\phi/dt=-{\delta S[\phi]}/{\delta\phi}$. Ref. [129] presents a particularly lucid side-by-side discussion of CL versus Lefschetz thimbles approaches for a simple quartic integral. There, it is shown that there are similarities in that the location of the distributions in the complex plane obtained from CL and thimbles follow each other closely. However, there are also important differences with regard to the weight distribution and the role of the residual phase across the thimble (see below), suggesting that, despite the structural similarity of the above equations, the methods behave quite differently in practice.

The crucial advantages of deforming the path integral to capture, effectively, a set of mean-field configurations and the corresponding fluctuations, are that $S_{I}[\phi]$ is locally constant and that the real part, which determines the weight $\exp\left(-S_{R}[\phi]\right)$, is maximally localized around the saddle point. In other words: Lefschetz thimbles are the best locations to carry out stochastic evaluations of path integrals. While the constant-$S_{I}[\phi]$ property is crucial, the method does not rely on finding the precise location of the critical points. Rather, it is based on finding a useful deformation which may or may not be close to the stationary phase contours attached to the critical points (wherever those may be), but where the variations in $\textrm{Im}S[\phi]$ are small [130]. Those ideas were crucial for the application of Refs. [131, 132] where they were used to calculate the properties of a low-dimensional field theory in real time. Furthermore, Ref. [133] found that, even if the holomorphic gradient flow of Eq. (96)) is used to push the deformation very close to the thimbles (which is naively the ideal situation), then the high barriers separating the thimbles would make the sampling very challenging in practice.

One of the difficulties of deforming a contour in a functional integral is the calculation of the associated Jacobian factor $J(t)$ due to the curvature of the thimble, which shows up explicitly upon parametrizing the complex-plane integral of Eq. (93) by a parameter $t\in(-\infty,\infty)$, namely

[TABLE]

While the above is a computational issue that is difficult but tractable [134], a more serious challenge lies in accounting for all the possible thimbles, which turns the above into a sum with unknown weights, i.e. in that case

[TABLE]

and

[TABLE]

where, crucially, the phases in the numerator and denominator of Eqs. (99) and (100) do not cancel. This issue is often referred to as a ‘global’ sign problem, as opposed to the ‘residual’ sign problem coming from the remaining curvature (i.e, variations in the imaginary part of the Jacobian across the thimble). In that regard, it is useful to note that (for a fixed set of input parameters), only a subset of thimbles contribute to the partition function. The global sign problem depends on the subset and the relative weights of each thimble, both of which are difficult to determine in practice (see e.g. Ref. [135]). However, the holomorphic flow bypasses that complication. It is then possible to track the global and the residual sign problems on the deformed surface by measuring the average phase.

Recently, contour deformation and Lefschetz thimbles have re-emerged as a method in quantum field theory and its applications have extended to condensed matter physics as well. References [136, 137, 138] were the first ones to propose that the sign problem in finite density QCD could be overcome using Lefschetz thimbles. Refs. [139, 140] proposed ways to overcome the residual sign problem across a thimble (see below). Reference [141] studied the structure of thimbles in fermionic theories. Ref. [142] used thimbles to avoid the sign problem in a mean-field analysis of the QCD partition function. while Ref. [143] interpreted the Silver Blaze problem of finite-density QCD in a one-site fermion model. Reference [144] used a generalization of the Lefschetz thimbles method to study the finite-density Thirring model in two spacetime dimensions, and Ref. [145] generalized that study to QED in two spacetime dimensions.

Other interesting connections have also been drawn. For instance, Ref. [146] showed that the steepest descent trajectories (i.e. the deformed integration contours mentioned above) can be interpreted as ground-state wave-functions of a supersymmetric Hamilton dynamics. Connections to the complex Langevin method and its convergence shortcomings (see also below) were pointed out in Ref. [147] and Ref. [148].

On the condensed matter side, Refs. [149, 150, 151] studied the Lefschetz thimbles representation of the Hubbard model and studied it in a hexagonal lattice away from half filling. Ref. [152] applied a variant of the Lefschetz thimbles method (so-called “tempered” Lefschetz thimbles method, developed in Ref. [153]) to the Hubbard model away from half filling.

In light of the excellent recent review article on Lefschetz thimbles and its applications [154], we limit ourselves to the above discussion.

2.9 Path optimization method

The path optimization method (POM) refers to a relatively novel complex plane approach, whose aim is to shift the integration away from the real line in order to minimize the sign fluctuations along the new integration path in the complex plane [155, 156, 157]. Although the sign problem will not be completely absent on the new integration contour, the hope is to ameliorate it enough such that the signal-to-noise ratio is manageable.

Similar to the above discussed method of Lefschetz thimbles, the path integral can be expressed in terms of a shifted path $\phi(t)$ parametrized by the real parameter $t$ and under consideration of the Jacobian $J[\phi(t)]$:

[TABLE]

Of course, in order for the Cauchy theorem to hold, the integration path should not enclose singular points of the Boltzmann factor. The method relies on a trial function which parameterizes the integration path in the complex plane by some parameters, collectively denoted as $\alpha$. Following Ref. [155], we may for instance expand the path as a sum over a complete set of polynomials

[TABLE]

although this is only one particular choice. The important features to consider (besides some technical aspects, see, e.g., [154]) is that the chosen parameterization allows contours with a more suitable sign-structure and that the Jacobian may be evaluated efficiently. To find the path with the mildest sign problem, the parameters $\alpha_{n}^{\rm R}$ and $\alpha_{n}^{\rm I}$ will be optimized by minimizing some cost function. In Ref. [155] the cost function

[TABLE]

was applied, but other cost functions can be used. Here $\theta_{\alpha}(t)$ is the complex phase of our parameterized integrand $J[\phi_{\alpha}(t)]{\rm e}^{-S[\phi_{\alpha}(t)]}$, as defined in Eq. (101), and $\theta_{0}$ is the complex phase of the original integrand. The cost function is then used to tune the parameters $\alpha_{n}^{R}$ and $\alpha_{n}^{I}$, such that $F[\phi_{\alpha}(t)]$ is minimized. This results in an enhanced phase factor, thus reducing or ideally eliminating the oscillations discussed in Sec. 2.1, which are the source of the sign problem.

For the optimization, several algorithms may be used, which mostly rely on the computation of the gradient in $\alpha$-space $\nabla_{\alpha}\left\langle{\rm e}^{\theta_{\alpha}(t)}\right\rangle$. A convenient (and important for practical implementations) feature is that it is typically enough to have a rough knowledge of the gradient such that efficient optimization is ensured. Naturally there is significant overlap with methods in Machine Learning (ML). A simple steepest descent method can be used, but other methods for minimizing the cost function have been explored in ML and neural network (NN) applications [158, 159].

While this method is very new, early work using shows it is able to overcome a severe sign problem in a one-dimensional toy model where CL fails [155], one dimensional Bose gases with chemical potential [160], as well as in a U(1) gauge theory and complex scalar field theory [161]. Work towards applications to full $3+1$ dimensional QCD has included $0+1$ dimensional QCD at finite density [162] and effective models for QCD [163, 164] and the Thirring model [165, 166].

3 The Langevin method for real and complex variables

Stochastic quantization as a method for treating Euclidean field theories has been around since Parisi and Wu first proposed the connection between the Euclidean field theories and statistical systems coupled to a heat bath [167]. It is now well-established as a successful tool for treating quantum many-body systems with a real Euclidean action [168]. This section examines the method for systems without a sign problem (referred to from here on as “real Langevin”) and its extension to systems with complex actions as a possible circumvention of the sign problem (“complex Langevin” or CL, the focus of this review). We present a pedagogical example to illustrate the real and complex Langevin methods using a simple toy model, and discuss some of the challenges that arise in using the complex Langevin method along with proposed solutions to those problems.

3.1 Complex Langevin: origins and modern re-emergence

Shortly after the introduction of the concept of real stochastic quantization, it was realized that the approach could be extended to the case of complex actions. Loosely speaking, using a Langevin equation rather than an importance sampling approach eliminates the restriction to real and positive semidefinite measures. This is due to the ability of the probability measure used in a Langevin method to be complexified – at least formally. In 1983 such a strategy was discussed independently by Klauder [169, 170] and Parisi [171] as an alternative to existing Monte Carlo methods and marked the first investigations of how the complex Langevin equation can be used to address the complex phase problem.

The elegant form of the approach as well as the potential to circumvent the sign problem drew considerable interest in the years after the initial proposals. Following the first successful numerical application for the quantum Hall effect by Klauder [172], the method was employed in many studies, albeit with mixed success. Unfortunately, the convergence of the method cannot be guaranteed a priori and even if convergence is achieved, spurious solutions with biased expectation values might be found [173, 174]. This is connected to subtle mathematical issues that arise in the case of complex weight and the structure of the associated complex Fokker-Planck equation. As a consequence, the initial flurry of interest stalled and progress on these matters slowed down over the years, despite early attempts to understand these shortcomings [175, 176, 177, 178]. Nevertheless, progress was made in several directions and the applicability of CL was investigated in a set of simplified models such as the chiral Schwinger model [179, 180] as well as toy problems for relativistic [181] and non-relativistic fermionic theories [182]. Interestingly, the method has also spread to the realm of physical chemistry and was used in simulations of polymeric fluids [183, 184, 185] as well as reaction simulations in the context of physical chemistry [186, 187, 188, 189, 190, 191] as a way to include beyond mean-field corrections.

As recently as the mid 2000s to early 2010s, the CL approach re-emerged as a method of interest in relativistic physics, particularly in the study of lattice QCD, when it was realized that some of the initially encountered problems are treatable with an improved integration strategy. In this new era, early applications to relativistic physics examined non-equilibrium QFT [192, 193], which can provide insights into high-energy physics, particularly heavy ion collisions. Moreover, in 2008, Aarts and Stamatescu demonstrated that CL could be applied to models of finite density QCD that exhibit a sign problem [194]. Shortly after that, Aarts demonstrated that CL could be used to circumvent the sign problem in the relativistic Bose gas with finite chemical potential [195, 196, 197]. This began a resurgence of interest in this method in the field of finite density Lattice QCD (LQCD), in which nonperturbative calculations of strongly interacting matter with finite baryon chemical potential are inhibited by the sign problem. This renewed interest led to work in the next few years on optimization of the method to prevent runaways and improve stability, using stochastic reweighting, gauge fixing, and adaptive step size algorithms [198, 199, 200], see also Fig. 7 for a history of the “CL activity on the arXiv”.

The successes of the method, and advances made in treating instabilities and singularities in the fermion determinant, have generated interest in applying CL to non-relativistic systems, particularly many-fermion systems, in which sign problems arise frequently. Work with the CL method in the context of non-relativistic systems is just beginning, but is already showing great promise. Discussion of these recent applications to relativistic and non-relativistic physical systems can be found in Secs. 4 and 5.

Despite the formal challenges of CL, it should be pointed out that, at least in principle, the existence of a well-defined probability measure on a complexified field space is guaranteed under certain conditions (often referred to as Weingarten’s theorem [201]). However, the challenge lies in constructing such a measure, which has been investigated by Salcedo and others  [202, 203, 204, 205, 206, 207, 208, 209]. While this is a very attractive area of research, we will not pursue it further in this review.

3.2 Stochastic quantization: path integrals and the Langevin process

The main ingredient of the CL method is the concept of stochastic quantization. The idea was introduced in the seminal 1981 paper of Parisi and Wu [167] and a few years later summarized nicely in the famous review article by Damgaard and Hüffel [168] 222See also Ref. [210] for a more formal introduction to stochastic quantization. Notably, stochastic quantization has played an important role both theoretically as well as computationally; in particular, it has been used extensively in field theory (see e.g. Ref. [211]) and condensed matter (see e.g. Ref. [212]), and is the precursor of the hybrid Monte Carlo algorithm [213, 214], which has been the workhorse of LQCD for decades. In the following we present a brief introduction to stochastic quantization in order to lay out the foundation for later considerations.

For a quantum field theory of a real field $\phi$ governed by a real action $S[\phi]$, stochastic quantization provides an intuitive way to understand path integrals of the form

[TABLE]

As a first step, we introduce a purely fictitious time variable $t$, which represents the direction of the stochastic evolution. This fictitious time evolution is then governed by a stochastic differential equation, namely the Langevin equation:

[TABLE]

The first term on the right hand side is called the drift term, also sometimes referred to as the classical flow as it constitutes the deterministic part of the time-propagation of the fields. The second term on the left hand encodes the random nature of the equation and is given by a white-noise with zero autocorrelation, i.e. $\tilde{\eta}\sim\delta(t-t^{\prime})$.

For practical purposes it is convenient to rewrite this equation in a discrete form which can be done by integrating both sides over the time interval $\Delta t$. This leads to the discrete Langevin equation

[TABLE]

where $\eta$ is typically chosen to be a Gaussian random variable with $\langle\eta\rangle=0$ and $\langle\eta^{2}\rangle=2\Delta t$. The angle brackets denote an average over $\eta$. This random process will produce configurations $\phi$ that follow a certain probability distribution. The key ingredient of stochastic quantization is the realization that the equilibrium distribution (if it exists) of the $d+1$ dimensional random process in Eq. (106) corresponds to the probability measure in the $d$-dimensional path integral Eq. (104). The extra dimension is simply the fictitious time $t$.

It is instructive to consider the above stochastic differential equation without the noise term. In that case, Eq. (105) reduces to an ordinary differential equation and its form is nothing but that of a gradient descent. Starting out at a random (non-pathological) state, this implies that the solution will converge to a stationary point of the action, which is the “mean-field” or classical solution. The simple interpretation of the noise term is that it represents quantum fluctuations around this classical solution. In order to reproduce the correct physics, we have to “add the correct amount of fluctuations” which is set by the Fluctuation-Dissipation theorem. Thus, stochastic quantization can be viewed as a very explicit form of quantization.

The random process results in a sequence of time-dependent configurations $\phi_{\eta}(t)$ distributed according to a time-dependent probability measure $P[\phi,t]$. The expectation value of a given observable $\mathcal{O}[\phi]$ is then given by

[TABLE]

To establish the validity of such a Langevin average, we need to investigate the temporal behavior of the expectation value and show that the time-dependent probability distribution depends on $S[\phi]$ in the way dictated by Eq. (104), at least at large enough $t$. Such a property will justify the use of temporal averages along the Langevin evolution to estimate the true expectation values of the theory. An instructive discussion can be found in [215], which we will follow closely.

As a first step, we take the fictitious-time derivative of the expectation value

[TABLE]

Note, that only the probability distribution carries a fictitious-time dependence. Alternatively, we may perform the same fictitious-time derivative by expanding the observable to second order in its $\phi$ dependence

[TABLE]

According to Eq. (106) we may write

[TABLE]

where ${\rm d}w$ is the so-called Wiener increment with the property

[TABLE]

and vanishing mean $\langle{\rm d}w\rangle=0$. Substituting in Eq. (109) and using the properties of ${\rm d}w$ yields

[TABLE]

which allows us to write

[TABLE]

where we defined the so-called Langevin operator

[TABLE]

with the drift $K[\phi]=-\frac{\delta S[\phi]}{\delta\phi}$, according to Eq. (106). We may again write this as an integral over configurations

[TABLE]

where the second line results from partial integration. Here we made the important assumption that the probability vanishes at the boundaries (or decays fast enough if the integration region is non-compact). These assumptions are in fact crucial and will be discussed in more detail below.

Comparing the above equation with Eq. (108) yields the Fokker-Planck (FP) equation:

[TABLE]

with the formal adjoint of the above Langevin operator

[TABLE]

which is also referred to as the FP operator or FP Hamiltonian. To show that the stationary solution of this equation is indeed our desired probability distribution we perform similarity transformation

[TABLE]

to rewrite the FP equation

[TABLE]

with

[TABLE]

This last equation reveals that, with a real action $S[\phi]$, our modified FP Hamiltonian is a self-adjoint and positive semidefinite operator, with a unique FP ground state $\psi_{0}={\rm e}^{-S[\phi]/2}$ and vanishing FP energy $E_{0}=0$. We can therefore project our probability over the complete set of eigenfunctions and non-negative eigenvalues of $\widetilde{L}_{r}^{\rm T}$, and see that our probability collapses to the ground state in the long time limit:

[TABLE]

Upon performing the back-transformation according to Eq. (119) we obtain

[TABLE]

which shows that the Langevin equation produces field configurations distributed according to the Boltzmann weight ${\rm e}^{-S[\phi]}$ in the limit of large fictitious time 333We have assumed here that the spectrum of $\widetilde{L}_{r}^{\rm T}$ is discrete, which may not be true in practice. This assumption can be relaxed but it is important that the $E_{0}=0$ eigenvalue be non-degenerate..

The above justifies the use of temporal averages to estimate equilibrium expectation values. In fact, expectation values are obtained in practice by integrating over a time $T$:

[TABLE]

where $t_{\rm th}$ reflects the equilibration time that is needed to approach the stationary probability distribution.

3.3 A practical guide to real Langevin

The above procedure is a well-established method for real-valued fields $\phi$ on a real manifold $\mathcal{M}$. In the following, we connect the concept with conventional Monte Carlo approaches based on Markov chains and highlight the similarities of these approaches in a practical way.

Generally, we are interested in expectation values of a given (Euclidean) field theory of the form of Eq. (107):

[TABLE]

Typically, the evaluation of such high-dimensional path integrals is achieved by stochastic sampling, i.e. by producing a sequence of random states $\phi$ (interchangeably called samples or configurations). However, rather than producing a random state from scratch at every step (which might be expensive), new states are obtained by changing or updating an existing one. This can be written in the generic form

[TABLE]

where the sample $\phi$ may be any representation of a physical state444In fact, finding suitable and efficient updates is a crucial part in devising any useful Monte Carlo method.. If the next state is only dependent on the current one and not on any previous states the process is called “memoryless” and the sequence of random samples represents a so-called “Markov chain”. This is the basis of the vast majority of most conventional Monte Carlo methods (often dubbed Markov-Chain Monte Carlo (MCMC) methods).

The Langevin method can be understood as such an approach, as we can see by recasting Eq. (106) into the form of Eq. (126):

[TABLE]

By virtue of the discussion in the previous section, we know that in the long-time limit the samples $\{\phi_{n}\}$ follow the desired probability distribution $P[\phi]$ from Eq. (125). Therefore, the strategy to evaluate the random sequence is identical to the one in regular Monte Carlo approaches: after starting out from a randomly produced configuration, we let the sampling process equilibrate for a certain thermalization time (typically a few multiples of the autocorrelation time – see below) before we start to collect samples. An unbiased estimator of the expectation value with a total number of samples $N$ is then given by

[TABLE]

and the statistical uncertainty (i.e., the variance of the mean) is estimated by

[TABLE]

Here, $\sigma$ denotes the standard deviation of the set of values $\{{\mathcal{O}}[\phi_{i}]\}$ and $\tau_{a}$ is the integrated autocorrelation time, which reflects the (statistical) dependence of samples over the their Langevin-time history555For any Markov chain method, estimation of $\tau_{a}$ might be a challenging task depending on the number of accessible samples. Techniques such as bootstrap and jackknife can be used to obtain a reliable error estimate that considers the statistical dependence of the samples. An educational summary can be found in Ref. [216].. From this expression we can also deduce that the number of samples is proportional to the elapsed Langevin time and the same conclusion holds: a longer Langevin evolution will result in better statistics and hence a smaller error bar.

It is important to note that the statistical uncertainty is not the only source of error, as we have introduced a systematic bias through the discretization of the Langevin equation in Eq. (106). Results must then be extrapolated to the limit of $\Delta t\to 0$ 666Notably, the hybrid Monte Carlo algorithm, mentioned above as a close cousin of RL, avoids this extrapolation in fictitious time by using Metropolis accept/reject steps. This property, however, does not imply that hybrid Monte Carlo can operate at arbitrarily large step sizes.. The discretization of the Langevin equation presented above, which corresponds to Euler integration, leads to a linear dependence on the integration step $\Delta t$. Higher-order integrators can be designed if the noise term is appropriately accounted for, as investigated within RL in Refs. [217, 218, 219] and later extended to the complex case in Ref. [220]. Within the studied model, these extensions showed great success in reducing computational cost at equal systematic bias as well as to bring finite-step dependence below the statistical uncertainty. It is noted, however, that the majority of stochastic quantization studies still rely on the linear discretization as it is often sufficient to obtain useful results at modest computational effort.

3.3.1 Toy problem I: a pedagogical example of real Langevin

In order to illustrate the RL method in a concrete numerical setting, we consider a simple integral as a toy problem for a $0+0$-dimensional field theory (i.e., the field $\phi$ depends neither on space nor time). Note that here we are not interested in a detailed study of the specific model at hand but rather aim to investigate the behavior of the RL method in a straightforward case. Of course it is highly inefficient to use RL for the solution of this simple problem, however, the section serves as a basic example of an application of the method. Furthermore, conclusions that generalize to more involved problems can be drawn by these simple considerations.

In this case, we consider the action

[TABLE]

with real couplings $\mu$ and $\lambda$. In the spirit of a real field theory, we will keep $\lambda$ positive and thus end up with two distinct scenarios: one in which $\mu>0$ (the single well anharmonic potential) and one in which $\mu<0$ (the double well potential).

We readily derive the discrete Langevin equation for our toy problem, according to Eq. (127):

[TABLE]

where $\eta$ denotes a standard Gaussian white noise. In principle, this is everything needed to calculate expectation values of the form Eq. (128).

In Fig. 8, a detailed analysis of two simulations at fixed $\lambda=0.4$ and $\mu=\pm 1$ is presented (left and right columns, respectively). The second row from the top shows the histograms of the sampled field values, which should follow the distribution ${\rm e}^{-S[\phi]}$ (exact solution shown with dashed lines) in the limit of large Langevin time $t\rightarrow\infty$. While in the single-well system (left column) this is the case to a good approximation, it is apparent that the double-well scenario still suffers from a slight asymmetry. This can occur when the random process gets “stuck” in an area of configuration space and does not easily move to another high probability area of configuration space (i.e. the other well).

This behavior can be further elucidated by the measured field values as a function of $t$ [row (c)]: on the left we see white noise centered around the expected value of [math] while on the right we observe several correlated plateaus, corresponding to either the negative or the positive well. Ultimately, this behavior leads to a signal-to-noise issue in the calculation of the expectation value $\langle\phi\rangle$, reflected by large statistical uncertainties for the double well case [right column of row (d)]. Due to the symmetry of the problem, however, the running average for $\langle\phi^{2}\rangle$ converges to the exact value relatively smoothly in both cases. Thus, by investigating one observable no profound statements can be made about a different one. While the autocorrelation of observable $A$ may be small and its statistical errors under control, observable $B$ could display erratic behavior and suffer from extremely slow convergence.

By tuning the parameters of the model, one could even study the extreme case where the two wells are separated by a barrier that cannot be surmounted by the random walk (signaling a breakdown of ergodicity, a topic that we will return to below). In such a situation, the expectation value $\langle\phi\rangle$ would indicate that the discrete symmetry is broken, which certainly is not a physical result for our model. This reflects the problem of meta-stability of any Markov chain method, which is often very hard to detect a priori. Generally, one needs to address this issue carefully in real simulations, for example by sweeping numerical parameters in a systematic manner.

Finally, by inspecting the Langevin equation in Eq. (131) it becomes apparent that the autocorrelation time between samples $\tau_{a}$ should be inversely proportional to the Langevin time step $\Delta t$, i.e., statistically independent samples will be more expensive as the integration step decreases. On the other hand, a coarser integration step will yield a larger systematic error and thus a balance must be found where both the computational effort as well as the precision are within reasonable bounds. This behavior is illustrated in Fig. 9, where we show the dependence for $\langle\phi^{2}\rangle$ as well as $\tau_{a}$ on the integration stepsize $\Delta t$. Indeed, we observe a systematic linear behavior of $\langle\phi^{2}\rangle$, which is expected due to the order of the Langevin equation (see discussion in the previous section). Further, also the integrated autocorrelation $\tau_{a}$ shows the expected behavior: it increases as the integration step approaches [math].

3.4 A practical guide to complex Langevin

Thus far we have seen that the concept of stochastic quantization works well, given that the action $S[\phi]$ is real, i.e. when there is no sign problem. But what about the much more general and interesting case of complex-valued actions? In that case, the probability distribution in Eq. (107) becomes a complex distribution

[TABLE]

while the field $\phi$ is still a real quantity. With such an action, a single step in the Langevin process according to Eq. (127) would result in an imaginary component for $\phi$. At least from a practical perspective, as remarked by Parisi, “nothing forbids to write a Langevin equation also for complex actions” [167]. Formal aspects aside (we return to those below), the idea is appliedthat one may perform calculations by extending the real stochastic process described above to a complex one. For a theory governed by a complex action $S[\phi]$, CL extends the target manifold of the field $\phi(x)$ to the complex plane by setting

[TABLE]

and analytically extending the domain of the action functional:

[TABLE]

Naturally, a necessary condition for this extension to be valid is that $S[\phi]$ must be a holomorphic function of $\phi$.

The two Langevin methods – real and complexified – are compared side-by-side in Tab. 1. With such an extension, the CL method proceeds very much in the same way as the real Langevin method, but now with a double system of coupled stochastic differential equations:

[TABLE]

where the real and imaginary drift functions $K_{\rm R}$ and $K_{\rm I}$ are found by taking the real and imaginary parts of the functional derivative of the complex action:

[TABLE]

The real and imaginary noise obey the properties shown in Tab. 1.

It is important to note that, while the amplitudes of the real and imaginary noise terms are related, the two Wiener processes are completely independent. Additionally it should be pointed out that, beyond the complexification of each real degree of freedom, the above (coupled) Langevin processes are themselves real. In practice, the imaginary noise is usually set to zero, which satisfies the constraints shown in Tab. 1 and it has been found to have the best numerical properties [221].

Although a rigorous derivation of this empirical “rule” is not available, some intuition may be gained as follows: For overly broad distributions of the imaginary noise (i.e. large $N_{I}$) the CL trajectory may be pushed into regions of the complex plane which are far away from the physical configuration space (which necessarily needs to be somewhat local in the imaginary directions as discussed further below). If the trajectory wanders off too far in the imaginary direction, it could take a long time to return and thus a long time to restore the proper sampling ratios required to reproduce physical results (in that sense, the strategy may be viewed as a variance reduction technique in a similar spirit to importance sampling). Eventually such a behavior will interfere with the decay of the sampled histogram in the imaginary directions, with possible implications for the formal validity of the algorithm, as already pointed out in Ref. [221]. Therefore, in order to avoid these potential complications, the noise term is typically chosen to be real, which seems to be the most convenient choice within the formal requirements.

We return to the justification for CL and related challenges after giving an illustrative example.

3.4.1 Toy problem II: a pedagogical example of complex Langevin

In order to see the CL machinery at work, we build on the toy problem in Sec. 3.3.1. Indeed, from a computational standpoint, CL is largely just the Langevin process of stochastic quantization with complex variables. Of course we could just turn our toy problem into a CL problem by using complex noise but, as mentioned above, real noise has preferable numerical properties. Instead we turn our toy problem into a complex field theory which results in complex drift terms.

The most natural route towards a complex field theory would of course be to consider complex-valued couplings $\mu$ and $\lambda$ in Eq. (130), which in fact has been considered before, see e.g. [215]. However, this would amount to solving a different theory as the couplings necessarily take on different values. Alternatively, we may rewrite the above problem with a suitable Hubbard-Stratonovich transformation which merely amounts to an alternative representation of the same physical scenario (note, however, that the representation given below is by no means unique). Moreover, this is very much in the spirit of real-world problems, where a HS-transform is often used to integrate out fermionic degrees of freedom. To achieve this, we insert a suitable factor of $1$ into the partition function in terms of an auxiliary variable $\sigma$:

[TABLE]

A shift $\sigma\rightarrow\sigma+{\rm i}\phi^{2}$ allows us to write

[TABLE]

and subsequently to integrate out the dependence on the old field $\phi$ (note that this is only possible in the case $\textnormal{Re}[\mu]>0$). Ultimately we obtain the “bosonized” version

[TABLE]

where we defined the “bozonized action”

[TABLE]

which is, by construction, a complex quantity. According to the discussion in the previous section we can still evaluate expectation values stochastically by using the complex Langevin equation

[TABLE]

Specifically, we write for the second moment of the initial field $\phi$ in terms of the new field $\sigma$:

[TABLE]

where the subscript $\sigma$ denotes averaging over different realizations of $\sigma$.

From this point on, we proceed exactly as in the real case, with the exception that we now have to deal with a complex variable $\sigma$. The qualitative dependence of the numerical results on the integration step size $\Delta t$ should still be the same. This is indeed the case, as apparent from the lower left panel of Fig. 10, where we show CL results for the model given by Eq. (130) with parameters $\mu=1.0$ and $\lambda=0.4$. We observe that the CL correctly reproduces the exact result. Interestingly the CL values show a much milder dependence on the integration step, which is likely a consequence of the specific representation and in any way should not to be interpreted as a general result.

It is instructive to investigate the sampled configurations by looking at the complex weight according to Eq. (142). This is shown in the upper left panel of Fig. 10 where it becomes clear that the Langevin process now samples a complex quantity which is necessary to represent the correct answer. The coloring highlights the sampled configurations of $\sigma$ corresponding to regions with ${\rm Re}[\sigma]>0$ and ${\rm Re}[\sigma]<0$. These two sets are connected through a point in configuration space which corresponds to a stationary point in the “classical flow”, i.e. the vanishing of the drift terms in Eq. (144). (The entire flow diagram is shown on the right panel of Fig. 10). The attractive fixed point (marked by the lower green dot) pulls the field towards the connection point, while the noise induces fluctuations around that point, which occur mainly in the real component, as we have chosen to use only real noise in Eq. (144) (such that the only way the imaginary part can change is through the complex part of the drift). Thus, the imaginary component stays approximately constant during the Langevin evolution.

The situation might change drastically if poles are encountered inside the domain of the distribution, which could lead to a breakdown of ergodicity, due in turn to a breakdown of holomorphicity of the action stated above (see also the discussion below). In the model considered here, we can find a pole at the point $\sigma=6\tfrac{\mu}{\lambda}$, which corresponds to the branch-point of the action (marked by the red cross). At first glance this looks very dangerous as there is also an attractive stationary point (upper green dot) right next to the pole which would suggest faulty behavior. However, the imaginary part of the drift points away from the pole and even if a trajectory approaches this area of configuration space, fluctuations (in the real direction) will kick the process back into a stable trajectory that decays towards the attractive fixed point below. In equilibrium, the distribution will thus be confined to the gray dashed line far away from the pole, ensuring correct behavior (i.e. it is approximately shifted from the real axis by a constant offset). In fact, the existence of such an attractive fixed point is a necessary condition for the existence of an equilibrium distribution of the Langevin process [222, 223]. Luckily, this appears to be the case for systems of physical interest.

3.5 Formal aspects and justification

The CL process defines a random walk in a complexified manifold, such that for a given configuration $\phi=\phi_{\rm R}+{\rm i}\phi_{\rm I}$ there is a well-defined probability $P[\phi_{\rm R},\phi_{\rm I},t]$ at time $t$. For a given observable $\mathcal{O}$, there will be an expectation value

[TABLE]

By virtue of the CL process, the real probability $P[\phi_{\rm R},\phi_{\rm I},t]$ obeys the FP equation:

[TABLE]

where

[TABLE]

It is not obvious a priori whether this process would reproduce the desired expectation values of the physical observables, i.e., whether $\langle\mathcal{O}\rangle_{P(t)}$ actually corresponds to the physical expectation value of the theory, at least in the large-$t$ limit. In fact, it is not even clear that the process would converge and if it does, whether it converges to the correct answer. Indeed, following the steps outlined above for the case of real Langevin, one finds that the resulting FP Hamiltonian is neither self-adjoint nor positive semidefinite, such that the proof of convergence to the desired probability distribution is spoiled.

The fundamental question underlying the validity of the CL approach is the relation between the CL distribution $P[\phi_{\rm R},\phi_{\rm I},t]$ and the desired complex distribution $\rho[\phi]$ in Eq. (132). The latter defines the physics of interest and is a fixed point of its own FP equation

[TABLE]

where

[TABLE]

which is obtained by temporal differentiation of

[TABLE]

Again, the absence of boundary terms at infinity when integrating by parts was assumed. More specifically, the crucial question is whether

[TABLE]

holds. Stated explicitly, the question arises whether the expectation value of a given observable with respect to the complex measure $\rho[\phi_{\rm R}]$ may be obtained by instead sampling the positive semi-definite probability measure $P[\phi_{\rm R}+{\rm i}\phi_{\rm I}]$.

In Refs. [221, 224] it was shown how the desired relationship Eq. (153) can be proven for holomorphic observables, as long as the action and the associated drift are holomorphic functions of $\phi$. The proof relies on analyzing the behavior of

[TABLE]

where $0\leq\tau\leq t$. The function $F(t,\tau)$ interpolates between the two expectation values of interest:

[TABLE]

where we have assumed that the initial conditions are chosen as

[TABLE]

We find that

[TABLE]

while, using the initial conditions of Eq. (156),

[TABLE]

where we have used Eq. (150) to shift the Langevin evolution operator from $\mathcal{O}$ to $\rho$ by transposition (which involves integration by parts).

Now the crucial observation is that if $F(t,\tau)$ is independent of $\tau$, then Eq. (153) holds, and the Langevin method is formally shown to be valid for complex-valued variables, i.e. to converge to the correct physical answers (assuming it converges). Naturally, this statement assumes that the expectation values in Eq. (153) agree at $t=0$, which can be ensured by choosing the initial condition of the Langevin process as in Eq. (156). The $\tau$ derivative of $F(t,\tau)$ again involves an integration by parts:

[TABLE]

where $L$ is the Langevin operator and $L^{T}$ its adjoint, as defined in Eq. (149). If integration by parts is carried out and – importantly – the boundary terms are zero, then $\frac{\partial}{\partial\tau}F(t,\tau)=0$. If the decay of

[TABLE]

and its derivatives is not fast enough to ensure that the boundary terms will vanish, then it cannot be guaranteed that the expectation values of the quantities of interest obtained via a Langevin process will converge to the correct values [221, 224]. It is this property that ultimately determines the applicability of CL to a given theory.

While the condition of fast decay of Eq. (160) was recognized in [221, 224], the precise rate was not immediately clear. In Ref. [225], the above arguments were reviewed by considering a finite step-size in Langevin time. It was then found that the above integration by parts is valid if the probability distribution of the drift term falls off faster than any power at large drift magnitude. In practice, it is very difficult to establish the behavior of Eq. (3.5), but it is perfectly possible to study the probability distribution of the drift and establish whether the decay is exponential.

3.5.1 Practical monitoring of the boundary terms

In practice, it would be useful to have an accessible diagnostic device which tells us whether the results of the CL evolution are trustworthy. Such a tool may be derived by re-examining the interpolating function in Eq. (154), which must be independent of the interpolating time $\tau$ in order for the boundary terms to vanish. It has been found [224] that the maximal violation, if there is any, of this condition occurs at $\tau=0$. Thus, we focus on the more straightforwardly accessible quantity

[TABLE]

which is weaker than the expression without the restriction to $\tau=0$, however, still a sufficient criterion to ensure correctness (up to some technical conditions, as shown in Ref. [224]).

Evaluating Eq. (3.5) at $\tau=0$ we observe that the first term vanishes, since $L^{T}P[\phi_{\rm R},\phi_{\rm I},t=\infty]=0$ by virtue of the Fokker-Planck equation. Therefore, only the second term contributes and we finally arrive at the so-called consistency conditions

[TABLE]

with the corresponding Langevin operator777 Note that this operator agrees with $L^{T}$ defined in Eq. (149) for holomorphic observables because of the Cauchy-Riemann equations.

[TABLE]

These conditions need to be fulfilled in order for the complex Langevin process to be correct. Technically, Eq. (162) represents an infinite tower of identities, since the condition would need to hold for all observables. However, it was argued that the check for a finite number of observables suffices in order to establish the validity of the results [221].

In fact, for practical purposes, it was realized that a slight re-definition of the interpolating function is more useful. Following Ref. [226] we introduce a cutoff $Y$ in the imaginary direction, such that Eq. (154) reads

[TABLE]

and the original interpolating function $F(t,\tau)$ is recovered in the limit $Y\to\infty$. Inserting this expression into the condition Eq. (161) finally yields

[TABLE]

where the temporal evolution of the observables $\tfrac{\partial}{\partial\tau}{\mathcal{O}}[\phi,\tau]=L_{0}O[\phi,\tau]$ has been inserted.888For the sake of simplicity, it was assumed that the integration along the real direction is benign and boundary terms are absent. In order to accommodate potential issues, a cutoff in the real direction may also be introduced, see, e.g., [227]. An explicit derivation may be found in the appendix of [228]. The form of this quantity also makes immediately apparent that a correct CL simulation requires a sufficient decay of the product $K_{\rm I}P{\mathcal{O}}$ in the imaginary direction.

Just like with regular observables, we must evaluate Eq. (3.5.1) in the limit $t\to\infty$, such that the relevant boundary term is given by $B_{{\mathcal{O}}}(Y)\equiv\lim_{t\to\infty}B_{{\mathcal{O}}}(Y,t,0)$. Note, that this evaluation itself does not require separate simulations at fixed values of the cutoff, but merely a certain post-processing of the sampled data. In practice, however, it may be challenging to assess the value of $B_{{\mathcal{O}}}(Y)$ for large cutoffs due to excessive noise. Therefore, it is convenient to investigate the boundary term as a function of $Y$, which typically allows one to identify a converged plateau region before the growing influence of the noise makes an extrapolation $Y\to\infty$ impossible. The value of the plateau may then be interpreted as the value of the boundary terms, which, if consistent with zero, marks the validtiy of the CL results. For non-zero values, the results suffer from incorrect convergence. Interestingly, the magnitude of the boundary term may be used to assess the systematic error of the CL values, as recently exploited in Refs. [229, 230].

Finally, the issue of correct convergence and boundary terms was very recently independently revisited in Ref. [231], where the understanding of these shortcomings was put on more solid mathematical footing.

3.6 Challenges and some solutions

Although our understanding of CL has progressed considerably during the past few years, the approach still faces some challenges that remain to be fully understood. These can be roughly divided into two kinds: mathematical and practical, which naturally have some overlap. In this subsection we attempt to summarize the current understanding of these issues.

3.6.1 Mathematical shortcomings

Without a doubt, the biggest challenge for CL is the lack of general mathematical proofs, the previous section notwithstanding. More specifically, as pointed out most recently in Ref. [223], it remains unknown whether the Langevin operators defined above as $L$, $L_{0}$, $L^{T}$, $L_{0}^{T}$ properly define unique stochastic processes, although in practice this is not typically an issue. Crucially, it remains unknown under what conditions the positive measure $P[\phi_{R},\phi_{I},t]$ converges to an equilibrium measure, although (again) there is substantial numerical evidence that such an equilibrium measure exists in many cases of interest.

3.6.2 Instabilities

One of the issues recognized early on (in fact, since the 1980s; see e.g. [173, 174]) is the appearance of instabilities in the form of runaway trajectories along the CL evolution. These can become frequent enough to completely spoil a calculation performed at fixed step size. Only with the availability of increased computational resources was it realized that too coarse an integration step would yield problematic trajectories. Finally, in Ref. [199], the need for adaptive step size integration of the complex Langevin equations was identified. It was found that such an approach provides a full solution to the problem of instabilities arising from the coarse integration step, by moderating the change between subsequently sampled configurations. In fact, this strategy has become standard in the field and virtually all modern applications of CL in both relativistic and non-relativistic physics implement the adaptive step (which is possible in several ways, see, e.g., [199]).

3.6.3 Convergence to incorrect values

The above issues aside, CL has been shown to fail in certain cases (due either to failure to converge or convergence to the wrong answer), but also appears to work in scenarios that lie outside the holomorphic-action regime mentioned above. In cases of failure, the behavior has been traced back to insufficient decay at infinity, or to a breakdown of ergodicity due to poles in the action (which should be expected in fermionic systems as the fermion determinant will vanish at specific points, thus leading to meromorphic drifts).

Boundaries at infinity

As outlined above, the behavior of boundaries at infinity is a relevant question, and a possible failure to meet this condition plagues models in both relativistic and non-relativistic physics. In particular, for gauge theories the complexification of the link variables leads to non-compact groups, e.g. SU(3) becomes SL(3,C). As we explain in Sec. 5, a similar effect is seen in nonrelativistic physics when using compact HS transformations. In either case, merely assuming that the derivative of $F(t,\tau)$ in Eq. (3.5) vanishes is a bad idea. For that to happen, the solutions to the FP equation must fall off sufficiently quickly along non-compact directions in the (complexified) space of field configurations (see in particular Refs. [232, 226] for a recent and insightful discussion of an exactly solvable case). That property is very difficult to determine a priori, but can be checked a posteriori following the arguments of Ref. [225]. Case studies show that in many instances, while the solutions fall off faster than exponentially in the real directions, the decay in the imaginary directions may be insufficient [221]. The prime source of the non-vanishing boundary terms are unphysical excursions far away from the real line, dubbed “the excursion problem.” Below we present several practical ways to mitigate such a problematic behavior.

Poles, ergodicity and singular drift

Besides the behavior near the boundary at large imaginary direction, another issue which could potentially spoil the formal justification of the CL method is the non-analyticity of the action. This is intimately connected to the occurrence of zeros of the measure, due to the relation

[TABLE]

such that the action diverges when $P[\phi]$ vanishes. Such potentially problematic points appear, for instance, for fermionic actions as zeroes of the determinant of the fermion matrix. Naively, one could argue that the problematic point is never sampled since it has zero probability measure. However, as with the boundaries at infinity, the success of a CL simulation crucially depends on the behavior of the probability distribution as the poles are approached. Formally, this can be elucidated by manually cutting out a region around the singularity, thereby extending the boundary of the integration domain. This procedure is justified as long as the probability measure vanishes around those poles sufficiently fast; in other words, one has to address the boundaries around the poles potentially alongside those at infinity mentioned above [222, 233, 234].

A detailed study of incorrect convergence due to poles in the drift function showed that the location of these poles, the decay of the drift function, and the behavior of the observables in the region near the poles all played a role in whether the method would return correct results [222, 225]. To be precise, there are three distinct scenarios, depending on the location of the poles with respect to the support of the probability measure:

*Poles on the outside – * These are unproblematic, as the influence of these singularities will be lost after thermalization such that the equilibrium distribution is not affected. Note that this is also the case for the toy problem discussed in Sec. 3.3.1 (visible as the red dot in Fig. 10, whereas the equilibrium distribution lies near the gray dashed line)

*Poles at the edge – * This type of singularity may potentially spoil the correctness due to a singular drift and therefore is a violation of the above derived rules for the decay of the histograms. Indeed, a careful investigation of several toy problems has elucidated such a behavior [222].

*Poles on the inside – * These are potentially problematic in two ways: As in the above case, a singular drift probem could appear, however, the additional complication is given by the formation of bottlenecks, which arise because of the vanishing probability distribution at these points [223]. As a consequence, the equilibrium distribution separates into disjoint patches and the Langevin process may fail to tunnel between them, resulting in a correct but restricted sampling of the configuration space and therefore a failure of ergodicity. If the relative weight of the avoided regions is small, the associated bias might be indiscernible and useful values may be obtained for practical purposes. Conversely, if the contributions of the avoided configurations may not be neglected, the CL process converges to erroneous values. Such a behavior is closely related to the concept of meta-stability in conventional Monte Carlo approaches.

For simple toy models, some of the problems connected to the occurence of zeroes of the measure can be inferred from the study of so-called classical flow-patterns, as also shown in the toy model studied above. For real-world applications, on the other hand, these are impossible to generate such that other diagnostic tools are necessary.

As mentioned above, a re-examination of the conditions for correctness in Refs. [235, 236, 225] revealed that a failure of CL in some cases has been attributed to the excursion and singular drift problems are actually due to the drift function falling off too slowly. This motivates the drift histogram as a diagnostic tool and it has indeed been shown to be a useful indicator for wrongfully converged values [237]. In fact, this type of analysis covers both types of boundary terms, the ones at infinity and around poles, in the same manner. It should be noted, however, that the observation of the drift histogram alone may be insufficient to distinguish the validity of multiple observables, as recently discussed in the context of the XY-model in Ref. [229].

Finally, it was argued in Ref. [147], using a semiclassical analysis, that when more than one saddle point in the complex plane contributes to the ensemble averages, the CL method can lead to incorrect answers due to the different complex phases associated with each saddle point. The interference of these complex phases is essential in phenomena such as the Silver Blaze phenomenon and real-time dynamics.

3.6.4 Practical solutions

The preceding discussion lays out the possible issues of a CL simulation with the tentative takeaway that certain theories seem to be intractable with this method. On the practical side, several practical strategies have been devised to enable CL simulations in otherwise problematic models and parameter regimes. These ideas try to either prevent uncontrolled excursions in non-compact complex direction which are ultimately responsible for the violation of the criteria for correctness, try to circumvent the issues of a diverging drift term, or a combination of both. In the following we provide an overview on the most successful strategies in this regard.

*Gauge cooling (GC) – * By now, this is a standard technique to prevent long excursions of the CL trajectory in the complex plane when simulating gauge theories [238, 239, 240, 241, 242, 243, 244, 245]. The underlying idea is that at each Langevin step, one can make a gauge transformation to keep the link variables, which live in the unbounded complexified manifold, close to the compact subgroup without interfering with the correctness of the simulation. For full QCD simulations, for instance, gauge transformations are performed in the extended non-compact group $SL(3,C)$ such that the sampled configurations are closed to the original compact $SU(3)$ manifold, measured by the unitary norm. Among all the strategies to mitigate the excursion problem, GC represents the best understood approach from a mathematical perspective and is, in fact, mathematically exact [242]. Finally, it is worthwhile to note that GC has been found to be necessary, however, sometimes not sufficient to bring CL simulations of gauge theories under control.

*Dynamic stabilization (DS) – * This approach was developed to further aid with the excursion problem [246, 247, 248, 249, 250]. The essential idea is to add a term to the Langevin drift $K[\phi]$ in the schematic form

[TABLE]

where $\alpha_{\rm DS}$ is a control parameter and $M$ acts only in the non-$SU(3)$ directions. The latter is chosen such that it grows rapidly with the distance from the $SU(3)$ manifold, thereby suppressing contributions with large deviations from the target manifold. Although this modification to the Langevin evolution cannot be derived from the action itself, it can be shown to vanish in the continuum limit.

*Regulators – * The idea behind a regulator, similar to the above presented strategies, is to suppress terms far from the original configuration space. The approach is similar in nature to DS and first appeared in a nonrelativistic application, namely Ref. [251], and was further discussed more recently in Ref. [226].

As opposed to DS, the resulting modifications on the Langevin equations are directly derived from the action, which is modified by a term of the form

[TABLE]

which effectively adds a harmonic oscillator trapping potential, i.e. a restoring force that prevents the field from running away.999In a relativistic setting, the modification of the action may be understood as a a mass-like term for the auxiliary field $\phi$. The modification, strictly speaking, amounts to solving a different theory and therefore introduces a systematic bias in the computed expectation values. To remove this bias, it is necessary to either show that the results are insensitive to the regulator strength $\xi$ or, equivalently, extrapolate $\xi\to 0$. Stated differently, the idea is to simulate a series of well-behaved systems with decreasing deviation from the physical system of interest in order to extrapolate to the physically correct point.

While the advantages of such a practical solution are clear, it is by no means a full solution and in many cases – especially at strong coupling or low temperatures – it is not possible to make $\xi$ small and obtain a converging calculation. Nevertheless, it has been found that this strategy is essential for simulations of non-relativistic fermions via CL and it is by now a standard approach in this context (see Sec. 5.2.2 for a detailed discussion).

*Deformation technique – * In order to mitigate the occurrence of singular drift values, a viable option is to deform the action by adding a fermion bilinear term to the action [252, 253, 254], ideally in such a way that the vacuum of the system is minimally affected [252]. For a suitable choice of this term, the effect is to suppress (near-)zero eigenvalues of the fermion matrix (in this context also often referred to as the Dirac operator) which are ultimately responsible for exceedingly large drift values [255, 256]. In analogy to the regulator prescription discussed above, it is necessary to extrapolate the coefficient of this additional coupling to zero to obtain results for the physical theory of interest. For gauge theories, this strategy is typically used in combination with GC in order to achieve stable and correct CL simulations, most recently for a full QCD simulation [254]. For non-relativistic theories, this approach remains unexplored to the best of our knowledge.

3.6.5 Schwinger-Dyson representation

Quantum field theories can be cast in terms of differential equations rather than path integrals; this is called the Schwinger action principle. A specific set of differential equations resulting from this principle is the Schwinger-Dyson equations, and these involve variations with respect to the source fields. Of the many solutions to this set of differential equations, only one at most corresponds to the Feynman path integral. It appears that the CL method is able to converge to multiple stationary distributions, and while at most one of those distributions corresponds to the desired distribution, it has been shown that these stationary distributions of the complex Langevin equation satisfy the Schwinger-Dyson equations [193, 257].

The initial conditions appear to play a role in this behavior: CL converges to certain solutions of the Schwinger-Dyson equation with one set of initial conditions and to another with a different set, and for certain conditions CL does not converge at all  [258, 259]. The role of the Schwinger-Dyson equations in the presence of zeros of the complex density $\rho$ was analyzed in detail in Ref. [260]. Further study of the convergence of CL to incorrect solutions must be an integral part of the understanding and application of this method to systems with complex actions.

4 Applications in relativistic physics: from toy models to finite density QCD

As already mentioned above, the first applications of CL were undertaken in the context of lattice gauge theories in the 1980s. Following several successful applications in this regard [261, 262, 174, 181] the initial wave of excitement simmered down due to observation of incorrect convergence and a poor understanding of the causes. After a period of reduced activity, the CL method resurfaced in the 2000s in the context of relativistic physics (see Refs. [192, 193, 194, 196, 197, 263]). Motivated by the revived interest, much work has gone into the testing of the methodology, its capabilities, and its limitations. While CL has successfully circumvented the sign problem in a few key areas, it still suffers from some of the problems discussed in Sec. 3. As a result, recent emphasis has been on determining regions of applicability and adjusting the method in order to prevent uncontrolled excursions of the Langevin evolution into the complex plane and ensure convergence to the correct solution.

While the main focus of this review is on applications of CL to nonrelativistic quantum many-body problems, such an overview would be incomplete without accounting for the developments that enabled these advances. Therefore, we provide a brief summary of selected applications of CL to relativistic field theories. For a more in-depth overview on the use of CL in relativistic physics, we refer the reader to Refs. [263, 264, 223, 265] for various status reports.

4.1 QCD-inspired toy models

In Lattice QCD, the main promise of CL is its potential to explore regions of the QCD phase diagram which are currently inaccessible due to a sign problem. One of the primary goals is to reliably treat the region of non-zero quark chemical potential (see, e.g., Refs. [266, 267, 268, 269]), which eludes efficient treatment with conventional methods due to a severe complex phase problem. In addition, theories with coupling to a topological charge (see Ref. [223]) or computation in Minkowski space to obtain real-time dynamics [192, 193, 198, 270] are also of great interest.

While the CL method has not yet been able to produce detailed solutions for these problems (some attempts notwithstanding, see below), a number of simpler models, which partly contain the phenomenology of QCD in certain limits, have been successfully studied. The findings in these reduced models allow deep insights into the behavior of the CL method in these limits and can help elucidate where problems may arise with such a treatment. In the following we intend to highlight some of these advances.

4.1.1 XY model

The XY model is an extremely useful benchmark case for CL as the sign problem can be circumvented in other ways, for example by using imaginary asymmetries or dual variables in a world-line formulation.. Results from both of these alternate approaches have been compared to calculations done with CL for the three-dimensional XY model at non-zero chemical potential [271, 200]. The CL results turn out to be very promising for the ordered phase at low temperatures but fail to reproduce known answers for the disordered (high-temperature) phase. From this failure, first criteria for the correctness of CL results have been extracted, as already discussed above. Similar conclusions on the applicability of the CL method have been reached in a recent study for the $O(3)$ model [272].

Moreover, the model is very sensitive to instabilities in the algorithm – just as heavy-dense QCD is, see also below. The instabilities in the algorithm can be eliminated by the addition of adaptive stepsize to the Langevin evolution, which was first shown to work in the three-dimensional XY model at non-zero chemical potential. As remarked above, the adaptive step size mitigates the issue of unstable CL trajectories by preventing overly large changes in the Markov process. Another way to improve the convergence of CL in the XY model is by dynamic stabilization, which forces the Langevin trajectory to remain near the real axis by means of unitary transformations [248]. These transformations modify the drift function to prevent large excursions in the imaginary direction [223, 247].

4.1.2 Polyakov and SU(3) spin models

Polyakov models at non-zero density were among the first theories under investigation with CL, and their exploration has led to methodological breakthroughs. Most notably, it was realized that the complexification of the group elements led to an enlarged manifold that needs to be constrained via gauge cooling in order to prevent unstable trajectories [238]. This strategy is key to studying more sophisticated theories.

The three-dimensional SU(3) spin model is an effective Polyakov loop model for QCD at nonzero temperature and density. It suffers from a sign problem at nonzero chemical potential and typically reweighting or the phase-quenched approximation are used to study its dynamics. Complex Langevin was used to successfully address this model, originally in Refs. [261, 273] and more recently in Refs. [220, 274]. The latter studies found remarkable accuracy across a phase transition. Some of these results are shown in Fig. 11 along with the average phase factor in the phase-quenched theory, which vanishes quickly with increasing system size. The success of CL within this model, despite the smallness of the average phase, underlines its ability to solve severe sign problems. Moreover, a comparison with the 3D XY model, for which CL fails for such a phase transition, sheds light on the justification of the approach and sets the stage for the derivation of further criteria of correctness.

Later, CL results for the SU(3) spin-model were also compared with the relative weights method across a wide range of chemical potentials [275, 276]. While values for the densities agreed perfectly in all cases, results for the Polyakov loop displayed some discrepancies at large chemical potentials. The shortcoming that was brought in connection with the appearance of a logarithmic term (and hence a gauge cooling) in the action, as already encountered in [277].

branch cut

4.1.3 Random matrix theory (RMT)

Random matrix theory has been proven generally useful for the analysis of QCD [278]. Perhaps most importantly from a CL perspective, just like full QCD, RMT models are “constructed” with complex fermion determinants and both exhibit a severe sign problem for small masses. As opposed to QCD, however, analytic solutions for RMT models are available, making them a valuable playground to assess the ability of CL to deal with complex-valued fermion determinants.

Some work has been done for models without a finite-density transition [279], where it was found that CL works well for large quark masses but incorrect convergence occurs for smaller values of the mass. The failure was linked to CL paths that wind around the origin of the complex plane, therefore encountering the ambiguity in the logarithmic part of the action, i.e., a branch cut. The latter criterion was initially believed to be a general marker of incorrectly converged CL results and indeed had some success in explaining the failure encountered in other models. However, a later study found that CL also fails to deliver reliable results for cases where such a point is on the edge of the domain of the distribution function where no full winding around the pathological point is possible [222]. In a follow-up study, CL was shown to be able to fully capture the chiral limit of RMT which was achieved by shifting the integration variables and changing to a polar representation [277].

In addition, an extended gauge cooling procedure was proposed and shown to deliver exact results in the regime of light quark masses [243, 280]. Here, the general idea is to avoid small eigenvalues of the Dirac operator which could possibly spoil the correctnes of CL [255, 256].

Moreover, models with finite-density phase transition, where many problems with CL may arise, has been addressed. A recent work on one of these models suggests that with suitably designed reweighting methods, CL is able to fully reproduce the known analytical solution of the model at the finite-density phase transition [281].

4.1.4 Thirring model

The Thirring model model constitutes one of the rare cases of an exactly solvable quantum field theory. CL studies of this model have been conducted in the $0+1$ dimensional case where it was found that CL is able to reproduce the exact results in the weak coupling regime [282]. At larger couplings the method is only able to reproduce exact results at small and large values of the chemical potential but fails for intermediate $\mu$. Similar findings have been reported for the $2+1$ dimensional scenario [283]. A later comparison (for both dimensions) with the fermion-bag method (see Sec. 2.5), underscored this behavior and showed that CL does not perform sufficiently well in this parameter regime [284, 285].

In an independent study of the $0+1$ dimensional Thirring model, essentially the same behavior was observed and it was found that the CL trajectories were sampled in the vicinity of zeroes of the fermion determinant [286]. Consequently, large values of the drift occur in such studies and the distributions of the drift term are no longer localized, thereby violating certain criteria for correctness. A solution in terms of reweighting and deformation was achieved. However, the scaling of these “remedies” is exponenphase diagramtial in system size.

4.1.5 Heavy-dense QCD (HDQCD) and hopping parameter expansion

A very prominent simplification of QCD is the limit of large quark masses and large densities, called heavy-dense QCD (HDQCD). In this limit, the kinetic term of the quarks is suppressed to the point that they become essentially static. In this sense, the limit may also be used as a starting point for an expansion in the hopping parameter, which approaches the dynamic limit with increasing order, see below.

Studies of HDQCD have been among the first objectives since CL was revived [194] and have led to profound insights and methodological advances, most notably in the light of adaptive step-size control [199] and the gauge cooling procedure [238]. These advances made it possible to map out the phase diagram in the full $T-\mu$ plane at several values of the inverse coupling $\beta$, and excellent agreement with known results was found in the large-$\beta$ regime [287] (this includes, for instance, the Silver Blaze behavior that was observed to sharpen with decreasing temperature). More recent studies have refined these findings [288, 289, 290, 291, 292]. However, it was also found that for $\beta\lesssim 5.9$, instabilities emerged where the histograms of the drift term developed skirts and sometimes even shifted the mean. These shortcomings were accompanied by an increasing unitary norm (i.e., the trajectories were far away from the $SU(3)$ submanifold), despite the use of gauge cooling. Close monitoring was found sufficient to get reliable ab initio insights in the entire $T-\mu$ plane, albeit with reduced precision in different parameter regimes. With advances in the assessment of systematic errors of CL, the failure was linked to boundary terms which significantly increase as $\beta$ is decreased [230]. In fact, such boundary terms are also present in the “correct” parameter regime, but they are sufficiently small such that the estimated systematic error is well below the statistical accuracy.

In addition to the phase diagrams in the heavy limit, CL studies have also been performed in combination with a hopping parameter expansion. It has been mentioned [293] that the effective theories resulting from such an expansion exhibit a much milder sign problem, sometimes allowing a comparison with reweighting. Low-order expansions have been used in Refs. [295, 296], where excellent agreement between Metropolis-based algorithms and CL has been observed, and it was reported that CL is vastly superior at larger lattices. Moreover, all-order hopping expansions have been simulated with CL [297, 293, 298] and were found to agree well with CL results for the full theory [294] (see below), which is considered to be a validation. Results for the density are shown in Fig. 12. It is worthwhile to note that two different expansions – one that includes poles in the measure and one that does not – agree, indicating the absence of the influence of poles in the drift term in this case, providing us with an optimistic outlook for the application of CL to full QCD.

4.1.6 Low-dimensional QCD at non-zero chemical potential

Another strategy to benchmark the applicability of CL in the context of QCD is to reduce the spatial dimensionality of the lattice, thereby considerably reducing the numerical effort. Applications include $0+1$ dimensional QCD, where the quark density and chiral condensate have been computed using CL, with gauge cooling used to prevent the simulation from deviating in the imaginary direction [299]. Simulations in $1+1$-dimensional QCD have also successfully provided us with results for the chiral condensate in the thermodynamic limit [300]. In recent years, work on combining reweighting techniques with CL have expanded the range of applications of the method to previously inaccessible parameter space in QCD, allowing for greater ranges of values for the mass, temperature, and chemical potential in $1+1$-dimensional QCD [301, 302, 303].

4.2 $3+1$ dimensional QCD at finite chemical potential

As mentioned already above, one of the big goals in the development of CL is the application to full QCD on the lattice at finite chemical potential, i.e., the consideration of dynamical quarks with physical masses. The earliest attempts to calculate the properties of full QCD at finite chemical potential is the work of Refs. [304, 305, 306, 294] which employed a staggered-fermion formulation. As shown in Fig. 12, the CL approach using gauge cooling agrees with the HQCD results in the valid region, i.e., for large quark masses. For values of $\beta\gtrsim 5.9$ the simulations are well under control and reliable results may be obtained. Similar to previous model studies, however, smaller values of $\beta$ exhibit large unitary norms despite gauge cooling and, as a consequence, develop skirts in the observable distributions. Another early application to the full theory was performed in Ref. [307], where an extensive comparison with data from studies based on reweighting techniques at $\mu=0$ was performed. Both methods agree in a wide parameter range, however, the conclusion was similar to the parallel studies of QCD, namely that lower values of $\beta$ remain problematic.

Later on, the results of Ref. [308] suggested CL with gauge cooling produces trustworthy solutions for $3+1$-dimensional QCD at zero temperature and finite chemical potential. Work on this same system using a combination of adaptive step size and gauge cooling to prevent runaways in the complex plane for weak coupling showed promise. In fact, the results share some important physical traits with the system. However, known results were not accurately recovered and CL results seem to be influenced by the phase-quenched theory. To be more specific, the observation of the transition from hadronic to nuclear matter (at $\mu\approx m_{N}/3$) is visible, as is evidence of saturation at large chemical potential [309]. The zero chemical potential limit disagreed with known results, but to a very small degree. With larger lattices and smaller coupling, the known results are still not accurately recovered. However, it was argued that the agreement becomes better in these cases [309]. This has been traced back to the fact that these calculations suffer from the appearance of zeroes in the fermion determinant in this regime [310] (see also [311, 312]). Very recent extensions of these investigations employed modified actions with additional irrelevant terms (for the continuum theory) that are designed to separate the physics between different mass scales [313]. This strategy seemed to improve the location of the phase transition, however, further computations, in particular towards the continuum limit, are required.

A parallel development started with investigations of full QCD at small $4^{3}\times 8$ lattices [314] where encouraging results have been obtained with gauge cooling alone. While problems associated with a singular drift have not been observed at these lattice sizes, somewhat increased autocorrelation times have been argued to be the harbinger of potential issues in that regard on larger lattices. A subsequent study on the same lattice size but with the deformation technique yielded well-controlled results and was able to extend the applicable parameter range [254, 315]. In fact, it was found that the deformation is a necessary tool to circumvent the increasingly problematic influence of the singular drift and it was argued that when $\mu$ is chosen within the region of validity of CL, the eigenvalue distribution of the Dirac operator exhibits a gap along the real axis [253]. Most strikingly, a delayed onset of the influence of the chemical potential on the quark density was observed, a signal that was argued to be in agreement with the Silver Blaze property. Later, the results were extended to larger lattices [316], where the formation of a plateau in the quark density has been found. This was interpreted as the onset of nuclear matter (i.e., the formation of a Fermi surface) which was observed to persist even for lattices up to $16^{3}\times 32$, which to date is the largest calculation of full QCD with CL [317, 318].

Although progress has been made in the exploration of full QCD in $3+1$ dimensions, many aspects of the role of CL in this regard remain to be investigated. One of the pressing questions is the general applicability of CL in the confined phase, which was argued to be questionable at least due to the structure of the Dirac operator [319, 320]. Note that this was suspected, on more practical grounds, already in Ref. [307]. In addition, the issue of erroneous convergence and the absence of a systematic error estimate constitutes a potential roadblock. To this end, recent advances in the analysis of boundary terms may constitute a way to handle this situation and, in fact, have also been applied to full QCD recently [229, 230]. Finally, it remains to be seen if the use of certain improved operators, a strategy that was recently employed in a precise calculation of the equation of state in the quark gluon plasma [321], will help to circumvent some of these issues. A very recent overview on the current status of full QCD simulations with CL can be found in Ref. [265].

4.3 Other applications in relativistic field-theories

Although most of the applications and development of CL have been achieved in the context of lattice QCD and related models, several other attempts have been made in the past. Here, we only provide a short overview of such studies.

4.3.1 Real-time dynamics

As mentioned above, a very intriguing application of the CL method is the direct simulation of QFTs in Minkowski spacetime, which suffers from a so-called “dynamical” sign-problem, originating from the complex exponent in the path integral (and typically avoided by going to the imaginary-time formulation). This potential was realized immediately after the CL prescription was originally proposed [322] which led to several applications in the context of high-energy physics and related models [322, 323, 324, 325, 326] (with early applications even in the field of gravity [327]), see Ref. [168] for an overview on the formalism and an early review.

After a period of relatively low activity, stochastic quantization once again was employed to study non-equilibrium QFT in the context of non-relativistic physics in the mid 2000s. Due to the non-perturbative nature of these non-equilibrium systems, standard approximation techniques fail. heavy ionulations may potentially be applied here but a Minkowski formulation suffers from a sign problem. In 2005, Berges and Stamatescu demonstrated the viability of CL to treat non-equilibrium QFT using first-principles simulations [192, 198], albeit, only close to the equilibrium case as the long-time limit is prone to erroneous convergence. Later progress built on these results to examine how CL could lead to breakthroughs in our understanding of QCD plasmas in heavy ion collisions, early thermalization, and other open questions in quantum field theory [193].

More recently, CL in combination with non-Markovian Langevin equations was used to address non-equilibrium physics for simple non-interacting field theories [270], leading to first sucesses and a potential route to future applications for interacting theories. Moreover, non-equilibrium aspects of scalar $0+1$ dimensional field theories (or, equivalently, quantum mechanics) have been studied within the CL algorithm and found promising results but also potentially troublesome, unphysical fixed points in the Langevin propagation [328]. To deal with the latter issue, a restriction in phase space was introduced, yielding correctly converged, yet potentially approximate, results. Finally, Ref. [328] also contains an elucidating comparison of real-time stochastic quantization with several state-of-the-art methods for the determination of dynamical properties.

4.3.2 Theories with topological terms

To (potentially) address the question of the strong CP problem, theories need to contain a so-called $\theta$-term, generically written as $S_{\theta}=-{\rm i}\theta Q$ with the topological charge $Q$. Such a term is purely imaginary, thereby introducing a phase problem. Besides results obtained by reweighting [329], a few CL studies have been conducted to examine the effect of such terms.

In Ref. [239, 330] the $SU(3)$ Yang-Mills theory was investigated, for the first time directly with a purely real $\theta$ - for imaginary $\theta$ the sign problem is absent and HMC may therefore be used as a benchmark in this regime. It was found that gauge cooling is essential for an efficient implementation and that CL is well controlled above an inverse coupling of $\beta\gtrsim 5.9$, below which the conditions for correctness were observed to be violated.

Additionally, the two-dimensional $U(1)$ gauge theory in the presence of a $\theta$-term was investigated [331]. It was shown that a naive implementation fails because of the crossing between different “topology sectors” (which implies the crossing of a branch cut). Reassuringly, the failure was observed to be accompanied by the violation of correctness markers. By removing a plaquette from the torus, i.e. “piercing” the manifold, the CL model was observed to deliver reliable results even at large values of $\theta$, previously inaccessible with methods based on reweighting.

4.3.3 Matrix models for superstring theory

While the overwhelming application of CL in the context of relativistic field theories has been in the context of QCD and related simplified models, several studies have explored quantum field theories beyond this scope. Recent progress was specifically made in the realm of matrix models, which are conjectured to provide a non-perturbative formulation of certain superstring theories [332] (in a similar fashion as lattice gauge theories constitute a numerically accessible non-perturbative formulation of QCD).

Amongst others [333, 334], particular interest was paid to the type IIB matrix model, which is conjectured to be a nonperturbative definition of type IIB superstring theory [335]. In its Euclidean version, a severe sign problem arises upon integrating out the fermions, leaving a complex-valued Pfaffian. The complex-valued measure was found to be responsible for the central phenomenon in these theories, namely a dynamical compactification of spacetime from ten to four dimensions, which rules out phase-quenched approaches to investigate these models [336]. Initial CL studies focused on the low-dimensional sector [252] and found that a singular drift problem arises which is treatable with the deformation technique. Subsequently, the study of this model was extended to six dimensions [337, 338], where the a breakdown from $SO(6)$ to $SO(3)$ was observed, marking a major improvement over reweighting-based methods. Most recently, these investigations have been extended to the original ten dimensions of the model in the Lorentzian [335] as well as the Euclidean version [339, 340] of the model, where clear signals of the spontaneous symmetry breaking have been observed.

5 Applications in non-relativistic matter: ultracold atomic gases

In contrast to relativistic theories, applications of CL to non-relativistic matter remain largely in their infancy (notwithstanding notable early work such as Ref. [182]), perhaps not least due to an early finding that conventional MC approaches are superior to CL in this regard because of uncontrolled runaway trajectories [341]. Recent advances concerning these issues made largely in the high-energy community, however, induced progress also for non-relativistic models. As described below, there have been several attempts to characterize bosons as well as fermions in a variety of situations [342], but as of this writing only one calculation has been successfully carried out for a strongly coupled fermionic system in 3D [343] (see the discussion in Sec. 5.3). Non-relativistic systems have difficulties of their own in the form of phase transitions and strongly coupled regimes, but they do not feature quintessential (and technically challenging) QCD elements such as non-abelian gauge fields. Nevertheless, as we shall see, some of the challenges faced by CL calculations are universal, as are the ideas to diagnose and tackle them. The different viewpoints taken in relativistic and non-relativistic applications leads to an important interplay between these two subfields and ultimately progresses the exploration of CL in both communities (see Sec. 3).

Among the vast array of non-relativistic systems that remain to be explored, some pressing candidates remain across different areas. In the condensed matter context, for example, the leading candidate is undoubtedly the repulsive 2D Hubbard model which faces a severe sign problem upon introducing a finite doping; a single exploratory study was conducted for the repulsive 2D Hubbard-model at half-filling [344]). Other highly interesting candidates include the spin-polarized electron gas, spin-isospin polarized nuclei, and neutron and nuclear matter, to name only a few. The list of applications outlined in this section may potentially also lead to further advances in those areas.

5.1 Non-relativistic bosons

One of the pioneering applications of CL to non-relativistic systems in the modern era (roughly within the last decade) concerned the study of a 3+1 dimensional bosonic quantum field theory in a rotating frame [344, 345]. The number density and condensate fraction with no rotation computed via CL was found to agree with mean-field calculations, and when rotation was introduced, the CL results showed quantized circulation for high condensate fraction. This is consistent with vortex formation in rotating superfluids, a phenomenon which has been directly observed using rotating ultracold bosonic atoms. The results for the circulation using CL do not agree with mean-field calculations for small condensate fraction, as expected due to the breakdown of mean-field theory to describe a system with strong quantum fluctuations. Besides this study in three spatial dimensions, an exploration in two spatial dimensions, where quantum fluctuations exhibit an even greater impact, was conducted [346]. However, in the latter study, the parameter regime where a quantized vorticity is expected remained out of reach, requiring further CL studies on the topic.

In addition to these single component systems, the effect of spin-orbit coupling in a two-dimensional Bose gas with two pseudospin components interacting via a delta function potential was explored [347]. In this case, the sign problem arises due to a first-order derivative in the imaginary time which was well controlled with the CL approach. Mean-field predictions were found to overestimate thermodynamic equations of state, in particular at large spin-orbit coupling.

Aside from these brief explorations, the application of CL to circumvent the sign problem in systems of non-relativistic bosons remains scarce.

5.2 Non-relativistic fermions in one dimension

One-dimensional (1D) Fermi gases have drawn considerable interest in the past, when exact solutions for various models have been derived by the use of the so-called Bethe ansatz (see e.g., Ref. [349] for an extensive review). The availability of these exact solutions, paired with a relatively modest computational cost, make these systems excellent benchmark scenarios for any newly devised method, such as the CL approach for non-relativistic quantum matter. Moreover, 1D systems have become accessible in experiments in recent years, which provides yet another motivation to study these exotic and intrinsically strongly-interacting systems.

It is worthwhile to note here that in the special case of 1D it is often possible to re-write relativistic and non-relativistic models using dual variables in a way that avoids sign problem. One would be able to then compute quantities by conventional Monte Carlo methods. While this is an option in 1D, those approaches do not generalize to higher dimensions, in contrast to the CL method. In fact, the dimensionality of the problem is mainly a question of computational effort, and all insights gained on the numerical behavior of CL simulations carry over to higher dimensions.

In this section, we will review recent advances that have been made in applying CL to one-dimensional fermionic systems that suffer from a sign problem. Specifically, we will discuss systems at zero temperature featuring attractive and repulsive interactions, finite spin-polarization as well as asymmetry in the masses of the fermions in Sec. 5.2.1. Furthermore, we recall results obtained at finite temperature in Sec. 5.2.2, where repulsively interacting fermions are shown as well as systems at asymmetric chemical potential.

5.2.1 1D fermions in the ground state

In the following, we will consider a system of $N=N_{\uparrow}+N_{\downarrow}$ fermions at fixed linear lattices of $N_{x}$ sites, which corresponds to working in the canonical ensemble. The physics in 1D is set by the dimensionless interaction parameter $\gamma=g/n$ where $g$ is the bare interaction and $n$ is the linear particle density.

Balanced systems

As a first benchmark of the CL method in the ground state, the spin- and mass-balanced Fermi gas was investigated [118, 351]. A comparison to various other methods is shown in Fig. 13, where generally good agreement is apparent across the entire range of interactions shown. While the CL results shown are at $N=5+5$ particles and a finite lattice of $N_{x}=40$, the curves for the Bethe ansatz expansions (strong [352] and weak [353] coupling) correspond to the thermodynamic limit (TL) of taking particle number and volume to infinity at constant density. The agreement on the repulsive sides indicates sufficiently large box sizes. On the attractive side agreement is observed; however, larger spatial lattices are needed as is suggested by the slight discrepancy between CL and the volume-extrapolated HMC values of Ref. [354], as well as results obtained by a DFT-RG approach in Ref. [348]. Another possible source of systematic bias is the finite (adaptive) integration step of $\Delta t=0.01$, which was used throughout this study. Although the $\Delta t$ dependence was checked in Ref. [118], where it was found that $\Delta t$ was sufficiently small (within the statistical uncertainty), the influence of $\Delta t$ could vary in other areas of parameter space, in particular with varying coupling strength.

While the agreement among the methods looks excellent, a closer look reveals interesting technical subtleties. During stochastic calculations it is instrumental to study histograms of all measured quantities to ensure correct behavior. As can be appreciated in the outer panels of Fig. 13, two distinct behaviors were observed in different regimes: the strongly attractive case ($\gamma=-3.0$) displays well-behaved, localized distributions; the strongly repulsive case ($\gamma=3.0$), on the other hand, exhibits a large amount of outliers with respect to an assumed Gaussian. The existence of so-called “fat tails” renders the simulation problematic, as this violates the conditions for correct behavior and thus undermines the validity of the CL approach. Therefore, it was conjectured that the CL values are possibly faulty in this parameter regime, which was later confirmed by a different few-body method [350] based on dual variables and the worm algorithm. A more precise understanding of the behavior at strong repulsive interaction, and eventually its resolution, would be very useful.

Mass-imbalanced Fermi-Fermi mixtures

Mixtures of two different fermion species, i.e., systems where the constituents have unequal masses, are challenging to address theoretically. While many 1D models are integrable via the Bethe ansatz, there is no such solution in the mass-imbalanced case. Nevertheless, these systems are of great interest as there are a number of experimental realizations available for these mixtures. The CL method can be straightforwardly extended to cover fermion systems with unequal masses by the use of spin-dependent dispersion relations, i.e. a spin-dependent mass $m_{\sigma}$ (or, equivalently, anisotropic hopping parameters). To quantify the mass asymmetry between two species, it is convenient to introduce the dimensionless relative mass imbalance $\bar{m}={(m_{\uparrow}-m_{\downarrow})/(m_{\uparrow}+m_{\downarrow})}$, which is restricted to the interval $[0,1]$. While the use of spin-dependent dispersion relations is unproblematic from a conceptual viewpoint, numerical difficulties might occur due to an increasing separation of scales. Thus, one may use different parametrizations of the mass imbalance depending on the specific implementation of the algorithm.

In Refs. [118, 351] the mass-imbalanced few-body problem was studied and the average mass was fixed to $m_{0}=1$. In this way, a comparison is made possible to results obtained with the iHMC method (imaginary asymmetry – introduced in Sec. 2.7). A comparison is shown in the left panel of Fig. 14 for fixed values of the interaction strength $\gamma$. Excellent agreement of CL and iHMC is reported across a wide range of mass imbalances. Starting at $\bar{m}\sim 0.6$ the iHMC results start to deviate as a consequence of the instability of the analytic continuation. The CL results, on the other hand, remain smooth and precise up to high mass imbalances. The agreement between these different methods provides confidence that the CL method is suitable to study this problem.

Recently, a worldline algorithm was adapted to study mass-imbalanced few-body systems without a sign problem [350, 355] in 1D and its results have been compared to the CL values obtained in [118] (shown in the right panel of Fig. 14). On the attractive side, excellent agreement is observed for all mass imbalances, which further validates the CL method in this setting. A comparison on the repulsive side, however, revealed that the CL values deviate from the wordline results, which is connected to the occurrence of heavy tails mentioned above and in Ref. [118]. It remains to be shown whether mass-imbalanced Fermi mixtures with repulsive interactions in the ground state are accessible with CL. This question is also intimately related to the applicability to the repulsive Hubbard model.

Spin- and mass-imbalanced systems

Very recently, the investigation of mass-imbalanced fermionic systems was extended to study pairing phenomena based on two-body correlation functions [356] which, to date, is the first determination of correlation functions with the CL approach. Concretely, the effect of unequal masses on the behavior of a spin-polarized system of two-component fermions was addressed. While it is known (theoretically) that such a system exhibits Cooper pairing at finite pair-momentum (referred to as Fulde-Ferell-Larkin-Ovchinnikov or FFLO pairing), the effect of mass imbalance in this situation remained relatively little explored. Notably, no exact solutions for these mass-imbalanced system are known to date whereas exact solutions from the Bethe ansatz exist for the mass-balanced case, even for finite spin imbalance [352], see Ref. [349] for a review on the Bethe ansatz.

With CL, it was shown that the appropriate correlation functions can be extracted accurately. Interestingly, in addition to the expected FFLO-type pairing correlations, unprecedented correlations in the density-density correlation function in momentum space (also referred to as shot-noise correlations) have been observed in case of finite mass imbalance, see Fig. 15. The extracted results for the mass-balanced case (leftmost panel in Fig. 15) agree with related DMRG studies of the spin-imbalanced Hubbard model at low filling [357]. Specifically, dominant correlations are observed at the points $(\pm k_{{\rm F},\uparrow},\mp k_{{\rm F},\downarrow})$ which indicates the expected pairing at the opposite Fermi momenta. With increasing mass imbalance, on the other hand, an additional correlation maximum emerges, suggesting a previously unknown scattering of fermions far below (heavy majority) or far above (heavy minority) the Fermi point of the heavier spin species. It is worthwhile to note that the underlying mechanism is still associated with pairing of fermions with a pair-momentum also assumed in case of conventional FFLO correlations, i.e. $q=|k_{{\rm F},\uparrow}-k_{{\rm F},\downarrow}|$. However, the relative momentum of the bound fermions is shifted.

The two-body quantities obtained from CL provide an exceptionally clean signal of the pairing instability. Moreover, the shot noise is of direct experimental interest, as it may be extracted with readily available time-of-flight measurements.

5.2.2 1D fermions at finite temperature

One of the first attempts to use CL for non-relativistic fermions was made in Refs. [251, 359], where the thermodynamics of a repulsive101010In fermionic systems with a repulsive coupling, the sign problem is present even for spin-balanced systems, as the Hubbard-Stratonovich transformation necessarily requires the determinants (and hence the probability measure) to be complex quantities. 1D Fermi gas was calculated using that method as well as third-order lattice perturbation theory (see Fig. 16). The density and pressure were computed for a range of interaction strengths and temperatures. The calculations were also compared with prior hybrid Monte Carlo results for attractive interactions, where the sign problem is absent and the convergence of the perturbative expansion can be assessed. On the repulsive side, CL was found to agree with the third-order perturbative expansion at weak coupling and away from the virial region (where the fugacity is small), beyond which perturbation theory was expected to break down. In the virial regime, a related study employed particle projection in combination with CL to obtain high-order virial coefficients in 1D [360]. This approach, while so far only applied in 1D, is applicable in any spatial dimension and could therefore be beneficial in extending the range of accessible temperatures of the virial expansion.

At the same time, the aforementioned work [251] explored for the first time the use of a regulator to avoid uncontrolled excursions of the auxiliary field into the complex plane. In those calculations, the excursions into the complex plane were highly problematic because the dependence of the action and the drift on the auxiliary field $\sigma$ involved hyperbolic functions. The HS transformation used in that work depends on $\sin\sigma$, and the complexification $\sigma\to\sigma_{R}+i\sigma_{I}$ results in

[TABLE]

As a consequence, a growing magnitude of $\sigma_{I}$ amounted to increasing the fermion-fermion coupling at an exponential rate. In such a situation, the calculation would either stall or result in a converged but incorrect answer. This exponential growth is similar to the problem found in gauge theories, as the complexified link variables representing the gauge field become unbounded in the same fashion in those theories. For lattice calculations of gauge theories, Refs. [247, 248] tackled the problem of large excursions by dynamic stabilization (see Sec. 3.6). Following a similar idea, Ref. [251] added a regulating term to the CL dynamics controlled by a parameter $\xi$, such that the new CL equations became

[TABLE]

As remarked above, the new term in the action can be understood as a harmonic oscillator trapping potential, i.e. a restoring force that prevents the field from wandering too far in the imaginary direction. That modification introduces a systematic effect that needs to be studied for each quantity of interest as a function of $\xi$ as $\xi\to 0$. Fig. 17 shows the running average of the density as a function of the Langevin time $t$ for several values of $\xi$ in the neighborhood of [math]. As is evident in that figure, there is a sizable window of small values of $\xi$ where CL converges.

Spin-polarized systems

Encouraged by the success of CL in calculating 1D systems with repulsive interactions, the follow-up study of Refs. [361, 359] extended the results of Ref. [251] to the spin-polarized case by introducing a non-vanishing asymmetry $\beta h=\beta(\mu_{\uparrow}-\mu_{\downarrow})/2$. Such an asymmetry leads to a sign problem for the case of attractive interactions and to a phase problem for repulsive interactions. As mentioned above, such a sign or phase problem can be avoided in 1D using specific methods which, unfortunately, do not generalize to higher dimensions.

In the study of Ref. [361], the density and magnetization equations of state were calculated with CL and compared with perturbation theory, imaginary polarization approaches (iHMC), and the virial expansion. A sample of those results is reproduced in Fig. 18 where the remarkable agreement among all of these approaches is apparent. However, at the largest chemical potential asymmetry (i.e. $\beta h=2.0$ in the figure), the results obtained with imaginary asymmetry start to deviate from the ones obtained by other means, including CL. This trend is not surprising, as the iHMC method is expected to be less robust in this limit (see also the discussion above). Reassuringly, as in the mass-imbalanced case discussed above, the CL values continue to agree with benchmark data in this regime, confirming the robustness of CL for spin-polarized systems. Interestingly, a study based on Lefschetz thimbles (see Sec. 2.8) finds a better agreement with iHMC in the virial regime for this value of the imbalance [362]. At larger chemical potential, on the other hand, excellent agreement between thimbles and CL was reported in the same study.

5.3 Non-relativistic fermions in three dimensions: the polarized unitary Fermi gas

One of the most studied many-body systems in recent years is the so-called unitary Fermi gas (UFG) which is associated with a diverging s-wave scattering length. Due to this dominant length scale, the microscopic information about the interaction between the spin-up and spin-down particles becomes effectively irrelevant, and all observed quantities may be written as universal functions of the fermion density and temperature (as these are the only scales left in the system). Its intricate behavior is linked to various many-body phenomena, and moreover it is realized to an excellent precision in numerous cold atoms experiments, having led to a plethora of precise measurements of its unique behavior (see, e.g., [363] for a comprehensive review).

For SU(2) symmetric fermions, i.e. for equal numbers of spin-up and spin-down particles, it is possible – albeit challenging – to study the UFG with conventional Monte Carlo approaches. Past stochastic studies of the UFG have been conducted using methods such as hybrid Monte Carlo (HMC) [364, 365, 366] and auxiliary field methods (AFQMC) on the lattice [367, 368, 369] as well as various flavors of diagrammatic Monte Carlo methods [370, 371] in the continuum. Building on these (and many more) theoretical and experimental advances, it was found that the UFG undergoes a phase transition from a strongly correlated normal fluid to a BCS superfluid at the unusually high transition temperature $T/T_{\rm F}\approx 0.17$. Moreover, the Bertsch parameter, which relates the non-interacting energy to the ground-state energy of the UFG, was determined accurately to be $\xi=0.37(2)$ in excellent agreement with high-precision experiments [372, 373].

For a broken SU(2) symmetry, in which the particle numbers of the up and down components differ, experiments have shed light on the phase diagram in the interesting regime of intermediate polarizations up to $p=(N_{\uparrow}-N_{\downarrow})/(N_{\uparrow}+N_{\downarrow})\approx 0.4$, above which a first-order phase transition (at, and very close to, zero temperature) from the superfluid to the normal state takes place [374, 375, 376, 377, 378]. On the theory side, progress was made in the highly-imbalanced limit of polarons through a variational ansatz for the $N+1$ wavefunction. However, for the full many-body problem with a finite spin-imbalance no sufficiently good ansatz is known. Diffusion MC calculations in the ground state, however, show good agreement with the experimental findings despite the fixed-node approximation [379]. This simplification, which enables such calculations, unfortunately limits the predictive power in terms of possibly exotic pairing states (such as the FFLO mechanism) and was therefore mainly limited to the normal phase.

Nevertheless, several studies based on mean-field theory (see, e.g., [380, 381] for a broad overview on such advances) and (semi-)analytic perturbative [382] as well as non-perturbative studies based on the functional renormalization group [383, 384, 385, 386] and the Luttinger-Ward (LW) formalism [387] have led to insights into the phase diagram spanned by the temperature and polarization. Whereas these studies do not suffer from a sign problem, their application requires the construction of truncation schemes (e.g., in terms of a set of correlation functions that is taken into account). Therefore, complementary studies based on stochastic approaches allow us to enhance our understanding of the spin-imbalanced UFG.

As stated above, a finite spin polarization introduces a sign problem in stochastic approaches which is especially severe at low temperatures, where the interesting physics takes place. To circumvent this sign problem, the CL method was applied to the spin-polarized unitary Fermi gas for the first time in Ref. [343], finding excellent agreement with available benchmark results in the balanced limit. In the left panel of Fig. 19, the density equation of state is shown for a balanced gas in comparison to previously obtained results from DHMC [364] and bold diagrammatic Monte Carlo (BDMC) [371] studies, as well as experimental values from the MIT group [371]. Across a wide range of $\beta\mu$ values, excellent agreement was reported down to low temperatures. In the vicinity of the phase transition to superfluidity, however, the CL results, along with all other theoretical values, start to deviate from the experimental values, which is attributed to finite lattice effects and can, in principle, be mitigated by extending the box sizes beyond $N_{x}=11$.

In the central and right panels of Fig. 19 first predictions for the density and magnetic equations of state are shown, respectively. It is important to note here that these systems require essentially the same computational effort as in the spin-balanced case, whereas a reweighting approach would suffer from an exponential increase in computational effort. As can be appreciated from the figure, excellent agreement with the virial expansion (see e.g. [388]) is achieved at high temperatures, which gives confidence on the reliability of the CL results in that regime. Moving towards lower temperatures, the virial expansion is expected to break down, whereas the CL results continue to evolve smoothly across the entire parameter range studied. With respect to the latter regime, a computation of the pressure equation of state of the polarized UFG (beyond the mean-field approximation) based on the LW formalism [387] is found to agree well with the CL results up to small positive values of $\beta\mu$. Above this regime, the LW study appears to overestimate the density already in the spin-balanced case, as opposed to the CL values.

This precise treatment of the spin-polarized UFG at finite temperature may considered one of the biggest successes of the CL method in the context of non-relativistic physics so far and sets the stage for the search of possibly inhomogeneous superfluid sectors in the phase diagram both at finite and zero temperatures. Such phenomena, most notably FFLO phases with a spatially inhomogeneous order parameter, are not only of great current experimental interest but also reminiscent of the dynamics underlying the QCD phase diagram at finite baryon chemical potential and low temperature, where such inhomogeneous phases may also exist [389, 390]. This highlights a phenomenological similarity between high-energy physics and cold-dilute quantum matter with the perspective that success in one of these fields likely will trigger progress in the other.

6 Summary and Outlook

In this report we have discussed the origin and methods to address the sign problem as it appears across a wide range of systems in relativistic and nonrelativistic quantum many-body physics. We have proposed a broad classification of the various methods into new variables, statistical, and complex plane approaches. After a brief overview of each of those categories, we focused on the complex form of stochastic quantization, namely CL.

Complex Langevin first appeared in the 1980s as a natural generalization of stochastic quantization for cases with a complex action. Operationally, stochastic quantization does not require a probability to be defined (there is no Metropolis accept/reject step), which made field complexification (itself needed in the presence of a sign or complex phase problem) appear more like a feature than a bug. From the mathematical standpoint, however, the challenges seemed daunting. It took decades for the community to begin to understand the properties and behavior of CL, to clarify the origin of its problems and limitations, and to propose solutions. Some of that progress was enabled by advances in hardware, as modern personal computers are powerful enough to run small but useful quantum field theory calculations with little wait; that was certainly not the case in the 1980s, nor in the 1990s when initial explorations of this method were underway.

As computer power continues to grow, and given the overall progress made during the last decade, there is good reason for optimism. Furthermore, there are now several groups and international collaborations around the world applying CL methods to many systems, which generate a wider range of situations than ever before from which insight can be gained. The most remarkable step forward, in the latest chapter in the history of CL, is the derivation of conditions for correctness and how they relate to behavior at the boundaries of the integration region (at infinity and at zeros of the complex weight). This new understanding has spawned new practical solutions such as gauge cooling, dynamic stabilization, and modified actions (regulators), which have enabled more applications than previously thought possible. Continued studies of CL will help illuminate when the method is reliable and when it is not. A more detailed understanding of the structure of the problems in CL might help develop new methods to ameliorate or solve those problems.

Acknowledgements

We are especially grateful to G. Aarts, F. Attanasio, G. Basar, S. Chandrasekharan, P. de Forcrand, C. Gattringer, E. Huffman, K. Langfeld, and E. Seiler for their extremely useful comments on earlier versions of this manuscript. We would also like to acknowledge G. Aarts, C. Gattringer, C. Ratti, D. Sexty, and H. Singh for kindly allowing us to reproduce their figures.

This material is based upon work supported by the National Science Foundation under Grant No. PHY1452635 (Computational Physics Program). C.E.B. acknowledges support from the United States Department of Energy through the Computational Science Graduate Fellowship (DOE CSGF) under grant number DE-FG02-97ER25308. J.B. acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant BR 4005/4-1 (Heisenberg program). J.B. and F.E. acknowledge support by the DFG grant BR 4005/5-1. J.B. and L.R. acknowledge support by HIC for FAIR within the LOEWE program of the State of Hesse. J.B. acknowledges support by the DFG – Projektnummer 279384907 – SFB 1245.