Auxiliary Field Quantum Monte Carlo for Multiband Hubbard Models: Controlling the Sign and Phase Problems to Capture Hund's Physics

TL;DR

This paper develops an auxiliary field quantum Monte Carlo method to accurately model multiband Hubbard models with Hund's coupling, effectively controlling sign and phase problems to predict magnetic and phase behaviors in strongly correlated materials.

Contribution

It introduces a complete AFQMC framework for the Hubbard Kanamori model, controlling sign and phase issues to achieve high accuracy in modeling Hund's physics.

Findings

Achieves less than 1% energy error for moderate to large Hund's coupling

Accurately predicts magnetic ordering and phase transitions

Enables predictive modeling of Hund's physics in complex materials

Abstract

In the study of strongly-correlated, many-electron systems, the Hubbard Kanamori (HK) model has emerged as one of the prototypes for transition metal oxide physics. The model is multi-band in nature and contains Hund's coupling terms, which have pronounced effects on metal-insulator transitions, high-temperature superconductivity, and other physical properties. In the following, we present a complete theoretical framework for treating the HK model using the ground state Auxiliary Field Quantum Monte Carlo (AFQMC) method and analyze its performance on few-band models whose parameters approximate those observed in ruthenate, rhodates, and other materials exhibiting Hund's physics. Unlike previous studies, the constrained path and phaseless approximations are used to respectively control the sign and phase problems, which enables high accuracy modeling of the HK model's ground state properties within parameter regimes of experimental interest. We demonstrate that, after careful consideration of the Hubbard-Stratonovich transformations and trial wave functions employed, relative errors in the energy of less than 1% can routinely be achieved for moderate to large values of the Hund's coupling constant. Crucially, our methodology also accurately predicts magnetic ordering and phase transitions. The results presented open the door to more predictive modeling of Hund's physics within a wide range of strongly-correlated materials using AFQMC.

Figures

Tables

Table 1. Table 1 : The ground state energy of the two-band, 6 × \times 1 HK model with N ↑ = N ↓ = 6 subscript 𝑁 ↑ subscript 𝑁 ↓ 6 N_{\uparrow}=N_{\downarrow}=6 over a range of parameters using ED and AFQMC. All energies and parameters are reported in units of t 𝑡 t .

U	U1	U2	J	ED	AFQMC
2.0	1.5	1.0	0.5	-3.773268	-3.774(3)
2.0	1.5	1.0	1.0	-4.234037	-4.230(6)
2.0	1.5	3.0	0.5	0.758540	0.755(4)
3.0	5.0	1.0	0.5	2.460374	2.466(5)
6.0	1.5	1.0	0.5	1.496509	1.503(6)

Table 2. Table 2 : The charge gaps of the three-band model at Δ = 0.2 Δ 0.2 \Delta=0.2 and Δ = 1.5 Δ 1.5 \Delta=1.5 for different system sizes calculated using AFQMC. GHF trial wave functions with FM order and AFM order are used at Δ = 0.2 Δ 0.2 \Delta=0.2 and Δ = 1.5 Δ 1.5 \Delta=1.5 , respectively. All of the other Hamiltonian parameters are held fixed at t = 1 𝑡 1 t=1 , U = 6 𝑈 6 U=6 , U 1 = U − 2 ​ J subscript 𝑈 1 𝑈 2 𝐽 U_{1}=U-2J , U 2 = U − 3 ​ J subscript 𝑈 2 𝑈 3 𝐽 U_{2}=U-3J , and J = 0.15 ​ U 𝐽 0.15 𝑈 J=0.15U . The electron density per band is 4 / 3 4 3 4/3 .

# of bands	Charge Gap ($\Delta =0.2$)	Charge Gap ($\Delta =1.5$)
4x4x3	0.222(29)	1.311(32)
6x6x3	0.103(27)	1.268(35)
8x8x3	0.015(72)	1.225(36)
$\infty$	-0.006(47)	1.201(41)

Equations

$$\hat{H}\equiv {\hat{H}}_1+{\hat{H}}_2\equiv {\hat{H}}_1+{\hat{H}}_U+{\hat{H}}_J,$$

$${\hat{H}}_1=\sum _{im\sigma }\sum _{j{m}^{\prime }{\sigma }^{\prime }}{t}_{im,j{m}^{\prime }}^{\sigma {\sigma }^{\prime }}{\hat{c}}_{im\sigma }^{†}{\hat{c}}_{j{m}^{\prime }{\sigma }^{\prime }},$$

$${\hat{H}}_U$$

$$\begin{array}{rl}{\hat{H}}_J=\sum _{i,m\ne m^{\prime }}{J}_{im{m}^{\prime }}& ({\hat{c}}_{im\uparrow }^{†}{\hat{c}}_{i{m}^{\prime }\downarrow }^{†}{\hat{c}}_{im\downarrow }{\hat{c}}_{i{m}^{\prime }\uparrow }\\ & +{\hat{c}}_{im\uparrow }^{†}{\hat{c}}_{im\downarrow }^{†}{\hat{c}}_{i{m}^{\prime }\downarrow }{\hat{c}}_{i{m}^{\prime }\uparrow }+H.c.).\end{array}$$

$$\begin{array}{rl}\hat{H}& ={\hat{H}}_1+{\hat{H}}_2\\ & =\sum _{ij,\sigma {\sigma }^{\prime }}{t}_{ij}^{\sigma {\sigma }^{\prime }}{\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j{\sigma }^{\prime }}\\ & +\sum _i{U}^i{\hat{n}}_{i\uparrow }{\hat{n}}_{i\downarrow }\\ & +\sum _{i<j}{U}_1^{ij}({\hat{n}}_{i\uparrow }{\hat{n}}_{j\downarrow }+{\hat{n}}_{i\downarrow }{\hat{n}}_{j\uparrow })\\ & +\sum _{i<j}{U}_2^{ij}({\hat{n}}_{i\uparrow }{\hat{n}}_{j\uparrow }+{\hat{n}}_{i\downarrow }{\hat{n}}_{j\downarrow })\\ & +\sum _{i<j}{J}^{ij}({\hat{c}}_{i\uparrow }^{†}{\hat{c}}_{j\downarrow }^{†}{\hat{c}}_{i\downarrow }{\hat{c}}_{j\uparrow }+{\hat{c}}_{i\uparrow }^{†}{\hat{c}}_{i\downarrow }^{†}{\hat{c}}_{j\downarrow }{\hat{c}}_{j\uparrow }\\ & +{\hat{c}}_{j\uparrow }^{†}{\hat{c}}_{i\downarrow }^{†}{\hat{c}}_{j\downarrow }{\hat{c}}_{i\uparrow }+{\hat{c}}_{j\uparrow }^{†}{\hat{c}}_{j\downarrow }^{†}{\hat{c}}_{i\downarrow }{\hat{c}}_{i\uparrow }).\end{array}$$

$$|{\Psi }_0⟩\propto \underset{\beta \to \infty }{\lim }\left(e^{-\beta \hat{H}}\right)|{\Psi }_I⟩.$$

$$|{\Psi }_0⟩\propto \underset{n\to \infty }{\lim }{\left(e^{-\Delta \tau \hat{H}}\right)}^n|{\Psi }_I⟩,$$

$$|{\Psi }^{(n+1)}⟩=e^{-\Delta \tau \hat{H}}|{\Psi }^{(n)}⟩.$$

$$e^{-\Delta \tau \hat{H}}\approx e^{-\Delta \tau {\hat{H}}_1/2}{e}^{-\Delta \tau {\hat{H}}_2}{e}^{-\Delta \tau {\hat{H}}_1/2}.$$

$$\begin{array}{rl}{e}^{-\Delta \tau {\hat{H}}_2}& \approx e^{-\Delta \tau {\hat{H}}_U}{e}^{-\Delta \tau {\hat{H}}_{U_1}}{e}^{-\Delta \tau {\hat{H}}_{U_2}}{e}^{-\Delta \tau {\hat{H}}_J}\\ & =e^{-\Delta \tau \sum _i{U}_i{\hat{n}}_{i\uparrow }{\hat{n}}_{i\downarrow }}{e}^{-\Delta \tau \sum _{i<j}{U}_1^{ij}({\hat{n}}_{i\uparrow }{\hat{n}}_{j\downarrow }+{\hat{n}}_{i\downarrow }{\hat{n}}_{j\uparrow })}\\ & e^{-\Delta \tau \sum _{i<j}{U}_2^{ij}({\hat{n}}_{i\uparrow }{\hat{n}}_{j\uparrow }+{\hat{n}}_{i\downarrow }{\hat{n}}_{j\downarrow })}\\ & e^{-\Delta \tau \sum _{i<j}{J}^{ij}({\hat{c}}_{i\uparrow }^{†}{\hat{c}}_{j\downarrow }^{†}{\hat{c}}_{i\downarrow }{\hat{c}}_{j\uparrow }+{\hat{c}}_{i\uparrow }^{†}{\hat{c}}_{i\downarrow }^{†}{\hat{c}}_{j\downarrow }{\hat{c}}_{j\uparrow }+H.c.)}.\end{array}$$

$$e^{-\Delta \tau \alpha {\hat{n}}_1{\hat{n}}_2}=e^{-\Delta \tau \alpha ({\hat{n}}_1+{\hat{n}}_2-1)/2}\sum _{x=\pm 1}\frac{1}{2}{e}^{\gamma x({\hat{n}}_1+{\hat{n}}_2-1)},$$

$$e^{-\Delta \tau \alpha {\hat{n}}_1{\hat{n}}_2}=e^{-\Delta \tau \alpha ({\hat{n}}_1+{\hat{n}}_2)/2}\sum _{x=\pm 1}\frac{1}{2}{e}^{\gamma x({\hat{n}}_1-{\hat{n}}_2)},$$

$${\hat{\rho }}_{ij}\equiv \sum _{\sigma }({\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j\sigma }+{\hat{c}}_{j\sigma }^{†}{\hat{c}}_{i\sigma }).$$

$${\hat{\rho }}_{ij}^2=\sum _{\sigma {\sigma }^{\prime }}({\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j\sigma }+{\hat{c}}_{j\sigma }^{†}{\hat{c}}_{i\sigma })({\hat{c}}_{i{\sigma }^{\prime }}^{†}{\hat{c}}_{j{\sigma }^{\prime }}+{\hat{c}}_{j{\sigma }^{\prime }}^{†}{\hat{c}}_{i{\sigma }^{\prime }}),$$

$$\begin{array}{rl}{\hat{H}}_J& =\sum _{i<j}{J}^{ij}({\hat{c}}_{i\uparrow }^{†}{\hat{c}}_{j\downarrow }^{†}{\hat{c}}_{i\downarrow }{\hat{c}}_{j\uparrow }+{\hat{c}}_{i\uparrow }^{†}{\hat{c}}_{i\downarrow }^{†}{\hat{c}}_{j\downarrow }{\hat{c}}_{j\uparrow }+H.c.)\\ & =\sum _{i<j}\frac{J^{ij}}{2}[{\hat{\rho }}_{ij}^2-\sum _{\sigma }({\hat{n}}_{i\sigma }+{\hat{n}}_{j\sigma }-{\hat{n}}_{i\sigma }{\hat{n}}_{j\sigma }-{\hat{n}}_{j\sigma }{\hat{n}}_{i\sigma })]\\ & =\sum _{i<j}\frac{J^{ij}}{2}{\hat{\rho }}_{ij}^2-\sum _{i<j,\sigma }\frac{J^{ij}}{2}({\hat{n}}_{i\sigma }+{\hat{n}}_{j\sigma })+\sum _{i<j,\sigma }{J}^{ij}{\hat{n}}_{i\sigma }{\hat{n}}_{j\sigma }.\end{array}$$

$$e^{-\Delta \tau {\hat{A}}^2/2}=\int dx\frac{1}{\sqrt{2\pi }}{e}^{-x^2/2}{e}^{x\sqrt{-\Delta \tau }\hat{A}},$$

$$\begin{array}{rl}& e^{-\Delta \tau \sum _{i<j}\frac{J^{ij}}{2}[\sum _{\sigma }({\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j\sigma }+{\hat{c}}_{j\sigma }^{†}{\hat{c}}_{i\sigma }){]}^2}\\ & =\prod _{i<j}\int d{x}_{ij}\frac{1}{\sqrt{2\pi }}{e}^{-x_{ij}^2/2}{e}^{x_{ij}\sqrt{-\Delta \tau J^{ij}}[\sum _{\sigma }({\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j\sigma }+{\hat{c}}_{j\sigma }^{†}{\hat{c}}_{i\sigma })]}.\end{array}$$

$${\hat{\rho }}_{ij}=\sum _{\sigma }{\delta }_{\sigma }({\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j\sigma }+{\hat{c}}_{j\sigma }^{†}{\hat{c}}_{i\sigma }),$$

$${\hat{\rho }}_{ij}^2=\sum _{\sigma {\sigma }^{\prime }}{\delta }_{\sigma }{\delta }_{{\sigma }^{\prime }}({\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j\sigma }+{\hat{c}}_{j\sigma }^{†}{\hat{c}}_{i\sigma })({\hat{c}}_{i{\sigma }^{\prime }}^{†}{\hat{c}}_{j{\sigma }^{\prime }}+{\hat{c}}_{j{\sigma }^{\prime }}^{†}{\hat{c}}_{i{\sigma }^{\prime }}).$$

$$\begin{array}{rl}{\hat{H}}_J& =\sum _{i<j}{J}^{ij}({\hat{c}}_{i\uparrow }^{†}{\hat{c}}_{j\downarrow }^{†}{\hat{c}}_{i\downarrow }{\hat{c}}_{j\uparrow }+{\hat{c}}_{i\uparrow }^{†}{\hat{c}}_{i\downarrow }^{†}{\hat{c}}_{j\downarrow }{\hat{c}}_{j\uparrow }+H.c.)\\ & =\sum _{i<j}-\frac{J^{ij}}{2}[{\hat{\rho }}_{ij}^2-\sum _{\sigma }({\hat{n}}_{i\sigma }+{\hat{n}}_{j\sigma }-{\hat{n}}_{i\sigma }{\hat{n}}_{j\sigma }-{\hat{n}}_{j\sigma }{\hat{n}}_{i\sigma })]\\ & =\sum _{i<j}-\frac{J^{ij}}{2}{\hat{\rho }}_{ij}^2+\sum _{i<j,\sigma }\frac{J^{ij}}{2}({\hat{n}}_{i\sigma }+{\hat{n}}_{j\sigma })-\sum _{i<j,\sigma }{J}^{ij}{\hat{n}}_{i\sigma }{\hat{n}}_{j\sigma }.\end{array}$$

$$\begin{array}{rl}& e^{\Delta \tau \sum _{i<j}\frac{J^{ij}}{2}[\sum _{\sigma }{\delta }_{\sigma }({\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j\sigma }+{\hat{c}}_{j\sigma }^{†}{\hat{c}}_{i\sigma }){]}^2}\\ & =\prod _{i<j}\int d{x}_{ij}\frac{1}{\sqrt{2\pi }}{e}^{-x_{ij}^2/2}{e}^{x_{ij}\sqrt{\Delta \tau J^{ij}}[\sum _{\sigma }{\delta }_{\sigma }({\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j\sigma }+{\hat{c}}_{j\sigma }^{†}{\hat{c}}_{i\sigma })]},\end{array}$$

$$e^{-\Delta \tau \hat{H}}=\int d\mathbf{x}p(\mathbf{x})\hat{B}(\mathbf{x}),$$

$$|{\Psi }^{(n)}⟩=\sum _k{w}_k^{(n)}|{\psi }_k^{(n)}⟩.$$

$$⟨\hat{A}{⟩}_{mix}$$

$$E_L[{\Psi }_T,\Phi ]\equiv \frac{⟨{\Psi }_T|\hat{H}|\Phi ⟩}{⟨{\Psi }_T|\Phi ⟩},$$

$$⟨\hat{A}{⟩}_{mix}=\frac{{\sum }_k{w}_k^{(n)}⟨{\Psi }_T|{\psi }_k^{(n)}⟩E_L[{\Psi }_T,{\psi }_k^{(n)}]}{{\sum }_k{w}_k^{(n)}⟨{\Psi }_T|{\psi }_k^{(n)}⟩}.$$

$$W_k^{(n)}=w_k^{(n)}⟨{\Psi }_T|{\psi }_k^{(n)}⟩.$$

$$\hat{V}=-\frac{1}{2}\sum _i{\left({\hat{v}}_i-⟨{\hat{v}}_i⟩\right)}^2-\sum _i{\hat{v}}_i⟨{\hat{v}}_i⟩+\frac{1}{2}\sum _i⟨{\hat{v}}_i{⟩}^2,$$

$$\sum _{i<j}\frac{J^{ij}}{2}{\hat{\rho }}_{ij}^2$$

$$\sum _{i<j}-\frac{J^{ij}}{2}{\hat{\rho }}_{ij}^2$$

Full text

I Introduction

Since the unanticipated discovery of high-temperature superconductivity in the cuprates, the single-band Hubbard model Hubbard (1963) has been the focus of an unparalleled level of theoretical scrutiny and associated algorithmic development.LeBlanc et al. (2015); Lieb and Wu (1968); Scalapino (2007); Gukelberger et al. (2015) Nevertheless, most materials exhibiting strong correlation, including most transition metal oxidesKuneš et al. (2008); Laad et al. (2006); Maeno et al. (1999) as well as the pnictides,Si et al. (2016); Yin et al. (2011); Haule and Kotliar (2009) fullerides,Nomura et al. (2016, 2015) and chalcogenidesSun et al. (2012); Si et al. (2016); Yin et al. (2011) possess multiple bands that cross their Fermi levels and are therefore fundamentally multi-band in nature.Georges et al. (2013) In recent years, it has become increasingly evident that some of the most significant effects in such multi-band materials stem from Hund’s coupling.Yin et al. (2011); Haule and Kotliar (2009); Johannes and Mazin (2009) According to Hund’s rules, electrons favor maximizing their total spin by first occupying different, degenerate bands in the same shell with parallel spins; only after they fill all available bands do they then doubly occupy the same bands.Georges et al. (2013) As such, the effective Coulomb repulsion among electrons in a half-filled shell is increased due to Hund’s rules, while that at any other filling is decreased. Hund’s effects therefore drive half-filled $d$- and $f$-electron materials closer to a Mott transition for a given Coulomb repulsion, yet drive non-half-filled materials away from a Mott transition while also increasing the correlation within their metallic phases. The consequences of these effects are perhaps best illustrated in 4$d$ transition metal oxides that have more than a single electron or hole in their 4$d$ shells. Koster et al. (2012); Fradkin et al. (2010); Dang et al. (2015a); Han et al. (2016) Unlike their rhodate counterparts, which possess a single hole in their shells, many ruthenates and molybdenates exhibit substantial mass enhancements,Ikeda et al. (2000) unexpected Mott Insulator transitions,Liebsch and Ishida (2007); Gorelov et al. (2010); Sutter et al. (2017) novel quantum phase transitions,Grigera et al. (2001) and even superconducting phasesMackenzie and Maeno (2003) – all of which may be attributed to Hund’s physics.

Despite both the prevalence and importance of Hund’s effects, they remain a challenge to describe. Most analytical and numerical treatments revolve around solving a multi-band Hubbard model, most often the Hubbard-Kanamori (HK) model,Kanamori (1963) containing a mixture of kinetic, Coulomb $U$, and Hund’s $J$ terms. Although analytical studies have been performed,Roth (1966); Lyon-Caen and Cyrot (1975); Khomskii and Kugel (1973); K.I. Kugel (1982); Ishihara et al. (1997) just as in the case of the single-band Hubbard model containing a repulsive $U$ term,LeBlanc et al. (2015); Zheng et al. (2017) accurate treatments of these models necessitate methods capable of treating strong correlation non-perturbatively. However, because these models possess significantly larger state spaces and involve additional pair-hopping and Hund’s exchange terms, they are often even more difficult to treat than the Hubbard model.

Due to the complicated interactions involved, there is no general analytical solution for these problems. Thus, numerical treatments are in high demand. To date, most numerical studies of multi-band models have employed Dynamical Mean Field Theory (DMFT)Georges et al. (1996); Metzner and Vollhardt (1989) either on its own or in combination with Density Functional Theory (DFT)Anisimov et al. (1997) because of DMFT’s ability to treat band and atomic effects on equal footing by self-consistently solving an impurity problem within a larger bath. DMFT has been very successful at mapping out multi-band phase diagrams at finite temperatures.Gorelov et al. (2010); Han and Millis (2018); Han et al. (2016); Dang et al. (2015a, b); Werner et al. (2008); Inaba et al. (2005) Nevertheless, DMFT is fundamentally limited by the accuracy and scaling of its impurity model solver. Some DMFT studies rely upon exact diagonalization (ED) to solve their impurity models, yet the computational cost of ED grows exponentially with the number of bands involved, thus thwarting its application to many-band models. Some DMFT algorithms employ continuous-time quantum Monte Carlo (CTQMC)Rubtsov et al. (2005) to solve their impurity models. CTQMC can solve larger impurity models than ED, but is still hampered by the sign problem, an exponential decrease in the signal to noise ratio observed in stochastic simulations,Loh et al. (1990) in certain parameter regimes and low temperature calculations remain difficult.Gull et al. (2011) A method that can accurately simulate larger system sizes at lower temperatures is thus in need.

One suite of techniques particularly well-suited for studying the large state spaces inherent to multi-band models are quantum Monte Carlo (QMC) techniques.Foulkes et al. (2001); Motta and Zhang (2018) Both finite temperature QMC methods, including CTQMCGull et al. (2011) and Hirsch-Fye QMCHirsch and Fye (1986) algorithms that have been employed as impurity solvers within DMFT, and ground stateMotome and Imada (1997, 1998) QMC algorithms have been developed and applied to the multi-band Hubbard model. Nonetheless, the Hund’s terms of the HK Hamiltonian have posed challenges for all of these methods. This is because Hund’s terms are not readily expressed as products of density operators and are therefore not readily amenable to standard QMC transformations. Straightforward decoupling of the exchange and pair hopping terms leads to a severe sign problem.Held, K. and Vollhardt, D. (1998) Attempts have therefore been made to simplify the Hund’s contribution to the Hamiltonian to make it more palatable to QMC methods by constraining its direction to the z-axis,Held, K. and Vollhardt, D. (1998); Han et al. (1998) but such treatments sometimes fail to properly capture the model’s expected physics. Several Hund’s-specific transformations have been proposed, including a discrete transformation by Aoki Sakai et al. (2004, 2006) and a continuous transformation by Imada. Motome and Imada (1997, 1998) Nevertheless, these transformations ultimately do not eliminate the sign problem and are limited to parameter regimes with only high signal to noise ratios. These parameter constraints obscure our fundamental understanding of multi-band physics.

In this paper, we present an Auxiliary Field Quantum Monte Carlo (AFQMC) framework especially suited for the study of ground state multi-band Hubbard models and demonstrate its accuracy over a range of realistic parameters using different signal-preserving approximations and trial wave functions. Key to our approach is the strategic use of two forms of both the continuous and discrete Hubbard-Stratonovich (HS) Transformations to decouple the Hund’s term: a charge decomposition for negative values of the Hund’s coupling parameter, and a spin decomposition for positive values of the Hund’s coupling parameter. We also employ an unconventional form of importance sampling in which we shift propagators instead of auxiliary fields so as to enable importance sampling of discrete transformed propagators. Unlike previous works, we furthermore utilize flexible Generalized Hartree-Fock (GHF) trial wave functions combined with the constrained path and phaseless approximations to tame the sign and phase problems, respectively. Altogether, we find that these improvements yield promising results for a variety of HK model benchmarks. Although the algorithm presented is designed for the ground state, it can easily be adapted for use in finite temperature methods.Liu et al. (2018); Zhang (1999) Our algorithm therefore paves the way to the high accuracy modeling of the low temperature physics of a wide range of multi-band models and materials over a dramatically larger portion of the phase diagram.

The remainder of the paper is organized as follows. In Section II, we outline the HK model, summarize the key features of the AFQMC method, and describe how the conventional AFQMC technique may be modified to best accommodate the HK Hamiltonian. In Section III, we then present benchmarks of our method’s performance within different parameter regimes, using different trial wave functions, and employing different approximations on two- and three-band HK models for which ED results may be obtained. Towards the end of this section, we also demonstrate the accuracy with which our techniques can predict the charge gaps and magnetic ordering of two-dimensional lattice models far beyond the reach of most other techniques. We conclude with a discussion of the broader implications of this work and future directions in Section IV.

II Methods

II.1 Hubbard Kanamori Model Hamiltonian

The HK model is a multi-band version of the Hubbard model designed to account for the competition between the spin and band degrees of freedom observed in the physics of $d$- and $f$- electron material.Kanamori (1963); Georges et al. (2013) In order to accomplish this, the model includes not only standard Hubbard on-site density-density interactions, but also inter-band density, exchange, and pair hopping terms. The full HK Hamiltonian, written as general as possible, reads

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

In the above, $\hat{c}_{im\sigma}^{\dagger}$($\hat{c}_{im\sigma}$) creates (annihilates) an electron with spin $\sigma$ in band $m$ at site $i$. $\hat{n}$ denotes the number operator and $\hat{n}_{im\uparrow}$, for example, represents the number of spin-up electrons at site $i$ in band $m$. $\hat{H}_{1}$ contains all one-body contributions to the Hamiltonian, including terms parameterized by the constants $t_{im,jm^{\prime}}^{\sigma\sigma^{\prime}}$ that describe spin-orbit coupling and the hopping of electrons in different bands between sites $i$ and $j$. $\hat{H}_{2}$ denotes the collection of all two-body operators. $\hat{H}_{U}$ contains all density-density interactions, including the intraband ($U$) and interband ($U^{\prime}$) Coulomb interactions, and the $z$- (or Ising) component of the Hund’s coupling. In contrast, $\hat{H}_{J}$ contains all of the terms that cannot be written as density-density interactions, which consist of the $x$- and $y$- components (spin-exchange) of the Hund’s coupling, ($\hat{c}_{im\uparrow}^{\dagger}\hat{c}_{im^{\prime}\downarrow}^{\dagger}\hat{c}_{im\downarrow}\hat{c}_{im^{\prime}\uparrow}+H.c.$), as well as the pair-hopping interaction ($\hat{c}_{im\uparrow}^{\dagger}\hat{c}_{im\downarrow}^{\dagger}\hat{c}_{im^{\prime}\downarrow}\hat{c}_{im^{\prime}\uparrow}+H.c.$), in which two electrons in a given band transfer as a pair to other bands. $J$ denotes the Hund’s coupling constant. Note that our formalism is general and allows for band- and site-dependent $U$, $U$’, and $J$ constants.

II.2 Modified Hubbard Kanamori Model Hamiltonian

In order to facilitate programming and the generalization of this HK Hamiltonian into a form in which all coupling constants are independent, we map the Hamiltonian given by Equations (1)-(4) into a one-band model whose terms only depend upon their band indices. If we now let $i$ and $j$ denote superindices that combine both lattice site and band information, then

[TABLE]

$t_{ij}^{\sigma\sigma^{\prime}}$ describes the hopping and spin-orbit coupling between different sites and bands. In keeping with the $\sum_{i,m<m^{\prime}}$ and $\sum_{i,m\neq m^{\prime}}$ summations in Equations (3) and (4), $\sum_{i<j}$ only sums over index combinations that reference different bands on the same site. In this modified HK Model, the $U$ term describes density-density interactions only between electrons with opposite spins in the same band, the $U_{1}$ term describes interactions between electrons with opposite spins in different bands on the same site, the $U_{2}$ term describes interactions between electrons with parallel spins in different bands on the same site, and the $J$ term describes spin-exchange and pair hopping interactions on the same site. Thus, in going from Equations (3) and (4) to Equation (5), the original $U^{\prime}$ term has become the $U_{1}$ term, the original $(U^{\prime}-J)$ term has become the $U_{2}$ term, and the $J$ term has been re-expressed. Using Equation (5), we map a multi-band model into a single-band model in which the number of lattice sites has been enlarged into the number of bands. Since there is no explicit index in our model, we can deal with any number of bands as long as the mapping is done correctly.

II.3 Overview of AFQMC

In the remainder of this work, AFQMC will be employed to obtain accurate numerical solutions to the HK Model. AFQMC is a quantum many-body method that solves the ground state Schrodinger Equation by randomly sampling an overcomplete space of non-orthogonal Slater determinantsZhang et al. (1997); Zhang (2003, 2013) and has consistently been demonstrated to be among the most accurate of modern many-body methods for modeling the Hubbard model over a wide range of parameter regimes. LeBlanc et al. (2015); Zheng et al. (2017); Chang and Zhang (2010, 2008); Shi and Zhang (2013) At its heart, AFQMC is an imaginary-time projection quantum Monte Carlo technique that applies a projection operator, $e^{-\beta\hat{H}}$, onto an initial wave function, $|\Psi_{I}\rangle$,

[TABLE]

In the limit of infinite imaginary projection time $(\beta\rightarrow\infty)$, it converges to the ground state wave function, $|\Psi_{0}\rangle$, as long as the initial wave function is not orthogonal to the ground state wave function. Because the projection operator cannot be evaluated for large values of $\beta$, it is discretized into $n=\beta/\Delta\tau$ smaller time slices for which it can be evaluated

[TABLE]

and the projection is carried out iteratively as follows

[TABLE]

For sufficiently small $\Delta\tau$, the projection operator may be factored into one- and two-body pieces via Suzuki-Trotter FactorizationSuzuki (1976); Trotter (1959)

[TABLE]

The two-body propagator may be further decomposed into the four terms given in Equation (5)

[TABLE]

A time step extrapolation is needed to make sure the Trotter error is negligible in the Monte Carlo simulation.

II.4 Hubbard-Stratonovich Transformation of the Modified Hubbard Kanamori Hamiltonian

According to Thouless’s Theorem, Thouless (1960) acting the exponential of a one-body operator on a determinant results in another determinant, reducing the process of projecting a one-body operator onto the wave function into standard matrix multiplication. Nevertheless, no such theorem applies to exponentials of two-body operators, which necessitates re-expressing these operators into integrals over one-body operators using the so-called Hubbard-Stratonovich transformation. Hubbard (1959)

In order to transform the two-body propagator given by Equation (10), both discreteHirsch (1983); Gubernatis et al. (2016) and continuousBuendia (1986) HS transformations need to be performed. The $U$, $U_{1}$, and $U_{2}$ terms are products of density operators, much like the conventional Hubbard $U$ term, and may therefore be decomposed using discrete transformations. For $\alpha<0$, where $\alpha$ may be denote $U$, $U^{1}$, or $U^{2}$, it is usually better to use the discrete charge decomposition

[TABLE]

where $\cosh(\gamma)=e^{-\Delta\tau\alpha/2}$, while for $\alpha>0$, it is usually better to use the spin decomposition

[TABLE]

where $\cosh(\gamma)=e^{\Delta\tau\alpha/2}$. In both Equations (11) and (12), $x$ represents the namesake auxiliary field that may assume the discrete values of $+1$ or $-1$. For the subsequent discussion, note that the charge decomposition is so named because it produces a one-body propagator involving the sum of $\hat{n}_{1}+\hat{n}_{2}$, which would be equivalent to the charge on a site if $1$ represented an up and $2$ a down spin on that site. Along similar lines, the spin decomposition is so named because it involves the difference between $\hat{n}_{1}$ and $\hat{n}_{2}$, which would represent the spin on a site under the same assumptions.

Because $\hat{H}_{J}$ contains terms that are not simple products of density operators, decomposing it is a much more challenging task. Past attempts have either neglected or simplified $\hat{H}_{J}$.Held, K. and Vollhardt, D. (1998); Han et al. (1998) Several techniques have employed exact decompositions,Motome and Imada (1997); Sakai et al. (2006, 2004) but all such decompositions are accompanied by a sign problem that thwarts explorations of wide swaths of the phase diagram. Unlike these past attempts, in the following, we define a unique decomposition that can be employed in both continuous and discrete transformations, and accompany it by importance sampling that first mitigates and the constrained path and phaseless approximations that eliminate the sign and phase problems. As part of our decomposition of $e^{-\Delta\tau\hat{H}_{J}}$, we first re-expressed $\hat{H}_{J}$ in terms of squares of one-body operators. Let

[TABLE]

Then,

[TABLE]

and $\hat{H}_{J}$ may be re-expressed as (see the Supplemental Materials for more details)

[TABLE]

The second term of Equation (15) consists of one-body operators and can be combined with the other one-body operators into $\hat{H}_{1}$. The third term consists of a product of density operators and can therefore be transformed according to either Equations (11) or (12). The first term, however, consists of a square that cannot be resolved into products of density operators. In order to decouple this two-body term, a continuous HS transformation must be employed. In general, the continuous HS transformation may be written as

[TABLE]

where $\hat{A}$ represents any one-body operator and $x$ denotes an auxiliary field, as before. Letting $\hat{A}\equiv\hat{\rho}_{ij}$, it follows that the most obvious way to transform the exponential formed from the first term of Equation (15) is using the charge decomposition

[TABLE]

As long as $J^{ij}<0$ for all $i,j$, all of the propagators produced by this transformation will be real, as is desirable within AFQMC simulations. However, if any of the $J^{ij}$ are greater than 0, $\sqrt{-\Delta\tau J^{ij}}$ will be complex resulting in a complex propagator that immediately introduces a complex phase into simulations. To prevent complexity from being introduced into the operators, in certain cases, we take a cue from the discrete case and define a continuous spin decomposition that involves the difference between spin up and down operators. Let

[TABLE]

where $\delta_{\uparrow}=1$ and $\delta_{\downarrow}=-1$, then (see the Supplemental Materials for further details)

[TABLE]

Using this to re-express $\hat{H}_{J}$, we have

[TABLE]

Employing this form for the decomposition, the exponential that stems from the first term of Equation (20) may now be transformed to yield

[TABLE]

which is real for $J^{ij}>0$. Using the charge decomposition (Equation (17)) when $J^{ij}<0$ and the spin decomposition (Equation (21)) when $J^{ij}>0$ thus completely eliminates complex propagators, easing simulation. In Section III.1, we compare the merits of using this mixed decomposition approach to exclusively relying upon the complex charge decomposition on the accuracy of our overall results.

Inserting the HS transformations defined by Equations (11), (12), (17), and (21) into Equations (9) and (10) and combining terms, one arrives at the final AFQMC expression for the projection operator

[TABLE]

where ${\bf x}=\{x_{1},x_{2},...,x_{N_{F}}\}$ denotes the set of $N_{F}$ total normally distributed auxiliary fields sampled at a given time slice, $\hat{B}(x)$ represents the amalgamation of all one-body operators, and $p(x)$ is a combination of all scalar functions of the fields. Example expressions for $\hat{B}(x)$ and $p(x)$ are given in the Supplemental Information. As is clear from Equation (22), the series of HS Transformations described ultimately maps the original two-body propagator into a weighted integral over one-body propagators that are functions of external auxiliary fields.

II.5 Sampling in AFQMC

II.5.1 The Sampling Process

One of the most computationally efficient ways of evaluating many dimensional integrals such as that given by Equation (22) is to use Monte Carlo sampling techniques. As described in more detail in previous publications,Hirsch (1985); Zhang et al. (1997); Zhang (2003, 2013) if $|\Psi_{I}\rangle$ is represented by a single Slater determinant, after each application of the projection operator, a new Slater determinant will be produced. Thus, if $k$ instances (so-called “walkers”) are initialized to $|\Psi_{I}\rangle$ and the projection operation given by Equation (22) is applied to each of them by independently sampling sets of fields, then a random walk through the space of non-orthogonal determinants is realized in which the overall wave function at time slice $n$, $|\Psi^{(n)}\rangle$, is represented by an ensemble of $k$ wave functions $|\psi_{k}^{(n)}\rangle$ with weights $w_{k}^{(n)}$

[TABLE]

Here, the $w_{k}^{(n)}$ consist of the products of numbers accumulated over all time slices by walker $k$, which can be a complex number.

Ground state observables at each time slice, such as the energies reported below, may then be computed by evaluating the mixed estimatorFoulkes et al. (2001) over the ensemble

[TABLE]

where $|\Psi_{T}\rangle$ denotes a trial wave function that approximates the true ground state wave function. To facilitate the evaluation of the mixed estimator, it is common to introduce the local energy

[TABLE]

such that Equation (24) may be simplified to

[TABLE]

After a sufficiently large number of time slices such that $|\Psi^{(n)}\rangle$ approaches the ground state, final estimates of $\langle\hat{A}\rangle$ may be obtained by averaging over each time slice expectation value.

A population control procedure Calandra Buonaura and Sorella (1998) is needed during the random walk. During this procedure, walkers with larger weights are replicated and those with smaller weights are eliminated probabilistically. The weight used in population control is

[TABLE]

When there is a sign or phase problem, $W^{(n)}_{k}$ may become negative or complex. As described in Section (II.5.5) and Section (II.5.6), $W^{(n)}_{k}$ is always positive or zero if the constrained path or phaseless approximations are employed.

II.5.2 The Sign and Phase Problems

Unfortunately, the “free” projection process just described is typically beset by either the signLoh et al. (1990); Ceperley (1991) or phase problems.Zhang and Krakauer (2003) These problems fundamentally stem from the fact that observables computed using a single Slater determinant, $|\Psi\rangle$, remain invariant to arbitrary rotations, $e^{i\theta}|\Psi\rangle$, of that determinant, where $\theta$ is a phase angle. Consequently, during the course of an AFQMC simulation involving complex propagators, walkers may accumulate infinitely many possible phases (as there are infinitely many possible phase angles, $\theta\in[0,2\pi)$), resulting in infinitely many possible determinants. Since these phases are directly multiplied into the walker weights of Equations (24) and (26), after many iterations, the walker weights end up populating the entire complex plane and many of the terms summed to compute weighted averages of observables cancel one another out. This cancellation leads to an exponential decline in observable signal to noise ratios that manifests as infinite variancesShi and Zhang (2013) called the phase problem. If transformations that preclude propagators from becoming complex are employed as described above, positive and negative versions of each determinant may still be generated, resulting in a somewhat less pernicious cancellation of positive and negative weights termed the sign problem. If left unchecked, the sign and phase problems render obtaining meaningful observable averages nearly impossible, thwarting AFQMC simulations. We therefore mitigate these problems using a combination of background subtraction, importance sampling, and either the constrained path (for the sign problem) or phaseless (for the phase problem) approximations.

II.5.3 Background Subtraction

One of the simplest ways of reducing variances within AFQMC is via background subtraction.Purwanto and Zhang (2005) As part of background subtraction, the two-body portion of a Hamiltonian is rewritten so that a mean field average is subtracted from each one-body operator. Thus, if the original two-body operator may be written as a square such that $\hat{V}=-\frac{1}{2}\sum_{i}\hat{v}_{i}^{2}$ to make it amenable to a HS Transformation, as part of background subtraction, it would be re-expressed as

[TABLE]

where $\langle\hat{v}_{i}\rangle$ denotes the mean field average of the operator $\hat{v}_{i}$ (see the Supplemental Materials for more details on how this mean field average is obtained). Because the modified $\hat{v}_{i}-\langle\hat{v}_{i}\rangle$ operator will be smaller in magnitude than the bare $\hat{v}_{i}$ operator, background subtraction reduces the variance involved in AFQMC simulations. In this work, we perform background subtraction on the only term in the Hamiltonian that is not a product of on-site densities, the $\frac{J^{ij}}{2}\hat{\rho}_{ij}^{2}$ term of Equation (15) or the $-\frac{J^{ij}}{2}\hat{\rho}_{ij}^{2}$ term of Equation (20), yielding

[TABLE]

and

[TABLE]

respectively.

II.5.4 Importance Sampling

In order to further reduce the variance of walker weights and to make our simulations more amenable to the constrained path and phaseless approximations, we additionally perform importance sampling, which aims to shift the center of the distribution from which we sample our auxiliary fields so that the most important fields are sampled more frequently. The conventional way of performing importance sampling in AFQMC simulations is by introducing a force bias that shifts each sampled field by an amount dependent upon the operator being transformed and the current walker wave function.Purwanto and Zhang (2004); Zhang and Krakauer (2003); Rom et al. (1997); Shi et al. (2015) Because we utilize a mixture of discrete and continuous transformations and force bias importance sampling is only applicable to continuous transformations, in this work, we employ a formally equivalent strategy in which we shift the propagators instead of the auxiliary fields.

For continuous HS Transformations, this may be accomplished by shifting the operator $\hat{A}$ by $\langle\hat{A}\rangle$ in Equation (16)

[TABLE]

where $\langle\hat{A}\rangle$ is the mixed estimator of $\hat{A}$

[TABLE]

If we define the dynamic force as $F\equiv\sqrt{-\Delta\tau}\langle\hat{A}\rangle$, then Equation (31) may be re-expressed as

[TABLE]

In order to realize this transformation, fields are sampled from the shifted Gaussian probability density function, $\frac{1}{\sqrt{2\pi}}e^{-(x-F)^{2}/2}$, and the propagator $e^{x\sqrt{-\Delta\tau}\hat{A}}$ is applied with weight $e^{\frac{1}{2}F^{2}-xF}$. The field distributions are now centered around the dynamic force, which can be shown to minimize the variance. If the dynamic force $F$ is a complex number, our auxiliary fields will have the same imaginary part to ensure the $x-F$ is real. Then the probability function $\frac{1}{\sqrt{2\pi}}e^{-(x-F)^{2}/2}$ stays in the real axis, which can be sampled by Monte Carlo.

Shifting the propagator within a discrete transformation proceeds in exactly the same fashion. Comparing Equations (16) and (12), the dynamic force needed to shift the propagator in Equation (12), for example, would be $F\equiv\gamma(\langle\hat{n}_{1}\rangle-\langle\hat{n}_{2}\rangle)$, resulting in the transformation

[TABLE]

As in the continuous case, in order to realize this transformation, fields are now sampled from a shifted probability density function, $e^{xF}/W$, where $W$ is the normalization factor, $W=e^{xF}+e^{-xF}$, and the propagator $e^{(-\Delta\tau\alpha/2+\gamma x)\hat{n}_{1}}e^{(-\Delta\tau\alpha/2-\gamma x)\hat{n}_{2}}$ is applied with weight $\frac{1}{2}We^{-xF}$. A shifted transformation may similarly be constructed for the discrete charge decomposition given by Equation (11). Propagators that include background subtraction may be shifted by simply replacing $\hat{A}$ with $\hat{A}-\langle\hat{A}\rangle$ in Equations (31) and (33) above (see the Supplemental Materials).

It can readily be proven that shifting auxiliary fields is equivalent to shifting propagators.Rom et al. (1997, 1998); Shi et al. (2015) Shifting propagators therefore entails a convenient way of introducing importance sampling when discrete transformations are involved. Overall, the importance sampled propagation produces the same observable averages as free propagation, but favors the sampling of determinants with larger overlaps with the trial wave function and suppressing the sampling of determinants with no overlap.

II.5.5 Constrained Path Approximation

In order to address the sign problem that may emerge when our propagators, $\hat{B}(\vec{x})$, are real, we employ the constrained path approximation. Zhang et al. (1997) Here, we impose this approximation by requiring that all walkers maintain a positive overlap with the trial wave function after each propagation step

[TABLE]

As in typical constrained path implementations, walkers with negative overlaps with the trial wave function will be killed (have their weights set equal to zero), preventing them from being propagated further. This condition will select for only walkers with positive determinants, eliminating the sign problem. It can be shown that if the trial wave function is the exact ground state wave function, this condition will be exact;Carlson et al. (1999) however, since the trial wave function is typically unknown, constraining the propagation path in this way results in a small, but consequential approximation.Chang et al. (2016); Shi and Zhang (2013)

II.5.6 Phaseless Approximation

In cases in which our propagators are complex, instead of employing the constrained path approximation, we employ the more general phaseless approximation.Zhang and Krakauer (2003); Purwanto and Zhang (2005) The phaseless approximation controls the phase problem by projecting complex walker weights onto the positive real axis according to the equation

[TABLE]

where $W^{(n)}_{k}$ is defined in Equation (27) and $\Delta\theta$, the phase angle, is defined as

[TABLE]

The use of the cosine function to project also ensures that the density of the walkers will vanish at the origin. Because this cosine projection does not affect walkers with real weights, in practical implementations, we apply Equation 36 to realize both the constrained path and phaseless Approximations.

II.6 Trial and Initial Wave Functions

Although AFQMC can readily accommodate multi-determinant trial wave functions, we restrict ourselves to employing single determinant trial wave functions that satisfy certain symmetriesShi and Zhang (2013) such as the free electron (FE), restricted Hartree-Fock (RHF), unrestricted Hartree-Fock (UHF), and generalized Hartree-Fock (GHF) wave functions. RHF wave functions preserve spin symmetry. While RHF and UHF wave functions separately conserve the number spin up and down electrons, GHF only fix the total number of electrons. Details about how these wave functions are generated may be found in the Supplemental Materials.

As illustrated in what follows, because GHF wave functions do not impose any spin symmetries and are therefore the most flexible of these wave function ansatzes, they enable the fastest AFQMC wave function relaxation to the global energy minimum. Nevertheless, when the number of up and down electrons must be fixed, UHF/RHF wave functions were employed instead. Even though our formalism permits our initial wave functions to differ from our trial wave functions, we take our initial and trial wave functions to be the same, except where otherwise noted.

III Results and Discussion

III.1 Two-Band Hubbard Kanamori Model Benchmarks

In order to test the accuracy of our theoretical framework, we began by benchmarking our method against ED results for the one-dimensional, two-band HK Model on $5\times 1$ and $6\times 1$ lattices with periodic boundary conditions small enough to diagonalize. For these benchmarks, we simplify the Hamiltonian given by Equations (1) and (2) so that hopping can only occur between adjacent sites within the same bands and may be described by a single site- and spin-invariant constant $t$, such that

[TABLE]

We moreover assume that the parameters are site-invariant, such that Ui=U, U${}_{1}^{ij}$ = U1, U${}_{2}^{ij}$ = U2, and Jij = J.

Table 1 presents our results for a $6\times 1$ HK model over a representative sampling of parameters at half filling. All of the calculations presented were initialized using 560 walkers and employed FE trial and initial wave functions, except for the $U$=3.0,$U_{1}$=5.0,$U_{2}$=1.0,$J$=0.5 case. In this case, it was found that an RHF trial wave function yielded a lower trial energy and manifested a different spin order (antiferromagnetic (AFM) order between two bands) than the FE solution. Thus, an RHF trial wave function was employed instead. This demonstrates that trial wave functions should first be analyzed to determine whether their global minima exhibit the correct order before using them to guide propagation within AFQMC. Unless otherwise noted, all of the results presented in this section were obtained using a charge decomposition for $J$ and the phaseless approximation to tame the related phase problem that emerges.

As is clear from the table, AFQMC results, are within $0.01t$ or less of exact results, with the smallest discrepancy occurring for the $U=2.0$, $U_{1}=1.5$ case and the largest occurring for the $U=6.0$ case. In all of these cases, exact results are within two standard derivations of the Monte Carlo results, despite the use of the phaseless approximation.

To pinpoint AFQMC systematic bias, as well as to better understand which regions of the phase diagram are the most challenging for AFQMC, we independently scanned through each of the $U$, $U_{1}$, $U_{2}$, and $J$ parameters holding the others fixed for a $5\times 1$ HK model. In Figures 1 and 2, we present our scans over $U$ and $J$; figures of our $U_{1}$ and $U_{2}$ scans are presented in the Supplemental Materials.

As shown in Figure 1, although the magnitude of the error bars grows with $U$, the relative error remains within 0.1% to 1% throughout this range. Similar trends are observed for $U_{1}$ and $U_{2}$. This gives us reason to believe that our method can readily accommodate some of the even larger $U$ values used in studies of strongly correlated materials. Nevertheless, much larger relative errors are observed as $J$ is varied, as depicted in Figure 2. This is consistent with previous work, which also implicates the $J$ terms as being most conducive to QMC errors.Held, K. and Vollhardt, D. (1998) Fortunately, for most real materials, $J$ is usually a small fraction of $U$. For small $J$ values, the relative errors are observed to remain less than 1% and are therefore controllable.

What may also be gleaned from Figure 2 is that the quality of the $J>1.5$ energies depends upon the type of trial wave function employed. While free propagation calculations yield results that are independent of the trial wave function, the quality of the constrained path and phaseless approximations fundamentally depend on the accuracy of the trial wave function employed. As depicted in Figure 2, the relative errors in the energies produced by FE trial wave functions surpass 10% and increase with increasing $J$; in contrast, the relative errors produced by RHF trial wave functions not only remain less than 10%, but plateau as a function of $J$. As $J$ increases, the RHF electron density becomes non-uniform, yielding a lower variational energy than the FE wave function. Figure 2 thus demonstrates that AFQMC becomes more accurate as trial wave functions better describe the ground state. Note that we also tested UHF and GHF wave functions, which all converged to the same states as RHF wave functions.

The accuracy of AFQMC predictions are also influenced by the constrained path and phaseless approximations employed. In Figure 3, we compare the errors produced by these approximations. As discussed in Section II.4, for $J>0$, the spin decomposition will yield real propagators that we constrain using the constrained path approximation, while for $J<0$, the spin decomposition will yield complex propagators that we constrain using the phaseless approximation. The charge decomposition behaves in the opposite fashion with respect to $J$. As shown in Figure 3, the constrained path approximation behaves significantly better than the phaseless approximation, which appreciably differs from the exact results for $|J|>1.5$. Indeed, the constrained path approximation nearly reproduces the exact results for $J<0$, only manifesting a slight deviation for larger positive values of $J$. These results attest to the fact that using the transformations we describe to prevent the phase problem from emerging is key to maintaining AFQMC accuracy. They also underscore that our method is capable of simulating -$J$ values, which have been unattainable in previous QMC simulations. We expect these trends in accuracy to generalize to models with more bands and higher dimensionality.

III.2 Application to Three-Band Hubbard Kanamori Models

In order to understand how our techniques generalize to models that approximate more realistic materials and their magnetic phase transitions, we constructed a three-band model with an adjustable band gap. As illustrated in Figure 4, in this model, three bands are located at each site, one band of which is lower in energy by a ‘band gap’ parameter, $\Delta$, than the other two degenerate bands. When $\Delta=0$, all three bands are completely degenerate. Similar to the two-band model, the hopping occurs between adjacent sites within the same bands, with hopping constant $t_{ij}=1$. While the band gap would be fixed in any given material, creating a separate $\Delta$ parameter enables us to sample a range of band gaps and, by extension, to drive magnetic ordering transitions. We moreover assume that $U^{i}=U$ and $J^{ij}=J=0.15U$ with $U_{1}^{ij}=U_{1}=U-2J$ and $U_{2}^{ij}=U_{2}=U-3J$, which are appropriate for the description of transition-metal oxides with a partially occupied $t_{2g}$ shell.Sugano et al. (1970) In the following discussion, we fix our filling such that an average of four electrons occupy the three bands at each lattice site.

As an initial step, we benchmarked our AFQMC method against ED results. Diverging from our previous two-band analysis, as part of our three-band benchmarks, we studied our model on two-dimensional lattices with periodic boundary conditions, only varying $\Delta$ and $U$ while keeping the other parameter relationships fixed in order to preserve realism. Our simulations were initialized with 560 walkers and GHF initial and trial wave functions for all of the benchmarks described below. The charge decomposition with the phaseless approximation was employed throughout this section.

In Figure 5, we illustrate how the energy and relative errors change as $\Delta$ is varied from 0 to 1 with $U=6$ on a 2$\times$2 lattice. At fixed $U$, the relative error remains fairly stable and less than 0.1% throughout this range. This may be anticipated since the band gap only modifies the magnitude of the one-body terms and does not change the phase of the model, which do not directly contribute to our method’s stochastic errors.

In Figure 6, instead of scanning $\Delta$, we scan $U$ with $\Delta=0.8$. As shown in Figure 6, the relative errors are larger in this case, but still range from 0.1% for $U<6$ to 1% for $U>6$. Errors would be expected to grow in this manner as the system becomes more correlated. Overall, the magnitudes of these relative errors suggest that AFQMC’s performance is promising.

The rationale for introducing the band gap $\Delta$ parameter is to enable tuning of the magnetic order of the model system. Intuitively, when the band gap is small, the three bands are nearly degenerate and the four electrons have the largest freedom to move among the bands. Such a situation would favor ferromagnetic (FM) order. However, when the band gap becomes sufficiently large, two electrons will populate the lower band, forcing the other two electrons to reside on the higher energy bands. Such a situation would favor AFM order.

This intuition was confirmed by comparing the AFQMC energies attained using trial wave functions with FM and AFM order, respectively (see Figure 7). Typically, GHF calculations converge to the lowest state with the same magnetic order as the initial state. Thus, in order to construct wave functions with FM order, a randomly initialized density matrix was supplied to the GHF self-consistent equations; to construct wave functions with AFM order, an AFM-ordered initial density matrix was supplied. Several independent GHF calculations were conducted for each system studied to guarantee that the final GHF wave functions produced attained their global minima. At large $\Delta$ ($\Delta\gtrsim 1.1$) values at which ferromagnetic order is disfavored, GHF calculations initialized with random density matrices often developed order. In these situations, FM wave functions produced at smaller values of $\Delta$ were used as trial wave functions in “FM” AFQMC calculations performed at larger $\Delta$ values. Figure 7 depicts the energies of AFQMC simulations performed with AFM and FM trial wave functions, respectively, as a function of band gap. All of the AFQMC energies presented here are the lowest energies we can obtain at each $\Delta$. At smaller $\Delta$s, trial wave functions with FM order led to the lowest AFQMC energies, while at larger $\Delta$s, AFM trial wave functions did so. This confirms that our model undergoes a ferromagnetic to antiferromagnetic transition at roughly $\Delta=1.1$. In contrast, Hartree-Fock theory predicts the transition at $\Delta=0.5$, which is reasonable since Hartree-Fock theory tends to fall into AFM order sector. An illustration of the AFM order exhibited by our model is depicted in Figure 4.

To further corroborate the phase transition we observe, we extrapolated the magnitude of the charge gap at $\Delta=0.2$ and $\Delta=1.5$. To do so, we computed the ground state AFQMC energies of 4$\times$4, 6$\times$6, and 8$\times$8-site systems, with three bands filled with four electrons situated at each site. The charge gap may be determined by computing $E_{N-1}+E_{N+1}-2E_{N}$, where $N$ denotes the total number of electrons in the system. To determine the charge gap in the thermodynamic limit, we fit a $1/L$ form, where $L$ denotes the total number of lattice sites, to the energies and extrapolated to the infinite $L$ limit (see the Supplemental Materials for more details). The energies produced using FM initial trial wave functions were used to ascertain the $\Delta=0.2$ charge gap, while those produced using AFM wave functions were employed to ascertain the $\Delta=1.5$ charge gap. The charge gaps obtained are presented in Table 2. After extrapolations, the $\Delta=0.2$ charge gap converged to -0.006(47) and the $\Delta=1.5$ charge gap converged to 1.201(41). As one would expect antiferromagnetic, not ferromagnetic, order to be accompanied by a charge gap, these extrapolations support our previous conclusions.

The successful determination of the magnetic order and charge gaps in this model system illustrate our method’s promise for accurately modeling realistic materials.

IV Conclusions

In summary, we have presented a ground state AFQMC framework suited for the study of the HK model, a multi-band model designed to capture the Hund’s physics of many $d$- and $f$-electron materials. Diverging with past QMC studies of the HK model, we employ a novel set of HS transformations to decouple the Hund’s coupling term while preserving the term’s essential physics. We find that by carefully combining these transformations with a form of importance sampling that shifts our propagators, well-optimized GHF wave functions, and the constrained path and phaseless approximations, we can accurately predict the energetics of benchmark lattice models and the magnetic order of much larger models that approximate realistic materials. Overall, we find that the phaseless version of our method produces nearly exact energies for small models for $-3<J<3$, a range of $J$ values which contains those commonly observed in experiment. This bodes well for the generalization of our method to other systems.

Our method may readily be extended to include spin-orbit coupling effects and negative $J$ values, which opens the doors to the highly accurate study of exotic, -$J$ fulleride physics.Nomura et al. (2016, 2015) In order to describe superconducting physics, our method can be adapted to use superconducting trial wave function forms, including Bardeen-Cooper-SchriefferShi et al. (2015); Carlson et al. (2011) and Hartree-Fock-BogoliubovShi and Zhang (2017) wave functions. We foresee our method having the most immediate impact as a way to delineate low-temperature phase diagrams currently beyond the reach of DMFT methods. As the same transformations and importance sampling techniques may readily be adapted into finite temperature AFQMC formalisms,Hirsch and Fye (1986); Liu et al. (2018); Zhang (1999) the same methods may be used to develop lower scaling, sign and phase problem free impurity solvers. We look forward to employing our methods to more accurately elucidate the complex many body physics of 4$d$ transition metal oxides such as the ruthenates, rhodates, and molybedenates in the near future.

Acknowledgements.

H.H., B.R., and H.S. thank Andrew Millis, Antoine Georges, and Shiwei Zhang for stimulating discussions, and Qiang Han for providing data and insight. H.H. and B.R. acknowledge support from NSF grant DMR-1726213 and DOE grant DE-SC0019441. H.S. thanks the Flatiron Institute for research support. The Flatiron Institute is a division of the Simons Foundation. This work was conducted using computational resources and services at the Brown University Center for Computation and Visualization, the Flatiron Institute, and the Extreme Science and Engineering Discovery Environment (XSEDE).