Lefschetz thimbles decomposition for the Hubbard model on the hexagonal lattice

TL;DR

This paper develops algorithms to analyze Lefschetz thimbles in lattice fermion models, specifically applying them to the Hubbard model on a hexagonal lattice, improving understanding of the sign problem and thimbles structure in large systems.

Contribution

It introduces efficient algorithms for Lefschetz thimbles analysis in lattice fermion models and applies them to the Hubbard model, revealing insights into thimbles complexity and sign problem mitigation.

Findings

Existence of an optimal Hubbard-Stratonovich transformation regime.

Agreement between Monte Carlo simulations and exact diagonalization on small volumes.

Higher average sign achieved with Lefschetz thimbles compared to traditional methods.

Abstract

We propose a framework to study the properties of the Lefschetz thimbles decomposition for lattice fermion models approaching the thermodynamic limit. The proposed set of algorithms includes the Schur complement solver and the exact computation of the derivatives of the fermion determinant. It allows us to solve the gradient flow (GF) equations taking into account the fermion determinant exactly, with high performance. We can find both real and complex saddle points and describe the structure of the Lefschetz thimbles decomposition for large enough lattices to extrapolate the results to the thermodynamic limit. The algorithms are described for a general lattice fermion model, with emphasis on two types of lattice discretizations for relativistic fermions (staggered and Wilson), as well as on interacting tight-binding models for condensed matter systems. We apply these algorithms to the Hubbard model on a hexagonal lattice, dealing with lattice volumes as large as 12x12 with 256 steps in Euclidean time, in order to capture the properties of the thimbles decomposition as the thermodynamic, low-temperature, and continuum limits are approached. The complexity of the thimbles decomposition appears to be dependent on the form of the Hubbard-Stratonovich (HS) transformation. In particular, the evidence is provided for the existence of an optimal regime for the hexagonal lattice Hubbard model, with a reduced number of thimbles being important in the overall sum. We have performed quantum Monte Carlo (QMC) simulations using GF to deform the integration contour into the complex plane and demonstrated the agreement with exact diagonalization on small volumes (8 sites in space) if optimal setup for HS transformation is used. The residual sign problem is compared with the state-of-the-art BSS-QMC. We show that the average sign can be kept substantially higher using the Lefschetz thimbles approach.

Tables

Table 1. Table 1: Comparison of observables for exact diagonalization, BSS-QMC (ALF) and two variants of HMC with gradient flow for a 2 × 2 2 2 2\times 2 lattice with N τ = 256 subscript 𝑁 𝜏 256 N_{\tau}=256 , U = 2.0 ​ κ 𝑈 2.0 𝜅 U=2.0\kappa and μ = κ 𝜇 𝜅 \mu=\kappa .

	$⟨\hat{K}⟩$	$⟨{\hat{S}}_x^{\left(1\right)}{\hat{S}}_y^{\left(1\right)}⟩$
ED	19.5781	-0.14624
BSS-QMC	19.587$\pm$0.002	-0.1466$\pm$0.0008
HMC, $\alpha =1.0$	19.65$\pm$0.31	-0.112$\pm$0.0069
HMC, $\alpha =0.8$	19.52$\pm$0.17	-0.142$\pm$0.0062

Table 2. Table 2: Comparison of the sign problem for BSS-QMC (ALF) and two variants of HMC with gradient flow for a 2 × 2 2 2 2\times 2 lattice with N τ = 256 subscript 𝑁 𝜏 256 N_{\tau}=256 , U = 2.0 ​ κ 𝑈 2.0 𝜅 U=2.0\kappa and μ = κ 𝜇 𝜅 \mu=\kappa .

	$⟨\cos \text{Im}S⟩$	$⟨\cos \arg J⟩$	$⟨{\Sigma }_G⟩$
BSS-QMC	0.2363$\pm$0.0032		0.2363$\pm$0.0032
HMC, $\alpha =1.0$	0.9627$\pm$0.0038	0.427$\pm$0.014	0.351$\pm$0.015
HMC, $\alpha =0.8$	0.797$\pm$0.022	0.915$\pm$0.008	0.644$\pm$0.028

Equations

$$⟨\mathcal{O}⟩=\frac{1}{\mathcal{Z}}\int \mathcal{D}\Phi \mathcal{O}[\Phi ]e^{-S[\Phi ]}=\frac{\int \mathcal{D}\Phi \mathcal{O}[\Phi ]e^{-S[\Phi ]}}{\int \mathcal{D}\Phi e^{-S[\Phi ]}}$$

$$=\frac{\frac{1}{{\mathcal{Z}}_{\mathrm{pq}}}\int \mathcal{D}\Phi \mathcal{O}[\Phi ]\frac{e^{-S[\Phi ]}}{e^{-S_R[\Phi ]}}{e}^{-S_R[\Phi ]}}{\frac{1}{{\mathcal{Z}}_{\mathrm{pq}}}\int \mathcal{D}\Phi \frac{e^{-S[\Phi ]}}{e^{-S_R[\Phi ]}}{e}^{-S_R[\Phi ]}}=\frac{⟨\mathcal{O}{e}^{-i{S}_I}{⟩}_{S_R}}{⟨e^{-i{S}_I}{⟩}_{S_R}},$$

$${\mathcal{Z}}_{\mathrm{pq}}=\int \mathcal{D}\Phi e^{-S_R[\Phi ]}$$

$$\frac{\mathcal{Z}}{{\mathcal{Z}}_{\mathrm{pq}}}=e^{-\beta V\Delta f}.$$

$$\mathcal{Z}={\int }_{{\mathbb{R}}^N}\mathcal{D}\Phi e^{-S[\Phi ]}=\sum _{\sigma }{k}_{\sigma }{\mathcal{Z}}_{\sigma },$$

$${\mathcal{Z}}_{\sigma }={\int }_{{\mathcal{I}}_{\sigma }}\mathcal{D}\Phi e^{-S[\Phi ]},$$

$${\frac{\partial S}{\partial \Phi }|}_{\Phi =z_{\sigma }}=0.$$

$$\frac{d\Phi }{dt}=\overline{\frac{\partial S}{\partial \Phi }},$$

$$\mathcal{Z}=\sum _{\sigma }{k}_{\sigma }{e}^{-i\mathrm{Im}S(z_{\sigma })}{\int }_{{\mathcal{I}}_{\sigma }}\mathcal{D}\Phi e^{-\mathrm{Re}S(\Phi )},$$

$$S=S_b+\ln \det M,$$

$$\frac{\partial \ln \det M}{\partial \Phi }=\mathrm{Tr}\left(M^{-1}\frac{\partial M}{\partial \Phi }\right).$$

$$M_{i,j}^{st}=2am{\delta }_{i,j}+$$

$$({\eta }_{i,1}{e}^{\mu a}{U}_{i,1}{\delta }_{i+\hat{1},j}-{\eta }_{j,1}{U}_{j,1}^{†}{e}^{-\mu a}{\delta }_{i-\hat{1},j})+$$

$$\sum _{\nu =2}^4({\eta }_{i,\nu }{U}_{i,\nu }{\delta }_{i+\hat{\nu },j}-{\eta }_{j,\nu }{U}_{j,\nu }^{†}{\delta }_{i-\hat{\nu },j})$$

$$M^{st}(U)=$$

$$\displaystyle\left(\begin{array}[]{ccccccc}B_{1}&e^{\mu a}[U_{(1,x),1}]&0&0&0&\ldots&e^{-\mu a}[U^{\dagger}_{(N_{\tau},x),1}]\\ -e^{-\mu a}[U^{\dagger}_{(1,x),1}]&B_{2}&e^{\mu a}[U_{(2,x),1}]&0&0&\ldots&0\\ 0&-e^{-\mu a}[U^{\dagger}_{(2,x),1}]&B_{3}&\ddots&0&\ldots&0\\ 0&0&\ddots&\ddots&\ddots&\ldots&0\\ \vdots&&&\ddots&\ddots&\ddots&\vdots\\ 0&0&\ldots&0&-e^{-\mu a}[U^{\dagger}_{(N_{\tau-2},x),1}]&B_{N_{\tau}-1}&e^{\mu a}[U_{(N_{\tau-1},x),1}]\\ -e^{\mu a}[U_{(N_{\tau},x),1}]&0&\ldots&0&0&-e^{-\mu a}[U^{\dagger}_{(N_{\tau-1},x),1}]&B_{N_{\tau}}\\ \end{array}\right).$$

$${\overline{M}}^{st}(U)=$$

$$\displaystyle\left(\begin{array}[]{cccccc}1&D_{1}&0&0&0&\ldots\\ 0&1&D_{2}&0&0&\ldots\\ 0&0&1&D_{3}&0&\ldots\\ 0&0&0&1&D_{4}&\ldots\\ \vdots&&&&\ddots&\\ -D_{2N_{\tau}}&0&0&&\ldots&1\\ \end{array}\right).$$

$$D_{2k}$$

$$D_{2k-1}$$

$${{\overline{M}}^{st}}^{-1}(U)=$$

$$\displaystyle\left(\begin{array}[]{cccccc}g_{1}&\ldots&\ldots&\ldots&\ldots&\bar{g}_{2N_{\tau}}\\ \bar{g}_{1}&g_{2}&\ldots&\ldots&\ldots&\ldots\\ \ldots&\bar{g}_{2}&g_{3}&\ldots&\ldots&\ldots\\ \ldots&\ldots&\bar{g}_{3}&g_{4}&\ldots&\ldots\\ \vdots&&&&\ddots&\\ \ldots&\ldots&\ldots&&\ldots&g_{2N_{\tau}}\\ \end{array}\right).$$

$${\bar{g}}_{i+1}=D_{i+1}^{-1}{\bar{g}}_i{D}_i.$$

$${\overline{M}}^W=$$

$$\displaystyle\left(\begin{array}[]{cccccc}A_{1}&C_{1}&0&0&\ldots&0\\ 0&A_{2}&C_{2}&0&\ldots&0\\ 0&0&A_{3}&\ddots&\ldots&0\\ \vdots&&&\ddots&&\vdots\\ 0&0&\ldots&&A_{N_{\tau}-1}&C_{N_{\tau}-1}\\ -C_{N_{\tau}}&0&\ldots&&0&A_{N_{\tau}}\\ \end{array}\right).$$

$$\det {\overline{M}}^W=\det Q^W\prod _i\det A_i,$$

$$Q^W=$$

$$\displaystyle\left(\begin{array}[]{cccccc}I&{A_{1}}^{-1}&0&0&\ldots&0\\ 0&I&C_{1}&0&\ldots&0\\ 0&0&I&{A_{2}}^{-1}&\ldots&0\\ \vdots&&&\ddots&\ddots&\vdots\\ 0&0&\ldots&&I&A_{N_{\tau}}^{-1}\\ -C_{N_{\tau}}&0&\ldots&&0&I\\ \end{array}\right).$$

$$D_{2k}=C_k,D_{2k-1}={A_k}^{-1},k=\mathrm{1...}{N}_{\tau }.$$

$$\frac{\partial \ln \det {\overline{M}}^W}{\partial U_{(t,x),\nu }}=\mathrm{Tr}\left({\bar{g}}_{2t}\frac{\partial C_t}{\partial U_{(t,x),\nu }}\right)-$$

$$\mathrm{Tr}\left({\bar{g}}_{2t-1}{A_t}^{-1}\frac{\partial A_t}{\partial U_{(t,x),\nu }}{A_t}^{-1}\right)+\mathrm{Tr}\left({A_t}^{-1}\frac{\partial A_t}{\partial U_{(t,x),\nu }}\right).$$

$$\frac{\partial \ln \det {\overline{M}}^W}{\partial U_{(t,x),1}}=\mathrm{Tr}\left({\bar{g}}_{2t}\frac{\partial C_t}{\partial U_{(t,x),1}}\right).$$

$$\hat{H}=-\kappa \sum _{⟨x,y⟩}({\hat{a}}_x^{†}{\hat{a}}_y+{\hat{b}}_x^{†}{\hat{b}}_y+\text{mbox}h.c)+\frac{U}{2}\sum _x{\hat{q}}_x^2+\mu \sum _x{\hat{q}}_x,$$

$$e^{-\beta \hat{H}}\approx ...e^{-\delta \hat{K}}{e}^{-\delta {\hat{H}}_U}{e}^{-\delta \hat{K}}{e}^{-\delta {\hat{H}}_U}...+O({\delta }^2)$$

$$\frac{U}{2}{\hat{q}}_x^2=\frac{\alpha U}{2}{\hat{q}}_x^2-\frac{(1-\alpha )U}{2}{\hat{s}}_x^2+(1-\alpha )U{\hat{s}}_x,$$

$$e^{-\frac{\delta }{2}\sum _{x,y}{U}_{x,y}{\hat{n}}_x{\hat{n}}_y}\cong \int D{\varphi }_x{e}^{-\frac{1}{2\delta }\sum _{x,y}{\varphi }_x{U}_{xy}^{-1}{\varphi }_y}{e}^{i\sum _x{\varphi }_x{\hat{n}}_x},$$

Full text

Hubbard model, sign problem, Lefschetz thimbles

I Introduction

One of the most widely used paths to study, non-perturbatively, the physics of strongly coupled quantum systems, in a fully ab-initio manner is through Monte Carlo simulations of the Feynman path integral. For many systems, the Euclidean formulation of the functional integral yields a real, positive-definite action whose evaluation can be performed using importance sampling. It often occurs, however, that the action is complex. Relevant examples of such systems abound in disparate branches of physics. The theory of the strong interactions, quantum chromodynamics (QCD), exhibits a complex action at finite baryon density Philipsen (2010); Aarts (2016). The study of the QCD phase diagram is essential to the understanding of the quark-gluon plasma, a strongly-coupled fluid which exists above the deconfinement transition, as well addressing deep questions in astrophysics and cosmology. Currently, these phases of matter are actively studied experimentally at several collider facilities, including RHIC and LHC. Additional examples from high-energy physics which are well-known to exhibit a complex action include gauge theories with the addition of a theta term, Chern-Simons theories, and matrix models which provide a non-perturbative definition of string theories. The sign problem also impedes ab-initio studies of many-body systems. The physics of Feshbach resonances in cold atomic Fermi gases, frequently realized in the laboratory and modeled with the unitary Fermi gas (UFG) is one prominent example Inguscio et al. (2007); Schäfer and Teaney (2009); Rammelmüller et al. (2018). A well-studied example from condensed matter physics is the Hubbard model. Despite its simplicity, the Hubbard model captures the physics of the Mott metal-insulator transition, and probably, high-temperature superconductors Baeriswyl et al. (1995); Lee et al. (2006). On a bipartite lattice at half-filling, the Hubbard model is free from the sign problem due to particle-hole symmetry. However, as soon as frustration or non-zero chemical potential appear in the Hamiltonian, one is faced with the sign problem. In fact, it could be argued that the majority of interesting systems are plagued by the sign problem. A key result, due to Troyer and Wiese Troyer and Wiese (2005), states that the sign problem is an NP-hard problem in a generic, Ising spin-glass system. It follows that a general solution to all sign problems is an unlikely proposition.

There are several approaches one can take to dealing with the problem. The most naive approach to deal with a system exhibiting a complex action is to absorb the imaginary part of the action in the observable and sample according to the real part of the action. This method is known as reweighting and is based on the following identity

[TABLE]

where $S=S_{R}+iS_{I}$, the ratio $\nicefrac{{e^{-S[\Phi]}}}{{e^{-S_{R}[\Phi]}}}$ is the reweighting factor, and

[TABLE]

is the phase quenched partition function. The angular brackets in (I) denote an ensemble average with respect to the measure $\mathcal{D}\Phi\;e^{-S_{R}}$. Although this sequence of expressions is nothing more than a rewriting of the standard thermal ensemble average, the practical calculation of observables using reweighting is exponentially difficult due to the sign problem. The last ratio in (I) is not well defined, since both the numerator and the denominator vanish exponentially as the spacetime volume is increased. This is a manifestation of what is coined as the overlap problem, i.e. the phase-quenched theory is different from the full theory which renders sampling from the former a highly ineffective approach to mimic sampling from the full theory. The technical issue at hand is the overlap of the ensemble sampled according to $S_{R}$ and the original ensemble that involves the entire action. A physical meaning can be attached to this difficulty by considering the average sign, $\langle e^{-iS_{I}}\rangle_{S_{R}}$, which can be understood as a ratio of two partition functions

[TABLE]

In (3) we have introduced the spatial volume $V$, inverse temperature $\beta$, and $\Delta f$, which is the free energy density difference between the two ensembles. Although $\Delta f$ is formulation dependent, one cannot cure the exponential scaling using naive reweighting. In any Monte Carlo calculation, the error on the mean scales with the computational time, $T_{CPU}$, as $1/\sqrt{T_{CPU}}$. Thus, in order to have the error on the average sign less than the value of the average sign itself, we must require that $T_{CPU}\gg e^{2\beta V\Delta f}$. We refer the interested reader to Aarts (2016) for a pedagogical and detailed presentation of the sign and overlap problems.

Recently, much progress has been made by complexifying the fields of systems suffering from the sign problem. This idea, which can easily be demonstrated in simple, one-dimensional integrals, has already been applied to a number of models, which we list below. Furthermore, in many of these cases, the initial severe sign problem was eliminated or substantially weakened. One successful approach along the previously mentioned lines is complex Langevin dynamicsAarts (2012); Aarts et al. (2013); Sexty (2014); Nagata et al. (2016); Aarts et al. (2017); Bloch et al. (2018); Anagnostopoulos et al. (2018); Scherzer et al. (2019). Another method, which we use in this study is the method of Lefschetz thimbles. Originally introduced in Witten (2010, 2011), it was not long after that lattice gauge theory practitioners sought to apply these methods to QCD at finite baryon density Cristoforetti et al. (2012). Pioneering studies using Lefschetz thimbles were performed on the relativistic Bose gas for lattices volumes up to $V=8^{4}$, showing good agreement with complex Langevin simulations Cristoforetti et al. (2013); Fujii et al. (2013); Cristoforetti et al. (2014). Several other studies have investigated a variety of other systems displaying a sign problem including $O(n)$ sigma models Tanizaki (2015), chiral random matrix ensembles Di Renzo and Eruzzi (2015), and the $U(1)$ one-link model using techniques borrowed from reweighting Bluecher et al. (2018). A significant hurdle was overcome as several groups extended the method of Lefschetz thimbles to interacting fermions in $0+1$ dimensions as well as at a single site Fujii et al. (2015); Tanizaki et al. (2016); Kanazawa and Tanizaki (2015); Alexandru et al. (2016a, b); Di Renzo and Eruzzi (2018). The successful application of these methods was then subsequently extended to field theories of strongly interacting fermions in both $1+1$ as well as $2+1$ dimensions Alexandru et al. (2017, 2018a, 2018b, 2018c). A short description of the results presented in this article and related to the Hubbard model originally appeared in our previous letter Ulybyshev et al. (2019a). In a subsequent preprint  Fukuma et al. (2019), an alternative approach to deal with ergodicity issues was applied to Hubbard model simulations within the thimbles formalism. Albeit, there is a big difference in the regimes studied by the two groups. In Ulybyshev et al. (2019a) and in this paper, the low temperature limit of the Hubbard model was studied, while the authors of Fukuma et al. (2019) reported results for significantly higher temperatures and thus of milder sign problem for which BSS-QMC has an average sign greater than 0.6 in the whole range of parameters studied in their paper. On the other hand, both in Ulybyshev et al. (2019a) as well as in this article, we are addressing the region of strong sign problem where BSS-QMC, even with an optimal setup, experiences an exponential decay of the average sign.

One of the main concerns about the efficacy of the Lefschetz thimbles method is the scaling of the number of thimbles with the system size and temperature. One might argue that the amount of important thimbles can grow exponentially once we approach the thermodynamic or low temperature limit. In this case, by construction, the method is unable to improve the sign problem. The problem is especially hard for fermionic models, since it is usually quite difficult to find the exact saddle points if the logarithm of the fermionic determinant is included in the action. Here, we propose a set of algorithms to address this issue and to elucidate, clearly and systematically, the non-trivial saddle point structure of the theory. Due to a more efficient calculation of the exact derivatives of the fermionic determinant, we are now able to reveal the construction of the Lefschetz thimble decomposition on large lattices and extrapolate our results to the thermodynamic limit. This also represents the main difference of our paper from earlier attempts to apply the Lefschetz thimbles decomposition to the Hubbard model Mukherjee and Cristoforetti (2014), where the thimbles decomposition was not optimised and only one thimble, out of many important ones, was taken into account. As a result, those simulations actually did not represent a full calculation of the functional integral, but rather represented only corrections to dynamical mean field theory (DMFT) results. Using a complete study of the saddle point structure of the Hubbard model, and identifying the advantageous regions in parameter space, one can safely proceed to address the sign problem using Lefschetz thimbles.

We start with a short introduction to the formalism, and proceed with the description of the method to solve the gradient flow equations for Wilson and staggered fermions. After this, we describe the application of the technique to the Hubbard model on the hexagonal lattice. First, we make a detailed study of the saddle points, which is an essential ingredient of the Lefschetz thimbles method. In particular, we explore the dependence of saddles on volume, the Hubbard coupling $U$, and chemical potential. Among other things, we discuss at length the algorithms used to search for saddle point configurations away from half-filling, when saddle points are shifted into complex space $\mathbb{C}^{N}$. Finally, in order to support our conclusions concerning the role of different saddle points, we perform Monte Carlo calculations over manifolds in complex space and compare results with exact diagonalization. In addition to that, we show that the average sign can be substantially increased even in comparison with BSS-QMC. This fact means that we can potentially construct a superior algorithm for dealing with the sign problem, if the additional computational costs associated with the gradient flow and integration over curved manifolds in complex space are improved upon.

II Lefschetz thimbles formalism

Let us first consider the complexification of the fields appearing in the functional integral (I), $\Phi\in\mathbb{C}^{N}$. This amounts to a shift of the contour of integration into complex space. We are allowed to do so, as Cauchy’s theorem tells us that one can choose any appropriate contour in complex space as long as the integral still converges and no poles of the integrand are crossed during this shift. As we will demonstrate, both of these conditions are satisfied. We now introduce one particularly useful representation, known as the Lefschetz thimble decomposition of the partition function Witten (2010, 2011),

[TABLE]

and $\sigma$ labels all complex saddle points $z_{\sigma}\in\mathbb{C}^{N}$ of the action, which are determined by the condition

[TABLE]

The integer-valued coefficients $k_{\sigma}$, are the intersection numbers and $\mathcal{I_{\sigma}}$ are the Lefschetz thimble manifolds attached to the saddle points $z_{\sigma}$. These manifolds, defined below, are the generalization of the contours of steepest descent in the theory of asymptotic expansions. We stress that if the saddle points are non-degenerate (${\left.\det\partial^{2}S/\partial\Phi^{\prime}\partial\Phi\right|}_{\Phi=z_{\sigma}}\neq 0$) and isolated, the relation (II) holds (for a generalization to the case of gauge theory see Witten (2011)).

The Lefschetz thimble manifold associated with a given saddle point is the union of all solutions of the following differential equation

[TABLE]

known as the gradient flow (GF) equations, which satisfy the following boundary condition: $~{}\Phi\in\mathcal{I}_{\sigma}:\Phi(t\rightarrow-\infty)\rightarrow z_{\sigma}$. Just as we made an analogy between the thimble and the contour of steepest descent, there is a second manifold associated with each saddle point which is analogous to the contour of steepest ascent. This manifold is known as the anti-thimble, $\mathcal{K}_{\sigma}$, and consists of all possible solutions of the GF equations (6) which end up at a given saddle point $z_{\sigma}$: $\Phi\in\mathcal{K}_{\sigma}:\Phi(t)=\Phi,\Phi(t\rightarrow+\infty)\rightarrow z_{\sigma}$. The intersection number $k_{\sigma}$ is defined by counting the number of intersections of $\mathcal{K}_{\sigma}$ with the original integration domain: $\mathbb{R}^{N}$, $k_{\sigma}=\langle\mathcal{K}_{\sigma},\mathbb{R}^{N}\rangle$. An example scheme of thimbles and anti-thimbles is drawn in the Fig. 1.

It is worth noting that thimbles and anti-thimbles are both real, $N$-dimensional manifolds embedded in $\mathbb{C}^{N}$. We now state two key properties of the thimbles, which follow from (6) coupled with the fact that the action, $S$, is regarded as a holomorphic function of the complex fields. These properties are that the real part of the action, $\mbox{Re}\,S$, monotonically increases along the thimble, starting from the saddle point and the imaginary part of the action, $\mbox{Im}\,S$, stays constant along it. The first property is essential in guaranteeing the convergence of the individual integrals in (II), while the second one obviously makes the method attractive with regards to the weakening of the sign problem. Using these crucial properties, it follows that neither thimbles nor anti-thimbles can intersect each other, no two saddle points can, in general, be connected by a thimble (with the very important exception which is discussed below), and all integrals on the r.h.s. of (II) are convergent.

As a result of the above discussion, it follows that (II) can be written as

[TABLE]

where we have explicitly written out the complex factors associated with different saddle points. Usually, thimbles can be classified as being either “relevant” or “irrelevant” using the intersection number. Relevant ones have their intersection number, $k_{\sigma}$, being nonzero and thus participate in the sum in (7). Conversely, a thimble is irrelevant if it has a zero intersection number. However, this classification can fail if the so-called Stokes phenomenon occurs for saddle points lying within $\mathbb{R}^{N}$. By definition, the Stokes phenomenon means that the saddle points are connected by a thimble. In this case, $k_{\sigma}$ is not well-defined and we need other tools in order to classify the saddle points. An example of such a situation will be demonstrated later on when we will study the actual saddle points for the Hubbard model.

As one can see, the initial sign problem has been split into two parts. The first part of the residual sign problem concerns the constant phase factors, $e^{-i\,\mathrm{Im}\,S(z_{\sigma})}$. The number of relevant thimbles, their weight, and the distribution of the imaginary part of the action at corresponding saddles define the remaining severity of the first part of the sign problem. An ideal situation arises when this sum only contains one dominant term. The second part of the residual sign problem relates to the fluctuations of the complex measure, $\mathcal{D}\Phi$, in the integration over the thimble. Potentially, there is the third source of the residual sign problem: residual fluctuations of $\mathrm{Im}\,S$, which appear if we are not following the thimble exactly. All these issues will be addressed below in our test calculations for the Hubbard model, where we perform a thorough classification of saddle points and then give an estimate for the fluctuations of the complex measure and the residual fluctuations of $\mathrm{Im}\,S$. We now present a description of our numerical methods.

III Algorithms

The GF equations (6) are the basis of the whole formalism. Here we present the set of algorithms, which allows us to solve them efficiently at least for lattices of moderate size. The main difficulty in implementing GF is the presence of the fermionic determinant in the action for a typical lattice field theory (or model -for the case of condensed matter systems) with fermions

[TABLE]

where $S_{b}$ is the bosonic part and the fermionic operator $M$ is more or less a sparse matrix with dimensionality $\sim N_{\tau}N_{s}$ (ignoring for the moment color and flavor indices). Here, $N_{\tau}$ is the Euclidean time extent of the lattice and $N_{s}$ is the number of degrees of freedom in space. The latter typically includes the number of sites in space (in the context of QCD one should also take into account the number of colors and flavors). The construction (8) is the same both for lattice field theories and interacting tight-binding models in condensed matter physics. The key element of our algorithms is the efficient calculation of the derivatives of the fermionic determinant with respect to the bosonic fields, which is essential for the solution of the GF equations. The derivatives of the logarithm of the fermionic determinant can be computed directly using the simple relation

[TABLE]

It turns out that this requires the knowledge of only a few elements of the fermion propagator $M^{-1}$, since the bosonic fields $\Phi$ enter the fermionic operator $M$ locally.

In the following considerations we rely on the special band structure of the fermionic operator. We start with unimproved staggered fermions, whose fermionic operator can be written as

[TABLE]

with the usual staggered phases $\eta_{i,\nu}=(-1)^{i_{1}+...+i_{\nu-1}}$ and gauge fields $U_{i,\nu}$. Here $\mu$ is the chemical potential and $m$ is the mass of fermions, and both are multiplied by the lattice spacing $a$. The four-dimensional index $i=(t,x)$ consists of both the temporal $t$ and the three-dimensional spatial part $x$. It is convenient to introduce the spatial part of the fermionic operator $B_{t}$, which contains all elements of the matrix (III), diagonal in Euclidean time direction for a given time slice $t$. After doing so, (III) can be rewritten as a block matrix consisting of blocks $N_{s}\times N_{s}$:

[TABLE]

$[U_{(t,x),1}]$ is a diagonal matrix which contains on the main diagonal all gauge field exponents in the Euclidean time direction for a given timeslice $t$. Following Hasenfratz and Toussaint (1992), the determinant of (III) is equivalent to the determinant of the following matrix:

[TABLE]

Now all blocks are of the size $2N_{s}\times 2N_{s}$:

[TABLE]

The same general form of the fermionic operator is also common for interacting tight-binding models in condensed matter physics. See Buividovich et al. (2019) and references therein for more details. The only difference is that the blocks would be of the size $N_{s}\times N_{s}$ in the case of an interacting tight-binding model, and their internal construction is also different. However, these details are not important for the present discussion.

The inverse fermionic matrix can also be written in terms of spatial $2N_{s}\times 2N_{s}$ blocks,

[TABLE]

The matrix ${\overline{M}^{st}}^{-1}(U)$ is dense, but here we explicitly show only those blocks which are needed for our calculations. In fact, in the trace in (9), only the off-diagonal blocks ,$\bar{g}_{n},n=1...N_{\tau}$, will contribute to the exact derivatives. Following Buividovich et al. (2019), an iterative procedure can be used to compute all needed elements of the fermionic propagator,

[TABLE]

Once we know one off-diagonal block $\bar{g}_{n}$ for some $n$, we can, in principle, reconstruct all of them. Of course, we need to invert the $D_{i}$ blocks after each update of the gauge fields, but taking into account their sparsity (31), it costs no more than $N_{s}^{3}$ operations for each block. In practice, these iterations typically can not last for more than $N_{\rm{Schur}}\sim 10$ time slices due to accumulation of round-off errors. Thus, we compute the Green’s functions $\bar{g}_{n}$ from scratch for each $n=kN_{\rm{Schur}},\ k=0,1,2...$ and use iterations (39) only in between for intermediate time slices. The Schur complement solver Ulybyshev et al. (2019b) is used for finding $\bar{g}_{n}$ from scratch, with additional simplifications described in Buividovich et al. (2019). The solver, including the iterations (39), scales as $N_{s}^{3}N_{\tau}$. However, despite the scaling being worse than that of iterative solvers, the method still gives substantial speedup in comparison with the calculation of fermionic determinant using stochastic estimators. The reason is twofold: 1) there is a very small prefactor in the scaling relation, which compensates for the $N_{s}^{3}$ term at least for lattices up to $N_{s}\sim 10^{3}$; 2) there is no need to repeatedly find solutions for multiple stochastic estimators, since we get all exact derivatives after one application of the Schur solver accompanied with the propagation through the entire Euclidean time extent of the lattice according to (39). A more careful analysis of the performance of the Schur solver and some benchmarks against an iterative solver were done in Ulybyshev et al. (2019b).

The situation is a bit more complicated for Wilson fermions. In this case, we used the derivation of the compressed form of the Wilson fermionic operator from Nagata and Nakamura (2010). Disregarding the constant determinant of the permutation matrix, the determinant of the Wilson fermionic operator is equivalent to the determinant of the following matrix,

[TABLE]

$A$ and $C$ are blocks of size $4N_{s}\times 4N_{s}$. A factor of four appeared due to the Dirac index of Wilson fermions, and $N_{s}$ includes both the number of sites in space and color degree of freedom as in the case of staggered fermions. We refer to Nagata and Nakamura (2010) for the derivation of this relation and for the exact form of the blocks $A$ and $C$. The gauge fields enter in both the diagonal and off-diagonal blocks, although the diagonal blocks $A_{i}$ include only spatial links, while the off-diagonal blocks $C_{i}$ include both spatial links and the links in Euclidean time direction. The determinant of the operator (III) can be simplified further as,

[TABLE]

where

[TABLE]

${Q}^{W}$ has exactly the form needed for the Schur solver (III), thus the blocks of $(Q^{W})^{-1}$ essential for the derivatives can be computed with exactly the same algorithm as described above with only one substitution,

[TABLE]

Equation (9) is also modified. For spatial links $\nu=2,3,4$ it can be written as

[TABLE]

For the links in the Euclidean time direction, the expression is a bit simpler, since they do not enter into the $A_{i}$ blocks,

[TABLE]

In both expressions above, $\bar{g}_{i}$ are the off-diagonal blocks of $(Q^{W})^{-1}$, enumerated according to the convention in (III). Despite the slightly more complicated expressions for Wilson fermions, one needs only to invert the spatial blocks and apply the Schur complement solver to compute all exact derivatives of the fermionic determinant. Thus, the scaling of the method is also $N_{s}^{3}N_{\tau}$.

IV The model

In order to demonstrate the performance of the described algorithms, we consider the Hubbard model on the hexagonal lattice at finite chemical potential. At half-filling, this model is known to exhibit a semimetal-to-insulator transition (Assaad and Herbut (2013); Sorella et al. (2012)), which gives us an appropriate scale for the interaction strength. Furthermore, the particle-hole symmetry at half-filling helps to identify and characterize the thimbles and saddle points before we increase the chemical potential. We start from the form of the Hamiltonian written in the particle-hole basis in order to have a manifestly positive-definite weight for the Hubbard field configurations at half filling

[TABLE]

where $\hat{a}^{\dagger}_{x}$ and $\hat{b}^{\dagger}_{x}$ are creation operators for electrons and holes, $\hat{q}_{x}=\hat{n}_{x,\text{el.}}-\hat{n}_{x,\text{h.}}={\hat{a}^{\dagger}_{x}}\hat{a}_{x}-{\hat{b}^{\dagger}_{x}}\hat{b}_{x}$ is the charge operator, $\kappa$ is the hopping parameter, $U>0$ is the Hubbard interaction, and $\mu$ is the chemical potential. Due to the van Hove singularity present in the density of states, one can clearly identify the scale where new physics is expected: $\mu=\kappa$ (see McChesney et al. (2010); Körner et al. (2017) and references therein). Special attention will be paid to this value of $\mu$ in the calculations that follow, since it is also the region with the most severe sign problem Ying and Wessel (2018), as can be seen from the test calculations made with BSS-QMC (see Fig. 2). Interaction strength is chosen in order to probe both sides of the antiferromagnetic (AFM) phase transition, which happens at $U\approx 3.8\kappa$ Assaad and Herbut (2013); Buividovich et al. (2018) on hexagonal lattice. Thus, we usually take tree values of $U$: $U=2\kappa$ of $U=3\kappa$ in semi-metallic phase, $U=3.8\kappa$ at the phase transition and $U=5\kappa$ in AFM phase. Temperature is chosen to be below the Neel temperature for $U=5\kappa$, so that we indeed are in ordered phase at that value of interaction strength Buividovich et al. (2018). Also, in this region of temperatures ($\beta\kappa=20...30$), the BSS-QMC method experiences the exponential decay of the average sign, reaching as small values as $10^{-2}...10^{-3}$ even on small $2\times 2$ lattices, as it is shown in Fig. 2. The last argument, why we think this temperature is small enough stems from comparison with other recent papers dealing with QMC calculations for doped Hubbard model: $\beta\kappa=5$ for BSS-QMC calculations in Huang et al. (2019) and $\beta\kappa=3$ in tempered Lefschetz thimbles method Fukuma et al. (2019). This comparison shows that in our saddle point analysis we could reach temperatures well below those achieved in the state-of-the-art QMC computations.

The next step in constructing the path integral formulation of the model is to introduce the Trotter decomposition:

[TABLE]

where $\hat{{K}}$ is the collection of all bilinear fermionic terms in $\hat{{H}}$, and $\hat{{H}}_{U}$ is the interaction part of the full Hamiltonian. Here we have introduced $\delta$, which specifies the discretization of Euclidean time, $N_{\tau}\delta=\beta$, where $N_{\tau}$ constitutes the Euclidean time extent of the lattice. Below, we will refer to $\beta$ in the units of inverse hopping.

One can obtain an additional, nonphysical, degree of freedom in the Hamiltonian, by applying the following identity to the interaction term

[TABLE]

where $\hat{s}_{x}=\hat{n}_{x,\text{el.}}+\hat{n}_{x,\text{h.}}$ is the spin operator. We can now simultaneously introduce two continuous auxiliary fields by applying the standard Hubbard-Stratonovich transformations to each four-fermion term in (60),

[TABLE]

The first four-fermionic term can be transformed into a bilinear using (61), and the second using (62). This is not the most general possible decomposition of four-fermionic terms into bilinear ones, but the one most commonly used in QMC algorithms with continuous auxiliary fields. This representation was first proposed in White et al. (1988) and was also used in several recent papers Beyl et al. (2018); Ulybyshev and Valgushev (2017). The parameter $\alpha\in[0,1]$ defines the balance between auxiliary fields coupled to the charge ($\hat{q}_{x}$) and spin ($\hat{s}_{x}$) density. This particular representation has an important advantage over discrete auxiliary fields in that it also works for non-local interactions, so that we do not need to introduce a new auxiliary field for every pair of interacting electrons.

The details of the construction of the path integral are straightforward and can be found in Ulybyshev et al. (2013); Smith and von Smekal (2014); Beyl et al. (2018). Here we simply state the explicit form of the partition function which we have used in our calculations:

[TABLE]

where the fermionic operators are given by

[TABLE]

This fermionic operator can be rewritten in the form (III), where

[TABLE]

Thus, all algorithms for the fast solution of the GF equations, described in the section III are fully applicable. We denote the field coupled to charge density as $\phi_{x,\tau}$, and the field coupled to spin density as $\chi_{x,\tau}$. The full action, which is used in Monte Carlo sampling, involves both the bosonic action of the auxiliary fields as well as the logarithm of the fermionic determinants, $S=S_{\alpha}-\ln(\det M_{\text{el.}}\det M_{\text{h.}})$. The total number of auxiliary fields is equal to $N=2N_{s}N_{\tau}$ if $\alpha\in(0,1)$, so that both fields participate, and $N=N_{s}N_{\tau}$ if $\alpha=0,1$, where only one type of field remains.

V Saddle points study

V.1 Saddle points at half-filling

Our goal is to study realistic lattice volumes in order to get a quantitative idea of what the thimbles decomposition (II) looks like as we approach both the thermodynamic limit in spatial volume and the continuum limit in Euclidean time. Unfortunately, at large lattice volumes, the fully analytical approach for finding saddle points (as was done in Ulybyshev and Valgushev (2017) on lattices with up to four sites) does not work. Thus, in this study we are using a completely different approach which is based on importance sampling and fast solutions of the GF equations, using the calculations of the derivatives of the fermionic determinant described in section III.

At half-filling, this method starts with the generation of lattice configurations using standard hybrid Monte Carlo (HMC) techniques. After this, we numerically integrate the GF equations for each field configuration for a finite flow time, in order to reach the local minimum of the action. At half filling, when thimbles are bounded within $\mathbb{R}^{N}$, the local minimum of the action always corresponds to a relevant saddle point. At the end of this sequence of steps, the distribution of lattice ensembles, taken after employing the GF procedure, gives an accurate characterization of the relevant saddle points at half-filling if the initial set of configurations was ergodic. An example of such a process is shown in Fig. 3. After generating configurations using HMC, one can observe the approach to the saddle point in our gradient flow routine. As noted, the real part of the action should monotonically decrease and eventually, at a certain flow time, converge to the value at the saddle. In general, the method scales similar to the Schur complement solver as $N_{s}^{3}N_{\tau}$.

A possible source of systematic error in our lattice set up is the discretization in Euclidean time that results from the Trotter decomposition. Thus, we first checked that we have already effectively arrived at the continuum limit in Euclidean time. In Fig. 4, the plot shows the histogram of the distribution of the action for the field configurations after GF. As the initial configurations were generated using HMC, the height of each bar corresponds to the exact weight of the thimble attached to the corresponding saddle point whose value of the action is denoted by the position of the bar. In Fig. 4 we display the histograms for two lattice spacings at fixed $\beta$. The results are almost identical, and thus we can claim that with $N_{\tau}=256$ at $\beta=20$, we are already close enough to the continuum limit in Euclidean time. This gives us confidence that our study of the features of the saddle points and thimbles is independent of the step size in Euclidean time. We will use the same style of plots to characterize the structure of the thimbles decomposition below.

We now proceed to study saddle points at different $\alpha$. One important thing to note is that at half-filling, we cannot faithfully sample the path integral at the extreme values $\alpha=1.0$ and $\alpha=0.0$. In both cases (see Beyl et al. (2018); Ulybyshev and Valgushev (2017); Wynen et al. (2018)), the product of fermionic determinants is equal to the square of some real-valued function

[TABLE]

Thus only one constraint, $F=0$, needs to be satisfied in order to have both fermionic determinants equal to zero. It follows that the dimensionality of the manifolds on which the determinant vanishes is equal to $N-1$ and therefore they cut $\mathbb{R}^{N}$ into disconnected regions. As a result, HMC can not penetrate through these domain walls Beyl et al. (2018); Ulybyshev and Valgushev (2017); Wynen et al. (2018), and we cannot rely on it to generate an ergodic set of configurations. However, as was shown in Fig. 5 (and one more example will also be shown below), even a small shift of $\alpha$ from these extreme values is enough to restore ergodicity. We explicitly check the value of squared spin per sublattice

[TABLE]

and compare the results from HMC with BSS-QMC, which does not have ergodicity issues due to the formulation in terms of discrete fields. We can thus safely use, e.g. $\alpha=0.01$ and $\alpha=0.99$ in order to gain an understanding of the thimbles decomposition when we have either a dominant spin-coupled field auxiliary field or a dominant charge-coupled auxiliary field, respectively.

We first study saddle points at $\alpha=0.01$, when the spin-coupled field $\chi_{x,\tau}$ is the most important and $\phi_{x,\tau}$ is always equal to zero at the saddle points. The results are shown in Fig. 6 at different values of the interaction strength, which correspond to the semi-metal (SM) phase ($U=3\kappa$), the region close to the phase transition ($U=3.8\kappa$), and the antiferromagnetic (AFM) phase ($U=5.0\kappa$). In all cases, there is a dominant saddle point corresponding to the smallest value of the action, but its dominance becomes less and and less pronounced as we move towards the AFM phase. A more detailed study of saddles is presented in Fig. 7 for the case of large interaction strength, $U=5\kappa$. The lowest saddle (Fig. 7(b)) is just a static solution which corresponds to $\chi_{x,\tau}=\pm\chi_{0}$, with the sign depending on the sublattice. Both these saddles are just two identical mean field solutions corresponding to antiferromagnetic ordering (Néel state). They appear as a consequence of the bipartite nature of the lattice. Since there are two stationary vacua, “instanton” solutions, which represent tunneling events between the two vacua, inevitably appear. We indeed observe these instantons which correspond to saddle points with larger actions, examples of which are shown in Fig. 7c-e. In (c), one sees how, at a given site, one can have virtually instantaneous tunneling of the value of $\chi$ between $+\chi_{0}$ and $-\chi_{0}$ and back again, where $\chi_{0}$ is the same value that appears in the mean-field configuration. In these configurations, the tunneling occurs locally both in space and in time. However, there are cases such as (d), where the tunneling from $\pm\chi_{0}$ to $\mp\chi_{0}$ occurs all across the lattice in space at some Euclidean time $\tau_{0}$, and at a later time $\tau_{1}$, the fields return to their original configuration. Thus, we have two “global instantons”. Finally, we have observed cases such as (e), where a similar pair of “global instantons” exists, with the caveat that the tunneling structure is violated locally in space. The identification of these examples with the action depicted in the histogram is described in the caption to Fig. 7.

Next, in Fig. 8 we display the results at a larger value of $\alpha$ where the auxiliary field which couples to charge-density starts to dominate and all saddle points are located at $\chi_{x,\tau}=0$. The latter fact automatically ensures that the structure of the saddle points is completely different from the one at small $\alpha$. The histograms for the same three couplings show that the situation improves and that the construction of the saddle points is now more regular since they are equally spaced in action. A comparison of the $6\times 6$ and $12\times 12$ lattices (Fig. 8(upper panel) and Fig. 8(lower panel)) shows that the number of saddles appearing in the histogram increases with the increasing volume, particularly at larger values of the coupling $U$. However, the general structure of saddle points remains essentially the same. This situation is demonstrated in Fig. 9. Here we have taken a $6\times 6$ lattice with $N_{\tau}=512$ at $U=5\kappa$, as an example (corresponding histogram is shown in the Fig. 4). However, the same field configurations were observed at saddle points for other $U$, $N_{\tau}$ and also at a volume of $12\times 12$. For all histograms, shown here for $\alpha=0.9$ (Fig. 8), the first bar corresponds to the vacuum saddle $\phi_{x,\tau}=\chi_{x,\tau}=0$. The next bar corresponds to the localized field configurations shown in Fig. 9-1(a). These localized features come in two types, differing only in the sign of the $\phi_{x,\tau}$ field: $\phi_{x,\tau}\rightarrow-\phi_{x,\tau}$. We will refer to these structures as “blob” and “anti-blob” in the subsequent discussion. The third bar in the histograms for $\alpha=0.9$ corresponds to the three combinations one can construct out of two of these localized objects: blob-blob, blob-anti-blob and two anti-blobs, where the objects are located at some spatial separation on the lattice. Two examples are shown in Fig. 9-2(a) and 9-3(a). All further saddle points consist of more complicated combinations of increasing number of blobs/anti-blobs that are localized somewhere within the lattice. The single blob shown in Fig. 9-1(a) has an action given by $S_{1}=S_{0}+\Delta S$, where $S_{0}$ is the action of the trivial vacuum. Both configurations in Fig. 9-2(a) and 3 have actions given by $S_{2}\approx S_{0}+2\Delta S$ to a very high precision. It follows that the actions of $n$-blob configurations should be concentrated around $S_{n}=S_{0}+n\Delta S$, with the width of the distribution slightly widening with increasing $n$. This is due to that fact that as the density of blobs increases, they are no longer well-separated and start to interact with each other.

These single and multi-blob configurations have consequences for the fermions, as we attempt to illustrate in Fig. 9. We first define the equal-time fermion Green’s function in position-time representation

[TABLE]

where we have written the expression for particles and an analogous expression exists for the holes. We compute this expression on a given saddle point configuration, for fixed spatial positions $x$ and $y$ as a function of $\tau$. This quantity forms a closed curve in the complex plane due to periodic boundary conditions for the auxiliary fields. Furthermore, for certain locations of the source and sink, this curve exhibits a non-trivial winding around the origin in the complex plane. We define the winding number of the propagator for a given source and sink location as follows

[TABLE]

where in the first equality we have used $z\equiv g(x,y,\tau)$ and $\gamma$ refers to the closed curve swept out by the propagator in the complex plane. For the one-blob configuration in Fig. 9, we have plotted the Green’s function contour for two different sinks, with the source fixed at the center of the blob. In Fig. 9-1(b), the sink is located on the opposite sublattice of the source and shows a non-trivial winding number of $+1$, while in Fig. 9-1(c) the sink is located on the same sublattice of the source and shows a trivial winding of [math]. We thus see that there exists a correlation between fermion winding number, saddle points, and sublattice symmetry.

We have observed that, for the multi-blob configurations, blobs with the same sign lie on the same sublattice while blobs with opposite signs lie on opposite sublattices. The latter is depicted in Fig. 9-2(a) where we have a configuration containing a blob-anti-blob pair, and in Fig. 9-2(b) and 9-2(c) we observe the same correlation between sublattice symmetry and fermion winding number that was observed for the one-blob configuration. However, in  9-2(d), we see a non-trivial winding number of $+2$ where the sink and source were taken to be the centers of the two blobs. A two-blob configuration is depicted in 9-3(a), where again, the winding number is trivial for source and sink on the same sublattice (Fig. 9-3(b)). The winding number is non-trivial and equal to $-1$ for source and sink on different sublattices (Fig. 9-3(c)). We note that winding number $\pm 2$ was not observed for the two-blob configuration. We assume that a similar correlation exists between the winding number and the construction of saddle point configurations with a larger number of (anti-)blobs, and thus the winding number can be used for the classification of saddle points. However, we have left the detailed study of this point to future work.

One expects that the dependence of the thimbles decomposition on the Hubbard coupling should reflect the changing physics in the strongly-coupled phase. The dependence of the real part of the action of the various saddles on the coupling $U$ at half-filling is shown in Fig. 10 for the case where the charge-coupled Hubbard field dominates ($\alpha=0.9$). In order to track the location of the saddles in a continuous manner we have used the GF in the downwards direction after small shifts of the on-site interaction $U$. This means that we start from saddle points at large $U$, then slightly decrease $U\rightarrow U-\delta U$ and search for the new locations of the local minima by starting GF from the old saddles. This procedure is repeated to cover the desired interval of $U$. We have found that sometimes the profiles obtained in this way experience sharp decays into the vacuum saddle. This behavior implies that the corresponding saddle point becomes irrelevant.

Before we proceed further, the last point needs to be clarified. Usually, a thimble and its corresponding saddle point are classified as “relevant” if their intersection number, $k_{\sigma}$, is nonzero. However, things can be different if the so-called Stokes phenomenon occurs. This situation implies that several saddles are now connected by one thimble. Here we consider this situation at half-filling, when there is no sign problem and all relevant thimbles and saddle points are confined within $\mathbb{R}^{N}$. In this case, all eigenvectors of the Hessian matrices, $\Gamma_{\sigma}$, for saddles located within $\mathbb{R}^{N}$ have their components either purely real or imaginary. At the local minimum of the action within $\mathbb{R}^{N}$, which is a relevant saddle point, all real eigenvectors of $\Gamma_{\sigma}$ correspond to positive eigenvalues. However, it can happen that some $\mathcal{N}_{\sigma}>0$ real eigenvectors correspond to negative eigenvalues of $\Gamma_{\sigma}$. This situation is illustrated in Fig. 11. Because the thimbles attached to local minima cover the entire $\mathbb{R}^{N}$, saddles which have at least one real eigenvector corresponding to a negative eigenvalue of $\Gamma_{\sigma}$, do not participate in the sum (7), and thus are irrelevant. Simply counting the intersection points is impossible in this case as $\dim(\mathbb{R}^{N}\cap\mathcal{K}_{\sigma})=\mathcal{N}_{\sigma}>0$ for such saddles. The decay of a saddle, if we start from a slightly shifted field configuration, means that a negative eigenvalue of $\Gamma_{\sigma}$ with a corresponding real eigenvector has appeared. It then follows that this situation indeed corresponds to the transition between relevant and irrelevant status for the saddle point.

Using these ideas we can interpret from Fig. 10 that at small coupling, the trivial vacuum is the only relevant saddle. As we move to larger coupling, multiple non-trivial relevant saddles appear and above $U\approx 3.4\kappa$, we see two of them which are evenly spaced. Thus, we should expect that, at fixed $\alpha$ and large $U$, more and more non-trivial multi-blob saddles become relevant once we approach AFM phase. This interpretation is also supported by the previous histograms cf. Fig. 8.

Finally, in Fig. 12(a) we show how the situation becomes worse as we further increase the parameter $\alpha$, thus suggesting that there exists a “sweet spot” which possesses an advantageous structure for the thimbles decomposition. This regime is illustrated in Fig. 12(b), where even for large lattices ($12\times 12$) at $U=3.8\kappa$ (which corresponds to the AFM phase transition), only the vacuum saddle contributes at $\alpha=0.8$. One can compare this situation with that depicted in the lower panel middle plot of Fig. 8. A more detailed study of this regime is made below, accompanied by the study of saddles points away from half-filling.

V.2 Saddle points at nonzero chemical potential

Away from half-filling one can not rely on the naive application of the GF equations in order to find the saddle points. This is due to the fact that the downward GF ends up on a saddle point only if the initial configuration was exactly on the corresponding thimble. Since we can not generate those configurations (at least without prior knowledge about the saddle points), another method should be employed. We use a procedure similar to Powell’s method to search for local minima. The algorithm is illustrated schematically in Fig. 13(a), for a single complex field. The minimization procedure consists of alternating GF steps for constant imaginary and real parts of the field. The even iterations consist of GF in the downward direction with fixed $\mbox{Re}\Phi_{j}=\Phi^{(R)}_{j}$, where $\Phi_{j}\equiv\Phi^{(R)}_{j}+i\Phi^{(I)}_{j}$ represents both complex auxiliary fields. The flow stops when it reaches the local minimum. The odd iterations consist of upward GF with fixed $\mbox{Im}\Phi_{j}=\Phi^{(I)}_{j}$ and terminate when a local maximum or zero of determinant has been reached, where $\mbox{Re}S\rightarrow\infty$. The convergence can be controlled by monitoring the quantity, $\Sigma_{D,Re/Im}\equiv\sum_{i}|\partial\mbox{Re}S/\partial\Phi^{(R/I)}_{i}|^{2}$ (with the sum running over all sites in the spatial and temporal directions) after each iteration, with $\Sigma_{D,Re}$ reaching the level of numerical precision (typically $10^{-10}$) during even iterations and $\Sigma_{D,Im}$ during odd iterations (assuming the flow did not collide with a zero of determinant). Some examples are shown in Fig. 13(b) and Fig. 13(c). In the former, one can see two examples of the iterations on a $2\times 2$ lattice. Here we see that one converges into the vacuum saddle, which is uniformly shifted into the complex plane ($\mbox{Re}\phi_{x,\tau}=\mbox{Re}\chi_{x,\tau}=\mbox{Im}\chi_{x,\tau}=0$, $\mbox{Im}\phi_{x,\tau}=\phi_{0}$), while the other converges into a non-trivial saddle, which is non-uniform both in space and Euclidean time. The latter figure demonstrates an example for a $6\times 6$ lattice, where the iterations collided with a zero of the determinant on the way, but nevertheless converged afterwards.

Away from half-filling, the initial configurations were generated using a phase quenched HMC, using the algorithms already described in Buividovich et al. (2019). Thus, only the absolute value of $\ln\det(M_{\text{el.}}M_{\text{h.}})$ was taken into account during the Monte Carlo procedure. Usually, the initial configurations are generated along some contour in $\mathbb{C}^{N}$, uniformly shifted from $\mathbb{R}^{N}$, in order to approach the thimble. This is not surprising as we have found that this constant shift into complex space applies to the vacuum saddle at $\mu\neq 0$. The procedure of using a constant shift was performed at $\alpha=0.8$ and $\alpha=0.9$, where the charge-coupled field dominates. If $\alpha=0$, the thimbles and saddles again lie within $\mathbb{R}^{N}$, since both fermionic determinants are real. However, as discussed previously, this property of the fermionic determinants leads to a loss of ergodicity for HMC. Thus, in order to explore the case when the spin-coupled field dominates, we use small $\alpha=10^{-4}$ and generate configurations without a shift into the complex plane. Even such a small, nonzero value of $\alpha$ is enough to restore ergodicity, as one can see in the inset in Fig. 14. This inset shows the history of $\arg S$ during one trajectory in HMC update. If $\alpha=0$, all thimbles have $\cos\arg S=\pm 1$ again due to the fact that $\det(M_{\text{el.}}M_{\text{h.}})\in\mathbb{R}$. Thimbles with different signs are separated by zeros of the determinants, since they are branch points of the logarithm. Here we have a small but nonzero $\alpha$, and thus the $\cos\arg S$ only approaches $\pm 1$. A sharp transition is observed in the inset in Fig. 14 which shows us that the algorithm still can tunnel between different thimbles. This tunneling was, in fact, quite frequent and was observed in more than half of the Monte Carlo updates. This is a further confirmation that the HMC is ergodic.

Another concern regarding our GF procedure is the question of convergence of the alternating iterations. Unfortunately, the procedure we have used does not converge for an arbitrary saddle. The criterion for the convergence of the procedure can be derived from the fact that the distance to the saddle point should decrease after a full round of four iterations. The exact formulation can be expressed in terms of the Hessian matrix $\Gamma=\begin{pmatrix}A&C\\ C&B\end{pmatrix}$. The Hessian is written in terms of $2N_{s}N_{\tau}\times 2N_{s}N_{\tau}$ blocks $A_{i,j}\equiv\partial^{2}\mbox{Re}S/\partial\Phi^{(R)}_{i}\partial\Phi^{(R)}_{j}$, $B_{i,j}\equiv\partial^{2}\mbox{Re}S/\partial\Phi^{(I)}_{i}\partial\Phi^{(I)}_{j}$, and $C_{i,j}\equiv\partial^{2}\mbox{Re}S/\partial\Phi^{(R)}_{i}\partial\Phi^{(I)}_{j}$. Using these matrices, the minimization procedure is guaranteed to converge if both $A$ and $-B$ are positive-definite, and each of the eigenvalues, $\lambda_{i}$, of the matrix $A^{-1}CB^{-1}C$, which characterizes the update of the fields after two subsequent iterations, satisfy

[TABLE]

The latter is actually a constraint on $|\arg\partial_{i}\partial_{j}S|$. If all of the second derivatives are real, $C=0$, and thus $|\lambda_{i}|=0$. If $|\arg\partial_{i}\partial_{j}S|$ increases, with $A$ and $(-B)$ still remaining positive-definite, the thimble in the vicinity of the saddle point is no longer parallel to $\mathbb{R}^{N}$, but starts to “rotate” in complex space. In the $1$D case illustrated in Fig.  13(a), $|\lambda|<1$ simply means that $|\arg\partial^{2}S|_{z_{\sigma}}|<\pi/4$.

We again start from small $\alpha$, which corresponds to a dominant spin-coupled field. In this case, all saddle points are located at $\phi_{x,\tau}=0$ and $\mbox{Im}\chi_{x,\tau}=0$ with their phases $\cos\arg S=\pm 1$. Results are shown in Fig. 14, separately for positive and negative saddles. In general, we observed a very large variety of saddles with non-uniform structures both in space and time. It is almost impossible to characterize them, since their actions form a quasi-continuum distribution. Furthermore, positive and negative saddles almost compensate each other in this case (see histograms in Ulybyshev et al. (2019a)), and thus the part of the residual sign problem stemming from the phase factors in eq. (II) is quite strong. The qualitative explanation for such behavior can be derived from the schematic illustrations in Fig. 15. At half filling, the fermionic determinants for electrons and holes are identical for $\alpha=0.0$, thus the sign problem is absent, but, according to the previously mentioned results (see Fig. 6, 7) we have many thimbles in $\mathbb{R}^{N}$, separated by zeros of the determinants. Once the chemical potential shifts from zero, the two determinants are no longer identical, the domain walls between thimbles are split, and “negative” thimbles immediately appear along the borders between “positive” thimbles. Since we observe a large variety of thimbles at half-filling, the situation can only become worse at $\mu\neq 0$.

Results for large $\alpha$ are shown in Fig. 16. At $\alpha=0.9$, the distribution of both $\mbox{Re}S$ and $\mbox{Im}S$ show the same characteristic behavior as was observed at half filling (Fig. 8) with saddle points equidistantly spaced in action. The difference is that the step in the action is now a complex number. More precisely, the properties of the saddle points can be understood from Fig. 17, which shows the evolution of the saddle points as one goes away from half-filling. We see a continuous evolution of the same system of blobs and anti-blobs. The difference between them in $\mbox{Re}S$ remains constant (Fig. 17(a)), while $\mbox{Im}S$ increases and blob and anti-blob configurations acquire opposite phases. As obvious from our previous discussion, the general rule for the approximate action of saddle point is now $S_{n_{1},n_{2}}=S_{0}+n_{1}\Delta S+n_{2}\overline{\Delta S}$, where $n_{1}$ is number of blobs and $n_{2}$ is the number of anti-blobs in the configuration. Thus, the actions of the saddle points form a triangle in complex space as can be observed in Fig. 16(a). Another interesting consequence is that not only the vacuum, but also configurations with an equal number of blobs and anti blobs have $\mbox{Im}S=0$, which effectively decreases the complexity of the sign problem.

As we go to finite chemical potential we can also attempt to visualize the saddle point configurations. Unlike Fig. 9, where the field is real, we plot the configuration in the complex plane with each point representing the value of the field at a given lattice site for $\alpha=0.9$ (Fig. 18). In these plots, the vacuum configuration is the trivial vacuum $\phi=\chi=0$ with an added constant, volume-independent shift of the imaginary part. In Fig. 18(a), we plot the configuration of a single anti-blob. The collection of points which extends furthest away from the vacuum all come from the localized region of space and time surrounding the center of the anti-blob. In Fig. 18(a), where we display a blob-anti-blob pair, we see that we have have two such collections of points extending away from the vacuum in opposite directions. Each collection comes from the region of space surrounding the centers of the blob and anti-blob respectively. An illustration of the effect of finite-volume on the non-trivial saddles is depicted in Fig. 18(c) for a $2\times 2$ lattice where, for a single blob, the structure of the distribution of the field values is distorted as compared to Fig. 18(a).

In the sweet spot regime at $\alpha=0.8$, we detect again only the vacuum saddle (Fig. 16(b)). In principle, such situation should be very beneficial for the thimbles decomposition, since the fluctuations of $\mbox{Im}S$ can be made arbitrarily small. Also, it should improve the ergodicity of the Monte Carlo process, since the integration manifold is no more divided into disconnected domains. However, we should stress that unlike the $\mu=0$ case, the distributions away from half-filling (Fig. 14 and Fig. 16) are exact only for $\alpha=10^{-4}$, since we are quite close to the thimble in this case. For $\alpha=0.8$ and $\alpha=0.9$, the histograms are only approximate as the initial configurations for the iterations approach the thimble, but do not lie exactly on it. Furthermore, “vertically oriented” saddles, which do not satisfy the convergence condition (70) can be missed. However, subsequent QMC calculations support the conclusion that the regime around $\alpha=0.8$ is indeed better for simulations than $\alpha\rightarrow 1.0$.

The optimal regime around $\alpha=0.8$ is studied in Fig. 19. We start from half-filling in Fig. 19(c), The lower boundary of this region in $\alpha$ corresponds to the splitting of the vacuum saddle into two mean-field saddles. This splitting is observed by launching GF from a slightly perturbed vacuum (Gaussian noise is added to $\phi_{x,\tau}$ and $\chi_{x,\tau}$). If $\chi_{x,\tau}$ returns to zero, the vacuum is stable, otherwise the final value of $\chi_{x,\tau}$ is non-zero, since the flow arrives at the mean-field saddle point. This is what we see in the $\chi$-profiles for the mean-field saddles in Fig. 19(c) both for $6\times 6$ and $12\times 12$ lattices. The jump upwards corresponds to the appearance of the mean-field saddle and marks the lower boundary of the optimal regime. The upper boundary is determined by the decay of the nontrivial saddles into vacuum, analogous to what was observed in Fig. 10. We use the symmetry, $S(\phi_{x,\tau},\chi_{x,\tau})=S(\phi_{x,\tau},-\chi_{x,\tau})$, and the fact that the saddles are located at $\chi_{x,\tau}=0$ for large $\alpha$. The Hessian matrix is block-diagonal in this case as $\partial^{2}S/\partial\chi_{x,\tau}\partial\phi_{x,\tau}=0$. Because it is enough to find at least one instability (negative eigenvalue of $\Gamma$, with real eigenvector), we can study the relevance of saddles separately for $\phi$- and $\chi$-directions. It can be done in two different ways, which should lead to identical results. The first is the calculation fully analogous to the one made for Fig. 10: we start from large $\alpha$, slightly decrease it and launch GF in downwards direction for both fields, $\phi_{x,\tau}$ and $\chi_{x,\tau}$. At sufficiently small $\alpha$ we see the decay of the saddle into vacuum (example of the flow history for the fields at one particular site is shown in Fig. 19(b)). Thus, we can plot the dependence of $\phi_{x,\tau}$ after the flow on $\alpha$ and a sharp drop to zero will mark the transition of a non-trivial saddle into an irrelevant one. This approach, however, does not allow us to understand, in which block of the Hessian matrix does the negative, “unstable” eigenvalue appear. Alternatively, one can first use GF, restricted to the $\phi$ fields in order to find the saddle after a small shift of $\alpha$. No instability was found in this case, and the non-trivial saddles can be found for all values of $\alpha$. Finally, we use these saddles, add noise to the $\chi$ fields, and launch the GF, restricted to $\chi$ fields for these configurations. In this case, the instability manifests itself in a finite value of the $\chi$ fields after the flow is performed (the configuration does not return to the initial saddle located at $\chi_{x,\tau}=0$). Thus, we can plot $\chi$ field after the flow and the sudden appearance of a nonzero value signals the transition of the non-trivial, relevant saddle into an irrelevant one. Both approaches are demonstrated in Fig. 19(c) for $6\times 6$ and $12\times 12$ lattice. At $\alpha\approx 0.89$, the final value of $\chi$ after the flow for $\chi$ fields jumps upwards. Simultaneously, the final value of $\phi$ fields in the full flow goes down to the level of numerical errors (typically around $10^{-10}$). This signals an instability in the $\chi$-channel, and thus the non-trivial saddle becomes irrelevant. Remarkably, the results depend neither on the type of saddle point nor on lattice size. We attribute this property to the localized nature of non-trivial saddle points at large $\alpha$. An important observation is that the width of the ”optimal regime” grows with increasing system size, since the lower border shifts to smaller $\alpha$ for the $12\times 12$ lattice. This lends strong support to the existence of this regime in the thermodynamic limit.

A similar set of calculations was performed for $\mu=\kappa$, where we have used GF restricted to $\mbox{Re}\chi$. The plot in Fig. 19(a) demonstrates how the flow switches from the stable regime to eventual decay. We essentially observe the same behavior for the non-vacuum saddles at large $\alpha$ (Fig. 19(d)): all non-trivial field configurations are unstable in the $\mbox{Re}\chi$ direction if $\alpha<0.89$. This suggests that at $\alpha>0.89$, the non-vacuum saddles shift into complex space with increasing $\mu$, remaining relevant, while at $\alpha<0.89$ they begin from irrelevant ones at $\mu=0$ and move into complex space for $\mu\neq 0$, remaining irrelevant. These situations are depicted in the schematic drawings in Fig. 20, where both relevant and irrelevant saddle points are shown at half-filling and at $\mu\neq 0$. Another possibility is that the saddles acquire a more “vertical” orientation with decreasing $\alpha$ (as shown in Fig. 21). It can happen that the GF along $\mbox{Re}\chi$ can take us away from zero also in this case. However, there are additional arguments against this based on the results of our HMC over manifolds in complex space at different $\alpha$ described in the next section.

VI HMC with gradient flow

In order to check the conclusions concerning the structure of the thimbles decomposition, made in the previous section, we performed several Monte Carlo calculations, with the integration manifolds shifted in complex space towards the thimbles. Following Alexandru et al. (2017), the sequence of deformations of the integration contour can be summarized by

[TABLE]

The general idea behind these deformations is shown in the schematic illustration of Fig. 21. First, we perform a uniform shift into the complex space, but only for the charge-coupled field: $\phi\rightarrow\phi+i\phi_{0}$. We work at large $\alpha$, and this uniform shift moves onto the thimble attached to the vacuum saddle (see Fig. 18) in the Gaussian approximation. Below we will denote this shift as $\Phi\rightarrow\Phi+i\,\Phi_{0}$. A further shift is made using the GF equations. The quantity $\tilde{\Phi}\in\mathbb{C}^{N}$ is the result of the evolution of the field determined by (6), starting from the Gaussian thimble $\Phi+i\,\Phi_{0},\Phi\in\mathbb{R}$ with flow time $\mathcal{T}$. We should stress that this algorithm does not follow thimble exactly, it only approaches it in the limit of infinitely large flow time. Thus, the third source of the residual sign problem, the residual fluctuations of $\mbox{Im}S$ appear in calculations. The complex-valued Jacobian of the transformation, $J=\mathcal{D}\tilde{\Phi}/\mathcal{D}\Phi$, also appears in the integral at this stage. The flow time plays a dual role in this transformation. First, it defines, how close we can approach the thimble, and thus it regulates the residual fluctuations of $\mbox{Im}S$. Second, if the flow time is too large, the flow lines can reach zeros of the determinant, which separate thimbles (see Fig. 21). In this case, the integration domain for the $\Phi$ fields is again split into separate regions and as a result the Monte Carlo process can hardly be expected to be ergodic. And finally, another contribution to the residual sign problem comes from the Jacobian, especially in the case of “vertically” oriented thimbles, as shown in Fig. 21.

We use the following strategy to sample the partition function (71): 1) The Jacobian is not taken into account in the Markov process employed to generate field configurations $\tilde{\Phi}$, and is left for the final reweighting; 2) The fields $\Phi$ are generated using HMC, according to the distribution $e^{-\mbox{Re}S[\tilde{\Phi}(\Phi+\Phi_{0})]}$; 3) The fields $\tilde{\Phi}$ are computed through the gradient flow evolution. Several examples from the second stage of the process are shown in Fig. 22 for one particular site of the lattice. The second stage requires an additional comment. HMC employs global updates of the fields, using molecular dynamics (MD) governed by $H=\frac{1}{2}\sum_{i}p_{i}^{2}+\mbox{Re}S[\tilde{\Phi}(\Phi)]$, where an artificial momentum, $p_{i}$, is introduced for each Hubbard field $\Phi_{i}$. In order to solve Hamilton’s equations, we need to compute the derivative $\partial\mbox{Re}S[\tilde{\Phi}(\Phi+\Phi_{0})]/\partial\Phi_{i}$. We calculate this quantity by shifting the initial fields $\Phi_{i}\rightarrow\Phi_{i}+\Delta\Phi$ and solving the GF equations for each shift. Examples of such calculations are shown in Fig. 23. These plots show that we can compute these derivatives to sufficient accuracy increasing the precision of the numerical integration for GF equations. These calculations also give us immediate access to all elements of the Jacobian. In practice, however, we need the $\det J$ only after the accept-reject step of the HMC procedure, while the trajectory typically consists of $\mathcal{O}(10^{2})$ steps. Thus, the calculation of the Jacobian plays only a subdominant role in computational efforts. The overall scaling of the method is $C_{1}\tilde{N}_{MD}\tilde{N}_{GF}N_{s}^{4}N_{\tau}^{2}+C_{2}N_{s}^{3}N_{\tau}^{3}$, where the first term corresponds to the HMC procedure used to generate field configurations and the second term accounts for the calculation of $\det J$ at the end of the trajectory. Here $\tilde{N}_{MD}$ refers to the number of steps in a MD trajectory which is typically $\mathcal{O}(10^{2})$, $\tilde{N}_{GF}$ refers to the number of steps in the integrator for GF equations which is typically $O(10^{1})$, and $C_{1}$ and $C_{2}$ denote volume-independent constants. In what follows, we will refer to this algorithm as HMC-GF. Several examples of configurations of the $\tilde{\Phi}$ fields, generated with this algorithm, are show in Fig. 24.

The Jacobian is left for the final reweighting, and thus the observables are computed using the following expression

[TABLE]

where the residual fluctuations of $\mbox{Im}S$ are also taken into account. The brackets $\langle\rangle$ denote the averaging over configurations generated with HMC-GF. We also take into account symmetries of the action in order to further improve the ergodicity of our set of field configurations, generated with HMC-GF

[TABLE]

The following metrics are used to estimate the severity of the sign problem: $\langle\cos(\mbox{Im}S)\rangle$ and $\langle\cos(\mbox{Im}\ln\det J)\rangle$ for configurations and the Jacobian respectively, and the joint sign $\langle\Sigma_{G}\rangle=\langle\cos(\mbox{Im}(-S+\ln\det J))\rangle$. The first metric characterizes the part of the residual sign problem which stems from the fact that the sequence of shifts (71) does not follow thimble exactly. The second metric characterizes the part of the residual sign problem which stems from the fluctuations of complex measure during integration over curved manifold in complex space. The last metric characterizes the entire residual sign problem. We also estimate the strength of the fluctuations of the Jacobian by computing $D_{J}$, the dispersion of $\mbox{Re}(\ln\det J)$.

The following choice is made for the parameters of the simulations: $2\times 2$ lattice ($N_{s}=8$), $N_{\tau}=256$, $U=2\kappa$, $\mu=\kappa$, $\beta=20$. This lattice is small enough to make a comparison with finite-temperature ED possible, but large enough to host non-trivial saddle points at large $\alpha$ (see Fig. 18(c)). Their form is only slightly different from the ones appearing at larger lattice sizes. These saddles also experience decay along the $\mbox{Re}\chi$ direction at $\alpha\approx 0.8$, similar to the $6\times 6$ and $12\times 12$ lattices studied above. Thus we can say that such a small lattice can in fact model the properties of the saddle points even at thermodynamic limit. On the other hand, we find that $N_{\tau}=256$ is large enough to probe both the low-temperature regime as well as the continuum limit in Euclidean time simultaneously. We further note that the state-of-the-art QMC algorithm for condensed matter systems, BSS-QMC, taken from the ALF package Bercx et al. (2017), experiences exponential decay of the average sign at these parameters, even in the optimal regime where the discrete auxiliary field is coupled to spin. It is thus apparent that the sign problem is already strong in this regime. We have also probed two different values of $\alpha$: $\alpha=1.0$, so that only the charge-coupled field $\phi_{x,\tau}$ participates in the integral, and $\alpha=0.8$ in order to probe the “optimal regime”, where only the vacuum saddle point was detected.

Our results for the computed observables are displayed in the Tab. 1 while the study of the sign problem is summarized in Tab. 2. We compute the kinetic energy, $\langle\,\hat{{K}}\,\rangle$, and the nearest-neighbor correlation function for the first component of spin $\langle\,\hat{S}^{(1)}_{x}\hat{S}^{(1)}_{y}\,\rangle$. Results at $\alpha=1.0$ substantially deviate from ED, while at $\alpha=0.8$ the results of HMC calculation are in agreement with ED. This seems to imply that at $\alpha=1.0$ ergodicity issues indeed appear as there are several relevant thimbles and thus GF collides with zeros of the determinant. Unfortunately, HMC can not tunnel through the barrier separating two thimbles in such situations. At $\alpha=0.8$, however, ergodicity is restored. Moreover, we do not observe the growth of the fluctuations of the Jacobian, which should appear if GF approaches “vertically” oriented thimbles (Fig. 21). This tells us that the thimbles attached to the non-vacuum saddles indeed become irrelevant or they are bypassed by the integration manifold constructed by GF. In both cases, these non-vacuum saddles are effectively unimportant at $\alpha\approx 0.8$. Noting that there exists a striking resemblance between the saddle point structure on the $2\times 2$ lattice and on larger lattices at large $\alpha$, table 1 lends very strong evidence that at $\alpha\approx 0.8$, there exists an optimal regime with only one important thimble surviving in the full decomposition even for large spatial volumes.

We also collected smaller statistics for the same lattice, but with $\beta=30$ with $N_{\tau}=384$ ($\alpha=0.8$). We increased the flow time from $\mathcal{T}=0.05$ for $\beta=20$ to $\mathcal{T}=0.08$ for $\beta=30$ to keep the residual sign problem roughly the same ($\langle\Sigma_{G}\rangle=0.785\pm 0.021$ for $\beta=30$ with $\mathcal{T}=0.08$ in comparison with $\langle\Sigma_{G}\rangle=0.644\pm 0.028$ for $\beta=20$ with $\mathcal{T}=0.05$). As we have already mentioned introducing the Hubbard model, in this case, where only local interactions are present, we have the special possibility to perform calculations with BSS-QMC, using discrete auxiliary fields coupled with spin density. The sign problem in such set up seems to be several orders of magnitude smaller than in the approach with continuous auxiliary fields. Thus, in the Lefschetz thimbles formalism, where we can use only continuous fields, we start from a much stronger sign problem. However, despite this fact we can keep the average sign much higher than in BSS-QMC (see the Tab. 2 and the Fig. 25b). Moreover, our calculations demonstrate another important property of the HMC-GF algorithm: despite the increased difficulty of the initial sign problem (average sign in BSS-QMC drops down in several times between $\beta=20$ and $\beta=30$), we can maintain the residual sign problem under control by increasing the flow time and moving the integration contour closer to the thimble: residual sign $\langle\Sigma_{G}\rangle$ even slightly increases for $\beta=30$ in HMC-GF. We expect similar behaviour for the dependence on lattice size $N_{S}$. The detailed study of the dependence on $N_{S}$ is left for upcoming papers, because we still need to improve our algorithms for constructions of the integration manifold for reasons listed below, and the estimates for the residual problem will likely change for modified algorithm.

Unfortunately, larger average sign is still not enough to beat the better scaling of BSS-QMC: its scaling is given by $N_{s}^{3}N_{\tau}$, which is substantially better than the dominant term $N_{s}^{4}N_{\tau}^{2}$ in the scaling of the HMC-GF. Thus, at these particular parameters, we can still use naive reweighting (generating $\sim 10^{5}$ configurations) to achieve error bars which are smaller than the ones computed in HMC-GF with $\sim 10^{3}$ configurations. The dispersion of the Jacobian is also noticeable but in general does not cause large problems. We simply need larger statistics to compensate for these additional fluctuations. We have determined $D_{J}=1.157$ for HMC-GF with $\alpha=1.0$ and $D_{J}=1.011$ for HMC-GF with $\alpha=0.8$. We noticed that the properties of the Jacobian become worse at $\beta=30$: $D_{J}=1.68$ and $\langle\cos\arg J\rangle=0.823\pm 0.018$ in this case (compare it with $\alpha=0.8$ case in Tab. 2). Thus, at very low temperatures and, possibly, at larger system sizes, fluctuations of the Jacobian might become a problem. Right now we leave this potential problems for further study.

As a final statement, we stress that from the point of view of lattice gauge theories, we should compare only with reweighting made with HMC running for continuous auxiliary fields, since there are no such cheap alternatives such as BSS-QMC for gauge theories or for the case of condensed matter systems with long-range interactions Hohenadler et al. (2014). The results of these benchmarks are shown in the Fig. 25(a). One can see that the improvement in comparison with plain HMC for continuous fields running inside $\mathbb{R}^{N}$ is quite huge. In fact, the van Hove singularity can not be reached with that method at all, as the average sign decays too quickly with increasing chemical potential. The plot essentially shows that even the present algorithm would allow us to go to the region completely inaccessible for standard HMC, thus it should already be the method of choice if we are simulating, for instance, the Hubbard-Coulomb model. Also, it suggests the feasibility of using the Lefschetz thimbles method for lattice gauge theories which suffer from the sign problem, since with HMC-GF in our studies of the Hubbard model we could reach the ratio $\mu/T=30$, at least for small lattices, while the sign problem typically becomes quite severe in lattice QCD at $\mu/T\sim 1$. However, the possibility of reaching such values of the parameters for QCD is at the moment purely speculative and should be checked explicitly in future dedicated simulations.

VII Conclusion

We proposed a set of new algorithms to study the properties of the thimbles decomposition for a lattice fermionic model approaching the thermodynamic limit. The method was tested for the case of the Hubbard model on the hexagonal lattice at various values of the interaction strength and chemical potential. Using this method, we could find the exact saddle points in the path integral formulation of the Hubbard model with different forms of the HS transformation, combining both the spin- and charge-coupled auxiliary fields with different weights. This represents the main physical result of the paper. At half-filling, in the case where the spin-coupled field dominates, we have observed a large number of different saddle points with an instanton-like structure. Away from half-filling, this regime develops a strong sign problem due to a large amount of thimbles with opposite phases which almost compensate each other. If the charge-coupled field is dominant, the structure of the saddle points is more regular, as they are built from two basic building blocks (localized field configurations) both at zero chemical potential and away from half-filling. We have also observed an intermediate, “optimal” regime, where our method detected only one relevant saddle, even for a volume of $12\times 12\times 256$. These results show that the thimble decomposition for the Hubbard model strongly depends on the form of the HS transformation and that a solution to the sign problem using Lefschetz thimbles must properly take into account the saddle point structure. Consequently, the residual sign problem also depends strongly on the aforementioned decomposition of the fields.

In order to check these findings, we have performed several benchmark HMC simulations for lattices with $N_{s}=8$ and $N_{\tau}$ as large as $384$. In the “optimal” regime, we were able to reach an agreement with ED, which suggests that the ergodicity issues and the sign problem can be simultaneously weakened. We were also able to obtain an average sign which is much larger than that obtained with conventional QMC techniques. At the moment the thimbles algorithms is still too costly (mainly due to the aforementioned $N_{s}^{4}N_{\tau}^{2}$ scaling with the system size), but the fact that the average sign can be kept higher than that of BSS-QMC gives us hope that we can construct a superior technique if we manage to decrease the computational costs associated with the integration over curved manifolds in complex space. Furthermore, the proposed set of algorithms can be helpful for the study of Lefshetz thimbles decomposition in strongly-interacting lattice models with fermions. The latter might be important for the better understanding of the corresponding physics through more accurate saddle point approximations of the path integral.

Acknowledgements.

MU would like to thank Prof. F. Assaad for enlightening discussions. SZ acknowledges stimulating discussions with Prof. K. Orginos. MU is supported by the DFG under grant AS120/15-1. CW is supported by the University of Kent, School of Physical Sciences. SZ acknowledges the support of the DFG Collaborative Research Centre SFB 1225 (ISOQUANT). This work was granted access to the HPC resources of CINES and IDRIS in France under the allocation 52271 made by GENCI. We are grateful to these computing centers for their constant help. We are grateful to the UK Materials and Molecular Modelling Hub for computational resources, which is partially funded by EPSRC (EP/P020194/1). This work was also partially supported by the HPC Center of Champagne-Ardenne ROMEO.