Testing density-functional approximations on a lattice and the applicability of the related Hohenberg-Kohn-like theorem

We present a metric-space approach to quantify the performance of approximations in lattice density-functional theory for interacting many-body systems and to explore the regimes where the Hohenberg-Kohn-type theorem on fermionic lattices is applicable. This theorem demonstrates the existence of one-to-one mappings between particle densities, wave functions and external potentials. We then focus on these quantities, and quantify how far apart in metric space the approximated and exact ones are. We apply our method to the one-dimensional Hubbard model for different types of external potentials, and assess the regimes where it is applicable to one of the most used approximations in density-functional theory, the local density approximation (LDA). We find that the potential distance may have a very different behaviour from the density and wave function distances, in some cases even providing the wrong assessments of the LDA performance trends. We attribute this to the systems reaching behaviours which are borderline for the applicability of the one-to-one correspondence between density and external potential. On the contrary the wave function and density distances behave similarly and are always sensitive to system variations. Our metric-based method correctly predicts the regimes where the LDA performs fairly well and the regimes where it fails. This suggests that our method could be a practical tool for testing the efficiency of density-functional approximations.

give better results when using more accurate densities than those found self-consistently giving rise to the idea of distinguishing between functional-driven and density-driven errors 4 .
Lattice models, such as the Heisenberg or Hubbard models 5,6 have been extremely important for the understanding of strongly-correlated many-particle systems. Despite their simplicity, when employed appropriately, these models are able to capture phenomena sufficiently accurately: the Hubbard model has been shown to reproduce the Mott metal-insulator transition, and has been recently associated to the behaviour of 'exotic' systems such as inhomogeneous superfluidity in spin-imbalanced systems 7,8 , chains of Bose-Einstein condensates 9 , or entanglement in nanostructures 10,11 . However, when interactions and inhomogeneities are included, even these lattice models can rapidly become computationally intractable as the size of the systems increases. This is where the powerful concept of DFT on a lattice 12-14 becomes useful. By using approximations within the approach of lattice DFT, specifically the local density approximation (LDA) [15][16][17] , the constraint on the maximum system size can be considerably loosened, and larger lattice systems can be accurately modelled. DFT on a lattice is built on foundations [18][19][20] similar to the original Hohenberg-Kohn theorem for standard DFT. However these foundations were completed only recently 21 when the last part of the Hohenberg-Kohn-type theorem−the one-to-one mapping between ground-state density and external potential, and ground-state wave function and external potential−was demonstrated, and their limitations discussed.
In this context, we evaluate the use of specifically designed metrics for density, wave function and potential to appraise lattice-DFT approximations beyond total energy arguments, and use the same tools to explore the applicability of the Hohenberg-Kohn-type theorem on a lattice. Our results show that indeed these two subjects are intimately connected. The method we propose applies to any fermionic single-band lattice system and could be used on any DFT approximation. Here we will concentrate on the LDA. As practical−and exactly solvable− test-bed examples, we will demonstrate our method on the one-dimensional Hubbard model considering short to medium size chains of varying particle number, interaction strength, and applied external potentials.

Theoretical and Computational Methods
DFT demonstrates that, for time-independent systems with a spin-independent external potential, ground-state particle density, the corresponding wave function, or the external potential are sufficient to describe a quantum system. We wish then to rigorously appraise DFT approximations by calculating how well they reproduce the exact values of these three quantities, and we will do so by using metric spaces.
We note that, to this aim, we cannot simply compare exact and Kohn-Sham systems, as nor the wave function, neither the external potential of the KS system are constructed to reproduce the exact ones. We then follow the route proposed in 22,23 : we will consider, together with the exact, that unique interacting system constructed to have the same density as the non-interacting Kohn-Sham system obtained by using the DFT approximation, and with the same particle-particle interaction operator as the exact interacting system. We will construct it via reverse engineering, by using the iterative scheme in ref. 21 to find the potential of the interacting lattice system given the approximated density. At this point, the distance between all quantities of interest can be assessed by using appropriate metrics.
In this paper, to illustrate our method, we focus on one of the most used density-functional approximations, the LDA; we will introduce below all the elements of the method outlined above.
The interacting-LDA system. The "interacting-LDA" system, the i-LDA introduced in ref. 22 , is the uniquely defined many-particle interacting system with the same many-body interaction operator as the original exact Hamiltonian, but whose ground-state density is the LDA density which is found by solving the Kohn-Sham equations within the LDA for the original many-body system.
The advantage of using the i-LDA system is that its many-body interacting wave function tends to the exact many-body wave function for all interaction regimes for which the LDA is a good approximation. This is a fundamental difference with the corresponding KS-LDA wave function which is by definition non-interacting, and therefore it is able to reproduce the exact many-body wave function only in the non-interacting limit. In addition, because the many-body interaction operator is the same for the exact and the i-LDA Hamiltonians, then the better the LDA performance, the closer the corresponding two external potentials v and v i − LDA should be.
The properties of the exact, KS-LDA, and i-LDA systems are summarized in Table 1.
Iterative scheme on a lattice. Consider a generic lattice Hamiltonian with an external potential v j at site j, and so described by Here =ˆ † n c c j j j the site j occupation operator and ˆ † c j ĉ ( ) j the fermionic creation (annihilation) operators. The ground state wave function for this system is |Ψ〉, such that where E is the ground state energy and we have used = 〈Ψ| |Ψ〉 n n i i . We may then use the iterative scheme of ref. 21 to find the potential which gives the density n i target from an initial trial potential v i (0) . The recursive formula is () is the external potential that reproduces the target density via the many-body Schrödinger equation. This approach is general and not restricted to |Ψ〉 being a ground-state. We numerically implement this scheme with 80% mixing of the previous potential to reduce the chance of instabilities.
Wave function and density metrics. To describe the closeness between the wave functions (or the densities) of the exact and i-LDA systems we consider the rigorous 'natural' metrics for wave functions (D ψ ) and densities (D ρ ) discussed in ref. 24 . For completeness we report below their expression where we have followed the convention in ref. 24 and normalised the wave functions to the particle number N.
To facilitate straightforward comparisons we use scaled metrics ψ D and ρ D that reside on [0, 1]: Potential metrics. For a finite lattice system of d sites and finite potential we may define a distance between potentials similar to the density distance created from the density norm in ref. 24 where v 1, j (v 2, j ) is the external potential of system 1 (2) at site j. As the sum of the potential is not constrained, unlike the density, then we have divided by the number of sites d to allow fair comparison between different system sizes. A metric d(x,y) must adhere to the following three conditions: . We see that Eq. (8) is symmetric on exchanging v 1 with v 2 , is zero if and only if v 1 and v 2 are equal, and as the absolute value satisfies the triangle inequality then so does Hence the potentials on a lattice of d sites with the distance Eq. (8) give rise to a metric space. The metric (Eq. 8) could be applied to any couple of systems whose Hamiltonians satisfy (Eq. 1); for the scope of this paper we will use it to compare systems whose Hamiltonians differ by the external potential term only, and in particular the two external potentials are obtained by using two different methods/approximations to solve the same physical problem. We note that the metric (Eq. 8) would not be suitable for a continuous variable system since the potential for many systems is unbounded unless a cut-off is used.
Physical potentials are only defined up to an additive constant c. This needs to be taken into account to prevent the unwanted situation where the wave function and density distances are both zero, but the potential distance is not. So, similarly to the wave function distance 24 that was created to be gauge-independent in ref. 24 , we wish to remove the arbitrarity of this constant for the potential. We define then and we find the ∈ R c that minimizes this sum using an iterative weighted least squares technique 25 , We consider Eq. (10) as a distance between classes of potentials such that each class contains all potentials which are equal up to constant. Eq. (10) is symmetric on exchanging 1 and v 2 belong to the same class. Then, to demonstrate that Eq. (10) is a metric between physically different potentials, it only remains to be shown that it satisfies the triangle inequality. Now for the c 1 and c 3 that minimize the sums we have obeys the triangle inequality then Combining the above equations leads to the triangle inequality We also consider a distance between potentials similar to the wave function distance defined in 24 using the sum of squares: is the Euclidean distance scaled by 1/(d) 1/2 and so is a metric for potentials on d sites. The same arguments as for D v A then can be used to show that D v B is also a metric. This time the minimization may be achieved analytically resulting in where μ is the mean value of Δv j = v 1, j − v 2, j . This leads to (19) which was used to quantify the match of two potential curves in ref. 26 . Here σ(Δv) is the standard deviation of Δv.
Both potential distances can be scaled in a standard way to a distance between 0 and 1 (see, e.g., 27 )

Results
We consider the one-dimensional Hubbard model (HM): where t is the hopping parameter, U the on-site interaction, and σ σ + † c c , i i , 1 , are creation and annihilation operators of fermionic particles with z-spin component σ = ±1/2 at site i. The metric-space analysis will be performed for three very distinct types of external potentials: homogeneous potential, harmonic potential and localized impurities. In all calculations we will set t = 1. For small chains (d ≤ 14 sites) we obtain exact data via Lanczos diagonalization with tolerance 10 −14 , while for larger chains we produce nearly-exact results using DMRG techniques whose parameters (finite-size algorithm, basis-size 80, truncation error 10 −5 ), for d = 10, produced deviations smaller than 0.001% from exact distances. Finally the inversion scheme Eq. (3) is performed with an average site error threshold of 10 −8 .
Metric-space analysis: homogeneous potential. We start by considering small Hubbard chains with no external potential, v i = 0 for all i, so the inhomogeneity comes from finite-size effects only. In the upper panel of Fig. 1 we show the distances between the i-LDA potential, wave function, and density with the corresponding First we observe that, for all distances, the chains with odd number of particles are systematically closer to the exact system than the ones with N even. The odd-even oscillation appears because for fixed n we have alternating non-magnetic (even chains, N ↑ = N ↓ = N/2) and magnetic systems (odd-N chains, N ↑ ≠ N ↓ ). For non-magnetic systems, the LDA will give the same results for both spin density components, so that any error in a spin component will be doubled-up when considering the total density. For magnetic systems, spin-up and spin-down densities are different and hence the (inhomogeneous) density oscillations will be different for each component. In this case, when the two components are summed up to yield the total density, cancellation of errors in the LDA results may occur. Indeed this is what systematically happens in this system, as illustrated in the intermediate and lower panels of Fig. 1.
We also find that, with the exception of ψ D , the distances for even (odd) number of particles behave qualitatively similarly: they decrease as the chain (and particle number) increases, which is consistent with the fact that the LDA becomes exact at the limit d → ∞. In particular, both ways of quantifying the potential distance capture the trend. As all plotted distances are scaled to a maximum of 1, we can deduce that, for homogeneous chains and all quantities analyzed, the LDA performance is closer than 10% to the exact results. However, when comparing the different quantities, we note that LDA results for densities are closer to the exact results, while external potentials are reproduced slightly less precisely by the LDA, no matter which of the potential metrics we use.
The wave function distance reproduces the odd/even N oscillation behaviour; however for chains with even (odd) number of particles the distance is increasing with number of sites. In addition the i-LDA wave function seems to perform, as the number of sites grows, increasingly worse than the other i-LDA quantities: it gives the overall closest distance for d = 4 but the farthest for d = 10; unfortunately the size of the Hilbert space becomes too big to perform the inversion scheme for larger values of d, so we cannot confirm the trend for longer chains. We attribute this worsening in performance to the scaling with N and d of the overall space: this only increases linearly with the sites for the other quantities, but the wave function (Hilbert) space exponentially grows from = ( )( )  they are normalized, and the numbers of spin-up and spin-down particles are fixed at N ↑ = N ↓ for even N and N ↑ = N ↓ + 1 for odd. This is achieved by assigning a pseudo-random number from [−1, 1] for each permissible Slater determinant then normalizing the resulting wave function. We see in Table 2 that picking a close wave function from this collection is much less likely than finding a wave function that gives a close density. This becomes more pronounced as the number of sites increases causing the size of the wave function configuration space to dramatically enlarge. The densities actually become closer to the exact on average as the sites increase, which we speculate as due to the relative homogeneous density of the exact system, while by 10 sites we find that a random wave function is almost always around the maximum distance from the exact.

Metric-space analysis: localized impurities.
Here the external potential is chosen to be a collection of localized repulsive impurities with the same strength V. We start by considering a chain of size d = 10 with open boundary conditions and with two impurities localized at the central sites. Figure 2 presents the distances for density, wave function and potential between the i-LDA and the exact system as a function of the impurities strength V. We see that the distances have a similar qualitative behavior: they are minimum at V = 0, show a peak for small/intermediate V's, and saturate for V≫ 1. The saturation occurs because the increasing repulsion V progressively devoids the impurity sites until their occupation becomes negligible, and then the density profile at the remaining sites remains unaffected by further increase of V. Interestingly the peak appears at significantly higher V value for the potential distances than for the wave function and density distances, which are in general quantitatively much more similar. At contrast with Fig. 1, here we are keeping N, and hence the size of the configuration space, fixed. The regime V ≫ 1 can be considered as equivalent to reducing it. Instead, what strongly varies is the potential strength. This is found to affect the potential distances differently from the other distances: for intermediate and large V's, the potential distances are substantially larger than the density and wave function distances.
In the regime V ≫ 1, impurities V mimic extra boundaries, and they make the LDA poorer, so one would expect this regime to correspond to the maximum distance; because of this the presence of a peak at small/intermediate V's is unexpected. In order to understand its appearance, in Fig. 3 we analyze the density distance for larger chains, d = 20 and d = 30, keeping the same impurity percentage and structure (centered impurities). Our results show that not only the peak persists, but also a minimum at V ≠ 0 appears (the presence of this minimum cannot be confirmed for d = 10 as both LDA and DMRG results do not converge accurately for 0 < V < 1). This  is surprising because the LDA is, by definition, exact at the spatially homogeneous limit: in our finite chains one would expect that V = 0 is the closest to this limit, as then only the boundaries induce inhomogeneities, so that V = 0 should correspond to the minimum distance from the exact system. We solve these conundra by analyzing the density profiles (Fig. 4) these show that both the dips and the peaks of Fig. 3 are generated by the finite-size of the chains considered. For small/intermediate V the peak reflects the fact that the impurities sites are still non-empty, so the corresponding strong density inhomogeneity contributes to make the LDA worse in comparison to the saturation regime where the impurities sites are empty. Instead the dip observed in Fig. 3 is related to the particular choice of locating the impurities symmetrically with respect to the boundaries: this symmetry somehow favors the LDA performance. When we simply displace the impurities to an asymmetric position, the dip disappears, as we can see in the upper panel of Fig. 5.
In general we find that the LDA's performance worsens for shorter chains and same impurity structure (Fig. 3), and for a symmetric but increasing spreading of the impurities with same chain length (lower panel of Fig. 5). It is easier to understand this behavior in the saturation regime V ≫ 1, where impurity sites are practically empty and then the chains get fragmented into increasingly smaller segments, bearing higher inhomogeneity. For intermediate values of V, when the impurity sites are just depleted but not empty, the LDA performance may even worsen, leading to the peaks visible in both Figs 3 and 5. This is because in this case the LDA has also to cope with simulating the highly inhomogeneous impurity density dips (see panels with V = 1 in Fig. 4).
Finally in Fig. 6 we consider the LDA density distances with respect to the impurity strength for distinct values of the interaction U and fixed chain length (d = 14). Impurities here are in smaller concentration with respect to Fig. 3 but still located symmetrically in the center of the chains, so maximizing homogeneous segments length with respect to the chain size. We find that the density distance properly describes the expected behavior that the LDA performs worse for higher interaction U. The results also show how sensitive the LDA is to finite-size effects: a symmetric but smaller concentration of impurities is here enough to destroy the dip-and-peak structure.
Metric-space analysis: harmonic potential. We now turn to consider systems confined by the harmonic potential kx 2 , centered between the two middle sites. We use open boundary conditions, d = 8 sites, interaction strength U = 2 and N = 2 particles. We analyze the LDA performance using the distance between the i-LDA and the exact systems as a function of k for potential, wave function and density.
Density and wave function distances are considered in the upper panel of Fig. 7: the LDA performance is very good for all k's, with less than 3% maximum error for k → 0. Also, as for the localized impurity case, density and wave function distances behave very similarly: they both present a local minimum at k ≈ 0.3, a local maximum around k = 1, and tend to zero as the confining potential increases further. We attribute this minimum-maximum structure to a competition between depletion of density in certain sites and effective reduction of the chain length, similarly to the mechanism we have described for Fig. 3. In the present case, as k starts to increase, the parabolic potential starts to deplete the chain ends' sites, excluding more sites as the potential strengthens. For large k values, all central sites become doubly occupied, the system freezes, and the distances between densities and between wave functions tend to zero. We see in the inset of the upper panel of Fig. 7 that the density can be very low away from the central sites for k = 8. In this case we find that the percentage error of our numerically found i-LDA density with the LDA density is 15% on the first site, but this is for a density around 10 −10 which is 100 times lower than our convergence check. On the second site the error is only 0.7% despite the density being around 10 −6 demonstrating that our inversion procedure is sufficiently robust. The intermediate panel of Fig. 7 displays the percentage error of the LDA relative energy with respect to the exact values. Also in this case the error is small (below 2%) and the curve presents a broad maximum feature, qualitatively similar to the one seen in the upper panel of Fig. 7. However the energy, an integrated function of the density, is not sensitive to local effects such as the competition between depletion and effective chain reduction, hence missing the minimum/maximum structure characterizing both density and wavefunction metrics.
Conversely, both potential distances display a behavior completely different from wave functions and densities (lower panel of Fig. 7): they are monotonically increasing with k after k = 0.8 with the minimum of D is at k ≈ 0.4. Notably the LDA performance, when measured in this way, worsens as k increases, up to more than 70% error. Also both potential distances do not display any obvious sensitivity to the freezing of the system: indeed the potential can arbitrarily increase with increasing k, even after saturation of the central sites (freezing) has occurred. In ref. 21 we discussed how, for DFT on a lattice, in the limiting cases of full-filling and/or potential tending to infinity, the one-to-one correspondence between density and external potential does not hold. The large-k case at hand comes close to these limiting behaviors, and this explains why the trend of the potential distances can become strikingly different − and even opposite, like in this case − with respect to the behavior of densities and wavefuctions distances.

Conclusion
Within the metric approach to quantum mechanics 24,28 , we used wave function, density and external potential distances to explore LDA performances for short and medium inhomogeneous one-dimensional Hubbard chains.
The analysis has been carried out by comparing the distance between exact many-body systems and their corresponding interacting-LDA systems, that is the interacting systems built to have the LDA ground-state densities. To derive the Hamiltonian of the interacting-LDA systems, we have applied the inversion scheme of ref. 21 to the LDA density of the one-dimensional Hubbard model, and obtained, to a high degree of accuracy, the potential of the corresponding interacting system.
We considered open boundary conditions and three distinct potentials: homogeneous, localized impurities, and harmonic confinement, for different chains' length, particle numbers and interaction strength.
The homogeneous and the localized-impurities analyses revealed several features about the LDA performance across different regimes of parameters. In particular we found interesting finite-size effects, with some counterintuitive behaviour by the LDA. Remarkably, our distances for wave functions and particle densities were sensitive to all of these effects: they correctly pointed out the regimes where the LDA performs fairly well and the regimes where it fails, as confirmed by the direct analysis of the systems' densities. Thus these metrics were proved to be useful for testing the pitfalls of the LDA and therefore could be a practical tool for testing the efficiency of any other density-functional approximation. In addition our work provided evidence for lattice systems that when the LDA density is close to its exact counterpart then so too is the interacting LDA wave function and vice-versa. Although this relationship was affected by the size of  the wave function's configuration space it was not strongly altered. We also found that the density distance is in general smaller than the others.
One important result from this work is that care must be taken when considering the distance between external potentials in appraising the performance of a density-functional approximation. At least for systems on a finite lattice, the fact that the one-to-one correspondence between density and potential fails in certain limits 21 ) between i-LDA and exact systems, versus confining potential strength k. In all panels: U = 2, d = 8 and N = 2, with confining potential kx 2 symmetric with respect to chain centre.
SCIentIFIC REPORTs | (2018) 8:664 | DOI:10.1038/s41598-017-19018-x seems to imply that the potential distance between exact and approximated systems may be misleading if used alone, even showing a behavior which is qualitatively opposite to the density and wave function distances for systems close to these limits. In these parameter regions, the potential distance suggests that the LDA behaves increasing poorly, failing to recognize that instead its performance is increasingly improving.