A v0-representability issue in lattice ensemble-DFT and its signature in lattice TDDFT

We study a small Anderson-impurity cluster using lattice density functional methods, and try to determine the exact exchange-correlation (XC) potential via reverse engineering. In doing so we find singlet–triplet degenerate interacting ground states which cannot be v0-represented in an ensemble-DFT sense. We also find that it is possible to represent a triplet ground state as a pure state, but not the singlet. We further investigate this behavior within time-dependent density-functional theory. Starting from a v0-representable ground state and entering the degeneracy region via a time-dependent perturbation, we determine via reverse-engineering the time-dependent XC potential for progressively slower perturbations. This analysis shows that, even in a constant external field, the XC potential must retain a time-oscillating pattern to ensure a density constant in time, a hint of the underlying representability issue in the ground state.


Introduction
Density functional theory (DFT) [1][2][3][4][5], and time dependent density functional theory (TDDFT) [6][7][8][9] are rather distinct theories: while both electing the density as basic variable, they address different problems and rest on structurally different proofs [2,6,[10][11][12][13][14], since one (DFT) but not the other (TDDFT) disposes of a minimum principle. Also when it comes to the Kohn-Sham (KS) potential v KS , which enables the mapping of an interacting system into a noninteracting one, the two theories face different issues in matters of v 0 -representability [12]. This point is easily appreciated if one considers that the time degree of freedom in TDDFT introduces a new level of complexity in v xc : a dependence on the initial state and on the history of the density (for a recent discussion, see Refs. [15,16]).
Motivation and aims of this work -Here, we wish to superficially touch upon some of these aspects, by studying a finite lattice model system, engineered to be non-v 0 -representable in the ground state, and examining how a representability issue induces signatures in v xc during the system's response to increasing slower perturbations. We will use DFT and TDDFT for lattice systems [17].
A brief survey of literature on lattice (TD)DFT -Model systems can contribute to clarify conceptual aspects of a theory, to address qualitative features of a physical phenomenon, or, when exactly solvable, to numerically benchmark approximate methods. This is also the case for DFT and TDDFT, which have also been applied to simple lattice model Hamiltonians [18]. The general idea of lattice DFT was introduced about 30 years ago [19] and then a local-density approximation LDA targeted for the 1D Hubbard [20] model was considered [21]. However it was only a few years later [22], that a practical use of an LDA was advocated for investigating inhomogeneous Hubbard models, and a parameterized analytical form of the exchange-correlation (XC) potential was proposed [23]. Since then, several works have appeared in lattice DFT, dealing with the proposal of new XC potentials [24][25][26][27], temperature effects [28], attractive interactions [29], and the quality of approximate KS potentials in metric spaces [30]. A recent review of lattice DFT, also about applications, can be found in reference [18].
Lattice TDDFT in the linear response regime was initially applied to the Hubbard dimer [31,32], but later also to the 3D homogeneous Hubbard model [33]. The Hubbard dimer has also been the object of a comprehensive study to gain insight into fundamental aspect of DFT [34]. For very recent work on linear response theory and the Hubbard dimer, we refer to reference [35].
Using TDDFT for real-time dynamics was proposed a decade ago [36][37][38][39], together with an adiabatic approximation for the 1D Hubbard model and the extension to spin-dependent systems. Already at this initial stage, part of the focus shifted towards basic issues such v 0 representability and the need for a proof of the Runge-Gross (RG) theorem on the lattice. An RG proof (a) Representation based on the site-orbitals in the square: in all cases, U = 1, t = 1, c = 0. This leaves with a and b as the only free parameters in the equilibrium Hamiltonian. (b) Representation based on the symmetry-adapted orbitals: as detailed in the main text, the system separates in a 3-orbital system (3-cluster) with interactions and a group of noninteracting clusters (6-clusters). (c) Same as in (a), but with the symmetry-breaking perturbation terms explicitly shown. In this case, the cluster is not separable as in (b).

The model
We study a 3 × 3-site Anderson cluster, with Hamiltonian where ij are nearest-neighbour sites, a † iσ (a iσ ) creates (destroys) an electron of spin projection σ at site i, n iσ = a † iσ a iσ , and t is the hopping parameter (see Fig. 1). A local (Hubbard-like) interaction of strength U between electrons with opposite spins is present at site 1. The onsite energies i are chosen fulfilling the natural C 4v symmetry of the lattice: i.e., according to Figure 1, 1 = a , 2−5 = b , 6−9 = c . Apart from a trivial overall energy shift, this choice leaves a two-dimensional parameter space for the on-site energies. Finally,V (τ ) represents a timedependent (τ labels time), one-body perturbation applied to the cluster, left unspecified for the moment. For later considerations, we denote by H 0 the static one-body part of the Hamiltonian, i.e. H 0 = iσ iniσ + t i,j ,σ a † iσ a jσ . Spatial symmetry permits to partition the 3×3 cluster in Figure 1a into separate sub-systems. To this end, one introduces the symmetry-adapted basis In terms of these states, where the original symmetryrelated orbitals are mixed (for example, the orbitals at the four sites with energy b in Fig. 1a), k|H 0 |k assumes a block-diagonal form, corresponding to one 3-orbital cluster and four non-interacting smaller clusters (Fig. 1b). The change of basis does not affect the diagonal orbitals energies, while transition amplitudes (the "hopping" parameters) are renormalized, as specified in Figure 1b. The block-diagonal structure persists for the matrix elements of H 0 among Slater determinants formed from the symmetry-adapted states. For the onsite interaction term in equation (1), we note that it only concerns the central site of the original cluster (which in the symmetry-adapted basis maps onto itself). Thus the interaction term only appears in the total-symmetric block (IRREP A 1 ), i.e. in the 3-cluster [26,73,74].
Since the energetics of each subpart can be in many ways established separately, this offers a suitable premise for the insurgence of issues about uniqueness The solid-black/dashed-grey arrows in the outermost occupied levels of the many-body-system schematically indicate that the singlet-triplet degeneracy is realized by mixing spin-up and down projections of single electrons in the 3-and 6-cluster, as shown in equation (3). KS panel: when the middle level of the 3-cluster and the third lowest level in the 6-cluster are aligned, this corresponds to the case a = 2 b , used to explore ensemble v0-representability of the many-body singlet-triplet degeneracy. When such two levels are not aligned, this provides a KS image of the triplet state. Right panel: reverse engineering of the KS potential along the line a = 2 b to determine an ensemble KS-image of the may body degenerate ground state. The vertical axis reports the absolute value of the maximum discrepancy between the many-body and the KS densities at the nine cluster sites. In the inset, a close-up of the region where the discrepancy reaches its minimum is provided.
of the potentials and KS v 0 -representability [41,66] in lattice (TD)DFT. In what follows, we consider states with four spin-up and four spin-down electrons.
A remark about reverse engineering -We will employ both (i) ground state and (ii) time-dependent reverse engineering schemes. In either case, the KS potentials are found by requiring that the KS density matches the exact one. However, for (i), the possible presence of a degenerate ground state poses an additional constraint. In the static KS picture, having or not the requirement of degeneracy corresponds to simply add or remove a constraint on the effect of the KS potential on the topmost occupied KS eigenvalues (namely, together with the usual condition of matching the density, the KS potential must induce or not induce an alignment of such topmost occupied levels). This second constraint plays no role for the timedependent analysis: for (ii), a standard reverse engineering procedure was used, as e.g. done in reference [39]. ensemble DFT. The engineered ground state degeneracy was obtained with the Hamiltonian parameters set to the values t = 1, U = 1, a = −2.69, b = −1 (in the following, c = 0 always, setting the scale for the onsite energies). In this way, in the ground level there are three electrons in the 3-site cluster and five electrons in the remaining 6-site system, thus realising a single-triplet degeneracy.
where A 3 σ and B 5 σ are the ground states of the 3-and 6-clusters with N A = 3 and N B = 5, and an unpaired spin of projection σ, respectively (see Fig. 2a, left panel).
Preliminary to the discussion of v 0 -representability, it is useful to consider the general energy level structure of the non-interacting system, when varying the on-site energies. We will try to represent a ground state with a C 4v -symmetric density, thus we can restrict our search to potentials with the same symmetry and, without loss of generality, scan a two-parameter space ( a , b , 0). As suggested by Figure 1b In the relevant range max(| a || b )| < 10, the 3-cluster always contributes the highest and the lowest level of the system. This result puts restrictions on the possible orders of the levels, and thus the only way to obtain a degenerate ground state is when the middle state in the 3-site cluster matches 0 or b (depending on the sign of the latter), with 3 electrons in the 3-cluster. This can be gathered looking at the KS panel of Figure 2, since, qualitatively, the energy structure of the energy levels shares common traits with the U = 0 case. Explicitly, degeneracy occurs for a = 0 (with the middle eigenvalue of the 3-site cluster equal to 0) and for a = 2 b (with the middle eigenvalue of the 3-site cluster equal to b ). Neither of these two possibilities offers a solution: on the one hand, when a = 0, the KS densities on varying b are always largely away from the many-body ones; on the other hand, when a = 2 b (i.e. the search is restricted to a one-parameter space), on varying (a negative) b one gets close to reproducing the many-body density (see plot in Fig. 2), but without truly reaching it (inset in the graph of Fig. 2).
These results show that the interacting density is not KSVR in the ensemble sense. 1 However, in a broader perspective, in DFT the main aim of the image system is to reproduce the many-body density, for the present system in the ensemble sense making use of the aufbau principle for the image system (i.e. a degeneracy of the image system could be a somewhat too strict requirement). We thus considered the effect of lifting the degeneracy between the fourth and fifth level of the KS image system. This is schematically shown in Figure 2: by removing the degeneracy, and with two electrons with the same spin projection respectively in the middle level of the 3-cluster and in the third lowest energy level of the 6-cluster (the third lower solid line in the 6-cluster of the KS panel), the interacting density could be reproduced by the total potential a = −1.90996 and b = −1.01479, with an error smaller than 10 −9 . Such occupation scheme is an acceptable image for the triplet S z = 1 state, since the density is reproduced and both the interacting system and the KS image are invariant under spin rotations. The situation was different for the singlet state: we considered a pure singlet representation where the fourth and fifth KS levels are not degenerate and the electrons, according to aufbau and Pauli principles, are in the fourth level. Still, no singlet image could be found, i.e. the system is not even KSVR as a pure state.

Approaching the degeneracy point via TDDFT
Can we find traces of the situation just discussed in a TDDFT treatment? More concretely, would ground-state degeneracy and non-KSVR reflect in the XC potential of TDDFT?
The strategy we choose here to answer is to move to another, slightly varied many-body ground state that is nondegenerate and KSVR, but still with the full symmetry of the square. From there, by exact time evolution, we drive slowly the many-body system to the degenerate, non-KSVR point, then we "stay" there some time and finally go back in the same way. While this is taking place, we extract the time-dependent XC potential via reverse engineering along the trajectory, according to the KS equations (T +v KS )ϕ κ = i∂ τ ϕ κ . As initial state to apply this protocol, we took the many-body singlet ground state that occurs when a = −1.8, b = −1. Such ground state |−1.8, −1 is nondegenerate and KSVR, and one can imagine to drive the system towards the point of interest by changing slowly a . However, using the C 4vadapted cluster partition, it results that for |−1.8, −1 there are N A = 2 electrons in the 3-site cluster, while for the ground state with a = −2.69, b = −1 (henceforth referred to as | − 2.69, −1 ), there are N A = 3 electrons in the 3-cluster. Since a square-symmetric system (including a time-varying perturbation a (τ )) conserves the number N A of particles in the 3-cluster, this implies a level crossing when moving along −1.8 → −2.69 (Fig. 3a). As we are interested in switching adiabatically from the initial singlet state to the critical singlet, we need to avoid such crossing. To smoothly drive the system from singlet |−1.8, −1; N A = 2 to singlet |−2.69, −1; N A = 3 we used a slow but symmetry-breaking perturbation as schematically shown in Figure 3b. The perturbation acts on sites 1, 5, 9 (see Fig. 1c The form of α(τ ) and β(τ ) is shown in Figure 3c. 2 In this way, we do not have to consider the triplet state because S 2 (and S z as well) are conserved along the induced trajectory. We concentrate on the interval τ ∈ [700, 1350], corresponding to the shaded area in Figure 3c. The choice of the time window is motivated by the fact that a smooth, featureless dynamical behavior was observed before τ = 700, but a significant change occurs at about τ = 750, before the singlet | − 2.69, −1; N A = 3 is reached at τ = 1000. Figure 4a shows the time-evolved densities of the interacting system. The black curve shows the central site density. Below that, one can see the four densities of sites 2-5 (blue curves). The four lowest curves correspond to the density at the corner sites 6-9 of the cluster (orange curves). In the middle time region τ ∈ [1000, 1050] where the Hamiltonian is kept constant, curves 2 to 5 as well as 6 to 9 join each other to high accuracy, reaching a value in very good agreement with the ground state densities for such Hamiltonian (we also monitored the time dependent average number of electrons in the 3-site cluster during the time evolution; this is 2 at τ = 0 and becomes 3 at good accuracy when 1000 ≥ τ ≥ 1050). Figures 4b and 4c concern the reverse-engineering procedure and the results for the KS potential; the accuracy of the procedure is displayed in panel b, while the actual potential is shown in panel c. We note that v KS starts 2 With reference to Figure 3c to exhibit oscillations at τ ≈ 750: this should be put in relation with the behavior of the KS instantaneous eigenvalues, shown in Figure 4d: it is apparent that the forth and fifth lowest KS eigenvalues also cross at where the oscillations of v KS begin. It is also apparent that driving the perturbation backwards (corresponding to t > 1050) does neither leads to a reverse KS crossing nor to the disappearance of the potential oscillations. We interpret such oscillations as a signature in TDDFT of the approaching of the representability failure of ground state ensemble DFT in our model. As already mentioned, the oscillating v KS reproduces the interacting density to high accuracy, also in the region where the Hamiltonian is kept constant (in this region, the densities are constant within 10 −4 ). The oscillation in v KS is of relatively high magnitude and thus definitely necessary to keep the density constant to this accuracy.

How slow is slow?
We now look more closely at how the oscillations of the KS potential depend on the degree of adiabaticity of the symmetry breaking perturbation. In Figure 5, we present results for three different perturbation speeds, where the highest speed (orange curve) is already considerably (eight times) less than the one in Figure 4. The remaining two cases, represented by green and blue curves respectively, correspond to half and quarter speed of the orange curve (the temporal profiles of the perturbation normalized to unity for these three cases are shown in panel c).
On slowing down the perturbation, the oscillation of the density at the central site (panel a) becomes constant within progressively higher accuracy (note the vertical scale), both at the degeneracy point of ground state degeneracy (i.e. in the interval where the value of the perturbation profile is 1 in panel c) and towards the end of the time evolution (where, in a final interval, the perturbation has returned to zero).
For each perturbation speed, the behavior of the KS potential (panel b) resembles what seen previously in Figure 4, i.e. with strong oscillations in both the central and final intervals. However, the oscillations diminish in magnitude for slower perturbations, and it would thus be interesting to determine if they tend to a limiting nonzero value for infinitely slow perturbations. We attempted an analysis of the frequency and amplitude of the oscillations (not shown) based on the three cases of Figure 5; however, this did not provide clear trends, because of the small data set available (three cases) and the non fully monochromatic character of the oscillations. At the same time, the robustness of the oscillations was confirmed by comparing, for the fastest of the three perturbations, results from different timesteps (as also done for the results of Figure 4).
Altogether, our results are strongly indicative that for any finite-speed perturbation such oscillations are fundamental in nature, while evidence of their existence at infinitely-slow speed remains inconclusive at this stage. Finally, it is perhaps worth mentioning that we also approached the degeneracy point from other initial values of the onsite parameters. In some instances, when such values were considerably "close" to degeneracy, our reverse-engineering procedure did not converge, hinting at a possible TDDFT representability issue for dynamics starting in a small neighbourhood of the degeneracy point.

Conclusions
By applying DFT and TDDFT to a simple, finite lattice model, and by engineering the system in such a way that it has a spin-degenerate ground level, we have shown that it may be not possible to obtain an image of a degenerate interacting triplet-singlet ground states within groundstate ensemble-DFT. In these situations a possibility is to resort to an approach where each degenerate state with a different symmetry (in our case triplet and singlet states have different spin-symmetry) is treated in its own symmetry subspace (i.e. so-called "DFT in subspaces" [75,76]). Our results show that, when going to a "DFT in subspaces" approach, we can obtain a KS image for the triplet case, but not for the singlet. Here, it is important to point out that the situation presented concerns the case of zero temperature, i.e. a limiting case of Mermin's formulation [4]. It may very well be that for finite temperatures (where, for both exact and image systems, thermal averages are involved) this failure of ground-state ensemble-DFT would disappear: in that case, beyond the aufbau principle, other higher energy states of the image system are included, and with more possibilities to represent the interacting density at thermal equilibrium.
Although rather peculiar in character, and with the limitations just addressed, the model at hand permits to examine how a lack of ground state v 0 -representability induces signatures in a TDDFT description. In particular, we found that starting from a non degenerate ground state, and approaching the degeneracy point via time evolution, a level crossing occurs in the KS system. This in turn unveils an interesting behavior of the XC potential, which retains an oscillating pattern (even when the external field is kept constant) to ensure a constant density.
We also tried to ascertain if such oscillations persist in the adiabatic limit (infinitely slow perturbation), but our results are not conclusive in this respect; at this stage, we can only safely conclude that, for perturbations of finite speed, the observed oscillations remain, a consequence of the level crossing, and an indication of the highly non-local nature of the XC potential.
From a conceptual point of view, a constant density can arise from an oscillating KS potential even for a constant external field, for example if the density at τ = 0 pertains to a system prepared in a non stationary state. This would make the behavior observed in our cluster to be understood as related to the initial state dependence of v KS . While this can play a role, it also not unreasonable to conjecture that our TDDFT results largely emerge as a manifestation of non-representability in the ground state; however, making this statement more precise is beyond the aims of this work.
Most likely, our results are not relevant to real systems where due to Coulomb interactions since a strict system separability does not occur (an approximate separability could be possibly argued/achieved for very distant subsystems). It is also possible that akin results would emerge in a system simpler than the one studied here. However, searching for such system(s) is also outside the scope of the present study. In any case, to our knowledge this specific behavior of the exact KS potential has not been explicitly reported before, and adds to the growing evidence of non-trivial features (spikes, discontinuities, steps, etc.) that v KS has come to be known to exhibit.
Author contribution statement C.V. and C.O.A. conceived the idea of the project and supervised T.R. while he performed all ground state and preliminary time-dependent calculations; C.V. performed the final time-dependent calculations. All authors contributed to the writing of the manuscript.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.