Approximate expression for the ground-state energy of the two- and three-dimensional Hubbard model at arbitrary filling obtained from dimensional scaling

TL;DR

This paper introduces a nonlinear dimensional scaling method to accurately approximate the ground-state energy of the Hubbard model in 2D and 3D, extending the applicability beyond existing QMC data.

Contribution

The authors develop a nonlinear scaling relation based on linear discrete dimensional scaling, providing accurate ground-state energy estimates for the Hubbard model at arbitrary fillings.

Findings

Accurate approximations of ground-state energy in 2D and 3D Hubbard models.

Predictions made for parameter regimes lacking QMC data.

Validation against auxiliary-field quantum Monte Carlo results.

Abstract

We generalize the linear discrete dimensional scaling approach for the repulsive Hubbard model to obtain a nonlinear scaling relation that yields accurate approximations to the ground-state energy in both two and three dimensions, as judged by comparison to auxiliary-field quantum Monte Carlo (QMC) data. Predictions are made for the per-site ground-state energies in two and three dimensions for $n$ (filling factor) and $U$ (Coulomb interaction) values for which QMC data are currently unavailable.

Figures

Tables

Table 1. Table 1: Expansion coefficients on the right-hand side of Eq. ( 28 ) [ 19 ] . For convenient reference, the last column shows the U = 0 𝑈 0 U=0 per-site kinetic energies, in units of t 𝑡 t .

$d$	$C_1\left(d\right)$	$C_3\left(d\right)$	$\tau (d,0)/t$
1	$-4\ln \left(2\right)$	$9\zeta \left(3\right)$	$-4/{\pi }^a$
2	$-4.631788$	$34.59816$	$-{1.6211}^b$
3	$-6.582948$	$56.61288$	$-{2.0048}^c$

Table 2. Table 2: Per site ground-state energy for the two-dimensional model resulting from the analyses in Secs.  2.2  and 3.2 , for several Coulomb repulsions U 𝑈 U and densities n 𝑛 n . The rightmost column shows the benchmarks set by the available Quantum Monte Carlo. Over the subset of values of U 𝑈 U and n 𝑛 n for which QMC data are available, the mean average error of the linear scaling approximation with respect to the QMC values is 10.5%. The nonlinear generalization, truncated after the cubic term, reduces this to 3.14%.

$U$	$n$	Linear scaling	Nonlinear scaling	QMC
3	0.25	-0.779	-0.793
3	0.50	-1.137	-1.181
3	0.75	-1.150	-1.225
3	0.90625	-1.011	-1.105
3	1.00	-0.879	-0.983
5	0.25	-0.771	-0.788
5	0.50	-1.082	-1.136
5	0.75	-1.003	-1.093
5	0.90625	-0.790	-0.896
5	1.00	-0.619	-0.725
6	0.25	-0.769	-0.786	-0.756a
6	0.50	-1.063	-1.120	-1.093a
6	0.75	-0.952	-1.046	-1.037a
6	0.90625	-0.717	-0.823	-0.817a
6	1.00	-0.535	-0.635	-0.651a
10	0.25	-0.764	-0.782
10	0.50	-1.018	-1.081
10	0.75	-0.829	-0.927
10	0.78	-0.783	-0.883	-0.843b
10	0.89	-0.578	-0.679	-0.644b
10	0.94	-0.472	-0.566	-0.554b
12	0.25	-0.762	-0.781
12	0.50	-1.005	-1.069
12	0.75	-0.794	-0.892
12	0.80	-0.709	-0.810	-0.772c
12	0.875	-0.562	-0.662	-0.628c
12	0.90625	-0.496	-0.591
12	1.00	-0.287	-0.356	-0.369c

Table 3. Table 3: Per site ground-state energy for the three-dimensional model resulting from the analyses in Secs.  2.2  and 3.2 , for several Coulomb repulsions U 𝑈 U and densities n 𝑛 n . Benchmarks set by the available Quantum Monte-Carlo data are tabulated. Over the subset of values of U 𝑈 U and n 𝑛 n for which QMC data are available, the mean average error of the linear scaling approximation with respect to the QMC values is 26%. The nonlinear generalization, truncated after the cubic term, reduces this to 1.65%.

$U$	$n$	linear scaling prediction	nonlinear scaling prediction	QMC
4.0	0.25	-0.954	-0.978
4.0	0.50	-1.320	-1.407
4.0	0.75	-1.299	-1.466
4.0	0.90625	-1.092	-1.303
4.0	1.00	-0.903	-1.124	-1.16a
6.0	0.25	-0.947	-0.974
6.0	0.50	-1.269	-1.369
6.0	0.75	-1.159	-1.346
6.0	0.90625	-0.883	-1.101
6.0	1.00	-0.661	-0.866	-0.89a
8.0	0.25	-0.943	-0.972
8.0	0.50	-1.237	-1.345
8.0	0.75	-1.069	-1.265
8.0	0.90625	-0.755	-0.966
8.0	1.00	-0.516	-0.696	-0.706a
10.0	0.25	-0.940	-0.970
10.0	0.50	-1.215	-1.328
10.0	0.75	-1.009	-1.209
10.0	0.90625	-0.670	-0.874
10.0	1.00	-0.421	-0.580	-0.58a
12.0	0.25	-0.9384	-0.969
12.0	0.50	-1.1996	-1.315
12.0	0.75	-0.9657	-1.167
12.0	0.90625	-0.6112	-0.807
12.0	1.00	-0.3544	-0.495	-0.49a

Equations

$$\frac{P}{P_u}=f_{univ}\left(\frac{T}{T_C}\right).$$

$$\frac{E(d)}{T(d)}=\frac{E(d^{\prime })}{T(d^{\prime })},$$

$$E(d)=\frac{E(d^{\prime })}{T(d^{\prime })}T(d)=\frac{T(d)}{T(d^{\prime })}E(d^{\prime })\approx \frac{T_s(d)}{T_s(d^{\prime })}E(d^{\prime }),$$

$$\frac{E(d)}{T_s(d)}=f_{d,d^{\prime }}\left(\frac{E(d^{\prime })}{T_s(d^{\prime })}\right).$$

$$H=-t\sum _{⟨i,j⟩,\sigma }\text{leavevmode}\left(c_{i,\sigma }^{†}{c}_{j,\sigma }+c_{j,\sigma }^{†}{c}_{i,\sigma }\right)+U\sum _i{\hat{n}}_{i,\uparrow }{\hat{n}}_{i,\downarrow },$$

$$\frac{ϵ(d,U)}{\tau (d,U)}\approx \frac{{ϵ}^{\mathrm{approx}}(d,U)}{\tau (d,U)}\equiv \frac{ϵ(d^{\prime },U)}{\tau (d^{\prime },U)},$$

$$\tau (d,U)=ϵ(d,U)-U\frac{\partial ϵ(d,U)}{\partial U}.$$

$$ϵ(1,U)=-4{\int }_0^{\infty }\frac{J_0(x)J_1(x)}{x\left[1+e^{Ux/2}\right]}\text{leavevmode}dx,$$

$$\frac{ϵ(d,U)}{\tau (d,0)}\approx \frac{ϵ(d^{\prime },U)}{\tau (d^{\prime },0)}.$$

$$ϵ(d,U)=\tau (d,0)f_d\left(x(U)\right),$$

$$x(U)\equiv \frac{ϵ(1,U)}{\tau (1,0)}.$$

$$\displaystyle\frac{\partial\epsilon(d,U)}{\partial U}=\dfrac{\epsilon(d,U)}{\tau(1,0)f_{d}\big{(}x(U)\big{)}}\dfrac{df_{d}}{dx}\frac{\partial\epsilon(1,U)}{\partial U}.$$

$$\displaystyle\tau(d,U)=\epsilon(d,U)\left(1-U\dfrac{1}{\tau(1,0)f_{d}\big{(}x(U)\big{)}}\dfrac{df_{d}}{dx}\frac{\partial\epsilon(1,U)}{\partial U}\right),$$

$$\frac{ϵ(d,U)}{\tau (d,U)}=\frac{ϵ(1,U)}{ϵ(1,U)-U\frac{\partial ϵ(1,U)}{\partial U}{y}_d(U)}$$

$$\displaystyle y_{d}(U)\equiv\dfrac{x(U)}{f_{d}\big{(}x(U)\big{)}}\dfrac{df_{d}}{dx}.$$

$$\displaystyle f_{d}\big{(}x(U)\big{)}\approx x(U),$$

$$\frac{ϵ(d,U)}{\tau (d,U)}\approx \frac{ϵ(1,U)}{\tau (1,U)}.$$

$$ϵ(d,U,n)\approx {ϵ}_{\text{mbox}L}(d,U,n)=\tau (d,0,n)x(U,n)(n\le 1),$$

$$x(U,n)=\frac{ϵ(1,U,n)}{\tau (1,0,n)}.$$

$$ϵ(d,U,n)=ϵ(d,U,2-n)+U(n-1).$$

$$ϵ(1,U\to \infty ,n)=\frac{ϵ(1,0,2n)}{2}=\frac{\tau (1,0,2n)}{2},$$

$$x(U\to \infty ,n)=\cos \left(\frac{\pi n}{2}\right),$$

$$ϵ(d,U)\approx {ϵ}_{\text{mbox}NL}(d,U)=\tau (d,0)\sum _{s=0}^N\text{leavevmode}{\alpha }_s{x}^s.$$

$$ϵ(d,U\to 0,n)=\tau (d,0,n)+\frac{U}{4}{n}^2$$

$$\frac{U}{4}{n}^2=(x-1)\tau (1,0).$$

$$ϵ(d,U\to 0)=\tau (d,0)+(x-1)\tau (1,0),$$

$${\frac{\partial ϵ(d,U)}{\partial x}|}_{x=1}=\tau (1,0),$$

$$\displaystyle\epsilon(d,U\gg t)=C_{1}(d)\frac{t}{U}+C_{3}(d)\frac{t^{3}}{U^{3}}+O\Big{(}\frac{t^{5}}{U^{5}}\Big{)},$$

$$ϵ(d,U\gg t,n=1)=\tau (1,0,n=1)\frac{C_1(d)}{C_1(1)}x+O(x^3),$$

$${\frac{\partial ϵ(d,U)}{\partial x}|}_{x=0}=\tau (1,0)\frac{C_1(d)}{C_1(1)}.$$

Full text

1 Introduction

Scaling relations pervade physics. In their simplest form, they arise, e.g., in dimensional analysis, as formalized by Buckingham’s theorem [1], according to which quantities pertaining to one of two geometrically and dynamically similar systems can be inferred by measuring or computing the corresponding quantity in the other, if the basic variables are suitably scaled.

Similar, but more sophisticated, scaling relations appear in statistical physics, for example in the context of universality. In critical phenomena [2], the dependence of some property $P$ on the temperature near the critical point, expressed in multiples of an appropriate unit $P_{u}$, is given by universal scaling functions of $T/T_{C}$, where $T_{C}$ is the critical temperature,

[TABLE]

Different physical systems in the same universality class can have different critical temperatures, but for the same ratio $T/T_{C}$, they present the same universal properties.

In Ref. [3], a different kind of scaling, denoted discrete dimensional scaling (DDS), was identified. DDS relates dimensionless ratios computed in systems differing in their dimensionality $d$.

Systems of different dimensionalities are normally considered to be completely distinct objects. In statistical mechanics, systems with different dimensionalities belong to different universality classes, characterized by different critical exponents [4, 5]. In condensed-matter physics, one-, two- and three-dimensional (1D,2D,3D) systems are known to display very different phenomenology. This behaviour is reflected by commonly used model Hamiltonians, such as the repulsive Hubbard model (HM) which in 1D displays a Luttinger-liquid and a Mott-insulator phase, in 3D follows the Fermi-liquid paradigm, and in 2D has a complex ground-state that is not fully understood, but characterized by strong antiferromagnetic correlations and possible superconducting and charge fluctuations [6, 7]. Little, if any, important physics remains the same on crossing the dimensional divide.

On the other hand, remarkable progress has been made in specific situations by exploiting relations between systems and quantities in different dimensionalities. In statistical physics, for example, the $\epsilon$-expansion [8] treats $d$ as a continuous variable and expands around $d=4$. In condensed-matter physics, dynamical-mean-field theory (DMFT) treats many-body models in physical dimensionalities by starting from the limit $d\to\infty$ [9].

The enormous success of the $\epsilon$-expansion and of DMFT in their respective domains notwithstanding, there is little work on systematically exploring connections between different dimensionalities. However, in addition to these two concrete examples, general mathematical considerations suggest that such interdimensional relations may be more ubiquitous than usually thought.

Nonrelativistic quantum systems are described by Schrödinger’s equation, a linear partial differential equation (PDE). The behaviour of PDEs and their solutions is determined by their dimensionality $d$, in conjunction with several other defining features, such as its characteristics (elliptic, parabolic, hyperbolic or mixed), order (highest derivative), degree (power of the highest derivative), linearity or nonlinearity (power of the solution function), and whether it has constant or variable coefficients, and a real or complex solution space, among others [10]. The Schrödinger equations defined by a Hamiltonian in dimensionalities $d$ and $d^{\prime}$ are the same in all these key aspects, except for $d$ itself. This strongly suggests that a system in $d$-dimensions cannot be arbitrarily different from its counterpart in $d^{\prime}$-dimensions and that some connection between their properties must exist.

A more intuitive way of thinking about this connection is to note that low-dimensional systems are the basic building blocks out of which higher-dimensional systems are composed: sites connect to form chains, parallel chains define a plane and stacked planes a three-dimensional lattice. Again, some features of the building blocks should be expected to be preserved in the constituting system.

The quantum-mechanical virial theorem for bound Coulomb systems [11] provides a beautiful example: since in any dimension $d$ the total energy is equal to the negative of the kinetic energy, $E(d)=-T(d)$, we have the dimensional scaling relation

[TABLE]

which is a first example of DDS, as it relates a ratio in dimension $d$ to the same ratio in dimension $d^{\prime}$.

For the Fermi liquid, relation (2) holds exactly, for any $d$, $d^{\prime}$ and density $n$, due to the virial theorem. Interestingly, the same relation also holds in two limiting cases of the Hubbard model, as it is by construction exact at $U=0$ and, due to the Hellmann-Feynman theorem, also at $U\to\infty$. Thus, Eq. (2) provides a concrete example of how physical or mathematical constraints (the virial theorem and the Hellmann-Feynman theorem, respectively) establish connections between quantities in different dimensionalities.

The basic linear relation (2) can be rewritten as

[TABLE]

where $T$ denotes interacting and $T_{s}$ noninteracting kinetic energies, which are easier to calculate. In this latter form, dimensional scaling becomes a tool for approximating the ground-state energy in $d$-dimensions in terms of the ground-state energy in $d^{\prime}$-dimensions by scaling the latter by a simple energy ratio.

What makes this observation potentially useful in studying ground-state properties, is that numerical calculations for the Hubbard model have shown that Eq. (3), even in its version employing noninteracting kinetic energies, provides fair approximations also for intermediate values of $U$ (including the physically realistic range), both at $n=1$ [3] and for $n\neq 1$ [12].

This suggests that we use dimensional scaling to build approximations for physical quantities in a dimensionality in which calculations are easier to perform and scale them to a different, harder-to-describe, dimensionality. Of course, the difference between the physics in the two involved dimensionalities must be approximately accounted for by the scaling factor, which cannot always be expected to be of the simple form in Eq. (3).

Motivated by these considerations, we develop, in the present paper, nonlinear generalizations of (2) and (3) for the Hubbard model, which are of a general form inspired by Eq. (1):

[TABLE]

We start by providing, in Sec. 2, a justification for the previously [3] ad hoc replacement of the interacting by the non-interacting kinetic energy on the right-hand side of Eq. (3). In section 3, we then construct the corresponding nonlinear dimensional scaling relation for the half-filled model and for arbitrary filling. Our main result is an explicit approximation for the scaling function $f_{d,d^{\prime}}\left(E(d^{\prime})/T_{s}(d^{\prime})\right)$, with $d^{\prime}=1$, which we use to calculate ground-state energies for $d=2$ and $d=3$. A detailed comparison with quantum Monte Carlo (QMC) data in the literature shows that the nonlinear generalization significantly improves the accuracy of scaling with noninteracting kinetic energies, both at $n=1$ and for $n\neq 1$, for physically realistic values of $U$. The nonlinear scaling relation is simple enough to allow us to make predictions for the per-site ground-state energy of the two- and three-dimensional Hubbard model in parameter ranges that have not yet been explored by QMC, and where Monte Carlo techniques may be hard to apply.

2 Linear dimensional scaling

The Hubbard model Hamiltonian reads

[TABLE]

where $c_{i,\sigma}^{\dagger}$ creates a fermion on site $i$ with spin $\sigma$ $(\uparrow,\downarrow)$. The first sum in (5) runs over all pairs of first neighbors and represents the kinetic energy. The hopping amplitude $t$, which sets the kinetic-energy scale, will define our unit of energy. For dimensionalities $d=2$ or 3, we will consider the model on a square or cubic lattice, respectively. In the second sum, $\hat{n}_{i,\sigma}=c_{i,\sigma}^{\dagger}c_{i,\sigma}$, so that any doubly occupied site adds the intra-site interaction $U$ to the system energy. We are interested in repulsive interactions ($U>0$).

In one dimension, this model is exactly solvable by the Bethe Ansatz [13, 14]. Away from half-filling, it is an example of a Luttinger liquid [15], with spin-charge separation. At half-filling, the system is a Mott insulator. An energy gap separates the ground-state from the charge excitations for any $U>0$ [13]. For $d>1$, no exact solution is known. The most reliable information about the properties of the model has been obtained numerically, mostly from QMC calculations [16].

2.1 Half-filled case

The dimensional scaling approach was originally introduced [3] as a procedure to obtain approximations to the ground-state energy of the half-filled Hubbard model in 2 and 3 dimensions by evaluating the approximate linear scaling relation

[TABLE]

for $d^{\prime}=1$. Here, $\epsilon(d,U)$ is the ground-state energy per site for the $d$-dimensional model with intra-site interaction $U$ at half-filling, $\epsilon^{\rm approx}(d,U)$ is an approximation to it defined by this equation, and $\tau(d,U)$ is the interacting kinetic energy per site.

As a direct consequence of the Hellmann-Feynman theorem, the kinetic energy can be calculated from the full ground-state energy as

[TABLE]

The energies in Eqs. (6) and (7) can be computed exactly for $d\to 0$ (approximated as the smallest possible Hubbard model, namely a Hubbard dimer with open boundary condition) and for $d=1$. In particular, for $d=1$ at half-filling [13]

[TABLE]

where $J_{0}$ and $J_{1}$ are zero and first-order Bessel functions.

Figure 1(a) displays the ratios $\epsilon(d,U)/\tau(d,U)$ for $d=0,1,2$ as functions of the ratio $U/t$, computed independently. For $d=0,1$, the calculation can be done analytically; for $d=2$ we have resorted to numerical results from auxiliary-field quantum Monte Carlo (AFQMC) calculations reported in Ref. [16]. From the weak dependence on the dimensionality $d$ in Fig. 1, we conclude that the simple scaling relation (6) is a fairly satisfactory approximation. Clearly, however, there is room for improvement.

Given $\epsilon(d^{\prime},U)$ for $d^{\prime}=0$ or $d^{\prime}=1$, and $\tau(d,U)$, Eq.  (6) yields $\epsilon(d,U)$ for $d>1$. Unfortunately, this approximation has no practical value, since $\tau(d,U)$ is only known in special limits. Computing the interacting kinetic energy proves as difficult as evaluating the total energy.

To overcome this limitation, Ref. [3] proposed a further approximation: to substitute the noninteracting kinetic energies $\tau(d,0)$ for $\tau(U,d)$ in the scaling relation. A simpler, albeit somewhat less accurate expression resulted:

[TABLE]

Compared to Fig. 1(a), the plots of the ratios $\epsilon(d,U)/\tau(d,0)$ for $d=0,1,2$ in Fig. 1(b) show more pronounced dependence on dimensionality. Nonetheless, the green balls representing the ratio for $d=2$, and even the red balls representing $d=3$ are still fairly close to the corresponding black open squares ($d=0$) 0 or orange rhombuses ($d=1$). The separations are smaller than one would expect from the apparently crude approximation converting Eq. (6) into Eq. (9). Moreover, as we shall see, that approximation warrants good agreement with Quantum Monte Carlo data describing the $U$ dependence of $\epsilon(U,d)$ for $d=2$ and $3$. Brief discussion of the connection between Eqs. (6) and (9) seems therefore appropriate.

One can easily see that Eq. (9), like Eq. (6), become exact in the $U=0$ and $U\to\infty$ limits. $U=0$ makes the numerator on each side of Eq. (9) equal to the denominator. For $U\to\infty$, the numerator on each side vanishes [18, 19], while the denominator is nonzero.

For finite $U$, a more elaborate argument is required. Without sacrificing generality, we may write the interacting energy in $d$-dimensions as a function of the ratio between the one-dimensional energy and the noninteracting one-dimensional energy:

[TABLE]

with the shorthand

[TABLE]

The subindex $d$ on the right-hand side of Eq. (10) is a reminder that the function $f_{d}$ depends on the dimensionality. Clearly, $f_{1}(x)\equiv x$.

Differentiate, then, both sides of Eq. (10) with respect to $U$, and divide each side of the resulting equality by the corresponding side of Eq. (10). The following expression results:

[TABLE]

Substitution of the right-hand side for $\partial\epsilon(d,U)/\partial U$ on the right-hand side of Eq. (7) then yields the equality

[TABLE]

and multiplication of both sides by $\epsilon(1,U)$ followed by straightforward manipulation shows that

[TABLE]

where

[TABLE]

Equation (14) is formally exact—formally, since the function $f_{d}(x(U))$ is unknown, except at $d=1$. To exploit it, we need an approximate form for $f_{d}(x)$ ($d>1$), which must satisfy two boundary conditions: $f_{d}(0)=0$, since $\epsilon(d,U\to\infty)=0$, and $f_{d}(1)=1$, since $\epsilon(d,U=0)=\tau(d,0)$. The simplest approximation consistent with those conditions amounts to assuming $f_{d}(x)$ independent of $d$:

[TABLE]

which makes Eq. (10) equivalent to Eq. (9).

Under this approximation, Eq. (15) reduces to $y_{d}(U)=1$, and the denominator on the right-hand side of Eq. (14) reproduces the right-hand side of Eq. (7) for $d=1$. In other words, the approximation $f_{d}(x)\approx x$ reduces Eq. (14) to the $d^{\prime}=1$ version of Eq. (6):

[TABLE]

Conversely, we can see that Eq. (17) implies $y_{d}(x)=1$. Equation (15) then becomes a first-order ordinary differential equation for $f_{d}(x)$, which can be easily solved under the boundary conditions $f(0)=0$ and $f(1)=1$ to show that $f_{d}(x)=x$.

In summary, Eq. (7) implies that Eq. (9), with $d^{\prime}=1$, holds if and only if Eq. (17) is valid. The physical arguments in Sec. 1 indicate that Eq. (6) is a good approximation. We therefore expect Eq. (9) to be a good approximation, which justifies the expedient adopted in Ref. [3], which substituted the latter equality for the former.

One might hope that Eq. (6) be exact, in which case Eq. (9) would also be exact. A nonetheless shows that the two equalities are inequivalent: for large $U$, the former is more reliable than the latter. Even without reference to numerical evidence, we can see that the linear approximation (16) is incompatible with Eq. (7)—a corollary of the Helmann-Feynman Theorem. The numerical data, more specifically the separation between the golden balls and the red lines in Figs. 2 and 3, corroborate this conclusion and offer an estimate of the deviations introduced by the linear approximation (16). This leads us to the discussion in Sec. 3, which considers nonlinear approximations. Before that, however, a brief digression is necessary, to recall that the linear approximation has been extended to arbitrary filling [12].

2.2 Arbitrary filling

To generalize Eq. (10), Ref. [12] let the density (site occupation) $n$ be an additional variable and wrote the expression

[TABLE]

where

[TABLE]

The restriction $n\leq 1$ on the right-hand side of Eq. (18) is of no practical consequence, for a particle-hole transformation maps the $n$-particle per site problem onto the $2-n$ problem [13]:

[TABLE]

For fixed density $n\leq 1$, in the non-interacting limit, Eq. (18) becomes exact, since $x(0,n)=1$. As $U\to\infty$, the per-site ground-state energy for $d=1$ is given by the expression

[TABLE]

which can be physically understood by considering that the limit $U\to\infty$ enforces a special version of the Pauli principle, one that allows single-particle states to be vacant or occupied by at most one electron, with either $s_{z}=1/2$ or $s_{z}=-1/2$. Except for this restriction, the single-particle states are the eigenstates of the $U=0$ Hamiltonian. In the ground state, each single-particle energy below the Fermi level is singly occupied. The Fermi level for $n$ electrons per site hence coincides with the Fermi level for a noninteracting system with $2n$ electrons per site. The ground-state energy is $\epsilon(1,U\to\infty,n)=-\frac{2}{\pi}\sin(\pi n)$, which implies that

[TABLE]

so that the interval $[0,\infty)$ for $U$ is mapped onto the interval $(\cos\left(\frac{\pi n}{2}\right),1]$ in the scaling variable $x(U,n)$.

The simple physical argument following Eq. (21) is valid (i) for any density $n$ at $d=1$ [20]111We consider here the thermodynamic limit. For finite chains, only with open boundary condition is Eq. (21) correct. To see how periodic boundary condition invalidates the argument, consider a trimer and a particular configuration comprising a spin-up electron at the first site, a vacancy at the second, and a spin-down electron at the third, which can be denoted $|A\rangle=|\uparrow,0,\downarrow\rangle$. We can now easily construct the sequence of states including $|A\rangle$ in which each state is coupled to the following one by the $U\to\infty$ Hamiltonian. Under open boundary condition, the sequence comprises three elements: $|0,\uparrow,\downarrow\rangle$, $|\uparrow,0,\downarrow\rangle$, $|\uparrow,\downarrow,0\rangle$, which match the sequence of spinless-electron states in which each state is coupled to the next by the $U=0$ Hamiltonian: $|0,1,1\rangle$, $|1,0,1\rangle$, $|1,1,0\rangle$. Under periodic boundary condition, the sequence of spinless-electron states is unchanged, but the sequence including state $|A\rangle$ now has six states: $|0,\uparrow,\downarrow\rangle$, $|\uparrow,0,\downarrow\rangle$, $|\uparrow,\downarrow,0\rangle$, $|0,\downarrow,\uparrow\rangle$, $|\downarrow,0,\uparrow\rangle$, $|\downarrow,\uparrow,0\rangle$. The same argument can be extended to lattices of arbitrary sizes and shows that, under periodic condition with $n\neq 1$, the spectrum of the $U\to\infty$ Hamiltonian cannot be mapped onto that of the $U=0$ Hamiltonian for spinless electrons. and (ii) at half-filling for $d>1$ [18, 21].

3 Nonlinear dimensional scaling

Although the linear dimensional scaling in Eq. (18) provides a simple approximation to $\epsilon(d,U,n)$ for $d=2$ and $d=3$, it is clear from Figs. 2, 3 and 4 that there is much room for improving accuracy. As in the previous section, we start by considering the half-filled case and then develop an approximation for arbitrary filling.

3.1 Half-filled case

Significant improvement results from substituting a more general, polynomial expression for the linear dependence in Eq. (9), that is, from writing the expression

[TABLE]

Since the linear approximation $\epsilon(d,U)\approx\tau(1,U)x(U)$ becomes exact at both the non-interacting ($U=0$, $x=1$) and strongly interacting ($U\to\infty$, $x=0$) limits, a Taylor expansion around either extreme would be unwieldy. It is preferable to treat the $U=0$ and the $U\to\infty$ limits on equal footing. To this end, we will take advantage of exact expressions for $\partial\epsilon/\partial U$ at the $U=0$ and $U\to\infty$ to determine the coefficients $\alpha_{2}$ and $\alpha_{3}$ and derive a cubic ($N=3$) form for $\epsilon_{\mbox{\tiny NL}}$.

As $U\to 0$, the Hartree-Fock approximation

[TABLE]

becomes asymptotically exact.

For $d=1$, Eq. (24) can be combined with the definition (11) to yield the expression

[TABLE]

Substitution of the right-hand side of Eq. (25) for $U/4$ on the right-hand side of Eq. (24) then yields the relation

[TABLE]

which determines the slope at $x=1$ ($U=0$):

[TABLE]

a result that is valid for all $n$.

To determine the large $U$ slope, we recall that, in this limit, the half-filled Hubbard model maps onto the $S=1/2$ antiferromagnetic Heisenberg model with effective interaction $J=4/U$. An accurate asymptotic expansion then describes the ground-state energy [19]:

[TABLE]

where $C_{1}(1)=-4\ln(2)$, $C_{3}(1)=9\zeta(3)$, where $\zeta$ is Riemann’s zeta function, $C_{1}(2)=-4.631788$, $C_{3}(2)=34.599816$, $C_{1}(3)=-6.582948$, and $C_{3}(3)=56.612880$.

Equation (28) can be rewritten in the form

[TABLE]

which implies that

[TABLE]

From the exact values of $\epsilon(d,U)$ and its slope at $x=0$ and $x=1$, we determine the coefficients $\alpha_{s}$ in Eq.(23):

[TABLE]

Figures 2 and 3 show that the cubic expression yields excellent agreement with QMC results for the half-filled Hubbard model in 2 and 3 dimensions, for all $U$. The mapping between the $[0,\infty)$ range of $U$ and the finite range $(0,1]$ of the scaling variable $x(U)$ has allowed us to exploit exact results at both extremes, to obtain a low-order polynomial yielding accurate results over the entire $U$ axis.

In the two-dimensional case, the deviations relative to the AFQMC data shrink from 8.6–22.4% for the linear to 1–3.6% for the cubic approximation. In three dimensions, the drop is from 22-27% to 0.1–3.2%. Moreover, the procedure we have described can be systematically improved: with higher order terms on the right-hand side of Eq. (28) as input, one can extend the algebra leading to Eqs. (31)-(34) to obtain higher order coefficients $\alpha_{s}$ ($s>3$) and increment the degree $N$ of the polynomial describing $\epsilon(d,U)$.

3.2 Arbitrary filling

The significant improvement afforded by nonlinear dimensional scaling at half-filling prompts us to generalize the approach in Sec. 3.1 to arbitrary filling. The following discussion leads to an approximate expression for the ground-state energy at any dimensionality. Unfortunately, due to the scarcity of accurate numerical data away from half-filling for the $3$-dimensional model, we can only check the accuracy of our expression for $d=2$.

To extend the analysis in Sec. 3 to $n<1$ we simply substitute $\tau(1,0,n)$ for $\tau(1,0)$, and $\tau(d,0,n)$ for $\tau(d,0)$ in Eqs. (26) and (29) and in Eqs. (31)-(34). Again to third order, we then have that

[TABLE]

where the $\alpha_{s}(n)$ ($s=0,\ldots,3$) are the coefficients given by Eqs. (31)–(34) after the aforementioned substitutions.

As explained in Sec. 2.2, for $n\leq 1$ the scaling variable $x$ is constrained to the interval $\cos\left(\frac{\pi n}{2}\right)\leq x(U,n)\leq 1$. By construction, both the per-site energy and its derivative $\partial\epsilon(d,U,n)/\partial x$ become exact at $x=1$, for all $n$. In contrast, at the minimum $x=\cos\left(\frac{\pi n}{2}\right)$ ($U\to\infty$), neither the per site energy, nor the derivative resulting from the same approximation are guaranteed to yield exact results. There are two exceptions, however: (i) at half-filling, we recover the results described in Sec. 3.1, and (ii) as $n\to 0$, the domain of $x$ reduces to a single point, $x=1$, for which the approximation $\epsilon(1,d,U)\approx\epsilon_{\mbox{\tiny NL}}(d,U,n)$ becomes exact.

Figure 4 ratifies these considerations. The agreement between the blue solid line representing the nonlinear approximation and the golden balls representing the AFQMC data [22] is excellent in the vicinity of the extremes, $n=0$ and $n=1$, while the small deviations around $n=0.5$ become more clearly visible as $U/t$ grows.

The same reasoning shows that, at low densities, the nonlinear and linear approximations become equivalent. At half-filling, by contrast, the corrections to the linear approximation defined by the terms proportional to $\alpha_{2}$ and $\alpha_{3}$ on the right-hand side of Eq. (35) are sizeable, even for $U$ as large as $8t$.

The AFQMC data taken as benchmarks, Fig. 4 shows that linear approximation overestimates the per-site ground-state energy, while the nonlinear approximation underestimates it. The same trend is apparent in Table 2, which displays more data for the two-dimensional model. A number of results without matching benchmarks are shown, to encourage additional QMC computations. With the same motivation, we present in Table 3 results for the three-dimensional model, for which we have been unable to find QMC results representative of the thermodynamic limit away from half-filling.

The deviations separating the blue lines from the golden balls in Fig. 4 leave room for further improvement of the nonlinear scaling procedure. Unfortunately, progress along that line of work calls for exact results for the $U\to\infty$ limits of the ground-state energy and its derivative with respect to $U$ for $d>1$ and $n\neq 1$, which are unavailable. We have tested a nonlinear scaling approximation for $d=2$ and $n<1$ built upon Eq. (21), which offers an (uncontrollable) approximation for the per-site ground-state energy. The resulting expression for $\epsilon_{\mbox{\tiny NL}}(2,U,n)$ is not shown because it is inferior to the one represented by the blue lines in Fig. 4.

4 Conclusion

We have presented a rationale justifying, for the first time, an approximation of capital importance in the dimensional scaling approaches introduced in Refs. [3] and [12]: the substitution of the noninteracting kinetic energy for the (unknown) interacting kinetic energy on the right-hand side of Eq. (6), which defines the scaling variable.

The same line of reasoning led to a nonlinear approximation for the per-site ground-state energy. Exact results for the per-site ground-state energy for the two- and three-dimensional Hubbard models at half-filling gave access to the energy and its derivative with respect to the Coulomb repulsion $U$ in the noninteracting ($U=0$) and strongly interacting ($U\to\infty$) limits. The resulting cubic approximations for the per-site ground-state energies, devoid of adjustable parameters, yielded substantial improvement over the agreement with accurate Quantum Monte Carlo data for the two- and three-dimensional models at half-filling.

Away from half-filling, a region of the parametric space of the two- and three-dimensional models for which exact results are unavailable in the $U\to\infty$ limit, the agreement between our results and the Quantum Monte-Carlo data for the two-dimensional model is still good. Tables 2 and 3 displays results generated by the linear and nonlinear scaling approaches, for a number of ($n$,$U$) pairs, with a view to motivating additional Monte-Carlo computations for the two- and three-dimensional models. In the absence of such QMC calculations, our results constitute predictions for the ground-state energy at the reported $(n,U)$ pairs. The simplicity of the cubic scaling expression permits to easily generate further data for other $(n,U)$ pairs in two and three dimensions.

In addition to providing predictions that can be checked by QMC, the availability of a relatively simple explicit expression for the ground-state energy as a function of $n$ and $U$ opens up the possibility of an occupation-number functional approach [26] to the inhomogeneous Hubbard model in $d>1$, within the local-density approximation, which could treat finite lattices of large sizes, subject to arbitrary external potentials. We plan to explore this application of the dimensional scaling approach in future work.

On the conceptual side, it remains a challenge to verify if the nonlinear expression for the Hubbard model can be further improved, by incorporating additional exact properties into higher-order terms, or exploiting the variational principle. More generally, it would be rewarding to search for further contraints, similar to the virial theorem for the Fermi liquid and the Hellmann-Feynman theorem for the Hubbard model, that enable interdimensional scaling laws for other systems and physical quantities.

Acknowledgments

This work was supported by the CNPq (grant number 312658/2013-3).

Appendix A Inequivalence between the linear approximations

To show that Eqs. (6) and (9) are inequivalent, we examine both equalities for large $U$. In that domain, the ground state energy is given by Eq. (28). The corresponding kinetic energy can be obtained from Eq. (13):

[TABLE]

and therefore

[TABLE]

while

[TABLE]

Substitution of the coefficients and kinetic energies in Table 1 on the right-hand side of Eq. (38) now yields the following expressions:

[TABLE]

[TABLE]

and

[TABLE]

We can now compare Eq. (37), which states that, to $O(t/U)$, the ratio $\epsilon(d,U\gg t)/\tau(d,U)$ is a constant, independent of $d$, with Eqs. (39)-(41), which show that the linear coefficient in the $t/U$ expansion varies with $d$. This distinction is sufficient to show that Eqs. (6) and (9) cannot be equivalent.