Dynamical exchange-correlation potential formalism for spin-$\frac{1}{2}$ Heisenberg and Hubbard chains: the antiferromagnetic/half-filled case

The exchange-correlation potential formalism previously introduced and applied to the one-dimensional Hubbard model has been extended to spin systems and applied to the case of the one-dimensional antiferromagnetic spin$-\frac{1}{2}$ Heisenberg model. Within the spin exchange-correlation potential formulation, a new sum rule for spin-systems is derived. The exchange-correlation potential for the Heisenberg model is extrapolated from exact diagonalization results of small antiferromagnetic Heisenberg clusters. This procedure is also employed to revisit and computationally improve the previous investigation of the exchange-correlation potential of the half-filled Hubbard model, which was based on the exchange-correlation potential of the dimer. Numerical comparisons with exact benchmark calculations for both the Heisenberg and the Hubbard models indicate that, starting from the exchange-correlation potential of a finite cluster, the extrapolation procedure yields a one-particle spectral function with favorable accuracy at a relatively low computational cost. In addition, a comparison between the ground state energies for the one-dimensional Hubbard and Heisenberg models displays how the well known similarity in behavior of the two models at large interactions manifests within the exchange-correlation potential formalism.


I. INTRODUCTION
Lattice models, in spite of their apparent simplicity, can be very valuable to reveal important features in low dimensional and highly correlated quantum systems.This certainly is the case of two highly paradigmatic models of condensed matter physics, namely the Hubbard [1] and spin- 1  2 quantum Heisenberg models [2].For several decades, these two models have been a testground for new theoretical and computational methods [3][4][5].Notably, they have been used to describe phenomena such as the Mott transition [6], high T c superconductivity [7], quantum spin liquids [8], and quantum entanglement [9,10].Furthermore, via suitable parameterization from first-principles ground-state calculations, they have also been used to describe the dynamical behavior of real materials, which is experimentally measurable via, e.g., neutron scattering and angle-resolved photoemission spectroscopy.This model approach is very useful when first-principles descriptions are too complicated to perform.(see e.g.[11][12][13][14]).
There are a number of approaches of increasing sophistication being continuosly developed to solve the Hubbard and Heisenberg models [15][16][17][18][19][20][21][22].Exact analytical solutions remain scarce.In one dimension (1D), both models are integrable and exactly solvable via Bethe ansatz [23,24].Yet, exact analytic treatments for higher dimensional or even extended 1D systems (e.g., with nextnearest-neighbor coupling) are in general not available.As it happens, already in 1D not all quantities of interest can be accessed: the Bethe ansatz provides information about the energy dispersion [25,26] but not, for example, the spectral weight, one of the more interesting quantities to consider when studying dynamical correlations, which are usually directly connected to experimental results.
On the numerical side, several approaches can be suitably employed for both models, such as Exact diagonalization (ED) [27], Quantum Monte Carlo (QMC) [28][29][30], and Density Matrix Renormalization Group (DMRG) [31][32][33], to name a few.ED gives exact and complete information about the system, but is restricted to small systems, thus unable to capture the thermodynamic limit features.DMRG and QMC are applicable to fairly large systems and with high accuracy in 1D [34][35][36], but for higher dimensions the computational cost increases rapidly [37][38][39].
The local-density approximation (LDA) and its extension to local-spin-density approximation (LSDA) are widely used in DFT [43,57,58].L(S)DA successfully describes many materials, but does not perform well in strongly correlated systems, and much effort has been devoted to improving it.With focus on model lattice systems, one way is to use the exact Bethe ansatz solution of the Hubbard model to approximate the correlation energy of an inhomogeneous lattice system [59].A similar employment of DFT has also been considered for the Heisenberg model [60].What is noteworthy about these L(S)DA approaches when applied to the Hubbard and Heisenberg models is that the exchange-correlation term has information about the lattice structure and dimensionality of the system.
From a different perspective, a formalism based on the dynamical exchange-correlation potential (Vxc) was recently introduced [61].The formalism is not limited by system size, system dimensionality, or type and range of the interaction, and it is thus useful to describe electronic and magnetic structures in general situations.A main feature of the dynamical Vxc formulation is that the coupling between the dynamical Vxc and the Green function occurs as a direct product in space and time.In contrast, the self-energy, which is traditionally used to calculate the Green function, acts on the Green function as a convolution in space and time.
As a first application of the framework, the lattice oneparticle Green function of the infinite 1D Hubbard chain was determined [61,62] using an extrapolation scheme, starting from the dynamical Vxc of the Hubbard dimer as input.In spite of the simplicity of the approximation used and the low computational load, the scheme provides estimates of the band gap and spectral function in favorable agreement with the results obtained from the Bethe ansatz and the Dynamical Density Matrix Renormalization Group (DDMRG) [63].One general conclusion from this investigation is that the Vxc formalism provides a simple picture of the one-electron spectrum: for a given momentum, a time-independent term in Vxc together with the kinetic energy term determine the main peak of the spectral function, while a time-dependent term in the form of an exponential couples the Green functions with different momenta and generates incoherent structures or satellite peaks.The energy variable appearing in the exponent can be understood as the main bosonic excitations of the system.
More recently, as a step towards the study of realistic systems, the Vxc of the homogeneous electron gas was calculated within the random-phase approximation [64] with the long-term aim of constructing the Vxc as a universal functional of the ground-state density within the local-density approximation.

II. THIS WORK, AND PLAN OF THE PAPER
In this work, the Vxc framework is extended to spin systems, more specifically to the 1D Heisenberg model.The Vxc-based equation of motion and the sum rule for the spin exchange-correlation hole are derived.Furthermore, the extrapolation scheme employed in the previous work for the 1D Hubbard chain is adopted.The essential idea of the extrapolation scheme is to start from the Vxc of a finite cluster (kernel), which can be calculated accurately using an exact diagonalization method or other methods such as the density-matrix renormalization group.By a suitable extrapolation, this is then used to determine the Green function of the corresponding lattice model.The spin Vxc framework within the extrapolation scheme is applied to calculate the spectral functions of the 1D spin− 1 2 antiferromagnetic (AFM) Heisenberg model in the thermodynamic limit, starting from the spin Vxc of small clusters.
In addition, the 1D Hubbard chain is revisited.In the previous work, the Hubbard dimer was the kernel, which was used to calculate the Green function of the 1D Hubbard chain.In this work, in order to improve the quality of the starting Vxc, the cluster size is enlarged so that additional information arising from interactions beyond nearest-neighbor is captured.The improved Vxc is then used to calculate the Green function of the halffilled 1D Hubbard chain.
To summarize, the main outcomes of the present work are: (i) derivation of the Vxc-based equation of motion and the sum rule of the spin exchange-correlation hole for the 1D Heisenberg model, which can be readily generalized to other spin systems; (ii) calculations of the spinon Green function for the 1D AFM Heisenberg lattice by extrapolating from a finite-cluster spinon Vxc; (iii) improved treatment of the Vxc of the half-filled 1D Hubbard lattice from the previous work by using as kernel a Vxc from a finite cluster; (iv) illustration on how in the Vxc formalism the well known large-U limit (where results from the Hubbard model match those from the AFM Heisenberg one) is recovered.
The plan of the paper is as follows: in Section III, we review briefly the general Vxc formalism.Then, in section III A and III B we extend and apply the approach to the 1D AFM Heisenberg model.Specifically, in section III C and III D, we derive an analytic expression for the spinon Vxc for a four-site chain, and compute the lattice dynamical structure factor by extrapolating the finite cluster Vxc to the infinite case.In section IV, we revisit the 1D Hubbard model and compute the exact Vxc of a finite cluster larger than the dimer, with which we improve previous results in the infinite chain limit.In Section V we discuss Vxc from a comparative perspective, addressing the ground-state energy for both the 1D AFM Heisenberg model and the half-filled 1D Hubbard model in the large U limit.Finally, in Section VI we provide some conclusive remarks and an outlook.

III. GENERAL FORMALISM AND APPLICATION TO THE HEISENBERG CHAIN
For a system with a one-body term and two-body interactions, the Hamiltonian reads Ĥ = dr ψ † (r)h 0 (r) ψ(r) where ψ(r) is the fermionic field operator and r = (r, σ) is a combined space and spin variable.The time-ordered Green function is defined in the Heisenberg picture as where the argument numbers label the space-time 1 := (r 1 , t 1 ), ⟨.⟩ denotes the zero-temperature ground-state expectation value, and T is the time-ordering symbol.The equation of motion in the Vxc formalism is given by [61] [ where the single-particle term contains the Hartree potential The Vxc reproduces the interaction term containing a special case of the two-particle Green function , i.e., For fermion field operators and in the presence of Coulomb interactions, the bare exchange part of Vxc can be obtained by considering the lowest order of the first term on the RHS of Eq. ( 6), A. Spin-spin interactions For systems with spin-spin interactions, an observable of central interest is the spin dynamical structure factor, whose longitudinal and transverse terms are and where Ŝz,+,− p (t) are the spin field operators in the Heisenberg picture.
For the Hubbard model, the spin dynamical structure factor can be obtained by solving a two-particle Green function but the equation of motion of the two-particle Green function contains three-particle Green function, and thus is generally difficult to solve.Simplification is however recovered for large repulsion, where charge transfer becomes less likely and spin correlations can be obtained by studying the AFM Heisenberg model.It is thus of fundamental and practical interest to discuss the Vxc formalism directly for the Heisenberg model.The isotropic 1D Heisenberg Hamiltonian with nearest neighbour (NN) exchange coupling is given by where for convenience we use an even total number of sites before taking the thermodynamic limit.We define the Green function with spin field operators in which the Heisenberg J is the analog of the twoparticle interaction in Eq. ( 1).From the Heisenberg equation of motion for the spin field operators, the equation of motion of the Green function reads where the interaction term is Here, and J pl = J(δ l,p+1 + δ l,p−1 ) for the 1D NN exchange coupling.One can define the spin exchange-correlation potential analogous to the charge case as follows: where the last two terms on the right-hand side , V H and V F , are the analog of the Hartree and exchange potentials, respectively: Consequently, a spin correlator g lpq (t) can be defined such that while the spin exchange-correlation hole ρ xc is defined as Denoting the total z-component of the spin by S z = l ⟨ Ŝz l ⟩, and observing that we can obtain a sum rule for general spin interactions: The detailed derivation is provided in Appendix A.
In this paper, we consider only the case of AFM coupling, i.e.J < 0 so that S z = 0.For a translationally invariant system, the Hartree and Fock terms (Eq.( 17), ( 18)) vanish, and thus the two-spinon Vxc is then with the corresponding exchange term given by

B. The infinite chain
We specialize now the description to the case of the homogeneous infinite Heisenberg chain, where ⟨S z p ⟩ ≡ s is site-independent due to translational symmetry.It is convenient to move to the momentum domain, with the Green function and V xc defined via the Fourier transform as and where the equation of motion for the Green function becomes In the momentum representation, the exchange term becomes where we have neglected the k ′ -dependence in the weight represented by G(k ′ , 0 − ), performed the sum over k ′ and subsumed all the constants into λ.To proceed further, the dynamical part of Vxc is separated from the static exchange term to finally arrive at the solution to the equation of motion Eq.( 26): In this expression, the (k-dependent) static exchange term V s determines the main peak of the spectral function, and the dynamical correlation term Z sp (k, t) produces the satellite structure.To attain an explicit solution, it is expedient to solve for a reference Green function by keeping only the static V s term in the equation of motion.This simplified solution contains the lower boundary of the two-spinon energy dispersion [15] and permits to determine the constant λ in Eq. ( 27) from the analytic form of the two-spinon spectrum.

C. A four-site spin chain
It is useful to start our discussion of Vxc in finite spinclusters by considering a four-site chain.This is the minimal cluster (with even number of sites) in which Vxc is nonzero.Furthermore, it is easy to obtain a compact analytical solution, that illustrates qualitatively several features present also in larger clusters (in which our solution is numerical in character).To illustrate the features of the four-site Vxc, we choose one of its diagonal elements as a representative case, namely where For the Green functions, the transformation reads such that the equation of motion is now where s µν = 2 pq U µp ⟨S z p ⟩δ pq U * qν .Comparing the equation of motion for the diagonal terms G µµ , to the infinite-chain equation of motion Eq. ( 26), we note the following: i) G µµ maps to G(k); ii) the contribution from fully off-diagonal terms V xc µν,δµ should be negligible; iii) V xc µγ,γµ , which maps to V (k), depends only on the difference of µ, γ; iv) the weights of the higher excitation term f 3 are relatively small.
According to i-iv) and ignoring the high energyexcitation contributions from f 3 , one thus arrives at an approximate expression for the matrix elements of V xc µγ,γµ : whereas (38) and as in (31), The analytic spinon Vxc in the bonding basis and its main excitation approximation are shown in Fig. 1.Ignoring the high-excitation factor f 3 reduces the fine-structure details in Vxc.Consequently, V xc BB,BB simplifies to a constant whereas V xc

BC,CB
oscillates with a single frequency and a constant magnitude, and all other components are negligible.36)-( 38) and related discussion).

D. Infinite chain from cluster extrapolation
In Fig. 2, we show ReZ sp (k, t) as obtained from the cluster Vxc discussed in the previous section.
It can be seen that ReZ sp (k, t) oscillates in time around a momentum-dependent term, a behavior that can be understood as due to a single quasiparticle-like main excitation.We therefore propose the following ansatz for Z sp in the infinite-chain case: where the amplitude A, the spinon excitation energy ω sp , and the shift term B all increase as k increases from 0 to π.The Green function is given by inserting the ansatz into Eq.(29): where the static potential is V s (k) = −Jπ| sin k|/2.
Expanding the last term on the RHS of Eq.( 40) to first order in e −iω sp (k)t , one gets an approximate Green function which in the frequency domain becomes From Eq. ( 42), it can be seen that the main peak position of the dynamical structure factor is given by V s + B.
The spinon excitation energy ω sp transfers weight from the main peak to higher-energy region resulting in satellite peaks at V s + B + ω sp .The relative weight between the main peak and the satellite is determined by the amplitude term A and the spinon energy ω sp .Specifically at k = π, the finite cluster solution gives nonzero B, which opens a spin gap that does not exist for the spin− 1 2 lattice.We attribute this to finite-size effects, and thus we adjust B to a smaller value in our extrapolation.
Based on our discussion so far, we now present the lattice case obtained by extrapolating the cluster Vxc.With Z sp obtained via ED from a twelve-site cluster, we estimate the parameters A, B and G(k, 0 + ) by linear interpolation.The spinon excitation energy is estimated by fitting the cluster ω sp to the two-spinon spectrum boundary, The longitudinal and transverse spin dynamical structure factors are then calculated from the spinon Green function.Since for a spin isotropic system S zz and S +− differ by a constant factor, we only calculate the spectral function of the Green function (Eq.( 40)), as shown in Fig. 3 (interestingly, approximating the term exp A ω sp exp(−iω sp t) − 1 by 1 + A ω sp exp(−iω sp t) − 1 gives no marked changes of the properties of G sp (k, ω)).A notable aspect in the behavior on the spin dynamical structure factor is that both the peak locations and the relative weights are close to the inelastic neutron scattering data from 1D compound KCuF 3 [21].Coming to more specific features, S(k, ω) is very small (i.e., close to zero) at small k, while, for a generic k, most of its spectral weight is concentrated around the main peak and the satellite peak.As k → π, the relative weight between the main peak and the satellite peak increases and the spectrum with broadening factor 0.1 is gapless.
While providing a good description of the dynamical structure factor for the 1D AFM Heisenberg model, the present implementation of the spin Vxc approach is also subject to some limitations.This can be seen by, e.g., comparing the dynamical structure factor from the Vxc approach with the two-spinon lower and upper boundaries (dashed line in Fig. 3).It is apparent that the main peak frequency ω = V s (k) + B is still slightly overestimated.To reduce the finite size effects due to a parameter B(π) originating from a twelve-site cluster, we set B(π) to be the same as the broadening factor, i.e. about 0.2 (see Fig. 2).However, the actual Bethe ansatz value of B(π) should be zero.The overall point is that, to obtain a more accurate dynamical structure factor, and to avoid the finite size effects inherent in the extrapolation from a small cluster, more powerful external methods need to be employed (e.g., the algebraic Bethe ansatz).
These considerations might reveal weaknesses of the extrapolation procedure.However, it must also be clearly stressed that this implementation of the Vxc approach captures most of the qualitative features of the 1D AFM Heisenberg model with a very low computational load and this central attractive feature of the method is expected to also apply in more challenging situations, e.g. in higher dimensions, where rigorous references like the Bethe ansatz are not available.

IV. IMPROVING THE TREATMENT OF THE 1D-HUBBARD LATTICE
Encouraged by the 1D AFM Heisenberg chain results obtained with a Vxc extrapolated from clusters, we now revisit the case of the 1D Hubbard Hamiltonian, using also in this case a Vxc obtained from small (Hubbard) clusters.In Eq. ( 44), p = 1, 2, • • • , N are the site labels (with N → ∞ eventually), σ =↑, ↓ is the spin label, ∆ is the hopping energy and U > 0 is the local repulsion.
In the site basis, the spin-up channel Green function is and the Vxc reads (46) where ρ p↓ is the spin-down particle density at site p.The exchange part of Vxc fulfils where G pp (0 − ) = i⟨ĉ † p↑ ĉp↑ ⟩ = iρ p↑ ; thus, the exchange part of Vxc of the Hubbard model is static and cancels the Hartree potential at half-filling, in contrast to the Heisenberg model in which the exchange part depends on the momentum.In general, the exchange part is timedependent [64].
Written in the momentum domain, the equation of motion for the Hubbard lattice takes the form where ε k = −2∆ cos k is the kinetic energy.Eq. (48) shows that the interaction term, expressed as the direct product of Vxc and Green function in space-time domain, is a convolution in the momentum domain.It has been shown [62] that the main peak position of the electron (hole) spectral functions can be described with V xc (k = 0), together with the kinetic energy, while V xc (k = π) plays an important role in determining the satellite peaks.One can also write the interaction term as a direct product in momentum domain, which gives explicitly the solution for the Green function: One can then use an N -site cluster with twisted boundary conditions [65] to parameterize G(k, 0 ± ) and thus the generalized Vxc in the momentum domain becomes

A. Extrapolation from finite clusters
The Vxc of clusters with 6 and 8 sites and with periodic boundary condition was computed using ED.The Hubbard U was chosen to be 7.74 with ∆ = 1, to allow for comparisons with previous work and the DDMRG results from the literature.In contrast to the dimer case, the cluster Vxc exhibits multiple sharp peaks as a function of time t.Time snapshots of Vxc as a function of k are shown in Fig. 4. For t ≃ 0, we have that ), but such behavior is unseen during the time evolution.The particle-hole symmetry leads to V xc (k, −t) = −V xc (k, t), and the increase of cluster size from N = 6 to 8 does not change qualitatively the characteristics of V xc as a function of k.
The dynamical properties of Vxc can be better illustrated through Z el , a generalisation of Vxc in the momentum basis defined in Eq. ( 51).Due to degeneracy, Z el (−k, t) = Z el (k, t) and, because of particle-hole symmetry, Z el (k, −t) = −Z el (π − k, t).To improve the simulation of Z el , we use a cluster with twisted boundary condition, that provides larger k-point sampling.The real part of Z el (k, t) with twisted boundary condition is shown in Fig. 5: for small k, it oscillates weakly in time (with small amplitude and long period).However, where the bandgap opens (k → π 2 ), the oscillation of ReZ el is more evident.For k → π, ReZ el exhibits sharp peaks at certain times.The peaks can be both positive and negative: mathematically, this means that some of the zeros of the Green function are located where the interaction term (Eq.( 46)) has nonzero finite (positive or negative) values.These spiky structures cannot be fitted into a weighted sum of several (but finite in number) oscillations, indicating that a model beyond the single-energy quasiparticle picture is necessary.
Provided with the numerically exact Vxc for N = 6, 8 clusters, we reconsider the approximate scheme proposed in the previous work based on the Hubbard dimer (N = 2) [62].The dimer admits two k-points (k = 0, π), with the corresponding approximate values for Vxc given by where η is a broadening factor.The spectral function of G h has a main peak at ω h k , determined by V xc (k = 0) and by the kinetic energy: gives rise to a continuous satellite region.Its relative weight to the main peak is V xc (π, 0), and its lower/upper boundaries are given by the minimum/maximum occupied state kinetic energy The dimer model [62] managed to capture the main structure of the hole spectra of the Hubbard lattice, but can be improved in several aspects: The main peak position given by the model is just the kinetic energy ε k = −2∆ cos k plus a constant determined by U , while the true k−dependence of ω h should be more complicated; the upper and lower boundaries of the satellite part given by the model are independent of k, which is also an oversimplification.Rewriting Eq. (50b) in the spirit of a factorization into a main-peak and a satellite term, where Z h,main k + Z h,sat (k, t) = Z el (k, t) for t < 0, one can see that: i) A momentum-dependent static term, Z h,main k , which is not present in the dimer model, together with ε k , determines the main peak; ii) the dispersion of ω h,lower and ω h,upper can be explained by the satellite term Z el,sat (k, t ′ ).Compared with Fig. 5, Z h,main k is seen to be the time-independent part around which Z el (k, t) oscillates; and Z h,sat (k, t) represents a series of excitation energies.The spike-like ReZ(k, t) for k → 0, t < 0 is a consequence of multiple excitation energies and large satellite peaks, while the less oscillatory ReZ(k, t) for k → π, t < 0 explains the lack of strong satellites of the hole spectral functions A h (k → π, ω).
Taking advantage of the physical picture given by the dimer model, we include the correction to the occupied k values by adding a set of momentum-dependent parameters, l 1,2,3 , such that i) α → αl 1 (k), ii) the main excitation determining the satellite boundaries (Eq.(55a), (55b)) becomes 2∆ → 2∆l 2 (k), and the effective kinetic energy in the summation of Eq. ( 53) becomes ε k ′ → −2∆ cos k ′ l 3 (k).The parameterized dispersion relations of the key frequencies are Thus, the hole part bandwidth for a given momentum, the satellite width, and the band gap respectively are This means that the main peak location, the bandwidth, and the satellite region width from cluster calculations can be used to determine the parameters l 1,2,3 , which are then used to calculate the lattice spectral functions A(k, ω) for k < π 2 .For k > π 2 , where the dimer model gives zero weight for the hole part spectrum, the cluster results show that the corresponding Vxc can be approximated with a single-energy excitation, where the parameters A, B and ω el k are estimated from cluster result, which is similar to the treatment for the spinon Vxc (Eq.( 39)).Combing the l 1,2,3 -involved occupied region and the A, B, ω el -involved unoccupied region, the hole part spectral function can now be calculated for the whole Brillouin zone.The spectral functions for selected k values are shown in Fig. 6.
Compared with the dimer model, the cluster Vxcbased parametrization improves the agreement with DDMRG in several aspects.Specifically, i) The missing weights for unoccupied k points appear when using as input a cluster Vxc.ii) The main peak positions (and thus the bandgap value as well) are more accurate.In fact, the bandgap value from the dimer model, αU , shows a discrepancy with the Bethe ansatz exact value at small U , do to the lack of long-range screening effects.Using a cluster Vxc, however, removes the disagreement.iii) Both boundaries and relative weight of the satellite structure are better described by the cluster Vxc and its momentum-basis generalization Z el ; iv) The total weight of the hole/electron part cannot be renormalized within the dimer model, because the non-interacting Green function used in the dimer model can only fix the total spec- tral weight: dωA h (k, ω) = θ(k F − k).With a cluster Vxc, using ⟨ĉ † k ĉk ⟩, we can rescale the total spectral weight for each k value.
Yet, the main peak ω h :s in Fig. 6 is in general lower than the one from DDMRG.This can be understood as due to the band gap narrowing upon increasing the numbers of site (the eight-site cluster we used leads to the overestimation of the gap and thus of the main peak position).
We conclude our discussion of the Hubbard chain, by considering its spectral functions in real space that we obtain starting from those in the momentum domain: where r = 0, 1, 2, • • • is in units of lattice parameter.A(r, ω) describes the correlation strength between two space points separated by r, at a given energy ω.The local case A(r = 0, ω) corresponds to the density of states.
Results for A(r, ω) with a eight-site kernel are shown in Fig. 7, whilst those from a six-site kernel with different U and r are reported in the appendix.The cluster Vxc result for A(r = 0, ω) shows better agreement with DDMRG than the dimer model.Also, the NN spectral weight at positive energy is predominantly negative, and for r ≥ 2, A(r, ω) exhibits nodal structures.Concerning the role of electronic correlations, spatial spectral functions with different U values become qualitatively alike at large repulsion (U > 4), but the band gap value keeps increasing with U .Finally, spectral functions calculated with eight-site and six-site kernels are qualitatively similar (see the appendix for the six-site case), with similarities in the overall shape and in the number of nodes.However, the estimated value of the band gap improves on increasing the cluster size.For both models, Vxc is extrapolated from a six-site kernel.

V. VXC FROM HUBBARD AND HEISENBERG MODELS: A COMPARATIVE DISCUSSION
It is well known that the 1D spin− 1 2 AFM Heisenberg model becomes equivalent to the 1D half-filled Hubbard model in the large U regime [67,68].After having discussed Vxc in the two models separately, it can be useful to look at both models together using as perspective the behavior of Vxc in such limit.Meanwhile, Z el and Z sp do not show a direct asymptotic behavior Z el U →∞ = Z sp , because they are coupled to the single-particle Green function (Eq.( 45)) and single spin-flipping Green function (Eq.( 12)) respectively.For the Hubbard model, the term corresponding to Z sp is coupled to the twoparticle Green function ⟨T ĉ † p↑ (t)ĉ p↓ (t) ĉ † q↓ ĉq↑ ⟩.Equation of motion of higher order Green function needs to be solved for the Hubbard model to calculate the 'higher order Vxc' that is comparable with the spinon Vxc under large repulsion.This means the Vxc formalism for Heisenberg model, having a similar sum rule (Eq.( 22)), reduces the difficulty in deriving the equation of motion and improves the interpretability via the quasiparticle picture.
Instead of solving the higher order Green function, we consider a more modest task of comparing the lattice ground state energies for the two models.In the large U limit [68], lim where E Hub 0 is the ground state energy of a N -site Hubbard ring with ∆ = 1, and E Heis 0 is the ground state energy of a N -site AFM Heisenberg ring with J = −1.Both energies can be calculated from the Green function via and In the frequency domain, To perform a comparison, we compute the ground state energy of the Hubbard lattice in two ways: i) by directly using the electron Vxc at different U values, and ii) by calculating E Heis 0 for a J = −1 Heisenberg lattice with the spinon Vxc, to be then used in the effective E Hub 0 of Eq. ( 61).The differences between the results from these two prescriptions and the exact Bethe ansatz solution are shown in Fig. 8.The E 0 results from ED for a six-site ring are also shown as a reference.
For U < 10, the repulsion strength is not large enough for Eq. ( 61) to be valid, leading to a discrepancy between the total energies for the two lattice models.However, in such region, E Hub,Vxc 0 (red dots) is already close to the exact Bethe ansatz value, and the difference gets smaller on increasing U .For U > 30, the ED results for the two models converge, meaning that the large repulsion limit is reached.The Vxc-based energies E 0 for the two models also converge to the exact Bethe ansatz value.
However, the effective Vxc-based Heisenberg result is rather accurate, with absolute error less than 10 −4 : this can be understood as a result of i) using the two-spinon upper and lower boundaries in the extrapolation, and ii) adjusting the B parameter from the cluster within the zero spin gap picture.In contrast, the Vxc-based Hubbard result is extrapolated without a good reference, and is more affected by the finite size effects.Thus, the difference with the Bethe ansatz result is larger.
As an overall remark, the comparative analysis of this Section shows the versatility of the Vxc approach across different lattice models, with results that are consistent with trends and benchmarks from other methods.

VI. CONCLUSION AND OUTLOOK
We have presented a novel exchange correlation potential (Vxc) formalism for the one-dimensional antiferromagnetic Heisenberg model, and derived a general new sum rule for spin systems.Our spin-formulation is a tailored extension of a previously introduced general framework for many-body systems that include both charge and spin degrees of freedom.Together with the new formulation, we have also devised a procedure to obtain, from a Vxc extracted from small finite clusters, an extrapolation to the thermodynamical limit.This procedure to access Vxc, originally devised for spin systems, has also permitted us to revisit and improve the treatment of the half-filled one-dimensional Hubbard model, a system already considered in earlier work within the Vxc approach.For both the 1D AFM Heisenberg model and the 1D Hubbard model, the static exchange term of Vxc was derived and shown to exhibit model distinctive properties.For the 1D AFM Heisenberg model, the static exchange term corresponds to the lower boundary of the two-spinon spectrum.For the Hubbard model, the local U leads to a constant V x , which cancels the Hartree potential.
For both models, the spectral functions calculated within the Vxc approach show favourable agreement with DDMRG and with experimental results.Furthermore, a single-energy quasiparticle picture can be used to explain the dynamics of the spinon Vxc for the 1D AFM Heisenberg model and the unoccupied/occupied part of the hole/electron Vxc for the 1D Hubbard model.Finally, we showed how the Vxc formalism captures the equivalence of the two models in the large U limit, by a comparative analysis via the lattice ground state energies.
In conclusion, our results indicate that the Vxc formalism provides an alternative way of calculating singleparticle Green function which is (computationally) costbeneficial but also physically well defined.Looking forward, we plan to apply such dimensionality-and interaction-insensitive scheme to models of increasing complication and higher dimensionality.At the same time, we intend to explore ways to devise Vxc approximations with the goal of improving over the local-(spin)density approximation of density functional theory, in a progression toward a first principle implementation for real materials.
and f i=1,2,3 are time oscillation factors determined by the difference between the spin excitation energies and the ground state energy.The full details and the explicit forms are given in Appendix B, together with other elements of Vxc.It is useful at this point to move from site orbitals {φ a } to bonding-like ones {ϕ µ }.In analogy to what is done with a Bloch basis in a Hubbard lattice, we set ϕ µ = U µa φ a , φ a = U aµ ϕ µ , in which µ = A, B, C, D, a = 1, 2, 3, 4 and the U matrix is

FIG. 2 .
FIG. 2. Real part of Z sp from a spin-1 2 AFM Heisenberg ring.Top (bottom) panel: results for a ring with 8 (12) sites.
52b) Here, the constant α depends only on U ∆ (the explicit dependence relation is shown in Appendix C together with a summary of the properties of the Vxc obtained from the Hubbard dimer) and 2∆ in the exponential represents the main excitation energy.In what follows, we use Eqs.(52a) and (52b) to compute the hole part of the spectral function, with the particle part obtainable via the particle-hole symmetry A e (k, ω) = A h (π − k, −ω).When |k| ≤ π 2 , the hole part of the Green function given by the dimer model is

FIG. 5 .
FIG. 5. Real part of Z el (k, t) of finite Hubbard chain with twisted boundary condition, U = 7.74, ∆ = 1, in the unit of U .Left and middle panel: with shorter time scale and fewer k-points.N = 8, 6, respectively.Right panel: N = 6, longer time and more k-points, peaks out of the color scale not shown.Z el (−k, t) = Z el (k, t) and Z el (k, −t) = −Z el (π − k, t).For a discussion of the negative peak in the middle panel, see the main text.

FIG. 7 .
FIG. 7. Spatial spectral functions, calculated with the eightsite Vxc as kernel (for the six-site case, see the Appendix) for U = 7.74, ∆ = 1 and broadening η = 0.1.64 k-points are used to approximate the k-integral, according to the Chadi-Cohen method [66].

FIG. 8 .
FIG. 8. Difference between the exact Bethe ansatz E Hub 0 N and the Vxc-based results for i) the 1D Hubbard model and ii) the 1D AFM Heisenberg model, and the ED results for iii) a six-site Hubbard cluster and iv) a six-site Heisenberg cluster.For both models, Vxc is extrapolated from a six-site kernel.