diff --git a/g.tex b/g.tex index 8fa7f34..58f9207 100644 --- a/g.tex +++ b/g.tex @@ -673,6 +673,8 @@ It is then a trivial matter to integrate out exactly the $n_k$ variables, leadin \\ + \sum_{p=1}^{\infty} \sum_{I_1 \notin \cD_0, \hdots , I_p \notin \cD_{p-1}} \qty[ \prod_{k=0}^{p-1} \mel{ I_k }{ P_k \qty( H-E \Id )^{-1} P_k (-H)(1-P_k) }{ I_{k+1} } ] + \\ + \times \mel{ I_p }{ P_p \qty( H-E \Id)^{-1} P_p }{ I_N } \end{multline} As an illustration, Appendix \ref{app:A} reports the exact derivation of this formula in the case of a two-state system. @@ -688,8 +690,9 @@ In fact, there is a more direct way to derive the same equation by resorting to where $H_0$ is some arbitrary reference Hamiltonian, we have the Dyson equation \be \label{eq:GE} - G^E_{ij}= \titou{G^E_{0,ij}} + \sum_{kl} G^{E}_{0,ik} (H_0-H)_{kl} G^E_{lj}. + G^E_{ij} = G^E_{0,ij} + \sum_{kl} G^{E}_{0,ik} (H_0-H)_{kl} G^E_{lj}, \ee +with $G^E_{0,ij} = \mel{i}{ \qty( H_0-E \Id )^{-1} }{j}$. Let us choose $H_0$ such that $\mel{ i }{ H_0 }{ j } = \mel{ i }{ P_i H P_i }{ j }$ for all $i$ and $j$. Then, The Dyson equation \eqref{eq:GE} becomes \begin{multline} @@ -794,7 +797,7 @@ Thus, from a practical point of view, a trade-off has to be found between the \t Let us consider the one-dimensional Hubbard Hamiltonian for a chain of $N$ sites \be - \hat{H}= -t \sum_{\expval{ i j } \sigma} \hat{a}^+_{i\sigma} \hat{a}_{j\sigma} + H= -t \sum_{\expval{ i j } \sigma} \hat{a}^+_{i\sigma} \hat{a}_{j\sigma} + U \sum_i \hat{n}_{i\uparrow} \hat{n}_{i\downarrow}, \ee where $\langle i j\rangle$ denotes the summation over two neighboring sites, $\hat{a}_{i\sigma}$ ($\hat{a}_{i\sigma}$) is the fermionic creation (annihilation) operator of an spin-$\sigma$ electron (with $\sigma$ = $\uparrow$ or $\downarrow$) on site $i$, $\hat{n}_{i\sigma} = \hat{a}^+_{i\sigma} \hat{a}_{i\sigma}$ the number operator, $t$ the hopping amplitude, and $U$ the on-site Coulomb repulsion. @@ -822,17 +825,19 @@ The parameters $\alpha$ and $\beta$ are optimized by minimizing the variational As discussed above, the efficiency of the method depends on the choice of states forming each domain. As a general guiding principle, it is advantageous to build domains associated with a large average trapping time in order to integrate out the most important part of the Green's matrix. + Here, as a first illustration of the method, we shall consider the large-$U$ regime of the Hubbard model where the construction of such domains is rather natural. -At large $U$ and half-filling, the Hubbard model approaches the Heisenberg limit where only the $2^N$ states with no double occupancy, $n_D(n)=0$, have a significant weight in the wave function. +Indeed, at large $U$ and half-filling, the Hubbard model approaches the Heisenberg limit where only the $2^N$ states with no double occupancy, $n_D(n)=0$, have a significant weight in the wave function. The contribution of the other states vanishes as $U$ increases with a rate increasing sharply with $n_D(n)$. -In addition, for a given number of double occupations, configurations with large values of $n_A(n)$ are favored due to their high kinetic energy (electrons move more easily). +In addition, for a given number of double occupations, configurations with large values of $n_A(n)$ are favored due to their high kinetic energy. Therefore, we build domains associated with small $n_D$ and large $n_A$ in a hierarchical way as described below. -For simplicity and decreasing the number of diagonalizations to perform, we shall consider only one non-trivial domain called here the main domain and denoted as $\cD$. -This domain will be chosen common to all states belonging to it, that is + +For simplicity and reducing the number of diagonalizations to perform, we shall consider only one non-trivial domain called here the main domain and denoted as $\cD$. +This domain will be chosen common to all states belonging to it, that is, \be \cD_i= \cD \qq{for} \ket{i} \in \cD. \ee -For the other states, we choose a single-state domain +For the other states, we choose a single-state domain as \be \cD_i= \qty{ \ket{i} } \qq{for} \ket{i} \notin \cD. \ee @@ -840,15 +845,15 @@ To define $\cD$, let us introduce the following set of states \be \titou{\cS_{ij} = \qty{ \ket{n} : n_D(n)=i \land n_A(n)=j }}. \ee -$\cD$ is defined as the set of states having up to $n_D^\text{max}$ double occupations and, for each state with a number of double occupations equal to $m$, a number of nearest-neighbor antiparallel pairs between $n_A^\text{min}(m)$ and $n_A^\text{max}(m)$. +which means that $\cD$ contains the set of states having up to $n_D^\text{max}$ double occupations and, for each state with a number of double occupations equal to $m$, a number of nearest-neighbor antiparallel pairs between $n_A^\text{min}(m)$ and $n_A^\text{max}(m)$. Here, $n_A^\text{max}(m)$ will not be varied and taken to be the maximum possible for a given $m$, $n_A^\text{max}(m)= \max(N-1-2m,0)$. Using these definitions, the main domain is taken as the union of some elementary domains \be - \cD = \cup_{n_D=0}^{n_D^\text{max}}\cD(n_D,n_A^{\rm min}(n_D)) + \cD = \bigcup_{n_D=0}^{n_D^\text{max}}\cD(n_D,n_A^{\rm min}(n_D)) \ee where the elementary domain is defined as \be - \cD(n_D,n_A^\text{min}(n_D))=\cup_{ n_A^\text{min}(n_D) \leq j \leq n_A^{\rm max}(n_D)} \cS_{n_D j} + \cD(n_D,n_A^\text{min}(n_D))=\bigcup_{ n_A^\text{min}(n_D) \leq j \leq n_A^{\rm max}(n_D)} \cS_{n_D j} \ee The two quantities defining the main domain are thus $n_D^\text{max}$ and $n_A^\text{min}(m)$. To give an illustrative example, let us consider the 4-site case. There are 6 possible elementary domains @@ -903,87 +908,90 @@ and Let us begin with a small chain of 4 sites with $U=12$. From now on, we shall take $t=1$. -The size of the linear space is ${\binom{4}{2}}^2= 36$ and the ground-state energy obtained by exact diagonalization is $E_0=-0.768068...$. The two variational parameters of the trial vector have been optimized and fixed at the values of $\alpha=1.292$, and $\beta=0.552$ with a variational energy of $E_\text{T}=-0.495361...$. In what follows -$|I_0\rangle$ will be systematically chosen as one of the two N\'eel states, {\it e.g.} $|I_0\rangle =|\uparrow,\downarrow, \uparrow,...\rangle$. +The size of the linear space is ${\binom{4}{2}}^2= 36$ and the ground-state energy obtained by exact diagonalization is $E_0=-0.768068...$. +The two variational parameters of the trial vector have been optimized and fixed at the values of $\alpha=1.292$, and $\beta=0.552$ with a variational energy of $E_\text{T}=-0.495361...$. +In what follows $|I_0\rangle$ will be systematically chosen as one of the two N\'eel states, {\it e.g.} $|I_0\rangle =|\uparrow,\downarrow, \uparrow,...\rangle$. Figure \ref{fig1} shows the convergence of $H_p$ as a function of $p$ for different values of the reference energy $E$. We consider the simplest case where a single-state domain is associated to each state. -Five different values of $E$ have been chosen, namely $E=-1.6,-1.2,-1,-0.9$, and -$E=-0.8$. Only $H_0$ is computed analytically ($p_{ex}=0$). At the scale of the figure error bars are too small to be perceptible. +Five different values of $E$ have been chosen, namely $E=-1.6,-1.2,-1,-0.9$, and $E=-0.8$. +Only $H_0$ is computed analytically ($p_{ex}=0$). At the scale of the figure error bars are too small to be perceptible. When $E$ is far from the exact value of $E_0=-0.768...$ the convergence is very rapid and only a few terms of the $p$-expansion are necessary. -In constrast, when $E$ approaches the exact energy, -a slower convergence is observed, as expected from the divergence of the matrix elements of the -Green's matrix at $E=E_0$ where the expansion does not converge at all. Note the oscillations of the curves as a function of $p$ due to a -parity effect specific to this system. In practice, it is not too much of a problem since -a smoothly convergent behavior is nevertheless observed for the even- and odd-parity curves. -The ratio, $\cE_{QMC}(E,p_{ex}=1,p_{max})$ as a function of $E$ is presented in figure \ref{fig2}. Here, $p_{max}$ is taken sufficiently large -so that the convergence at large $p$ is reached. The values of $E$ are $-0.780$, $-0.790$, $-0,785$, $-0,780$, and $-0.775$. For smaller $E$ the curve is extrapolated using -the two-component expression. The estimate of the energy obtained from $\cE(E)=E$ is $-0.76807(5)$ in full agreement with the exact value of $-0.768068...$. +In contrast, when $E$ approaches the exact energy, a slower convergence is observed, as expected from the divergence of the matrix elements of the Green's matrix at $E=E_0$ where the expansion does not converge at all. +Note the oscillations of the curves as a function of $p$ due to a parity effect specific to this system. +In practice, it is not too much of a problem since a smoothly convergent behavior is nevertheless observed for the even- and odd-parity curves. +The ratio, $\cE_{QMC}(E,p_{ex}=1,p_{max})$ as a function of $E$ is presented in figure \ref{fig2}. Here, $p_{max}$ is taken sufficiently large so that the convergence at large $p$ is reached. +The values of $E$ are $-0.780$, $-0.790$, $-0,785$, $-0,780$, and $-0.775$. For smaller $E$ the curve is extrapolated using the two-component expression. The estimate of the energy obtained from $\cE(E)=E$ is $-0.76807(5)$ in full agreement with the exact value of $-0.768068...$. +%%% FIG 1 %%% \begin{figure}[h!] \includegraphics[width=\columnwidth]{fig1} -\caption{1D-Hubbard model, $N=4$, $U=12$. $H_p$ as a function of $p$ for $E=-1.6,-1.2,-1.,-0.9,-0.8$. $H_0$ is -computed analytically and $H_p$ (p > 0) by Monte Carlo. Error bars are smaller than the symbol size.} +\caption{One-dimensional Hubbard model for $N=4$ and $U=12$. +$H_p$ as a function of $p$ for $E=-1.6$, $-1.2$, $-1$, $-0.9$, and $-0.8$. +$H_0$ is computed analytically and $H_p$ ($p > 0$) by Monte Carlo. +Error bars are smaller than the symbol size.} \label{fig1} \end{figure} - +%%% FIG 2 %%% \begin{figure}[h!] \includegraphics[width=\columnwidth]{fig2} -\caption{1D-Hubbard model, $N=4$ and $U=12$. $\mathcal{E}(E)$ as a function of $E$. -The horizontal and vertical lines are at $\mathcal{E}(E_0)=E_0$ and $E=E_0$, respectively. -$E_0$ is the exact energy of -0.768068.... The dotted line is the two-component extrapolation. +\caption{One-dimensional Hubbard model for $N=4$ and $U=12$. +$\cE(E)$ as a function of $E$. +The horizontal and vertical lines are at $\cE(E_0)=E_0$ and $E=E_0$, respectively. +$E_0 = -0.768068\ldots$ is the exact energy. +The dotted line is the \titou{two-component} extrapolation. Error bars are smaller than the symbol size.} \label{fig2} \end{figure} -Table \ref{tab1} illustrates the dependence of the Monte Carlo results upon the choice of the domain. The reference energy is $E=-1$. +Table \ref{tab1} illustrates the dependence of the Monte Carlo results upon the choice of the domain. +The reference energy is $E=-1$. The first column indicates the various domains consisting of the union of some elementary domains as explained above. -The first line of the table gives the results when using a minimal single-state domain for all states, and the last -one for the maximal domain containing the full linear space. The size of the various domains is given in column 2, the average trapping time -for the state $|I_0\rangle$ in the third column, and an estimate of the speed of convergence of the $p$-expansion for the energy in the fourth column. -To quantify the rate of convergence, we report the quantity, $p_{conv}$, defined as -the smallest value of $p$ for which the energy is converged with six decimal places. The smaller $p_{conv}$, the better the convergence is. +The first line of the table gives the results when using a minimal single-state domain for all states, and the last one for the maximal domain containing the full linear space. +The size of the various domains is given in column 2, the average trapping time for the state $|I_0\rangle$ in the third column, and an estimate of the speed of convergence of the $p$-expansion for the energy in the fourth column. +To quantify the rate of convergence, we report the quantity, $p_{conv}$, defined as the smallest value of $p$ for which the energy is converged with six decimal places. +The smaller $p_{conv}$, the better the convergence is. Although this is a rough estimate, it is sufficient here for our purpose. -As clearly seen, the speed of convergence is directly related to the magnitude of $\bar{t}_{I_0}$. The -longer the stochastic trajectories remain trapped within the domain, the better the convergence. Of course, when the domain is chosen to be the full space, -the average trapping time becomes infinite. Let us emphasize that the rate of convergence has no reason to be related to the size of the domain. -For example, the domain ${\cal D}(0,3) \cup {\cal D}(1,0)$ has a trapping time for the N\'eel state of 6.2, while -the domain ${\cal D}(0,3) \cup {\cal D}(1,1)$ having almost the same number of states (28 states), has an average trapping time about 6 times longer. Finally, the last column gives the energy obtained for -$E=-1$. The energy is expected to be independent of the domain and to converge to a common value, which is indeed the case here. +As clearly seen, the speed of convergence is directly related to the magnitude of $\bar{t}_{I_0}$. +The longer the stochastic trajectories remain trapped within the domain, the better the convergence. +Of course, when the domain is chosen to be the full space, the average trapping time becomes infinite. +Let us emphasize that the rate of convergence has no reason to be related to the size of the domain. +For example, the domain ${\cal D}(0,3) \cup {\cal D}(1,0)$ has a trapping time for the N\'eel state of 6.2, while the domain ${\cal D}(0,3) \cup {\cal D}(1,1)$ having almost the same number of states (28 states), has an average trapping time about 6 times longer. +Finally, the last column gives the energy obtained for $E=-1$. +The energy is expected to be independent of the domain and to converge to a common value, which is indeed the case here. The exact value, $\cE(E=-1)=-0.75272390...$, can be found at the last row of the Table for the case of a domain corresponding to the full space. -In sharp contrast, the statistical error depends strongly on the type of domains used. As expected, the largest error of $3 \times 10^{-5}$ is obtained in the case of -a single-state domain for all states. The smallest statistical error is obtained for the "best" domain having the largest average -trapping time. Using this domain leads to a reduction in the statistical error as large as about three orders of magnitude, nicely illustrating the -critical importance of the domains employed. +In sharp contrast, the statistical error depends strongly on the type of domains used. +As expected, the largest error of $3 \times 10^{-5}$ is obtained in the case of a single-state domain for all states. The smallest statistical error is obtained for the "best" domain having the largest average trapping time. +Using this domain leads to a reduction in the statistical error as large as about three orders of magnitude, nicely illustrating the critical importance of the domains employed. \begin{table}[h!] \centering -\caption{$N$=4, $U$=12, $E$=-1, $\alpha=1.292$, $\beta=0.552$,$p_{ex}=4$. Simulation with 20 independent blocks and $10^5$ stochastic paths -starting from the N\'eel state. $\bar{t}_{I_0}$ -is the average trapping time for the -N\'eel state. $p_{\rm conv}$ is a measure of the convergence of $\cE_{QMC}(p)$ as a function of $p$, see text.} +\caption{$N=4$, $U=12$, $E=-1$, $\alpha=1.292$, $\beta=0.552$, $p_{ex}=4$. +Simulation with 20 independent blocks and $10^5$ stochastic paths starting from the N\'eel state. +$\bar{t}_{I_0}$ is the average trapping time for the N\'eel state. +$p_{\rm conv}$ is a measure of the convergence of $\cE_{QMC}(p)$ as a function of $p$, see text.} \label{tab1} \begin{ruledtabular} -\begin{tabular}{lcccl} +\begin{tabular}{lrrrl} Domain & Size & $\bar{t}_{I_0}$ & $p_{\rm conv}$ & $\;\;\;\;\;\;\cE_{QMC}$ \\ \hline -Single & 1 & 0.026 & 88 &$\;\;\;\;$-0.75276(3)\\ -${\cal D}(0,3)$ & 2 & 2.1 & 110 &$\;\;\;\;$-0.75276(3)\\ -${\cal D}(0,2)$ & 4 & 2.1 & 106 &$\;\;\;\;$-0.75275(2)\\ -${\cal D}(0,1)$ & 6 & 2.1& 82 &$\;\;\;\;$-0.75274(3)\\ -${\cal D}(0,3)$ $\cup$ ${\cal D}(1,1)$ &14 &4.0& 60 & $\;\;\;\;$-0.75270(2)\\ -${\cal D}(0,3)$ $\cup$ ${\cal D}(1,0)$ &26 &6.2& 45 & $\;\;\;\;$-0.752730(7) \\ -${\cal D}(0,2)$ $\cup$ ${\cal D}(1,1)$ &16 &10.1 & 36 &$\;\;\;\;$-0.75269(1)\\ -${\cal D}(0,2)$ $\cup$ ${\cal D}(1,0)$ &28 &34.7 & 14&$\;\;\;\;$-0.7527240(6)\\ -${\cal D}(0,1)$ $\cup$ ${\cal D}(1,1)$ &18 & 10.1 & 28 &$\;\;\;\;$-0.75269(1)\\ -${\cal D}(0,1)$ $\cup$ ${\cal D}(1,0)$ &30 & 108.7 & 11&$\;\;\;\;$-0.75272400(5) \\ -${\cal D}(0,3)$ $\cup$ ${\cal D}(1,1)$ $\cup$ ${\cal D}$(2,0) &20 & 4.1 & 47 &$\;\;\;\;$-0.75271(2)\\ -${\cal D}(0,3)$ $\cup$ ${\cal D}(1,0)$ $\cup$ ${\cal D}$(2,0) &32 & 6.5 & 39 &$\;\;\;\;$-0.752725(3)\\ -${\cal D}(0,2)$ $\cup$ ${\cal D}(1,1)$ $\cup$ ${\cal D}$(2,0) &22 & 10.8 & 30 &$\;\;\;\;$-0.75270(1)\\ -${\cal D}(0,2)$ $\cup$ ${\cal D}(1,0)$ $\cup$ ${\cal D}$(2,0) &34 & 52.5 & 13&$\;\;\;\;$-0.7527236(2)\\ -${\cal D}(0,1)$ $\cup$ ${\cal D}(1,1)$ $\cup$ ${\cal D}$(2,0) & 24 & 10.8 & 26&$\;\;\;\;$-0.75270(1)\\ -${\cal D}(0,1)$ $\cup$ ${\cal D}(1,0)$ $\cup$ ${\cal D}$(2,0) & 36 & $\infty$&1&$\;\;\;\;$-0.75272390\\ +Single & 1 & 0.026 & 88 &$-0.75276(3)$\\ +$\cD(0,3)$ & 2 & 2.1 & 110 &$-0.75276(3)$\\ +$\cD(0,2)$ & 4 & 2.1 & 106 &$-0.75275(2)$\\ +$\cD(0,1)$ & 6 & 2.1& 82 &$-0.75274(3)$\\ +$\cD(0,3)\cup\cD(1,1)$ &14 &4.0& 60 & $-0.75270(2)$\\ +$\cD(0,3)\cup\cD(1,0)$ &26 &6.2& 45 & $-0.752730(7)$ \\ +$\cD(0,2)\cup\cD(1,1)$ &16 &10.1 & 36 &$-0.75269(1)$\\ +$\cD(0,2)\cup\cD(1,0)$ &28 &34.7 & 14&$-0.7527240(6)$\\ +$\cD(0,1)\cup\cD(1,1)$ &18 & 10.1 & 28 &$-0.75269(1)$\\ +$\cD(0,1)\cup\cD(1,0)$ &30 & 108.7 & 11&$-0.75272400(5)$ \\ +$\cD(0,3)\cup\cD(1,1)\cup\cD$(2,0) &20 & 4.1 & 47 &$-0.75271(2)$\\ +$\cD(0,3)\cup\cD(1,0)\cup\cD$(2,0) &32 & 6.5 & 39 &$-0.752725(3)$\\ +$\cD(0,2)\cup\cD(1,1)\cup\cD$(2,0) &22 & 10.8 & 30 &$-0.75270(1)$\\ +$\cD(0,2)\cup\cD(1,0)\cup\cD$(2,0) &34 & 52.5 & 13&$-0.7527236(2)$\\ +$\cD(0,1)\cup\cD(1,1)\cup\cD$(2,0) & 24 & 10.8 & 26&$-0.75270(1)$\\ +$\cD(0,1)\cup\cD(1,0)\cup\cD$(2,0) & 36 & $\infty$&1&$-0.75272390$\\ \end{tabular} \end{ruledtabular} \end{table} @@ -1027,10 +1035,10 @@ laptop. Of course, it will also be particularly interesting to take advantage of All these aspects will be considered in a forthcoming work. \begin{table} -\caption{$N=4$, $U=12$, and $E=-1$. Dependence of the statistical error on the energy with the number of $p$-components calculated -analytically. Same simulation as for Table \ref{tab1}. Results are presented when a single-state domain -is used for all states and when -${\cal D}(0,1) \cup {\cal D}(1,0)$ is chosen as main domain.} +\caption{$N=4$, $U=12$, and $E=-1$. +Dependence of the statistical error on the energy with the number of $p$-components calculated analytically. +Same simulation as for Table \ref{tab1}. +Results are presented when a single-state domain is used for all states and when $\cD(0,1) \cup \cD(1,0)$ is chosen as main domain.} \label{tab2} \begin{ruledtabular} \begin{tabular}{lcc} @@ -1051,8 +1059,10 @@ $8$ & $2.2 \times10^{-5}$ &$ 0.05 \times 10^{-8}$\\ \begin{table} -\caption{$N=4$, $U=12$, $\alpha=1.292$, $\beta=0.552$. Main domain = ${\cal D}(0,1) \cup {\cal D}(1,0)$. Simulation with 20 independent blocks and $10^6$ paths. -$p_{ex}=4$. The various fits are done with the five values of $E$} +\caption{$N=4$, $U=12$, $\alpha=1.292$, $\beta=0.552$, and $p_{ex}=4$. +Main domain = $\cD(0,1) \cup \cD(1,0)$. +Simulation with 20 independent blocks and $10^6$ paths. +The various fits are done with the five values of $E$} \label{tab3} \begin{ruledtabular} \begin{tabular}{lc} @@ -1072,11 +1082,11 @@ $E_0$ exact & -0.768068...\\ \end{table} \begin{table} -\caption{$N$=4, Domain ${\cal D}(0,1) \cup {\cal D}(1,0)$} +\caption{$N=4$, Domain $\cD(0,1) \cup \cD(1,0)$} \label{tab4} \begin{ruledtabular} \begin{tabular}{cccccc} -$U$ & $\alpha,\beta$ & $E_{var}$ & $E_{ex}$ & $\bar{t}_{I_0}$ \\ +$U$ & $\alpha,\beta$ & $E_\text{T}$ & $E_{ex}$ & $\bar{t}_{I_0}$ \\ \hline 8 & 0.908,\;0.520 & -0.770342... &-1.117172... & 33.5\\ 10 & 1.116,\;0.539 & -0.604162... &-0.911497... & 63.3\\