Uniqueness of density-to-potential mapping for fermionic lattice systems

Mapping many-body systems to a lattice while neglecting interactions beyond a given order has proven a powerful aid to devise tractable models —such as the Hubbard [1] and Heisenberg (see, e.g., ref. [2]) models —and further the understanding of their properties. Lattice systems are of recent interest for quantum chemistry calculations when the one-particle wave functions can be appropriately mapped to sites in a one-dimensional lattice and DMRG techniques are used to model the system (see, e.g., refs. [3–5]). In addition, lattice systems, both bosonic and fermionic, may be simulated experimentally using an optical lattice (see, e.g., refs. [6–8]). The latter is of great interest to the quantum technology community for developing quantum simulators. Fermionic lattice systems have been used to model the entanglement of nanostructures [9] and, using time-dependent density-functional theory, have been employed to investigate the time evolution of the out-of-equilibrium Mott insulator [10]. Time-dependent density-functional theory has also been used to create a method to simulate driven lattice gases with interactions [11] and used with other classical lattice Hamiltonians to model particle transport against a bias [12]. Lattice Hamiltonians have been used to describe the process of excitation-energy transfer and the role of intramolecular vibrations in the photosynthetic process [13]. Lattice systems are also used to simulate molecular systems, as for example the Pariser-Parr-Pople (PPP) model [14,15], similar to the inhomogeneous Hubbard model, can describe molecular systems with delocalised electrons in π-orbitals and has recently been used to model armchair polyacenes in ref. [16].

In the last decades density-functional theory (DFT) [17] in all its flavors has been a very important tool to predict properties of materials and structures. DFT is in principle an exact theory and allows to map many-body problems into non-interacting systems subject to a different external potential which accounts for the effects of the many-body interactions and inherits somewhat their complexity. The important concept of DFT on a lattice (L-DFT) was introduced by Gunnarsson and Schönhammer to investigate fundamental problems in DFT [18–20]. To improve efficiency of modelling and prediction power it is now possible to use approximations in DFT to efficiently model lattice systems (see, e.g., ref. [21] and references therein).

However there are still fundamental unsolved issues related to the theory itself. The power of DFT rests on the Hohenberg-Kohn [17] and Hohenberg-Kohn–like theorems [22–24] where the uniqueness of the mapping between a set of basic physics functions (e.g., the ground-state wave functions, their particle densities, and external potentials) is proven. Crucially, with respect to Hohenberg-Kohn–like theorems for DFT on a lattice [25–27], it has never so far been proven —or disproven— that a unique mapping between ground-state particle density (or wave function) and the external potential actually exist, leaving the basis of the theory incomplete.

In this article we prove that for finite fermionic lattice systems, the ground-state density uniquely determines the external potential, with the exception of the sites corresponding to nodes of the wave function and the extreme case of fully filled levels (double occupancy at all sites). We relate this exception to accidental and systematic non-uniqueness according to the classification by Capelle and Vignale [28]. To our knowledge this fundamental finding has not been demonstrated for a lattice system before, and it completes, for this class of systems, the one-to-one correspondence between ground states, their particle densities, and the external potential at the base of the L-DFT Hohenberg-Kohn theorem. We note that our demonstration also shows that the mapping from wave function to potential is unique not just for the ground state, but also for excited states.

After the formal demonstration, we develop a practical inversion scheme to determine the external potential from a given density for a lattice system and apply it to the one-dimensional finite Hubbard model with N particles to illustrate our findings. Our scheme applies to any fermionic lattice system.

The L-DFT HK theorem on a lattice [25] proves that when dealing with non-degenerate ground states the particle density uniquely determines the wave function on a lattice. Starting from this we will demonstrate that the one-to-one correspondence between the external potential and related ground state indeed holds true when considering a system on a lattice. We will show that given a finite lattice of L sites with particle number 0 < N < 2LM where M is the number of bands and where, for each energy level, the allowed site states are empty, singly or doubly occupied, and the interaction and external potential are not infinite, then the external potential is uniquely determined by the wave function.

We first consider ground-state nodeless wave functions in a single-band lattice, and extend afterwards to the case in which nodes are present, as well as to multiband lattices. Let us assume that there are two different potentials $\sum_{j}v_j$ and $\sum_{j}v'_j$ that give the same wave function $|\Psi\rangle$. Then

$$\begin{equation} \left(\hat{H}_{0}+\sum_{j}v_j\hat{n}_{j}\right)|\Psi\rangle=E|\Psi\rangle \end{equation} \tag{ 1 }$$

and

$$\begin{equation} \left(\hat{H}_{0}+\sum_{j}v'_j\hat{n}_{j}\right)|\Psi\rangle=E'|\Psi\rangle. \end{equation} \tag{ 2 }$$

Subtracting eq. (2) from eq. (1) gives

$$\begin{equation} \sum_{j}(v_j-v'_j)\hat{n}_{j}|\Psi\rangle=(E-E')|\Psi\rangle, \end{equation} \tag{ 3 }$$

which we may write as

$$\begin{equation} \sum_{j}\tilde{v}_j\hat{n}_{j}|\Psi\rangle=\tilde{E}|\Psi\rangle. \end{equation} \tag{ 4 }$$

We first consider site 1 and site 2. The occurrence of any two types of occupation is sufficient for the following proof to hold. We first consider a wave function that contains terms, with non-zero coefficient, no matter how small, of the form ${\left|{\uparrow \downarrow,S,X}\right\rangle}$ and ${\left|{S,\uparrow\downarrow,X}\right\rangle}$ where S is single occupation and X represents some configuration of the remaining particles into the remaining sites. This would be, for example the case for the Hubbard model ground state with the hypotheses of N > 2, no full filling, no complete spin polarization and not-infinite potential and interaction.

We then have

$$\begin{equation} \sum_{j}\tilde{v}_j\hat{n}_{j}{\left| {\uparrow\downarrow,S,X}\right\rangle}=(2\tilde{v} _{1}+\tilde{v}_{2}+k){\left| {\uparrow\downarrow,S,X}\right\rangle}, \end{equation} \tag{ 5 }$$

where k includes the terms due to X and

$$\begin{equation} \sum_{j}\tilde{v}_j\hat{n}_{j}{\left| {S,\uparrow\downarrow,X}\right\rangle}=(\tilde{v}_{1} +2\tilde{v}_{2}+k){\left|{S,\uparrow\downarrow,X}\right\rangle}. \end{equation} \tag{ 6 }$$

But $|\Psi\rangle$ is an eigenfunction of $\sum_{j}\tilde{v}_j\hat{n}_{j}$ so

$$\begin{equation} 2\tilde{v}_{1}+\tilde{v}_{2}+k=\tilde{v}_{1}+2\tilde{v}_{2}+k, \end{equation} \tag{ 7 }$$

which gives

$$\begin{equation} \tilde{v}_{1}=\tilde{v}_{2}. \end{equation} \tag{ 8 }$$

A similar procedure for sites 2 and 3, then sites 3 and 4 etc. leads to all the $\tilde{v}_{j}$ being equal so

$$\begin{equation} \tilde{v}_{j}=c. \end{equation} \tag{ 9 }$$

This means that the v_j and $v'_{j}$ may differ at most by a constant.

A similar proof follows also by considering terms ${\left|{\uparrow\downarrow,0,X} \right\rangle}$ and ${\left|{0,\uparrow\downarrow,X}\right\rangle}$ or terms such as ${\left|{S,0,X}\right\rangle}$ and ${\left|{0,S,X}\right\rangle}$. These would apply, for example to the Hubbard model with two particles or one particle, respectively. If we have half-filling with either $U\rightarrow \infty$ or complete spin polarization then we only have single occupation and the potential is not unique as would be the case for the Heisenberg model¹ .

For a multiband system, for more than two particles, the proof may be extended by considering terms of the form ${\left|{(B_{a})\uparrow\downarrow,(B_{b})S,X} \right\rangle}$ and ${\left|{(B_{a})S,(B_{b})\uparrow\downarrow,X}\right\rangle}$ (and similarly for different types of occupations) where we have written the occupation of the highest occupied band for sites 1 and 2, and the B_k represent some occupation of the lower bands at a site consistent with the total number of particles. This leads to

$$\begin{equation} (2+n_{a})\tilde{v}_{1}+(1+n_{b})\tilde{v}_{2}+k=(1+n_{a})\tilde{v}_{1}+(2+n_{b}) \tilde{v}_{2}+k, \end{equation} \tag{ 10 }$$

where n_k is the density contribution due to B_k. Hence $\tilde{v}_{1}= \tilde{v}_{2}$, so, by iterating the procedure to the other sites, we can show that $\tilde{v}_{j}=c$ for all j. The proof for N = 1 and N = 2 follows in a similar way.

We wish now to extend the above to ground-state wave functions with nodes. For simplicity we will consider a single node at site 3: extension to multiple nodes is straightforward. Now the procedure is carried out for site 1 and site 2 then site 2 and site 4 etc. thereby bypassing the node and again resulting in v_j and $v'_{j}$ only differing by a constant except at the node sites where this proof does not give any constraint on the potential.

This concludes our demonstration that the system density indeed determines uniquely the potential almost everywhere, with the sole exception of the wave function node sites. As nodes need to be isolated sites on the lattice we can then claim that for all practical purposes (f.a.p.p.) the ground-state density uniquely determines the potential on the lattice.

We note that for the limiting case of full filling we also do not prove that there is a unique potential, but this is consistent as any potential or interaction can then be applied with no change to the site occupations.

We classify the exceptions to uniqueness due to the Pauli exclusion principle completely determining the occupation of all sites as "accidental non-uniqueness" according to ref. [28]. Here Capelle and Vignale use this to label non-uniqueness due to a peculiar feature of the wave function, such as complete spin polarization in the case of spin-density-functional theory. We class non-uniqueness due to nodes as both accidental and systematic, in that once a node is known to occur the potential may be varied, up to a point, at the node-site to systematically construct other potentials that give the same wave function.

The above demonstration can be straightforwardly extended to show that f.a.p.p. an excited-state wave function uniquely determines the potential by using that any wave function, including excited states, uniquely determines the corresponding density. Conversely, given the external potential (for fixed particle-particle interaction) the excited-state wave functions are uniquely determined, so a one-to-one correspondence between an excited-state wave function and the external potential is now demonstrated.

Hence with finite interaction and potential f.a.p.p. the ground-state density, ground-state or excited-state wave functions uniquely determine the external potential on a lattice of L sites with 0 < N < 2LM particles where M is the number of bands. To corroborate our findings we will now develop an inversion scheme with a simple recipe of taking the nodes into account. Our findings in particular mean that, when our scheme converges to a potential that reproduces a density, this is, at least for ground states and f.a.p.p., the unique potential. We will explicitly show that the inversion scheme achieves the same final potential from different initial trial potentials.

For a generic lattice Hamiltonian with an external potential v_i at each site we may write

$$\begin{equation} \hat{H}=\hat{H}_{0}+\sum_{j}v_j\hat{n}_{j}. \end{equation} \tag{ 11 }$$

Given a wave function $|\Psi\rangle$ for this system, we then have

$$\begin{equation} \langle\Psi| \hat{n}_{i}\hat{H} |\Psi\rangle=En_{i}, \end{equation} \tag{ 12 }$$

where we have used $n_{i}=\langle\Psi |\hat{n}_{i}|\Psi\rangle$. In terms of the constituents of the Hamiltonian this is

$$\begin{equation} \langle\Psi|\hat{n}_{i}\hat{H}_{0} |\Psi\rangle +v_{i}\langle\Psi |\hat{n}_{i}^{2} |\Psi\rangle +\sum_{j\neq i}v_{j} \langle\Psi |\hat{n}_{i}\hat{n}_{j}|\Psi\rangle =En_{i}. \end{equation} \tag{ 13 }$$

We then rearrange eq. (13) similarly to ref. [29] and use this to create an iterative scheme to find the potential which gives the density $n_{i}^{\textit{target}}$ from an initial trial potential $v_{i}^{(0)}$

$$\begin{eqnarray} v_{i}^{(k+1)}&=\frac{1}{\langle\Psi^{(k)} |\hat{n}_{i}^{2}|\Psi^{(k)}\rangle} \bigl(n_{i}^{\textit{target}}E^{(k)}-\langle\Psi^{(k)}|\hat{n}_{i}\hat{H}_{0} |\Psi^{(k)} \rangle \nonumber\\ &-\sum_{j\neq i}v_{j}^{(k)}\langle\Psi^{(k)} |\hat{n}_{i}\hat{n}_{j}| \Psi^{(k)} \rangle\bigr). \end{eqnarray} \tag{ 14 }$$

Using the identity

$$\begin{eqnarray} &\langle \Psi^{(k)} | \hat{n}_{i}\hat{H}_{0} |\Psi^{(k)}\rangle +\sum_{j\neq i}v_{j}^{(k)} \langle\Psi^{(k)}|\hat{n}_{i}\hat{n}_{j}|\Psi^{(k)}\rangle= \nonumber\\ &-\,v^{(k)}_{i}\langle\Psi^{(k)} |\hat{n}_{i}^{2}|\Psi^{(k)}\rangle+E^{(k)} n^{(k)}_{i}, \end{eqnarray} \tag{ 15 }$$

we may rewrite eq. (14) in the simpler form

$$\begin{equation} v_{i}^{(k+1)}=\frac{\left(n_{i}^{target}-n_{i}^{(k)}\right)E^{(k)}}{\langle\Psi^{(k)} |\hat{n}_{i}^{2}|\Psi^{(k)}\rangle} +v_{i}^{(k)}.\end{equation} \tag{ 16 }$$

Following arguments similar to ref. [29], we note that this is expected to converge for E < 0, and so define the general expression

$$\begin{equation} v_{i}^{(k+1)}=\frac{\left(n_{i}^{(k)}-n_{i}^{target}\right)|E^{(k)}|}{\langle\Psi^{(k)} |\hat{n}_{i}^{2}|\Psi^{(k)}\rangle} +v_{i}^{(k)}. \end{equation} \tag{ 17 }$$

At convergency $v_{i}^{(k+1)}=v_{i}^{(k)}$ is the external potential that reproduces the target density via the many-body Schrödinger equation. We note that this approach can be applied to any of the system eigenstates, including the ground state. Our result eq. (17) is similar to that derived for continuous spatial co-ordinates in ref. [29] except for the change in the denominator. We will numerically implement this scheme with 80% mixing of the previous potential to reduce the chance of instabilities.

We now apply our scheme to a finite lattice system with fixed and finite particle-particle interaction to find the unique potential, and hence the wave function, corresponding to a density. To explore the behavior of the scheme at the nodes, and for practical examples of its implementation, we consider the one-dimensional Hubbard model (HM)

$$\begin{eqnarray} H_{HM}&=-t\sum_{i,\sigma}\left(c_{i,\sigma}^{\dagger}c_{i+1,\sigma}+c_{i+1,\sigma} ^{\dagger}c_{i,\sigma} \right) \nonumber\\ &\quad +\,\tilde{U}\sum_{i}\hat{n}_{i,\uparrow}\hat{n}_{i,\downarrow}+\sum_{i}v_i\hat{n}_{i}. \end{eqnarray} \tag{ 18 }$$

We use open boundary conditions, L sites, interaction strength of $U=\tilde{U}/t$ and N particles. We first demonstrate that we can achieve the same final potential, up to an additive constant, independent of the starting potential. For the exact system we use a harmonic potential

$$\begin{equation} v_{i}=\frac{1}{2}k\left(i-\frac{(L+1)}{2}\right)^{2} \end{equation} \tag{ 19 }$$

with k = 1 to generate the target density. We then apply the iterative scheme (17) starting from a harmonic potential with k = 2 and then starting from $k=0.1$. Here we require a mean absolute difference of $10^{-8}$ between the trial and target site densities for convergence.

We see in the top panel of fig. 1 that the same final potential is found, up to an additive constant, regardless of the starting potential. Figure 1 (lower panel) shows that the rate of convergence, however, may be affected by the trial potential.

Zoom In Zoom Out Reset image size

To investigate the use of the scheme when nodes are present we consider excited states to produce nodes with a finite potential. Our argument does not show what values the potential takes at a site with a node, so in practice, we set the potential at these sites to zero at all iterations of eq. (17). We then attempt to use the scheme to find a finite potential, if any, that gives the required density. In the calculations we prevent the possibility of the potential becoming too large, and possibly tending towards an infinite potential solution, by creating a random potential at all sites with maximum value of 10 if any of the sites has a potential greater than 1000.

We first consider the case of a target density with a node at a single site (see fig. 2), where the potential is not known beforehand, and consider U = 0. Here our scheme indeed converges towards a potential such that the third excited state of the system reproduces the target density with a node at site 3, as demonstrated in fig. 2. We note that in this case the final potential at site 3 is not unique but nor is it arbitrary: by direct solution of the Schrödinger equation, we find that the potential can be lowered at site 3 to around −0.2 before the density changes.

Zoom In Zoom Out Reset image size

The approach can also be applied to a density with multiple nodes. To find a density of this form that could be created from a known finite potential we considered the exact solution for one particle with a constant potential. This system with seven sites has three nodes in its third excited state. Figure 3 demonstrates that we can indeed recreate this density by finding a constant potential, albeit a different constant, when using the scheme with the k = 1 harmonic potential as the initial guess.

Zoom In Zoom Out Reset image size

We proved that for a single or multiband fermionic lattice system without full filling nor fully polarised at half-filling the external potential is uniquely determined by the ground-state wave function except at sites with zero-density nodes and hence by the ground-state density due to the unique mapping between ground-state densities and wave functions [25]. We completed in this way the demonstration of the Hohenberg-Kohn–like theorem for Lattice-DFT. We also demonstrated a one-to-one correspondence between external potential and excited states. As a practical application, we have adapted the inversion scheme of ref. [29] for a generic lattice system and applied it to the Hubbard model. We have shown that with this scheme the external potential corresponding to a lattice density, and given the other constituents of the Hamiltonian, may be found. Our results will prove useful for further fundamental work in density-functional theory and also in the testing, construction and improvement of approximations in lattice density-functional theory with potential applications in all fields in which lattice Hamiltonians are used. Future work will consider extending these results to finite-temperature density-functional theory [30,31] for lattice systems [32].

IDA gratefully acknowledges support from a University of York-FAPESP combined grant. VVF acknowledges support from FAPESP and CNPq. We are thankful to Klaus Capelle for useful discussions.