Applying the tempered Lefschetz thimble method to the Hubbard model away from half-filling

TL;DR

This paper advances the tempered Lefschetz thimble method to effectively address the sign problem in quantum Monte Carlo simulations of the Hubbard model away from half-filling, achieving accurate results on small lattices.

Contribution

It develops a refined algorithm for precise expectation value estimation and demonstrates its effectiveness on the Hubbard model.

Findings

The method accurately reproduces exact results for the Hubbard model.

The algorithm effectively mitigates the sign problem in small lattice simulations.

Global equilibrium and sample sufficiency criteria improve estimation accuracy.

Abstract

The tempered Lefschetz thimble method is a parallel-tempering algorithm towards solving the numerical sign problem. It uses the flow time of the gradient flow as a tempering parameter and is expected to tame both the sign and multimodal problems simultaneously. In this paper, we further develop the algorithm so that the expectation values can be estimated precisely with a criterion ensuring global equilibrium and the sufficiency of the sample size. To demonstrate that this algorithm works well, we apply it to the quantum Monte Carlo simulation of the Hubbard model away from half-filling on a two-dimensional lattice of small size, and show that the numerical results agree nicely with exact values.

Tables

Table 1. Table 1: TLTM parameters and the results

$\beta \mu$	1	2	3	4	5	6	7	8
$T/\left(\beta \mu \right)$	1/10	1/10	1/10	1/10	1/10	1/11	1/11	1/12
$A$	8	9	10	11	11	11	11	11
$t_c/T$	0.7	0.6	0.6	0.6	0.6	0.5	0.5	0.5
$a_c$	5	5	6	6	6	5	5	5
$N_{\mathrm{swap}}$	10	10	12	12	12	12	12	12
$N_{\mathrm{conf}}$	5k	5k	10k	10k	15k	25k	15k	15k
$a_{\min }$	0	0	2	3	5	5	6	7
${\chi }^2/\mathrm{DOF}$	0.53	0.43	0.47	0.12	0.45	0.39	1.07	0.72

Equations

$$⟨\mathcal{O}(x){⟩}_S\equiv \frac{{\int }_{{\mathbb{R}}^N}dx{e}^{-S(x)}\mathcal{O}(x)}{{\int }_{{\mathbb{R}}^N}dx{e}^{-S(x)}}.$$

$${\dot{z}}_t^i$$

$$⟨\mathcal{O}(x){⟩}_S=\frac{{\int }_{{\Sigma }_t}dz{e}^{-S(z)}\mathcal{O}(z)}{{\int }_{{\Sigma }_t}dz{e}^{-S(z)}}({\Sigma }_t\equiv z_t({\mathbb{R}}^N)),$$

$$⟨\mathcal{O}(x){⟩}_S$$

$$\displaystyle=\frac{\bigl{\langle}e^{i\theta_{t}(x)}\mathcal{O}(z_{t}(x))\bigr{\rangle}_{S^{\rm eff}_{t}}}{\bigl{\langle}e^{i\theta_{t}(x)}\bigr{\rangle}_{S^{\rm eff}_{t}}}.$$

$$e^{-S_t^{\mathrm{eff}}(x)}$$

$$e^{i{\theta }_t(x)}$$

$${\dot{J}}_t$$

$$\prod _a\exp [-S_{t_a}^{\mathrm{eff}}(x_a)].$$

$$\displaystyle\frac{\bigl{\langle}e^{i\theta_{t_{a}}(x)}\mathcal{O}(z_{t_{a}}(x))\bigr{\rangle}_{S^{\rm eff}_{t_{a}}}}{\bigl{\langle}e^{i\theta_{t_{a}}(x)}\bigr{\rangle}_{S^{\rm eff}_{t_{a}}}}$$

$$\approx \frac{{\sum }_{k=1}^{N_{\mathrm{conf}}}\exp [i{\theta }_{t_a}(x_a^{(k)})]\mathcal{O}(z_{t_a}(x_a^{(k)}))}{{\sum }_{k=1}^{N_{\mathrm{conf}}}\exp [i{\theta }_{t_a}(x_a^{(k)})]}\equiv {\bar{\mathcal{O}}}_a.$$

$$H$$

$$+U\sum _{\text{x}}(n_{\text{x},\uparrow }-1/2)(n_{\text{x},\downarrow }-1/2).$$

$$⟨n{⟩}_S$$

$$e^{-S[\varphi ]}$$

$$M^{a/b}[\varphi ]$$

$$n[\varphi ]$$

$$\Delta S_a$$

$$\displaystyle=({\rm Re}\,S(z^{\prime}_{a})-\ln\,\bigl{|}\det\!J^{\prime}_{a}\bigr{|})-({\rm Re}\,S(z_{a})-\ln\,\bigl{|}\det\!J_{a}\bigr{|}).$$

$$\displaystyle w_{a}(x,y)={\rm min}\Bigl{(}1,\,e^{-S^{\rm eff}_{t_{a}}(y)-S^{\rm eff}_{t_{a+1}}(x)+S^{\rm eff}_{t_{a}}(x)+S^{\rm eff}_{t_{a+1}}(y)}\Bigr{)}.$$

$$w_a(x,y)e^{-S_{t_a}^{\mathrm{eff}}(x)-S_{t_{a+1}}^{\mathrm{eff}}(y)}=w_a(y,x)e^{-S_{t_a}^{\mathrm{eff}}(y)-S_{t_{a+1}}^{\mathrm{eff}}(x)}.$$

$$\{(x_a^{(k)},z_a^{(k)},J_a^{(k)}){\}}_{k=1,2,\dots ,N_{\mathrm{conf}}},$$

$$\displaystyle\frac{\bigl{\langle}e^{i\theta_{t_{a}}(x)}\mathcal{O}(z_{t_{a}}(x))\bigr{\rangle}_{S^{\rm eff}_{t_{a}}}}{\bigl{\langle}e^{i\theta_{t_{a}}(x)}\bigr{\rangle}_{S^{\rm eff}_{t_{a}}}}$$

$$[{\theta }_a^{(k)}\equiv {\theta }_{t_a}(x_a^{(k)})].$$

$$H_1$$

$$H_2$$

$$=(U/2)\sum _{\text{x}}(n_{\text{x}}^a-n_{\text{x}}^b{)}^2-N_{s}U/4.$$

$$n\equiv (1/N_s)\sum _{\text{x}}(n_{\text{x},\uparrow }+n_{\text{x},\downarrow }-1)=(1/N_s)\sum _{\text{x}}(n_{\text{x}}^a-n_{\text{x}}^b).$$

$$e^{-\beta H}=(e^{-ϵ(H_1+H_2)}{)}^{N_{\tau }}\simeq (e^{-ϵH_1}{e}^{-ϵH_2}{)}^{N_{\tau }},$$

$$e^{-ϵH_2}$$

$$=e^{N_sϵU/4}\prod _{\text{x}}\int \frac{d{\varphi }_{\ell ,\text{x}}}{\sqrt{2\pi }}{e}^{-(1/2){\varphi }_{\ell ,\text{x}}^2+i\sqrt{ϵU}{\varphi }_{\ell ,\text{x}}(n_{\text{x}}^a-n_{\text{x}}^b)}.$$

$$Z_{\mathrm{QMC}}$$

$$=(e^{ϵU/4}/\sqrt{2\pi }{)}^{N_{\tau }{N}_s}\int [d\varphi ]e^{-S[\varphi ]}.$$

$$e^{-S[\varphi ]}$$

$$\times {\mathrm{tr}}_b\prod _{\ell }{e}^{ϵ{\sum }_{\text{x},\text{y}}(\kappa K-\mu 1{)}_{\text{x}\text{y}}{b}_{\text{x}}^{†}{b}_{\text{y}}}{e}^{{\sum }_{\text{x}}(-i\sqrt{ϵU}{\varphi }_{\ell ,\text{x}})b_{\text{x}}^{†}{b}_{\text{x}}},$$

$$e^A{e}^B=e^C\Rightarrow e^{\hat{A}}{e}^{\hat{B}}=e^{\hat{C}},$$

$$\mathrm{tr}{e}^{\hat{A}}=\det (1+e^A).$$

$$e^{-S[\varphi ]}$$

Full text

sign problem, Hubbard model, parallel tempering

††preprint: KUNS-2764

I Introduction

The sign problem is one of the major obstacles when performing numerical calculations in various fields of physics. Typical examples include finite density QCD [1], quantum Monte Carlo (QMC) calculations of quantum statistical systems [2, 3, 4], and the numerical simulations of real-time quantum field theories.

Among a variety of approaches, two algorithms have taken attention as potential candidates to generically solve the sign problem for systems with complex action; one is the complex Langevin method [5], and the other is a class of algorithms utilizing the Lefschetz thimbles [6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. Although both the algorithms make use of complexification of variables and analytic continuation of integrands, their methodologies are fairly different; the former algorithm attempts to replace the complex Boltzmann weight by a real positive weight defined in the whole complex space, while the latter deforms the integration region in the complex space so as to reduce the phase oscillation. At this stage, each algorithm has its own advantage and disadvantage. The former is advantageous in that it is relatively fast with computational cost $O(N)$ ($N$: the degrees of freedom), but it suffers from the so-called wrong convergence problem [16, 17, 18, 19]. The latter is generally free from the wrong convergence problem if only a single thimble is relevant in evaluating the expectation values of physical observables of interest. The disadvantage is its expensive numerical cost, which is $O(N^{3})$ because of the need to calculate the Jacobian determinant. When multiple thimbles are relevant, one needs to take care of the multimodality of the distribution. The tempered Lefschetz thimble method (TLTM) was thus proposed in [14] to tame both the sign and multimodal problems simultaneously, where the system is tempered by the flow time of the antiholomorphic gradient flow (see also [15] for a similar idea).

In this paper, we further develop the TLTM, proposing an algorithm which allows the precise estimation of expectation values with a criterion ensuring global equilibrium and the sufficiency of the sample size. The key is the use of the fact that the expectation values should be the same for all flow times. To demonstrate that this algorithm works well, we apply it to the QMC simulation of the Hubbard model away from half-filling.

The application of Lefschetz thimble methods to the Hubbard model has already been considered by several groups [20, 21, 22] (see also [23, 24] for recent study), and the relevance of the contributions from multiple thimbles has been reported. In this paper, we consider a two-dimensional periodic square lattice of size $N_{s}=2\times 2$ with the inverse temperature decomposed to $N_{\tau}=5$ pieces, and numerically evaluate the expectation values of observables as functions of the chemical potential with other parameters fixed to some values. We show that the TLTM (the implementation of tempering combined with the above algorithm for precise estimation) give results that agree nicely with exact values, simultaneously resolving the sign and multimodal problems.

We comment that the extent of seriousness of the sign problem in the QMC simulation of the Hubbard model depends heavily on the choice of the Hubbard-Stratonovich variables. In this paper, in order to apply the Lefschetz thimble method, we exclusively consider a Gaussian Hubbard-Stratonovich variable that leads to a complex action. There the sign problem is actually severe as we will see below, and one needs to seriously consider a dilemma between the sign and multimodal problems, which can be solved by the TLTM as stated above. However, the temporal size considered here is still small ($N_{\tau}=5$), and for such a high temperature regime one can resort to other methods than the Lefschetz thimble methods with a different type of Hubbard-Stratonovich variables (see discussions in section III).

This paper is organized as follows. In section II after briefly reviewing the TLTM [14], we give a new algorithm which allows the precise estimation of expectation values with a criterion ensuring global equilibrium and the sufficiency of the sample size. This algorithm is applied to the Hubbard model in section III, and we discuss about the obtained numerical results. We there also make a comment on the sign averages obtained by other methods. Section IV is devoted to conclusion and outlook. Five appendices are given for more detailed discussions on various topics.

II Tempered Lefschetz thimble method

Let $x=(x^{i})\in{\mathbb{R}}^{N}$ be a real $N$-dimensional dynamical variable with action $S(x)$ which may take complex values. Our main concern is to estimate the expectation values

[TABLE]

We assume that $e^{-S(z)}$ and $e^{-S(z)}\,\mathcal{O}(z)$ are entire functions over ${\mathbb{C}}^{N}$ when $x$ is complexified to $z=(z^{i})\in{\mathbb{C}}^{N}$. Then, due to Cauchy’s theorem for higher dimensions, the right-hand side does not change under continuous deformations of the integration region as long as the boundary at infinity is kept fixed so that the integrals converge. The sign problem will get reduced if ${\rm Im}\,S(z)$ is almost constant on the new integration region.

In [11, 11, 12, 13, 14, 15] such a deformation $x\to z_{t}(x)$ ($t\geq 0$) is made according to the antiholomorphic gradient flow:

[TABLE]

Equation (1) can then be rewritten as

[TABLE]

which can be further rewritten as a ratio of reweighted integrals over ${\mathbb{R}}^{N}$ by using the Jacobian matrix J_{t}(x)\equiv\bigl{(}\partial z_{t}^{i}(x)/\partial x^{j}\bigr{)} [11]:

[TABLE]

Here, $S^{\rm eff}_{t}(x)$ and $\theta_{t}(x)$ are defined by

[TABLE]

and $J_{t}(x)$ obeys the following differential equation [11] (see also footnote 2 of [14]):

[TABLE]

with $H(z)\equiv(\partial_{i}\partial_{j}S(z))$. Under the flow (2), the action changes as (d/dt)S(z_{t}(x))=\bigl{|}\partial_{i}S(z_{t}(x))\bigr{|}^{2}\geq 0, and thus ${\rm Re}\,S(z_{t}(x))$ increases except at the critical points $z_{\ast}$ ($\partial_{i}S(z_{\ast})=0$), while ${\rm Im}\,S(z_{t}(x))$ is kept constant. In particular, in the limit $t\to\infty$, the deformed region will approach a union of $N$-dimensional submanifolds (Lefschetz thimbles) on each of which ${\rm Im}\,S(z)$ is constant, and thus the sign problem is expected to disappear there (except for a possible residual sign problem arising from the phase of the complex measure $dz$ and a possible global sign problem caused by phase cancellations among different thimbles). However, in the Monte Carlo calculation one cannot take the $t\to\infty$ limit naïvely, because the potential barriers between different thimbles become infinitely high so that the whole configuration space cannot be explored sufficiently. This multimodality of distribution makes the Monte Carlo calculation impractical, especially when contributions from more than one thimble are relevant to estimating expectation values. A key proposal in [12] is to use a finite value of flow time that is large enough to avoid the sign problem but simultaneously is not too large so that the exploration in the configuration space is still possible. However, it is a difficult task to find such value of flow time in a systematic way, as we will discuss at the end of section III and in Appendix E.

The TLTM [14] is a tempering algorithm that uses the flow time as a tempering parameter. There, the global relaxation of the multimodal distribution is prompted by enabling configurations around different modes to easily communicate through transitions in ensembles at smaller flow times. Among other possible tempering algorithms, the parallel tempering algorithm [25, 26] (also known as the replica exchange MCMC method; see [27] for a review) is adopted in the TLTM [14] because it does not need to introduce the probability weight factors of ensembles at various flow times and because most of relevant steps can be done in parallel processes.

In the TLTM (see Appendix A for the summary of the algorithm), we first fix the maximum flow time $T$ which should be sufficiently large such that the sign problem is reduced there. A possible criterion is that the sign average $|\langle e^{i\theta_{T}(x)}\rangle_{S^{\rm eff}_{T}}|$ is $O(1)$ in the absence of tempering. This process can be carried out by a test run with small statistics. We then enlarge the configuration space from ${\mathbb{R}}^{N}=\{x\}$ to the set of $A+1$ replicas, $({\mathbb{R}}^{N})^{A+1}=\{(x_{0},x_{1},\ldots,x_{A})\}$. We assign to replicas $a$ ($a=0,1,\ldots,A$) the flow times $t_{a}$ with $t_{0}=0<t_{1}<\cdots<t_{A}=T$. The action at replica $a$, $S^{\rm eff}_{t_{a}}(x_{a})$, is obtained by solving (2) and (7) with its own initial conditions $z_{t=0}^{i}=x_{a}^{i}$, $J_{t=0}=\boldsymbol{1}$. We set up an irreducible, aperiodic Markov chain for the enlarged configuration space such that the probability distribution for $\{(x_{0},x_{1},\ldots,x_{A})\}$ eventually approaches the equilibrium distribution proportional to

[TABLE]

This can be realized by combining (a) the Metropolis algorithm (or the Hybrid Monte Carlo algorithm) in the $x$ direction at each fixed flow time and (b) the swap of configurations at two adjacent replicas. Each of the steps (a) and (b) can be done in parallel processes. After the system is well relaxed to global equilibrium, we estimate the expectation value at flow time $t_{a}$ [see (4)] by using the subsample at replica $a$, $\{x_{a}^{(k)}\}_{k=1,2,\ldots,N_{\rm conf}}$, that is retrieved from the total sample $\{(x_{0}^{(k)},x_{1}^{(k)},\ldots,x_{A}^{(k)})\}_{k=1,2,\ldots,N_{\rm conf}}$:

[TABLE]

The original proposal in [14] is to use (9) at the maximum flow time, $\bar{\mathcal{O}}_{a=A}$, as an estimate of $\langle\mathcal{O}\rangle_{S}$.

Recall here that the left-hand side of (9) is independent of $a$ due to Cauchy’s theorem, and thus the ratio $\bar{\mathcal{O}}_{a}$ should yield the same value within the statistical error margin if the system is well in global equilibrium. In practice, this is not true for small $a$’s due to the sign problem, where the estimate of the sign average, \bigl{|}\overline{e^{i\theta_{t_{a}}}}\bigr{|}\equiv\bigl{|}(1/N_{\rm conf})\sum_{k}e^{i\theta_{t_{a}}(x^{(k)}_{a})}\bigr{|}, can be smaller than its statistical error ($\simeq 1/\sqrt{2N_{\rm conf}}$; the value for the uniform distribution of phases). In this case, the statistical error of the ratio $\bar{\mathcal{O}}_{a}$ cannot be trusted, which means that such $\bar{\mathcal{O}}_{a}$ should not be used as an estimate of $\langle\mathcal{O}\rangle_{S}$.

Based on the observation above, we now propose an algorithm which allows a precise estimation of $\langle\mathcal{O}\rangle_{S}$ with a criterion ensuring global equilibrium and the sufficiency of the sample size. First, we continue the sampling until we find some range of $a$ (to be denoted by $a=a_{\rm min},\ldots,a_{\rm max}(=A)$) in which \bigl{|}\overline{e^{i\theta_{t_{a}}}}\bigr{|} are well above $1/\sqrt{2N_{\rm conf}}$ and $\bar{\mathcal{O}}_{a}$ take the same value within the statistical error margin. We will require that the $1\sigma$ intervals around \bigl{|}\overline{e^{i\theta_{t_{a}}}}\bigr{|} be above $3/\sqrt{2N_{\rm conf}}$. Then, we estimate $\langle\mathcal{O}\rangle_{S}$ by using the $\chi^{2}$ fit for $\{\bar{\mathcal{O}}_{a}\}_{a=a_{\rm min},\ldots,a_{\rm max}}$ with a constant function of $a$. Global equilibrium and the sufficiency of the sample size are checked by looking at the optimized value of $\chi^{2}/{\rm DOF}=\chi^{2}/(a_{\max}-a_{\min})$. Note that the parameters determined by this procedure (such as $N_{\rm conf}$, $a_{\rm min}$, $a_{\rm max}=A$) can vary depending on the choice of observable $\mathcal{O}$.

We close this section with a few comments. First, in the TLTM a sufficient overlap of the distributions at adjacent replicas is expected even for large flow times as long as the spacings are not too large. This is because the distributions at large $a$’s ($\propto\exp[-S^{\rm eff}_{t_{a}}(x)]$) have peaks at the same points in ${\mathbb{R}}^{N}$ that flow to critical points in ${\mathbb{C}}^{N}$. This is in sharp contrast with the situation in other tempered systems, where the distribution often changes rapidly as a function of the tempering parameter so that an enough overlap cannot be achieved for realistically meaningful small spacings. Second, the optimal form of $t_{a}$ is a linear function of $a$ when flowed configurations are close to a critical point. This is because the optimal choice for the overall coefficients in tempering algorithms is exponential (see, e.g., [28, 29]) and because the real part of the action grows exponentially in flow time near critical points. Finally, the computational cost in the TLTM is expected to be $O(N^{3-4})$ due to the increase caused by the tempering algorithm (which will be $O(N^{0-1})$). Note that this growth of computational cost can be compensated by increasing the number of parallel processes.

III Application to the Hubbard model away from half-filling

Let $\Lambda=\{{\bf x}\}$ be a $d$-dimensional lattice with $N_{s}$ lattice points. The Hubbard model describes nonrelativistic lattice fermions of spin one-half, and is defined by the Hamiltonian (including the chemical potential)

[TABLE]

Here, $c_{{\textbf{x}},\sigma}$ and $c_{{\textbf{x}},\sigma}^{\dagger}$ are the annihilation and creation operators on site ${\textbf{x}}\in\Lambda$ with spin $\sigma$ $(=\uparrow,\downarrow)$ obeying the anticommutation relations $\{c_{{\textbf{x}},\sigma},c_{{\textbf{y}},\tau}^{\dagger}\}=\delta_{{\textbf{x}}{\textbf{y}}}\,\delta_{\sigma\tau}$ and $\{c_{{\textbf{x}},\sigma},c_{{\textbf{y}},\tau}\}=\{c_{{\textbf{x}},\sigma}^{\dagger},c_{{\textbf{y}},\tau}^{\dagger}\}=0$, and $n_{{\textbf{x}},\sigma}\equiv c_{{\textbf{x}},\sigma}^{\dagger}c_{{\textbf{x}},\sigma}$. $K_{{\textbf{x}}{\textbf{y}}}$ is the adjacency matrix that takes a nonvanishing value ($\equiv 1$) only for nearest neighbors, and we assume the lattice to be bipartite. $\kappa\,(>0)$ is the hopping parameter, $\mu$ is the chemical potential, and $U\,(>0)$ represents the strength of the on-site repulsive potential. We have shifted $n_{{\textbf{x}},\sigma}$ as $n_{{\textbf{x}},\sigma}\to n_{{\textbf{x}},\sigma}-1/2$ so that $\mu=0$ corresponds to the half-filling state, $\sum_{\sigma}\langle n_{{\textbf{x}},\sigma}-1/2\rangle=0$.

We approximate the grand partition function ${\rm tr}\,e^{-\beta H}$ by using the Trotter decomposition with equal spacing $\epsilon$ ($\beta=N_{\tau}\epsilon$), and rewrite it as a path integral over a Gaussian Hubbard-Stratonovich variable $\phi=(\phi_{\ell,{\textbf{x}}})$. Then the expectation value of the number density $n\equiv(1/N_{s})\sum_{{\textbf{x}}}(n_{{\textbf{x}},\uparrow}+n_{{\textbf{x}},\downarrow}-1)$ is expressed as (see Appendix B for the derivation)

[TABLE]

where $\phi_{\ell}\equiv(\phi_{\ell,{\textbf{x}}}\,\delta_{{\textbf{x}}{\textbf{y}}})$ and $\prod_{\ell}$ is a product in descending order. Note that $n[\phi]$ in (14) can be replaced by $(i\sqrt{\epsilon U}N_{\tau}N_{s})^{-1}\,\sum_{\ell,{\textbf{x}}}\phi_{\ell,{\textbf{x}}}$, which is more preferable in Monte Carlo calculations because statistical errors will be reduced due to the averaging over $\ell$. The charge-charge correlation, $\langle n_{\textbf{x}}\,n_{\textbf{y}}\rangle_{S}$ ($n_{\textbf{x}}\equiv n_{{\textbf{x}},\uparrow}+n_{{\textbf{x}},\downarrow}-1$), can also be evaluated as a path integral by simply replacing $n_{\textbf{x}}$ by $(i\sqrt{\epsilon U})^{-1}\,\phi_{\ell=0,{\textbf{x}}}$ when ${\textbf{x}}\neq{\textbf{y}}$. As for the observables that are not directly constructed from $n_{\textbf{x}}$, the expectation values can be evaluated by using the formula (42).

We now apply the TLTM to the Hubbard model on a two-dimensional periodic square lattice of size $2\times 2$ (thus $N_{s}=4$) with $N_{\tau}=5$. We first estimate $\langle n\rangle_{S}$ numerically by using the expressions (11)–(14) for various values of $\beta\mu$ with other parameters fixed to be $\beta\kappa=3$, $\beta U=13$. Note that the physical quantities depend only on the dimensionless parameters $\beta\mu$, $\beta\kappa$, $\beta U$ for fixed $N_{\tau}$.

The complex action (12) gives rise to a serious sign problem, as can be seen in the left panel of Fig. 1.

However, we should note that the extent of the seriousness of the sign problem heavily depends on the choice of the Hubbard-Stratonovich variables, and actually, the sign problem can be avoided for the above parameters within the BSS (Blankenbecler, Scalapino and Sugar)-QMC method [30]. In fact, the right panel of Fig. 1 shows the sign averages calculated by using a public code called ALF (Algorithms for Lattice Fermions) [31] that is based on the discrete variables introduced in [32, 33]. We see that the sign averages are above 0.98 for all the range of $\beta\mu$ studied here.111 We thank a referee for suggesting us to investigate this point.

Following the general prescription and writing $x=(x^{i})=(\phi_{\ell,{\textbf{x}}})$ $(i=1,\ldots,N)$ with $N=N_{\tau}N_{s}$, we introduce the enlarged configuration space $(\mathbb{R}^{N})^{A+1}=\{(x_{0},x_{1},\ldots,x_{A})\}$. We here brief the setup of the parameters relevant to the TLTM (see Appendix D for more details). We set $t_{a}$ to be piecewise linear in $a$ with a single breakpoint whose position will be tuned such that the acceptance rates of the swapping process at adjacent replicas are almost the same for all pairs (being roughly above 40%).222 This functional form of $t_{a}$ is best suited to the case where the deformed region reaches the vicinity of all the relevant Lefschetz thimbles at almost the same flow time and such a linear form is effective also for the transient period. For each value of $\beta\mu$, we make a test run with small statistics to adjust parameters. This gives the values $T/(\beta\mu)=1/12$–$1/10$, $A=8$–$12$, $N_{\rm conf}=5,000$–$25,000$, varying on the value of $\beta\mu$. We make a sampling after discarding 5,000 configurations, and from the obtained data $\{\bar{n}_{a}\}_{a=a_{\rm min},\ldots,a_{\rm max}}$ we estimate $\langle n\rangle_{S}$ by using the $\chi^{2}$ fit.

As an example, let us see Fig. 2, which shows \bigl{|}\overline{e^{i\theta_{t_{a}}}}\bigr{|} and $\bar{n}_{a}$ at various replicas for $\beta\mu=5$.

The left panel shows that the $1\sigma$ intervals around \bigl{|}\overline{e^{i\theta_{t_{a}}}}\bigr{|} are above $3/\sqrt{2N_{\rm conf}}$ for $a=5,\ldots,11$ (and thus we set $a_{\rm min}=5$ and $a_{\rm max}=11$). The right panel shows that the data $\{\bar{n}_{a}\}$ in this range give the same value within the statistical error margin. The $\chi^{2}$ fit gives the estimate $\langle n\rangle_{S}\approx 0.221\pm 0.012$ (exact value: 0.212) with $\chi^{2}/{\rm DOF}=0.45$.

Figure 3 shows the thus-obtained numerical estimates of $\langle n\rangle_{S}$ as a function of $\beta\mu$.

We also display the estimates obtained without tempering (at the same maximum flow times $T$) and those from the original reweighting method (i.e. $T=0$), together with the values obtained by the explicit evaluation of the trace under the Trotter decomposition with $N_{\tau}=5$ and for the continuum imaginary time (i.e. $N_{\tau}=\infty$) (see Appendix C). We see that the exact values are correctly reproduced when the tempering is implemented, while there are significant deviations when not implemented. As in the $(0+1)$-dimensional massive Thirring model [14], the deviation reflects the fact that the relevant thimbles are not sampled sufficiently. In fact, from Fig. 4, which shows the distribution of averaged flowed configurations $\hat{z}\equiv(1/N)\sum_{i}z_{T}^{i}$ at $T=0.5$ for $\beta\mu=5$,

we see that, although the flowed configurations are widely distributed over many thimbles when the tempering is implemented, they are restricted to only a small number of thimbles when not implemented.

Three comments are in order. First, a larger value of the sign average does not necessarily mean a better resolution of the sign problem, as can be seen from Fig. 5. In fact, when only a very few thimbles are sampled, the sign average can become larger than the value in the correct sampling due to the absence of phase mixtures among different thimbles.

Second, whether the multimodality can affect the estimates of expectation values depends on the choice of observables. In fact, from the discrepancies of the sign averages in Fig. 5, we see that the multimodality must be severe in the region $\beta\mu\leq 9$. However, the estimates of $\langle n\rangle_{S}$ almost agree between the two methods with and without tempering in the range $7\leq\beta\mu\leq 9$. This means that the operator $n$ is not sensitive to the multimodality in this range. To find an observable that is sensitive to the multimodality, we estimated the nearest-neighbor charge-charge correlation $\langle n_{\textbf{x}}\,n_{\textbf{y}}\rangle_{S}$ with the same sample.333We thank the referee for suggesting us to investigate the expectation values of observables other than the number density operator. The results are shown in Fig. 6, where we see a significant discrepancy at $\beta\mu=9$ between the two methods.

Such discrepancies become more manifest if we look at the observables that are not directly constructed from the number density operator $n_{\textbf{x}}$. As an example, we show in Fig. 7 the expectation values of the kinetic energy operator (without the factor “$-\kappa$”) $K\equiv\sum_{{\textbf{x}},{\textbf{y}}}\sum_{\sigma}K_{{\textbf{x}}{\textbf{y}}}\,c^{\dagger}_{{\textbf{x}},\sigma}c_{{\textbf{y}},\sigma}$, which are estimated for the same sample as above by using the formula (42).

We there notice two things. One is that the discrepancies between the two methods now become significant for all the range $7\leq\beta\mu\leq 9$. The other is that the precision of the TLTM becomes worse compared to the case for the observables that are constructed solely from $n_{\textbf{x}}$. In fact, those observables that are not directly constructed from $n_{\textbf{x}}$ (such as $K$) contain matrix elements of $M^{a/b}[\phi]^{-1}$, and may have divergently large values in the vicinity of zeros of the fermion determinants $\det M^{a/b}[\phi]$. In this case, precise estimation will require a larger sample size and a more accuracy in integrating flow equations compared with operators constructed solely from $n_{\textbf{x}}$. We expect that a similar attention must be paid when one applies the TLTM to finite density QCD. We leave a further investigation of this point as a future investigation.

Finally, from Fig. 8, we see that it should be a difficult task to find an intermediate flow time (without tempering) that avoids both the sign problem (severe at smaller flow times) and the multimodal problem (severe at larger flow times) (see Appendix E for more detailed discussions).

Generically, flowed configurations repeatedly experience the event at which they get trapped to fewer number of Lefschetz thimbles, so that there is a large ambiguity in distinguishing the larger and smaller flow times in the first place.

IV Conclusion and outlook

In this paper, we proposed an algorithm for the TLTM which allows a precise estimation of expectation values. We confirm the effectiveness by applying it to the two-dimensional Hubbard model away from half-filling.

We should stress that our study in this paper is still at an exploratory level. In fact, the lattice must be enlarged much more both in the spatial and imaginary time directions to claim the validity of our method for the sign problem in the Hubbard model, revealing the phase structure of the model. In doing this, it should be important to check whether the computational scaling is actually $O(N^{3-4})$ as expected.

More generally, we should keep developing the algorithm further so that it can be more easily applied to the three major problems listed in Introduction. There should also be other interesting branches of fields where the TLTM may shed new light on the theoretical understanding through a numerical analysis, such as the Chern-Simons theory [34] and matrix models that generate random volumes [35].

When we were preparing the first version of the manuscript, there appeared an interesting paper [23] (see also its detailed version [24]), where the sign and ergodicity problems are also studied for the Lefschetz thimble method applied to the Hubbard model away from half-filling. In our method (TLTM), the two problems are solved simultaneously by tempering the system with the flow time, where one does not need to know a detailed structure of thimbles. In contrast, in [23, 24] they redundantly introduce two continuous Gaussian Hubbard-Stratonovich variables with a parameter representing the mixture of the two variables (see also [22]). With knowledge of thimble structures, they tune the parameter in such a way that only a few number of thimbles become relevant to the evaluation, and obtain results for a $2\times 2$ hexagonal lattice $(N_{s}=8)$ with $N_{\tau}=384$ and $\beta=30$. It would be interesting to introduce such redundant integration variables also in the TLTM so as to reduce the global sign problem (possible cancellation of phases among different thimbles), which we observe also depends heavily on the choice of integration variables.

Acknowledgements.

The authors thank Yoshimasa Hidaka, Issaku Kanamori, Norio Kawakami, Yoshio Kikukawa, Jun Nishimura, Akira Ohnishi, Masaki Tezuka, Asato Tsuchiya and Urs Wenger for useful discussions. They also thank an anonymous referee of Physical Review D for giving us valuable comments, which were very helpful in improving the first version of the manuscript. This work was partially supported by JSPS KAKENHI (Grant Numbers 16K05321, 18J22698 and 17J08709) and by SPIRITS 2019 of Kyoto University (PI: M.F.).

Appendix A Summary of the algorithm

We summarize the algorithm of the TLTM (we partially repeat the presentation of [14]):

$\bullet$ Step 0. We fix the maximum flow time $T$ which should be sufficiently large such that the sign problem is reduced there. A possible criterion is that the sign average $|\langle e^{i\theta_{T}(x)}\rangle_{S^{\rm eff}_{T}}|$ is $O(1)$ in the absence of tempering. This can be carried out by a test run with small statistics. We then pick up flow times $\{t_{a}\}$ from the interval $[0,T]$ with $t_{0}=0<t_{1}<\cdots<t_{A}=T$. The values of $A$ and $t_{a}$ are determined manually or adaptively to optimize the acceptance rate in Step 3 below. Practically, once $A$ is determined, $t_{a}$ can be chosen to be a piecewise linear function of $a$ [see the argument for (74)].

$\bullet$ Step 1. For each replica $a$, we choose an initial value $x_{a}\in{\mathbb{R}}^{N}$ and numerically solve the differential equations (2) and (7) to obtain the triplet $(x_{a},z_{a}\equiv z_{t_{a}}(x_{a}),J_{a}\equiv J_{t_{a}}(x_{a}))$.

$\bullet$ Step 2. For each replica $a$, we use the Metropolis algorithm to update the value of $x_{a}$. To be explicit, we take a value $x^{\prime}_{a}$ from $x_{a}$ using a symmetric proposal distribution, and recalculate the triplet $(x^{\prime}_{a},z^{\prime}_{a},J^{\prime}_{a})$ using the $x^{\prime}_{a}$ as the initial value. We then update $x_{a}$ to $x^{\prime}_{a}$ with the probability ${\rm min}(1,e^{-\Delta S_{a}})$, where

[TABLE]

We repeat the process sufficiently many times such that local equilibrium is realized for each $a$. Step 1 and Step 2 can be performed in parallel processes.

$\bullet$ Step 3. We swap the configurations at two adjacent replicas $a$ and $a+1$ by updating $(x_{a},x_{a+1})=(x,y)$ to $(x^{\prime}_{a},x^{\prime}_{a+1})=(y,x)$ with the probability

[TABLE]

One can easily see that this satisfies the detailed balance condition with respect to the global equilibrium distribution (8) because

[TABLE]

We repeat the process several times so as to reduce autocorrelations. This procedure can also be performed in parallel processes by choosing a set of independent pairs.

$\bullet$ Step 4. By repeating Step 2 and Step 3, we obtain a sequence of triplets,

[TABLE]

for each $a$, with which we estimate the expectation value at flow time $t_{a}$:

[TABLE]

Here, $N_{\rm conf}$ is chosen to be large enough so that we find some range of $a$ (to be denoted by $a=a_{\rm min},\ldots,a_{\rm max}$ with $a_{\rm max}=A$) in which the $1\sigma$ intervals around \bigl{|}\overline{e^{i\theta_{t_{a}}}}\bigr{|}=\bigl{|}(1/N_{\rm conf})\sum_{k}e^{i\theta_{t_{a}}(x^{(k)}_{a})}\bigr{|} are above $3/\sqrt{2N_{\rm conf}}$ and $\bar{\mathcal{O}}_{a}$ take the same value within the statistical error margin.

$\bullet$ Step 5. The expectation value of $\langle\mathcal{O}\rangle_{S}$ is estimated by the $\chi^{2}$ fit from the data $\{\bar{\mathcal{O}}_{a}\}_{a=a_{\rm min},\ldots,a_{\rm max}}$ with a constant function of $a$. Global equilibrium and the sufficiency of the sample size $N_{\rm conf}$ is checked by looking at the optimized value of $\chi^{2}/{\rm DOF}=\chi^{2}/(a_{\max}-a_{\min})$.

In the above algorithm, we have implicitly assumed that the action at $t_{0}=0$ does not exhibit multimodality. If this is not the case, we further introduce other parameters (such as the overall coefficient of the action) as extra tempering parameters or prepare flow times $\{t_{a}\}$ with $t_{0}<0$ [14].

Appendix B Derivation of eqs. (11)–(14)

For a bipartite lattice, we specify which sublattice x belongs to by the sign $(-1)^{\textbf{x}}=\pm 1$. We first make the so-called particle-hole transformation, $c_{{\textbf{x}},\uparrow}=a_{\textbf{x}}$ and $c_{{\textbf{x}},\downarrow}=(-1)^{\textbf{x}}\,b_{\textbf{x}}^{\dagger}$. Then the one-body part $H_{1}$ and the two-body part $H_{2}$ of the Hamiltonian (10) are rewritten, respectively, as

[TABLE]

In the last equation, we have used the identity $n^{a}_{\textbf{x}}\,(\equiv a_{\textbf{x}}^{\dagger}a_{\textbf{x}})=(n^{a}_{\textbf{x}})^{2}$ and $n^{b}_{\textbf{x}}\,(\equiv b_{\textbf{x}}^{\dagger}b_{\textbf{x}})=(n^{b}_{\textbf{x}})^{2}$. Note that the number density operator is written as

[TABLE]

In order to perform a Monte Carlo simulation, we approximate $e^{-\beta H}$ in the grand partition function by using the Trotter decomposition with equal spacing $\epsilon$ ($\beta=N_{\tau}\epsilon$):

[TABLE]

and rewrite $e^{-\epsilon H_{2}}$ at the $\ell$-th position from the right to the exponential of a fermion bilinear by using a Gaussian Hubbard-Stratonovich variable $\phi_{\ell,{\textbf{x}}}$:

[TABLE]

Then, the approximated grand partition function takes the following path integral form:

[TABLE]

Here, $[d\phi]\equiv\prod_{\ell,{\textbf{x}}}d\phi_{\ell,{\textbf{x}}}$, and the action $S[\phi]$ is given by

[TABLE]

where $\prod_{\ell}$ is an ordered product ($\prod_{\ell}f_{\ell}\equiv f_{N_{\tau}-1}\cdots f_{1}f_{0}$), and ${\rm tr}_{a}$ (or ${\rm tr}_{b}$) represents the trace over the Fock space created by $a_{\textbf{x}}^{\dagger}$ (or by $b_{\textbf{x}}^{\dagger}$). The fermion trace in (26) can be evaluated explicitly by using the following formulas that hold for the operator $\hat{A}\equiv\sum_{{\textbf{x}},{\textbf{y}}}\,A_{{\textbf{x}}{\textbf{y}}}\,a_{\textbf{x}}^{\dagger}a_{\textbf{y}}$ constructed from an $N_{s}\times N_{s}$ matrix $A=(A_{{\textbf{x}}{\textbf{y}}})$:

[TABLE]

(The first equation can be readily proved by the fact that $A\mapsto\hat{A}$ is a Lie algebra homomorphism. The second equation can be easily understood by moving to a diagonalizing basis for $A$.) We thus find that the action becomes

[TABLE]

where $\phi_{\ell}$ is a diagonal matrix of the form $\phi_{\ell}=(\phi_{\ell,{\textbf{x}}}\,\delta_{{\textbf{x}}{\textbf{y}}})$. Note that, while the action is real-valued for the half-filling case $(\mu=0)$ due to the identity $M^{b}[\phi]|_{\mu=0}=(M^{a}[\phi]|_{\mu=0})^{\ast}$, it is generically complex-valued when $\mu\neq 0$.

The expectation values of such observables that are made solely from the number density operators $n_{\textbf{x}}\equiv n_{{\textbf{x}},\uparrow}+n_{{\textbf{x}},\downarrow}-1=n^{a}_{\textbf{x}}-n^{b}_{\textbf{x}}$ can be evaluated as a path integral over $\phi$ by simply replacing $n_{\textbf{x}}$ by $(i\sqrt{\epsilon U})^{-1}\,\phi_{\ell=0,{\textbf{x}}}$, as easily proved by using the operator identity

[TABLE]

For example, the expectation value of the number density operator, $\langle n\rangle_{S}$, can be rewritten to a path integral form as in (11).

As for observables of general form, one can resort to the Wick-Bloch-de Dominicis theorem,

[TABLE]

to obtain the following expression:

[TABLE]

where $\Delta^{a/b}[\phi]=M^{a/b}[\phi]^{-1}$.

Appendix C Evaluation of the trace under the Trotter decomposition

The Hilbert space $\mathbb{V}$ of the Hubbard model after the particle-hole transformation is the tensor product of two Fock spaces, $\mathbb{V}=\mathbb{V}^{a}\otimes\mathbb{V}^{b}$, each constructed by acting $a_{\textbf{x}}^{\dagger}$ or $b_{\textbf{x}}^{\dagger}$ on the Fock vacuum $|0\rangle$. In this appendix, we give the explicit forms of the matrix elements that appear in the trace under the Trotter decomposition:

[TABLE]

Here, the one-body part $H_{1}$ and the two-body part $H_{2}$ of the Hamiltonian are given by [see (20) and (21)]

[TABLE]

The number density operator is given by

[TABLE]

and we have introduced the transfer matrices corresponding to $H_{1}$ and $H_{2}$:

[TABLE]

We first introduce a one-dimensional ordering to the set of all spatial coordinates, $\Lambda=\{{\textbf{x}}\}$, and take a basis of $\mathbb{V}$ to be

[TABLE]

where the states

[TABLE]

are labeled by the subsets of ordered coordinates, $X=\{{\textbf{x}}_{1},{\textbf{x}}_{2},\ldots,{\textbf{x}}_{m}\}\subset\Lambda$ (with ${\textbf{x}}_{1}<{\textbf{x}}_{2}<\cdots<{\textbf{x}}_{m}$), $Y=\{{\textbf{y}}_{1},{\textbf{y}}_{2},\ldots,{\textbf{y}}_{n}\}\subset\Lambda$ (with ${\textbf{y}}_{1}<{\textbf{y}}_{2}<\cdots<{\textbf{y}}_{n}$). We will denote their sizes by $|X|=m$, $|Y|=n$.

The matrix elements of $T^{a}_{1}=e^{-\epsilon H^{a}_{1}}$ are then given by the following determinants:

[TABLE]

as can be easily proven by investigating the action of $T^{a}_{1}$ on the state $|X^{\prime}\rangle$:

[TABLE]

where the coefficients do not vanish only when $|X|=|X^{\prime}|\,(=m)$. The matrix elements of $T^{b}_{1}=e^{-\epsilon H^{b}_{1}}$ can also be given in the forms of determinant,

[TABLE]

We thus obtain the explicit forms of the matrix elements $(T_{1})_{XY,X^{\prime}Y^{\prime}}=(T^{a}_{1})_{XX^{\prime}}\,(T^{b}_{1})_{YY^{\prime}}$.

As for $T_{2}=e^{-\epsilon H_{2}}$, we note that $H_{2}$ acts on $|X\rangle\otimes|Y\rangle$ diagonally:

[TABLE]

The coefficients $(h_{2})_{XY}$ can be calculated easily to be

[TABLE]

where $\theta()$ is the logical step function and $\bar{X}$ stands for the complement of the set $X$, $\bar{X}=\Lambda\setminus X$. The matrix elements of $T_{2}$ is then given by $(T_{2})_{XY,X^{\prime}Y^{\prime}}=e^{-\epsilon(h_{2})_{XY}}\,\delta_{XX^{\prime}}\,\delta_{YY^{\prime}}$.

Finally, the matrix elements of $n$ are given by

[TABLE]

With the matrix elements given above, $\langle n\rangle_{S}$ can be expressed as

[TABLE]

Appendix D Summary of the parameters in the computation

We summarize the parameters relevant to the TLTM in the estimation of $\langle n\rangle_{S}$. We order the termination times $t_{a}$ for replicas $a$ as $t_{0}=0<t_{1}<\cdots<t_{A}=T$ ($T$: the largest flow time), and set $t_{a}$ to be a piecewise linear function of $a$ with a single breakpoint at $a=a_{c}$, by assuming that the deformed region reaches the vicinity of all the relevant Lefschetz thimbles at almost the same flow time and that the linear form is effective also for the transient period:

[TABLE]

Each Monte Carlo step consists of 50 Metropolis tests in the $x$ direction and $N_{\rm swap}$ swaps of configurations at adjacent replicas, and the flow equations (2) and (7) are integrated numerically by using the adaptive Runge-Kutta of 7-8th order. For each value of $\beta\mu$, we make a test run with small statistics and adjust the parameters $A$, $t_{c}$, $a_{c}$ in such a way that the acceptance rates of the swapping process at adjacent replicas are almost the same for all pairs (being roughly above 40%). After this, we make another test run of 1,000 data points to adjust the width of the Gaussian proposal in the Metropolis test in the $x$ direction so that the acceptance rate is in the range 50%–80%. This width varies depending on replicas $a$ and the values of $\beta\mu$. Using the adjusted parameters, we get a sample of size $N_{\rm conf}$ after discarding 5,000 configurations, and analyze the data by using the Jackknife method, with bins whose sizes are adjusted by taking account of autocorrelations. Finally, from the obtained data $\{\bar{n}_{a}\}$ $(a=a_{\rm min},\ldots,a_{\rm max}(=A))$ [see (9)], we estimate the expectation value $\langle n\rangle_{S}$ by using the $\chi^{2}$ fit with a constant function of $a$. We confirm that the system is in global equilibrium and the sample size is sufficient by looking at the optimized value of $\chi^{2}/{\rm DOF}$ with ${\rm DOF}=a_{\rm max}-a_{\rm min}$. The obtained results are summarized in Table 1.

Appendix E More on the fine-tuning of flow time without tempering

In order to understand the difficulty to find such an intermediate value of flow time that avoids both the sign and multimodal problems (without tempering), let us see the right panel of Fig. 8, which is the counterpart of Fig. 2 (with tempering) for the same $\beta\mu=5$. We see that the estimated values have large statistical errors at smaller flow times (due to the sign problem) while they have small statistical errors around incorrect values at larger flow times (due to the trapping of configurations at a small number of thimbles). The best flow time must be at the boundary between the two regions, but it should be a difficult task to find such value out of the set of flow times with finite spacings. In fact, if one takes a flow time from the smaller region, then, although the obtained estimate may happen to be close to the correct value, it must have a large statistical error. On the other hand, if a flow time is taken from the larger region, it will give an incorrect value (but with a small statistical error because only a small number of thimbles are sampled).

In order to understand Fig. 2 and Fig. 8 as reflecting the extent of the sign and multimodal problems, let us see Fig. 9, which depicts the normalized histograms of phases $\theta_{t_{a}}(x)$ for $\beta\mu=5$ with tempering (top) and without tempering (bottom).

We see that at smaller flow times the histograms are almost flat for the both cases (giving rise to the sign problem), but at larger flow times those without tempering become almost unimodal (reflecting the trapping at a small number of thimbles) while those with tempering correctly come to have various peaks (which may not be so obvious from the figure because there are many peaks and each peak is broadened by the Jacobian determinant).