Simple parametrization for the ground-state energy of the infinite Hubbard chain incorporating Mott physics, spin-dependent phenomena and spatial inhomogeneity

Simple analytical parametrizations for the ground-state energy of the one-dimensional repulsive Hubbard model are developed. The charge-dependence of the energy is parametrized using exact results extracted from the Bethe-Ansatz. The resulting parametrization is shown to be in better agreement with highly precise data obtained from fully numerical solution of the Bethe-Ansatz equations than previous expressions [Lima et al., Phys. Rev. Lett. 90, 146402 (2003)]. Unlike these earlier proposals, the present parametrization correctly predicts a positive Mott gap at half filling for any U>0. The construction is extended to spin-dependent phenomena by parametrizing the magnetization-dependence of the ground-state energy using further exact results and numerical benchmarking. Lastly, the parametrizations developed for the spatially uniform model are extended by means of a simple local-density-type approximation to spatially inhomogeneous models, e.g., in the presence of impurities, external fields or trapping potentials. Results are shown to be in excellent agreement with independent many-body calculations, at a fraction of the computational cost.


Introduction
The archetypical strong-correlation phenomenon is the Mott insulator [1]. It is well known that the first-principles description of such systems by means of densityfunctional theory (DFT) encounters severe difficulties. The single-particle (Kohn-Sham) gap calculated with standard DFT methodology is not the same as the many-body gap, even in principle and if no approximations at all are made during the calculation. In the particular case of the Mott insulator, it is known that a proper description of the Mott gap is obtained by adding to the single-particle Kohn-Sham gap a correction arising from the derivative discontinuity of the exchange-correlation functional [2][3][4][5][6][7].
Many modern density functionals do have an implicit derivative discontinuity due to their orbital dependence, but this affords at best an incomplete description of the Mott state and the Mott gap [3,6,7]. Only very few density functionals have an explicit discontinuity as a function of the density, among them the 2-electron functional of Mori-Sanchez, Cohen and Yang [3] which is based on the so-called flat-plane condition, devised by the same authors, [3][4][5] and the Bethe-Ansatz LDA (BALDA) of Lima et al. [8,9], which is based on an approximate analytical parametrization of the Bethe-Ansatz solution for the ground-state energy of the one-dimensional Hubbard model that becomes exact in several important limits. Both of these functionals do properly account for the Mott insulator. Only the BALDA, however, has been parametrized in a way that allows its application to a wide range of physical systems. As a consequence, BALDA and variations thereof has been applied to inhomogeneous correlated manyelectron systems as well as to correlated many-atom systems in optical lattices and traps [9][10][11][12][13][14][15][16][17][18].
However, in these applications it has become clear that the parametrization by Lima et al., also referred to as LSOC parametrization, after the initials of its developers, still does not provide a fully correct description of the Mott gap. In particular, at small U, the LSOC expression for E c has a derivative discontinuity, but the resulting Mott gap is negative, effectively predicting the Mott state to be a (strange) metal instead of an insulator. It is not trivial to correct this behaviour, as any change to the LSOC expression must preserve the exact limits and properties already built into it. Thus, instead of merely algebraic adjustments, it becomes necessary to understand and improve the physics missing from the LSOC parametrization in this regime.
A key aspect of the Mott insulator is that it is a nonmagnetic state of matter, i.e., the insulating nature is not the result of antiferromagnetism. Therefore, the solution to the problem just described must be expressed, within DFT, in terms of the charge-density only, and cannot make use of spin densities and spin-density-functional theory (SDFT). On the other hand, many correlated systems do have magnetic phases. Therefore, a more complete description of strong correlations must properly account for both, the insulating state and various types of magnetic states, as well as their possible coexistence. While such description is possible within Bethe-Ansatz based SDFT by performing a fully numerical solution of the spin-resolved Bethe-Ansatz equations and interpolating between the resulting data points every time the exchange-correlation energy needs to be evaluated, this procedure is very inconvenient for Kohn-Sham calculations, where the functional must be evaluated hundreds or thousands of times during the iterations towards selfconsistency. A simple parametrization of the spin dependence would allow straightforward application of BA-DFT methodology to spindependent phenomena in electronic systems and hyperfine-label-dependent phenomena in optical lattices.
The present paper reports progress along both of these lines. In Section 2 we identify a shortcoming of the standard parametrization used in BALDA and propose a simple, ad hoc but physically motivated, variation of it that is more accurate and whose Mott gap properly remains positive for all values of the interaction U. The revised parametrization also considerably improves the description of the metallic (e.g., Luttinger liquid) phases. In Section 3 we employ exact analytical results extracted from the Bethe Ansatz to further generalize this revised parametrization to spin-dependent situations. In Section 4 we use a simple local-density type approximation to extend our results to spatially nonuniform systems. Density-Matrix Renormalization Group (DMRG) and Lanczos calculations are performed to test and validate this approximation. Section 5 contains a brief summary.

Improved description of the Mott gap
The task at hand is to obtain a simple and accurate analytical approximation for the ground-state energy E 0 of the inhomogeneous one-dimensional Hubbard model (1DHM), wheren iσ =ĉ † iσĉ iσ is the spin-resolved particle-density operator at site i,ĉ † iσ andĉ iσ are fermionic creation and annihilation operators, t is the hopping parameter (taken to be the unit of energy), U the on-site interaction and V i an on-site potential that makes the system spatially inhomogeneous. In applications to electrons in crystal lattices, V i can describe inequivalent atoms in the lattice, while in applications to cold atoms in optical lattices it accounts for the trapping potential.
The commonly used expression for the per-site ground-state energy (e 0 ) of the 1DHM is the LSOC parametrization, given by [8] where n = n ↑ + n ↓ is the charge density and the interaction U enters e 0 (n, U) through the interaction function β(U), which is determined from By construction, this expression becomes exact for U → 0 and any n (where β = 2), for U → ∞ and any n (where β = 1), and for n = 1 and any U (where 0 ≤ β ≤ 1), and provides a reasonable approximation to the full Bethe-Ansatz solution inbetween.
In the LSOC parametrization, the interaction function β(U) is independent of the particle density and of spin. This independence is algebraically very convenient, as it allows one to determine β(U) outside the selfconsistency cycle of DFT instead of having to recalculate it any time the charge (or spin) density changes. However, it is physically incorrect, as the relation between the bare interaction parameter U and the correlation energy E c must depend on the charge density, e.g. through screening. As a consequence of the density-independence of β, LSOC retains the sinusoidal density-dependence of the energy, which is correct only at U = 0 and U → ∞, also for all intermediate values of U.
Thus, it becomes necessary to modify the interaction function in a way that allows it to change with the density. The additional density dependence, however must not spoil the behaviour in the three limits exactly recovered by LSOC. This condition severely restricts the modifications that could possibly be made to the LSOC parametrization. The form we adopt is where β(n, U) = β(U) α(n,U ) and α(n, U) = n 3 √ U /8 . We stress that this particular form has not been derived from first principles but is a physically motivated ad hoc modification designed to restore the density-dependence of the interaction function through the replacement β(U) → β(n, U), while preserving all exact limits obeyed by the LSOC expression. The specific form chosen for the exponent α(n, U) is a consequence of the later generalization to spin-dependent phenomena, as explained in Sec. 3. Figure 1 compares the present parametrization (4) to data obtained from a fully numerical (FN) solution of the Bethe-Ansatz integral equations. For comparison purposes, the earlier LSOC parametrization (2) is also included. To distinguish the present from the LSOC expression, we label the curves corresponding to the former by FVC. As the insets show, the relative deviation of FVC data from FN data is typically less than 2%, and at most ∼ 4%. This is the same size of error of the local-density approximation itself (quantified by comparing BALDA/FN data for small Hubbard chains to results from exact diagonalization), so that to within the accuracy of the LDA the present parametrization is a faithful representation of the full BA solution for the entire parameter range, including values of U ≫ 6t that cannot be realized in solids but occur in systems of trapped cold atoms.
Our main motivation for developing Eq. (4), however, was the incorrectly negative Mott gap predicted by the LSOC expression between U = 0 and U = 2t. The Mott gap E gap can be evaluated analytically from the expression for e 0 (n, U), either by explicitly calculating the derivative discontinuity or by taking total-energy differences [9]. From the LSOC parametrization one obtains [9] E LSOC gap (U) = U + 4 cos π β(U) ,  whereas the present parametrization leads to As illustrated in Fig. 2, the additional terms arising from the present expression push the gap upward, avoiding that it becomes negative between U = 0 and U = 2t as was, incorrectly, the case within the LSOC expression [9]. For U > 2t both the LSOC and the present (FVC) gap are positive, the latter being significantly closer to the numerical BA results than the former, although neither reproduces the subtle nonperturbative behavior for U → 0, in spite of the logarithmic term in Eq. (6). Only within the latter, however, the Mott gap is everywhere positive.

Extension to spin-dependent phenomena
In spin-polarized situations, the energy depends on the spin density m = n ↑ − n ↓ , in addition to the charge density n = n ↑ + n ↓ and the interaction U. This dependence is not included in the LSOC parametrization, which can therefore not be applied to study spin-dependent phenomena for electrons in solids or hyperfinepolarization-dependent phenomena for atoms in optical lattices. We note that BALSDA calculations for the 1DHM have already been performed in the context of cold atoms in optical lattices [19][20][21]. In lieu of an analytical parametrization, these works resorted to a fully numerical (FN) solution of the BA integral equations. An analytical parametrization would substantially simplify this approach, making the BALSDA as easily implementable as the BALDA one. An additional advantage of the analytical approach over the numerical one is that in the solution of the BA integral equations one cannot specify from the outset the density and magnetization of the system one is interested in. Rather, one has to specify the upper and lower limits of the integrals, and obtains the densities as part of the solution. This is clearly inconvenient for DFT, where the energies have to be evaluated as functions of the densities. More generally, it is desirable to be able to specify the system under study in terms of physical observables, such as the densities n and m, instead of in terms of auxiliary quantities, as the limits of the BA integrals.
Analytical expressions also permit one to derive further analytical results for other quantities. Our present derivation of closed expressions for the Mott gap is one example, and the analytical derivation and solution of Euler equations determining the phase diagram of harmonically trapped fermions on optical lattices [17] is another. For all these reasons, a simple but reliable parametrization of e 0 (n, m, U) can be useful for various types of calculations, within SDFT and beyond. Therefore, we next present an analytical parametrization of e 0 (n, m, U) and use it to construct a Bethe-Ansatz local-spin-density approximation (BALSDA) for the 1DHM.
In order to generalize the LSOC and FVC expressions to spin-dependent situations, we one more time follow the basic philosophy of constructing an analytical interpolation function recovering exactly known limits as function of the physically relevant variables. Several such exact results for e 0 (n, m, U) are known [22][23][24]. For non-interacting systems (U = 0), For infinite interaction (U → ∞), For half-filled unpolarized systems (n = 1, m = 0), Finally, for maximum magnetization (m = n), e 0 (n, m = n, U) = − 2 π sin(πn).
The LSOC and FVC parametrizations take m = 0, and express e(n, U) as a controlled analytical interpolation between (8), (9) and the m = 0 limit of (7). The construction of a more complete interpolation, recovering all four limits as functions of n, m and U, is strongly constrained by these limits, but still not unique. We therefore impose five additional common-sense criteria: (i) avoid high-order polynomials, which can produce unphysical wiggles, (ii) avoid unusual special functions, and (iii) keep the form similar to the LSOC parametrization. This third condition is useful because by now the LSOC parametrization has been implemented and used by many groups [9][10][11][12][13][14], so it will be easier to update to a new parametrization which has a similar form to the old one. However, as we have argued above, the LSOC expression does not properly describe the density-dependence of the ground-state energy at intermediate U. Therefore, we also require, as condition (iv), our spin-dependent generalization to reduce to the present Eq. (4) for m = 0 instead of to Eq. (2). The final, fifth, additional condition is motivated by computational efficiency and is explained in Sec. 4 below. Even with these additional common-sense criteria, the form of the parametrization is not uniquely determined, and a very large function space can still be explored. The particular choice made below was obtained by starting from the LSOC form and then building in, one-by-one, the additional exact limits and criteria. Many different variations have been explored, but we stopped when arriving at one whose deviation from the fully numerical solution of the Bethe-Ansatz equations was less than the typical error bar of the local-density approximation for this type of system. At this point, further improvements in the form of the parametrization become indistinguishable, in applications to inhomogeneous systems, from the intrinsic error of the LDA.
All four exact limits and five supplementary conditions are incorporated by the choice where β(n, m, U) = β(U) α(n,m,U ) , and α(n, m, U) = n 2 − m 2 n 15/8 Here, β(U) is the same quantity employed in the LSOC parametrization. For zero magnetization, α(n, m, U) reduces to α(n, U) used in Eq. (4). Equation (11) is valid for U ≥ 0 and n ≤ 1, but it can be extended to n > 1 and to U < 0 by standard particle-hole transformations [9,17,22,24].   (10) is of particular interest because this condition has no counterpart in the spin-independent situation. Moreover, it establishes a coupling between the charge and the spin dependence. Physically, maximum spin means that the Pauli principle keeps all fermions maximally apart for any U, just as infinitely repulsive interactions (U → ∞) do for any degree of spin-polarization. For n = m the spin-dependent parametrization must thus reduce to the same limit as for U → ∞ and recover the value β = 1. On the other hand, for m = 0 the expression should reduce to the earlier form (4). This double requirement is the explanation of the particular form chosen for α(n, m, U) of Eq. (14) and, consequently, for α(n, U) of Eq. (4). Figure 3 presents a comparison of the spin-dependence of this parametrization with data obtained from a fully numerical solution of the BA integral equations for intermediate parameter values, where the expression is not exact already by construction. Clearly, the spin-dependence of the homogeneous system is recovered to within the same precision as the charge dependence.

Local approximation for inhomogeneous systems
The main use of DFT and SDFT is in calculations for spatially inhomogeneous systems, where translational symmetry is broken and the density becomes position dependent. In the case of lattice models, inhomogeneity means that not all sites are equivalent. By means of the L(S)DA prescription, expressions for the energy of the homogeneous system, such as those described in the preceding section, can be used on a site-by-site basis to approximate the corresponding energy of the inhomogeneous system.
Explicitly, the local-spin-density approximation to any energy component E is given by where L is the number of lattice sites, and e = lim L→∞ E(L)/L is the per-site energy of the homogeneous system. This expression approximates the energy of the inhomogeneous system, where n i and m i vary from site to site, by evaluating the energy density of the homogeneous system site by site at the densities of the inhomogeneous one.
In DFT, including model DFT, this prescription is usually applied to the correlation energy, which for the Hubbard model can be defined as E c = E 0 − E M F , where E M F is the mean-field approximation to the ground-state energy E 0 . Since E M F is simple to obtain, the task to approximate the correlation energy E c [n i , m i , U] of the inhomogeneous Hubbard model, is thus reduced to that of approximating the per-site ground-state energy of the homogeneous one, e 0 (n, m, U). This is the quantity that we extracted above from the Bethe-Ansatz equations. The minimization of the resulting energy functional is conveniently carried out via selfconsistent Kohn-Sham (KS) calculations of the charge and spin densities and of the resulting ground-state energies. In such KS calculations the correlation potentials (obtained by differentiating the correlation energy with respect to n and m, or, equivalently, to n ↑ and n ↓ ) are evaluated once in each of the iterations of the selfconsistency cycle. In these iterations, the densities n i and m i change, but the interaction U, being a parameter of the Hamiltonian, does not. We have therefore expressed the integral in Eq.(3) and the resulting transcendental equation for β exclusively in terms of U, in order to guarantee that the only slightly timeconsuming steps of the calculation take place only once in each calculation, outside the selfconsistency cycle. This is the fifth additional condition on the spin-dependent parametrization, alluded to above.
As a simple example of such KS calculations, which serves to illustrate all essential aspects, we consider open Hubbard chains, where the spatial inhomogeneity stems from the boundaries, which give rise to charge and spin-density oscillations in the bulk. Representative results for a chain with L = 100 sites are displayed in Fig. 4, where we compare the ground-state density profile ( Fig.4-a) and spin-density profile ( Fig.4b) obtained from DMRG to LSDA profiles obtained from using the fully numerical solution to the BA integral equations and from the present parametrization. The BALSDA/FVC and BALSDA/FN ground-state energies of the same system deviate from the DMRG energy by 0.01% and 0.64%, respectively. The local densities follow the same trend, and deviate from DMRG by 0.42% and 0.58%. For the local magnetization, the corresponding numbers are 2.20% for BALSDA/FVC and 5.86% for BALSDA/FN. Remarkably, for all three quantities the parametrized results are closer to the DMRG benchmark data than the numerically defined BALDA. This shows that the particular form chosen for the proposed parametrizations allows for considerable error cancellation. On 32 processors the DMRG calculation took approximately 17 hours, while the BALSDA calculations required approximately 40 seconds. Of course, for high-precison calculations, as well as for the calculation of quantities that are not easily extracted from densities and energies, DMRG is still essential.
A more complex case is depicted in Fig. 5, which shows density and magnetization profiles for parabolically confined systems in a periodic chain. For two different values (k = 0.05 and k = 0.5) of the curvature of the confining potential (whose form is schematically indicated by the dashed (green) curve), the data points show the siteresolved particle density and spin density obtained by exact (Lanczcos) diagonalization, fully numerical BA-LSDA and our presently proposed parametrization. To within the accuracy that can be expected from a local-density approximation for this type of system (a few percent) the agreement between all three sets of calculations is excelent, for both charge and spin distributions. The amplitude of the magnetization-density oscillations is overestimated by both flavours of local approximations, which is consistent with previous observations for similar approximations and systems [20].
Next, we compare, in Fig. 6, the numerically exact ground-state energy obtained from Lanczos diagonalization of a small open Hubbard chain with L = 15 sites to the LSDA ground-state energies obtained from using the fully numerical solution to the BA integral equations and from using the present parametrization. Up to U ∼ 4t the   BALSDA/FN data are almost identical to the exact data, which attests to the quality of the local approximation. For U larger than ∼ 5t, the present parametrization is, once again, better than the conceptually superior fully numerical LSDA, due to error cancellation. The inset shows that for larger systems the behaviour is qualitatively the same.
Finally, we point out that successful applications of our (then unpublished) spin-Energy parametrizations for the 1D Hubbard model: gaps, spin and inhomogeneity 12 dependent expression (11) to the study of spin-polarized transport across a correlated nanoconstriction [11] and to the calculation of occupation probabilities of exotic superfluids in spin-imbalanced systems [25] have already been reported.

Summary
In summary, we have constructed a simple and reliable parametrization for the groundstate energy of the homogeneous one-dimensional Hubbard model, for arbitrary fillings, spin-polarizations and interactions. For the first time, a qualitatively and quantitatively correct description of the Mott gap is obtained from a simple density functional with a proper explicit derivative discontinuity. This parametrization can be used in its own right, for the homogeneous model, whenever simple expressions for the ground-state energy and for the resulting Mott gap are required. However, its main application is as input for local-density and local-spin-density approximations, which allow one to efficiently minimize the energy and extract energies, density profiles and related quantities for spatially inhomogeneous models. Since in KS calculations one never diagonalizes the interacting Hamiltonian, but only the auxiliary noninteracting one, systems with thousands of sites can be dealt with, even in the absence of any simplifying symmetry.