{'fraction_non_alphanumeric': 0.03757234574148455, 'fraction_numerical': 0.012745501122742654, 'mean_word_length': 4.1980930461943124, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Solving lattice density functionals close to the Mott regime

We study a lattice version of the local density approximation (LDA) based on Behte ansatz (BALDA). Contrary to what happens in Density Functional Theory (DFT) in the continuum and despite its name, BALDA displays some very non-local features and it has a discontinuous functional derivative. The same features prevent the convergence of the self-consistent Kohn-Sham cycle thus hindering the study of BALDA solutions close to a Mott phase or in the Coulomb blockade regime. Here we propose a numerical approach which, differently from previous works, does not introduce ad hoc parameters to smear out the singularity. Our results are relevant for all lattice models where BALDA is applied ranging from Kondo systems to harmonically trapped Hubbard fermions. As an example we apply the method to the study of a one-dimensional lattice model with Hubbard interaction and a staggered potential which can be driven from an ionic to a Mott insulating state. In the Mott regime the presence of a"vacuum"allows us to calculate the different contribution to the gap and to highlight an ultranonlocality of BALDA.


I. INTRODUCTION
The extension of Density Functional Theory (DFT) to treat lattice fermions dates back to the 80's and it has been recently the subject of a revived interest.  One motivation to develop Lattice DFT (LDFT) is that it provides a "sand-box" environment where one can study the subtleties of DFT itself, clarify the origin of inaccuracies in approximate functionals [1][2][3][4][5][6][7][8][9][10][11][12][13] and test ideas on new functionals. 11 Another motivation is provided by the problem of solving lattice models in the presence of an inhomogeneous potential. Lattice models are at the basis of our understanding of the phenomenology of strongly correlated, magnetic and disordered systems. Their fundamental relevance has in recent years motivated a number of successful experiments with ultracold atomic gases in optical lattices [34][35][36] fueling at the same time the development and refinement of efficient theoretical tools (see e.g. [37]) among which LDFT has become particularly useful. 21,38,39 Static and time-dependent [14][15][16] lattice DFT were indeed used to investigate the physics of Hubbard models with on-site interaction 13,[17][18][19][20][21][22] , Kondo models [23][24][25][26][27][28] , disordered interacting lattice models 31,32 and spin liquids 33 .
The Local Density Approximation (LDA) is the commonest approximation to the exchange-correlation (xc) functional of DFT, already proposed by Hohenberg and Kohn 40 , it was first applied to a lattice model by Gunnarsson and Schönhammer 3,7 and subsequently by Lima et al. 10 In the case of the one-dimensional Hubbard model there is the advantage that the exact solution of the homogeneous reference model is known by Bethe Ansatz 41 (BA). One can thus base the LDA on the exact energy or use the approximate but accurate analytical expressions avaialble 9,10,17 . One intriguing aspect of such Bethe Ansatz LDA (BALDA) functional is that, differently from the standard continuum LDA, it has a discontinuous functional derivative. Such non-analytic behavior stems directly from electron correlation and it has important consequences on the relation between the charge gap and the gap in the Kohn-Sham (KS) spectrum. 13,[42][43][44][45][46][47][48][49][50] Furthermore, as it will be explained below, it is ultimately responsible for an ultranonlocality of BALDA, which is absent in standard continuum LDA.
Derivative discontinuities and the associated cusp singularities in the exchange-correlation energy undermine the convergence of self-consistent KS equations. So far various approaches have been proposed to solve this problem for the BALDA functional including smoothening of the cusp minimum by going to a finite temperature 19 or by adding an ad hoc parameter 20,24 or by relying on Thomas-Fermi approximation 22 . Here we present a solution which allows us to treat in a simple and clean manner the cusp singularity of the BALDA functional at zero temperature: instead of the site occupation we use the LDA local chemical potential as a variable and we develop self-consistent equations. As an example, we apply our method, which we call µ-BALDA, to a Hubbard model subject to a staggered spin-independent site potential, also known as the ionic Hubbard model 10,51-53 (IHM). At half-filling, by modifying the ratio between Hubbard interaction and the staggered potential it can be tuned continuously from an ionic to a Mott insulating regime 51,54 . In both regimes we calculate the xc potential and the charge gap by applying µ-BALDA.
The paper is organized as follows. In Section II we give a description of the model and of the basics of Lattice DFT. In Section III we present a brief review of BALDA, we explain how the cusp problem emerges and we introduce the µ-BALDA algorithm. We then apply our method to the IHM in Section IV. As a proxy to a solid with a surface we study a bulk system with high binding energy attached to a zero binding energy chain representing the vacuum. 13 This geometry allows to have a Typeset by REVT E X well defined Kohn-Sham potential at all sites, even those with integer density, and to highlight an ultra-nonlocal behavior of LDA in the lattice, computing also the different contributions to the charge gap in the discontinuous situation. We conclude in Sec. V.

II. LATTICE DENSITY FUNCTIONAL THEORY
Let us start by outlining of the basics of LDFT. We consider a Hubbard chain in an inhomogeneous static field v x : where c † xσ creates an electron with spin σ =↑, ↓ at site x and n xσ = c † xσ c xσ , while U and t are respectively the interaction constant and the hopping amplitude. Notice that we have included an on-site contribution in the definition of the "kinetic energy" i.e. in the first term on the r.h.s. of Eq. (1). With this choice and "external potential", v x = 0, the lowest energy one-particle state has zero energy in analogy with the continuum model, which will be useful below to simulate a "vacuum" region (Sect. IV).
As in standard Kohn-Sham continuum DFT, 55 also in lattice DFT, the total ground-state energy can be obtained by minimizing a functional written as the sum of three density dependent terms where T KS is Kohn-Sham kinetic energy functional, E Hxc is the Hartree-exchange-correlation (Hxc) functional and with ρ x denoting the density at site x. Notice that the "functional" is actually a multivariable function of the on-site densities ρ x . The functional E Hxc is "universal" in that it does not depend on the external potential v x . Minimization of the functional with the constraint wheret denotes the hopping operatortϕ κ,x = −t(ϕ κ,x−1 + ϕ κ,x+1 − 2ϕ κ,x ) and ε N κ indicates the κ-th eigenvalue of the N -particle system and κ = kσ includes the orbital k and spin components. The Kohn-Sham potential is defined as, Even if not explicitly indicated v s x and v Hxc x depend on the number of particles N . In the lattice formulation the functional derivative with respect to the density becomes a partial derivative with respect to the on-site density 2,29 leading to the following definition of the Hxc potential, v Hxc The ground state density of N particles is composed of all occupied KS orbitals, Since the v Hxc x is a functional of the total density, Eqs. (3,6) have to be solved self-consistently.
Due to the constraint on the total number of particles the Hxc potential is defined up to a constant both in the continuum and in the lattice. One can extend DFT by considering ensemble densities 42,49,50,56 . In this case even the constant term in v Hxc x is determined. Intrinsic to the ensemble formulation of DFT are derivatives discontinuities of the exchange-correlation energy functional when the total density crosses an integer N . As first discussed by Perdew et al. 42 these lead to a discontinuous uniform change, ∆ xc , in the KS potential v s x when an integer filling is approached from the left or from the right, namely, where N ± = N ±η with η = 0 + . It can be shown that the charge gap of the system ∆ c ≡ E 0 (N + 1) + E 0 (N − 1) − 2E 0 (N ), with E 0 (N ) the N -particle ground state energy, satisfies, where ∆ KS = ε N N +1 −ε N N is the single particle Kohn-Sham gap.
Another important result is the DFT version of Koopman's theorem which is valid both in the continuum 43,57 and in the lattice 13 and it relates the highest occupied Kohn-Sham eigenvalue ε N N to the ionization energy, where One can also show 58 that if particles are bound in a finite region of space around the origin, v Hxc where e hom (U, ρ) is the energy per-site of the standard Hubbard model defined as, We do not include here the on-site term in the kinetic energy which does not affect Eq. (10) but it affects the zero of the homogeneous chemical potential defined below.
In the one-dimensional case the energy of the uniform system can be calculated exactly for all fillings by BA 41 and one can easily relate the appearance of a finite ∆ xc to the presence of a cusp singularity in the BA energy density, e hom (U, ρ), at ρ = 1. The physical consequences of this non-analytic behavior and the solutions of the related technical difficulties in the implementation of KS-DFT are the subject of the following sections. Most of the results presented below can be generalized to higher dimension using numerical solutions A.

Local-Density Approximation
Most difficulties in the implementation of BALDA are related to the fact that BA energy, e hom (U, ρ) has a cusp According to Ogata and Shiba 60 in the U/t → ∞ the system can be mapped to a spin-less fermion model. Therefore the chemical potential for 0 < n < 1 has a simple sine behavior. At n = 1 the chemical potential jumps by the Mott gap. In the picture we took U = 10 to have a finite gap. For n > 1 the chemical potential has again a sine form. For comparison we also plot the chemical potential of noninteracting electrons and the HF chemical potential.
at ρ = 1. The implications for BALDA functional are evident in Fig. 1(a) where we show two energy curves describing the typical dependence of BALDA energy on the occupation, ρ x , of a site: when the energy minimum with respect to ρ x is located away or at ρ x = 1, the functional has respectively a round or a cusp behavior at equilibrium. In the latter case, the Hxc potential is a discontinuous functional of the density and, as recently discussed in Ref. [19], the convergence of KS self-consistent cycle is not guaranteed. More precisely, when the density of a site is away from half-filling, ρ x = 1, the Hxc potential is obtained as usual as the derivative of the Hxc-energy [Eqs. (5), (10)], v Hxc where the local chemical potential,μ hom (U, ρ), coincides with the chemical potential of a homogeneous system with density, ρ and interaction U , On the contrary, for ρ x = 1, the standard definition given in Eq. (12) becomes ambiguous, indeed the derivative of e hom (U, ρ) has a jump whose amplitude equals the Mott gap ∆ hom Mott of the uniform system. Moreover, Eqs. (12) and (13) together with the definition of Eq. (11), imply that v Hxc andμ hom have the same limits when ρ → 1 ± , as shown in Fig.1(b). We can thus define: v ± Hxc =μ hom ± and we haveμ hom Fig. 1(c) shows v + Hxc , v − Hxc (full lines labeled 1 ± ) as a function of U/t.
Not only one value of the density (ρ x = 1) corresponds to a continuum of values of v Hxc but also certain values of v Hxc correspond to three values of the density as one can see in Fig. 1(d). This leads to numerical instabilities if one uses a standard approach to solve Kohn-Sham equations.
As mentioned in the introduction, various methods have been developed to treat these problems : some modify the BALDA functional relying on a sort of Thomas-Fermi approximation 20,24 , others entail the introduction of an ad hoc parameter 22 or a finite temperature to smoothen the discontinuity 19 . Here we propose a different route which allows to solve the KS equations selfconsistently without any additional parameter. Instead of ρ x , we take the chemical potential of the homogeneous system as a variable. Since, as shown in Fig. 1b, Eq. (13) is a monotonous growing function of ρ, it can be inverted.
We useμ x as the independent variable to avoid confusion with the functionμ hom thus the inverse of Eq. (13) is written as, Then we express the LDA Hxc potential in terms ofμ x , v Hxc KS equations with the Hxc potential given in Eq.(15) define the density as an implicit functional ofμ: where the square brackets abbreviate the multivariable functions, i.e. we set ρ At self-consistence Eqs. (14) and (16) should give the same result, we therefore have to impose Eqs. (17) can be solved iteratively starting from some initial guess for the chemical potentials µ (0) x . Here we expand F x to arrive at a set of linear recursive equations: (18) where we defined the vector µ = (µ 1 , ...µ L ).
For a reasonable initial setμ Eqs.(18) converge rapidly to the final solution while, in the presence of cusps, the conventional KS iteration scheme does not converge. We stress that our algorithm treats the cusp and normal sites in a unified way. One does not need to assume or guess beforehand which site will be a cusp. If the self-consistent solutionμ x falls inside (outside) the range [μ hom + ,μ hom − ], the site x is a cusp (normal) site (see Fig.4(a)). It is thus well-suited to describe systems where compressible and incompressible phases coexist. The gradient in Eq. (18) can be computed analytically using perturbation theory or numerically from finite differences. We found the latter option to be faster.
To understand the results discussed below concerning the inhomogeneous potential case, it is useful to discuss the behavior ofμ hom and v Hxc in the uniform case as a function of interaction strength and filling. At zero filling [curves labeled 0 in Fig. 1(c)] and independently of interaction strength it is easy to see that our definitions lead to v Hxc = 0 (solid blue curve) andμ hom = −2t (dashed line). For larger fillings at a given interaction strength, v Hxc increases due to the effective repulsion among particles. For small U , the Hxc potential is dominated by the Hartree part and it behaves as v Hxc ≈ U ρ x /2. This behavior is clear up to U/t ∼ 2 from the curves labeled 1 ± in Fig. 1c which represent v ± Hxc . For larger U/t the two curves clearly diverge. Such difference is actually present for any U/t (but exponentially small for small U/t) and represents the Mott gap. Furthermore, for large interaction, v − Hxc tends to a constant while v + Hxc increases linearly with U . The limiting value of v − Hxc can be understood from the fact that the charge sector of the uniform Hubbard model can be mapped for large U to a spinless fermion model with bandwidth 4t describing the lower Hubbard band for ρ x < 1. 60 This point is also elucidated by Fig. 2. There we show the filling dependence of the chemical potential of a Hubbard chain for large U/t, along with the Hartree-Fock and the U = 0 chemical potentials. We see that, as the filling increases, the "large U " chemical potential changes from −2t for zero filling to 2t for filling ρ = 1 − while the HF chemical potential increases linearly. Notice that the behavior of the exact BA chemical potentialμ hom for U = 6t shown in Fig.1b is well approximated by the "large U" chemical potential. We also see that v Hxc has the same limiting values as µ hom for ρ → 1 ± . For ρ x = 1 + both the chemical potential and v Hxc jump by and amount equal to the Hubbard gap which is of order ∆ Mott = U − 4t in this limit. The behavior for filling larger than 1 can be understood using particle hole symmetry.
The fact that v Hxc (ρ x < 1) saturates for large U (instead of having the naive mean-field behavior v Hxc = U ρ x /2) is typical of approaches where correlations are taken into account so electrons can avoid the large Coulomb cost. The present picture can be compared with similar results obtained using Gutzwiller approximation (see Fig. 1c of Ref. [61]). There the self-energy and the uniform chemical potential for filling ρ x < 1 saturate at large U and it jumps by the Mott gap on passing from ρ x = 1 − to ρ x = 1 + . One can interpret the limiting value of v Hxc (ρ x < 1) as the effect of a residual kinematic interaction between quasiparticles in the lower Hubbard band.

IV. APPLICATION TO THE IONIC HUBBARD MODEL
We now apply µ-BALDA to the study of the ionic Hubbard model (IHM), a Hubbard model with on-site interaction and a staggered potential.
In closed lattice models, the single particle potential that yields the ground state density of the interacting system is determined up to a constant since adding a constant to the potential leaves the charges invariant. In open systems which have a "vacuum" region one can fix the constant in such a way that the potential vanishes when both the density and the external potential vanish, far from the region where the particles are bound. In Ref. [13] we considered the open system shown in Fig.  3, with a vacuum region to completely determine the exact Hxc potential of the "bulk" Hubbard chain from a knowledge of the exact density.
In BALDA, in principle, the zero of the Hxc potential is determined by the explicit expression (12). However, as explained above, the potential is not determined by this expression when ρ x = 1 thus to remove any ambiguity on the potential we find it useful to consider the same geometry as in Ref. [13]. It will become clear below that in reality the vacuum region is not essential but it is enough to have just one site in the system in which ρ x = 1 to completely determine the potential. Still we use the geometry of Ref. [13] which is conceptually simple and analogous to the real situation in a solid. In addition this geometry (or any geometry with a confining potential) presents a challenge to conventional algorithms because, as we shall see, cusp sites coexist with non-cusp sites.
Specifically, we consider a Hubbard chain of L B sites with a large binding energy, called "the bulk", followed by a chain of L V sites with zero binding energy, termed "the vacuum", with open boundary conditions. The external potential v x includes both a step-like potential which accounts for the "work function" of the "solid" and the staggered field, where w 0 > 0 denotes the well-depth and V > 0 is the amplitude of the staggered field. Here we study the model around half-filling and we consider both the band-insulating (BI) and the Mott-insulating phase (MI) the latter appearing when the Hubbbard interaction, U , dominates over the staggered potential, U 2V (see e.g. Refs. [13,51,54]). A. Shift in the exchange correlation potential and ultranonlocality of the lattice local density approximation Figure 4 shows the Hxc potential (panels a and b) and the density (panels c and d) when the bulk chain is halffilled (N = L B ) and when one particle is added or subtracted respect to half-filling (N = L B ± 1).
At half-filling and for large on-site repulsion the bulk is in the Mott phase and the system becomes nearly incompressible. It therefore does not react to the staggered potential and all bulk sites have charge ρ x ∼ 1 (orange line in Fig. 4(c-d)). When this occurs the KS potential becomes nearly constant in the bulk and v Hxc acquires a staggered component behaving as V (−1) x which "screens" out the external potential, as shown in panels (a) and (b) of Fig. 4. This is a strong correlation effect captured by BALDA, as opposed to LDA in the continuum which would not be able to describe a similar situation in an heteroatomic chain of atoms.
What determines the value of v Hxc in the bulk with respect to the vacuum? On a closer look at the density in Fig. 4(d) one sees that a small amount of charge leaks to the vacuum from the site with x = 39, which is the site with small ionization energy closer to the surface. Having a density smaller than 1, this site is not affected by the cusp problem and it has a well-defined Hxc potential, v Hxc 39 ≈ v − Hxc ≈ 2t. Moreover, since the density has to be homogeneous and all the sites must have nearly the same effective potential, the Hxc potential of the other odd sites will be also very close to v − Hxc while the Hxc potential of the even sites has to satisfy the relation, v Hxc Notice that all these sites are cusp sites with ρ x = 1, so their Hxc potential is not determined by the usual local relation v Hxc = v Hxc (ρ x ) but by the ultranonlocal condition that all sites need to have the same density. Our numerical algorithm correctly converges to this solution.
For N = L B + 1 the role of odd and even sites is reversed. Indeed, as shown by the squares in Fig. 4(c) in this case, some of the even sites have a density slightly above 1 and, having their Hxc potential determined by Eq. (12), they play the role of reference sites while the odd sites have all unitary occupations and their potential is fixed by the even sites potential. In this case we have: v Hxc 2x = v + Hxc and v Hxc 2x+1 = v + Hxc − 2V . The case N = L B − 1 is very similar to the case N = L B , in this case, as explained above, some of the odd sites have a density slightly below 1 and they play the role of reference sites.
From the above discussion we conclude that the external staggered potential is screened in the even sites for N = L B , L B − 1 and in the odd sites when N = L B + 1. The net result is that v Hxc jumps by a constant quantity of order of ∆ xc ∼ ∆ hom Mott − 2V , with ∆ hom Mott ∼ U − 4t and where we used the fact that v + Hxc −v − Hxc = ∆ hom Mott . The relation with the charge gap will be discussed in Sec. IV B.
It is easy to compute the exact ionization energies in the atomic limit directly from Eq. (1) setting t = 0. In order to satisfy Koopmans theorem and in the case of a bulk uniform system (as found in the Mott regime) the exact KS potential in the bulk should satisfy v s x = −I. Therefore in this limit we can obtain the exact xc potential quite easily. Interestingly, the KS potentials obtained by µ-BALDA converge to the same results as shown next.
At half-filling ionization occurs from odd sites leading to I ≃ w 0 − V and v s ∼ −w 0 + V in agreement with the µ-BALDA results . In the case N = L B + 1 (shown with blue lines in Fig. 4) the added charge will reside mainly on the even sites (v x = −w 0 − V ) so that ionization will occur from these sites. In the atomic limit in the Mott regime one obtains I = w 0 + V − U . This leads to v s ∼ −w 0 − V + U which implies v Hxc 2x+1 ∼ U − 2V ∼ v + Hxc − 2V and v Hxc 2x ∼ v + Hxc again in agreement with the µ-BALDA results.
Let us conclude the discussion of Fig. 4 with a technical remark. In the exact many-body solution only the uniform compressibility is zero but the compressibility is small but finite at finite momentum.Therefore, differently from what predicted by BALDA, in the exact solution, for all V = 0 the charge imbalance between odd and even sites is non-vanishing. For large U/t, however, the deviations from the BALDA ρ x = 1 solution are expected to be small, the compressibility being trivially zero at all momenta in the atomic limit. Figure 5 shows how the xc potential of odd and even sites evolves from the band insulating regime at small interaction to the Mott insulating regime at large interaction. In the band insulator regime the xc potential tends to screen the external potential but it lies outside the interval [v − Hxc , v + Hxc ] marked by the dot-dashed red lines. In this regime the occupation of the sites is larger or smaller than 1. As the interaction increases the system becomes less compressible and the charges tend to approach one. When the xc potential of the different sites (blue and orange curves) hit the dot-dashed red lines, (or more precisely the µ x fall inside [μ hom + ,μ hom − ]), the system enters the Mott insulating regime. The density becomes homogeneous and the difference between the odd and even site potential equals 2V . Depending on whether the reference sites have a density smaller or larger than 1, even sites align to v − Hxc (a) or odd sites align to v + Hxc (b). We stress again that the above results shows how, within BALDA, strong correlation leads to a ultranonlocality of the Hxc potential. The potential of the whole system is indeed fixed by the presence of just one or few sites whose density is slightly below or above half-filling.
As mentioned above, it is enough to have a single site in the system with density different from 1 to obtain the absolute value of the xc potential in all the sites. Leakage to the vacuum is a natural way to obtain such a reference site but an impurity would work as well. For a large periodic closed system the absolute value of the xc potential can be determined by a limiting procedure considering small deviations from half filling from below or from above.

B. Charge Gap
As stated in Eq. (8) the fundamental charge gap, ∆ c can be written as the sum of two terms, the KS gap, ∆ s , and a contribution coming from the discontinuity of v xc upon the addition of an infinitesimal amount of charge to the system [Eq. (7)]. The latter is site-independent and it can be calculated, in some simple cases even analytically 13,62 . A rigorous definition of the discontinuity in v xc requires using ensemble DFT 49,50,56 . Alternatively ∆ xc can be estimated employing a formula first derived by Sham and Schlüter in Ref. [44], based on finite differences, which reads (19) where ρ N N +1 is the N + 1'th KS orbital density for N particles and v xc x (N ) denotes the xc potential of the Nparticle system at site x. The difference between v xc (N ) and v xc (N + 1) on the r.h.s. of the above equation is in general site-dependent due to the fact that the charge added is finite.
The estimate of Eq. (19), which we adopt here, converges to the exact result of Eq. (7) when the local change in the density upon addition of an electron becomes infinitesimal small, as in extended systems in the thermodynamic limit. Eventually we remark that a simple expression of the error ∆ N,N +1 xc −∆ xc can be obtained along the lines of the Supplementary Material of Ref. [13] .
We now come to the discussion of µ-BALDA results for the charge gap. As mentioned in the Introduction, the ionic Hubbard model displays a transition between an ionic and a Mott insulating regime. As thoroughly discussed in a number of papers (see e.g. Refs. [54] and References therein), the Mott regime appears when the Hubbard interaction dominates over the staggered potential and the charge gap has a minimum at the transition point.
In Fig. 6 we plot the charge gap (obtained as total energy differences) as a function of V (a) and U (b) obtained with µ-BALDA and compared to practically exact numerical results.
In the Mott regime the charge gap equals the xc discontinuity and the µ-BALDA KS gap ∆ KS vanishes as the density becomes homogeneous. On the contrary, in the ionic regime the charge gap approximately coincides with the KS gap while the xc discontinuity becomes vanishingly small, since in this regime there are no half-filled sites.
Taking into account that the µ-BALDA solution is homogeneous in the Mott phase we can give an analytical expression for the behavior of the µ-BALDA charge gap in the Mott regime. In this case the charge gap in Eq. (8) is exhausted by the discontinuity in the potential obtained in previous Section, where the last expression is valid for large U/t. In the large U limit the critical value of the potential is given by, Clearly the transition from the Mott regime to the band insulating regime occurs when V ≈ 3.0t in µ-BALDA and V ≈ 4.3t in Lanczos.
These expressions explain the linear behavior observed in the explicit solution of µ-BALDA in the Mott regime [small V in panel (a) and large U in panel (b)] and also the critical values V c ≈ 3t in Fig. 6(a) and U c = 2V +4t ≈ 6t. The latter result overestimates the exact critical U on panel (b) which we attribute to the inaccuracy of the large U approximation at the critical value.
In general BALDA underestimates the charge gap in the Mott regime [small V (a) and large U in (b)] and it fails to describe the non-linear dependence of the gap for small V (a). These failures can be related to the fact that in the Mott regime BALDA tends to generate a ground state which is more homogeneous than the true one as we remarked in the previous Section. The small inhomogeneity in the exact density will be associated to a small Kohn-Sham gap in the exact Kohn-Sham spectrum which will tend to diminish the discrepancy. For large V in the band insulating regime the µ-BALDA charge gap rapidly converges to the exact result which approximately coincides with the KS gap. This shows that in the band insulating regime, despite the strong interaction, BALDA works remarkably well.
It is interesting to notice that Eqs. (20), (21) converges to the exact result in the atomic limit t = 0. In this case the charge gap is given by ∆ c ≃ |U − 2V | and the transition occurs at U ≃ 2V .

V. CONCLUSIONS AND DISCUSSIONS
Probably one of the most popular LDFT approaches called BALDA has till now required ad hoc regularizations of the xc energy to describe the Mott phase. In the present work, we have solved this long-standing problem, developing a new method to find the exchange-correlation potential of a lattice system in the Mott regime with a fully self-consistent procedure. Differently from previous works, our algorithm, which we call µ-BALDA uses the local chemical potentials as variables. As an example we apply the method to the study of the transition between Mott and band insulating regimes of the ionic Hubbard model. Beside the general methodological progress, we obtain several results: (i) we have shown that in the Mott regime the external potential is completely screened by the xc potential (ii) we highlight an ultra-nonlocality of LDA in the lattice, i.e. we show that in the presence of the discontinuity, one or few "reference sites" in the system are capable of fixing the whole xc potential; (iii) we calculate separately the different contribution to the gap in the different regimes providing, in particular, an analytical understanding of the behavior of the µ-BALDA charge gap in the Mott regime. iv) In the Mott phase, due to correlation, ionization occurs from different sites for a system slightly below or above half-filling. This gives a contribution to the ionization energy which can not be captured in a single particle picture but which is captured by µ-BALDA. Beside these nice qualitative features we have discussed the quantitative errors of BALDA in the Mott regime and the high-accuracy in the band insulating regime even when the interaction is not small.
In this paper we extract the energy of the homogeneous system from exact BA 41 , but the method can be clearly also applied when the approximate analytical parameterization introduced by Lima et al. 9,10 is used. Our results are therefore relevant for all 1D lattice models where BALDA is applied, including Kondo systems 23,25-28 , dynamical Coulomb blockade treated with time-dependent DFT 24 , harmonically trapped Hubbard electrons 19,22,33 and spinless Fermions with neighboring interaction 18 . µ-BALDA is actually a rather general approach and it could be in principle generalized also to higher dimension or to treat other discontinuous functionals.