Determinantal (or Auxiliary Field) Monte Carlo

In this section, we will briefly describe an application of this general algorithm to the quantum mechanical many-body problem of interacting electrons on a lattice, working in the grand-canonical ensemble. The basic idea of this approach was presented some time ago by Blankenbecler, Scalapino and Sugar. [13]

Quantum Monte Carlo relies on the fact that $d$-dimensional quantum problems can be intepreted as classical problems in $(d+1)$ dimensions through Feymann's Path Integral representation. The task consists in identifying Ising-like fields that would allow us to evaluate the partition function and mean values using the same Metropolis algorithm we have used before in the equivalent classical model. Suzuki [14] was the first to apply this concepts after generalizing an idea by Trotter. [15]

Suppose we want to evaluate the expectation value of a physical observable ${% \hat{O}}$, at some finite temperature $T=1/\beta$. If ${\hat{H}}$ is the Hamiltonian of the model, this expectation value is defined as,

$$\langle {\hat{O}}\rangle _{T}={\frac{{Tr({\hat{O}}e^{-\beta {\hat{H}}})}}{{% Tr(e^{-\beta {\hat{H}}})}}},$$ (307)

where the notation is the standard. From now on, let us concentrate on the particular case of the one band Hubbard model which was defined previously. The Hamiltonian of this model, with the addition of a chemical potential, can be naturally separated into two terms as,

$${\hat{K}}=-\mathrm{t}\sum_{\langle \mathbf{ij}\rangle ,\sigm... ...thbf{i}}(n_{\mathbf{i}\uparrow }+n_{\mathbf{% i}\downarrow }),$$ (308)

$${\hat{V}}=\mathrm{U}\sum_{\mathbf{i}}(n_{\mathbf{i}\uparrow ... ...c{{1}}{{2}% }})(n_{\mathbf{i}\downarrow }-{\frac{{1}}{{2}}}).$$ (309)

Discretizing the inverse temperature interval as $\beta =\Delta \tau L$, where $\Delta \tau$ is a small number, and $L$ is the total number of time slices, we can apply the well-known Trotter's formula to rewrite the partition function as,

$$Z=Tr(e^{-\Delta \tau L{\hat{H}}})\sim Tr(e^{-\Delta \tau {\hat{V}}% }e^{-\Delta \tau {\hat{K}}})^{L},$$ (310)

where a systematic error of order $(\Delta \tau )^{2}$ has been introduced, since $[{\hat{K}},{\hat{V}}]\neq 0$. In order to integrate out the fermionic fields the interaction term ${\hat{V}}$ has to be made quadratic in the fermionic creation and annihilation operators by introducing a decoupling Hubbard-Stratonovich transformation. At this stage, we can select from a wide variety of possibilities to carry out this decoupling i.e. we can choose continuous or discrete, real or complex fields, belonging to different groups. In particular, and for illustration purposes, here we use a simple transformation using a discrete ``spin-like'' field [16],

$$e^{-\Delta \tau \mathrm{U}(n_{\mathbf{i}\uparrow }-{\frac{{1... ...\lambda (n_{\mathbf{i}\uparrow }-n_{\mathbf{i}\downarrow })},$$ (311)

which is carried out at each lattice site $\mathbf{i}$, and for each temperature (or imaginary-time) slice $l$. The constant $\lambda$ is defined through the relation $\cosh (\Delta \tau \lambda )=\exp (\Delta \tau \mathrm{U/2})$. The transformation Eq.(314) reduces the four-fermion self-interaction of the Hubbard model to a quadratic term in the fermions coupled to the new spin-like field $s_{\mathbf{i},l}$. Thus, in this formalism the interactions between electrons are mediated by the spin field. Now we can carry out the integration of the fermions. While this is conceptually straightforward, and for a finite lattice of $\mathrm{N}\times \mathrm{N}$ sites it gives determinants of well-defined matrices, arriving to the actual form of these matrices is somewhat involved, and beyond the scope of this review. Then, here we will simply present the result of the integration (more details can be found in Refs. [17] and [18]. The partition function can be exactly written as,

$$Z=\sum_{\{s_{\mathbf{i},l}=\pm 1\}}\det M^{+}(s)\det M^{-}(s),$$ (312)

where

$$M^{\sigma }=I+B_{L}^{\sigma }B_{L-1}^{\sigma }...B_{1}^{\sigma },$$ (313)

and

$$B_{l}^{\pm }=e^{\mp \Delta \tau \nu (l)}e^{-\Delta \tau {\hat{K}}}.$$ (314)

$I$ is the unit matrix, and $\nu (l)_{\mathbf{ij}}=\delta _{% \mathbf{ij}}s_{\mathbf{i},l}$. Usually the physical observable ${\hat{O}}$, can be expressed in terms of Green's functions for the electrons moving in the spin field. Then, expressions similar to Eq.(315) can be derived for the numerator in Eq.(310). Once the partition function is written only in terms of the spin fields, we can use standard Monte Carlo techniques (such as Metropolis or heat bath methods) to perform a simulation of the complicated sums over $s_{\mathbf{i},l}$ that remain to be done. The probability distribution of a given spin configuration is given in principle by ${\frac{{1}}{{Z}}}\det M^{+}\det M^{-}$ (unless it becomes negative, see next section).