Low-and high-energy localization landscapes for tight-binding Hamiltonians in two-dimensional lattices

Localization of electronic wave functions in modern two-dimensional (2D) materials such as graphene can impact drastically their transport and magnetic properties. The recent localization landscape (LL) theory has brought many tools and theoretical results to understand such localization phenomena in the continuous setting, but with very few extensions so far to the discrete realm or to tight-binding Hamiltonians. In this paper, we show how this approach can be extended to almost all known 2D lattices and propose a systematic way of designing LL even for higher dimensions. We demonstrate in detail how this LL theory works and predicts accurately not only the locations, but also the energies of localized eigenfunctions in the low-and high-energy regimes for the honeycomb and hexagonal lattices, making it a highly promising tool for investigating the role of disorder in these materials.


I. INTRODUCTION
The concept of Anderson localization, initially introduced in tight-binding models in the context of condensed matter physics [1], has been applied in the continuous setting to all types of waves-quantum [2], classical [3][4][5][6][7], or even gravitational [8]. In this setting, the recent theory of the localization landscape (LL) [9] has brought new insights and methods to address the wave localization properties in systems such as gases of ultracold atoms [10], disordered semiconductor alloys [11], or enzymes [12], and has been successfully extended to Dirac fermions [13] and nonscalar field theory [14]. In this work we show that the whole machinery of the LL can be generalized to tight-binding systems for most known one-dimensional (1D) and 2D lattices, allowing us to broaden the range of predictive scope of this fruitful approach to a large class of new 2D materials [15] (see also [16] and the related collection of papers).

II. TIGHT-BINDING AND LOCALIZATION LANDSCAPES
Tight-binding models are commonly used to study perfect [17] as well as disordered lattices [18][19][20]. The tightbinding HamiltonianĤ with on-site disorder and nearestneighbor coupling is defined as where ψ ≡ (ψ n ) n∈[ [1,N]] is the wave function defined on the sites of the lattice (numbered from 1 to N here), V n is the on-site potential at site n, t is the coupling constant between neighboring sites (assumed to be constant here), n indicates the ensemble of nearest neighbors of site n, and b n is its cardinal. In the following, one will assume that t has value 1, thus setting the energy unit, and that V n = W ν n , where ν n is an i.i.d. random variable with uniform law in [−0.5, 0.5], with W being therefore the disorder strength for a given lattice. Let us start our study by presenting results on the honeycomb lattice. Figures 1(a) and 1(e) show this lattice and its celebrated dispersion relation in the tight-binding approximation, respectively [21]. We solve the Schrödinger equation for the Hamiltonian defined in Eq. (1) on the honeycomb lattice, with the on-site potential depicted in Fig. 1(b). In Fig. 1(c) are displayed the first four eigenstates which, as expected, exhibit a finite spatial extension typical of Anderson-localized modes. At the other end of the spectrum, Fig. 1(g) illustrates a feature that has no continuous counterpart: the existence of high-energy localized modes (the last four eigenstates are displayed in the example). This phenomenon is well known for instance in the case of 3D Anderson localization on a cubic lattice at low disorder strength, in which the spectrum of the Hamiltonian is symmetric in the range [−6 − W/2; 6 + W/2] and exhibits a transition (the mobility edge) between localized and delocalized states at both ends [22].
In the following, we show how to build the two discrete localization landscapes displayed in Figs. 1(d) and 1(i). These landscapes allow us to accurately predict the location of the localized modes near the two band edges (low and high energy), as well as their energies, without solving Eq. (1). We then generalize this method to a wide class of 2D lattices.
Let us first summarize the salient features of the LL theory in the continuous setting. For any positive definite Hamiltonian H in such a setting, the localization landscape u is defined as the solution to One of the main results of the LL theory is that the quantity 1/u-which has the dimension of an energy-acts as an effective potential confining in its wells the localized states at low energy [23]. Moreover, the energy of the local fundamental state inside each well was found to be almost proportional to the value of the potential 1/u at its minimum inside the well [24], where d is the embedding dimension of the system. 1/u can be understood as a renormalization of the disorder at the selfadapted local scale.
In the tight-binding setting, Lyra et al. [25] have studied a 1D linear chain with nearest-neighbor coupling and have shown that the positions of the localized modes are given by two different localization landscapes. The low-energy LL is obtained by solving the analog of Eq. (2) in the discrete setting, i.e.,Ĥu = 1 [u ≡ (u n ) n∈[ [1,N]] ; 1 is a vector filled with 1] with the same boundary conditions as the eigenvalue problem. Another LL, called the dual localization landscape (DLL), gives the position of the envelope of the highly oscillating, high-energy, wave functions. More recently, Wang and Zhang [26] have proved mathematically that, in any dimension, the reciprocals of these discrete LL and DLL act indeed as effective confining potentials in a tight-binding system at both low-and high-energy regimes, respectively. Figure 1(d) shows the reciprocal of the LL, 1/u ≡ (1/u n ) n∈[ [1,N]] , computed on the honeycomb lattice with the on-site disorder depicted in Fig. 1(b). Note that a shift H →Ĥ + V shift has been performed in (1) to ensure a positive definite Hamiltonian; see Appendix A. As already observed for continuous systems, the role of the effective confining potential played by 1/u is revealed through its basins, labeled following their depth in Fig. 1(d). Indeed, one can observe the correspondence between the deepest wells of 1/u and the positions of the first eigenmodes plotted in Fig. 1(c). As analyzed by Arnold et al. [24] in the continuous setting, two almostequal eigenvalues can lead to a different ordering in the values of the minima of 1/u, thus inducing a mismatch in the correspondence. This effect, which does not affect the ability of the LL to predict the position of localized modes, is visible in Figs. 1(c) and 1(d) with the first and fourth eigenstates and minima and has been analyzed in detail in Appendix B.
Finally, we have quantitatively tested that, regardless of the lattice, the tight-binding LL efficiently pinpoints the localized modes (see Appendix C).

III. HIGH-ENERGY LANDSCAPE FOR SYMMETRICAL DOS
On the other hand, the high-energy situation is much more intricate. The symmetry of the honeycomb lattice allows us a straightforward derivation of the landscape governing the high-energy localized states, namely the DLL. Indeed, the tight-binding Hamiltonian (1) can be decomposed intô H =Ĥ 0 +V, whereĤ 0 stands for the uniform honeycomb lattice with zero on-site energy andV accounts for the disordered on-site potential. The unperturbed part of the Hamiltonian displays the usual chiral symmetry for a bipartite lattice, zĤ0 z = −Ĥ 0 , where the Pauli-like matrix z acts on the sublattice degree of freedom: it keeps the amplitudes on the A sites fixed but inverts those on the B sites ( z = P A − P B , which is the difference between the respective projectors on the two sublattices). Due to the diagonal nature of the disordered potential, the complete Hamiltonian obeys the symmetry z (Ĥ 0 +V ) z = −(Ĥ 0 −V ). The latter property is exemplified in Figs. 1(e) and 1(f): unlike the DOS of the uniform lattice, the DOS of a given realization of the disordered system is not symmetric with respect to the origin, but the DOS obtained by inverting the sign of all on-site energies is the exact symmetric of the original situation.
Let us call φ ≡ (φ n ) n∈[ [1,N]] the eigenstates of the inverted Hamiltonian ordered by increasing eigenvalues. The lowenergy states of the inverted Hamiltonian now correspond to the high-energy states of the original Hamiltonian through φ = z ψ. Since the high-energy eigenstates oscillate with a period equal to the nearest-neighbor distance, the new low-energy states appear as "demodulated" versions of their high-energy counterparts; see Figs. 1(g) and 1(h). We can now therefore use the localization landscape for the inverted system, but with u now being the solution toĤ u = 1 with In the example of Fig. 1(i), one can clearly see how the deepest wells (which are different from the low-energy wells) pinpoint the locations of the localized states. Beyond the honeycomb lattice, this spectrum inversion strategy can be deployed for other lattices with symmetric band structure, as the 1D dimer chain or the 2D Lieb lattice (see Appendix C, Table II).
As mentioned in the Introduction, the localization landscape also provides accurate estimates of the localized eigenvalues in the continuous setting [24]. However, the generalization of the simple Eq. (3) to tight-binding Hamiltonians has never been studied beyond the simple d-cubic lattice, nor its extension to the higher part of the spectrum. We plot on Fig. 2(a) [2(b)] the lowest (highest) eigenvalues of Eq. (1) versus the local minimum values of the effective potential 1/u (1/u ) at the position of localized eigenstates for a honeycomb lattice with N = 2135 sites and for a given disorder W = 3. Each scatter plot corresponds to 100 realizations of the disordered potential. In both cases, a direct proportionality is clearly observed for the lowest part of the plots, with Pearson coefficients of the linear regression close to 0.99. A general study of the quality of the proportionality is presented in Appendix B. Note also that in order to obtain this proportionality (which is more than a simple linear dependency), one has to choose V shift so that the shifted potential has a minimum value close to zero (see Appendix A).
These observations indicate that the discrete low-energy localization landscape performs as well as its continuous analog in predicting energy and spatial distribution of localized modes without resolving an eigenvalue equation. Moreover, the high-energy DLL also exhibits the same properties. For both ranges of energy, LL and DLL provide a good estimate of the integrated density of states, as shown in Appendix D. All these results are not restricted to the honeycomb lattice. Figures 2(c) and 2(d) show that the proportionality is obtained for a large variety of "canonical" lattices (1D: chain, dimer chain, chain with second-neighbor coupling; 2D: square, honeycomb, Lieb, hexagonal, Kagome, tts [27]) and in a wide range of strength disorder. Note that to quantify the latter unequivocally for different lattices, the parameter W is not the best suited. Indeed, for a given value of W , the relative weight of the potential term in (1) compared to the kinetic term depends on the connectivity of the discrete Laplacian m (ψ m − ψ n ). The number of edges of the graph on which this operator is relying is given by the number b n of nearest-neighbor couplings, which itself depends on the elementary motif of each given lattice. Therefore, we use a less contingent quantity, namely the inverse participation ratio (IPR) defined for a given eigenvector ψ ( j) = n ψ ( j) n |n as IPR j = n |ψ ( j) n | 4 /( n |ψ ( j) n | 2 ) 2 . More precisely, in Fig. 2, we consider the first (c) and last (d) 3% of the eigenstates to compute the slopes that are plotted versus the mean IPR corresponding to the same 3% of the eigenstates.
With one noticeable exception for the Lieb lattice, the values of the slope s appear to evolve continuously between s = 1 + d/4 and s = 1, both for the lowest and the highest eigenvalues [see Figs. 2(c) and 2(d)]. Moreover, all the curves bunch into two smooth master curves, one for each space dimension. In the weak disorder limit, i.e., IPR → 0, that is to say when the influence of the disordered potential is small compared to the Laplacian term, one can reasonably expect that the continuous result of Eq. (3) still holds for both the lowest and the highest part of the spectrum. This is indeed observed: the slopes fall on the (1 + d/4) limit for IPR → 0. In the other limit case, when the disorder is so strong that an eigenstate is localized on a single site ( IPR → 1), the LL is essentially supported locally on the same site. This means that the eigenstate and the LL are locally proportional, a feature already observed in the continuous setting. Equation (2) at the only site n supporting the wave function therefore becomes Hu n = 1 ≈ Eu n ; hence E ≈ 1/u n and a slope s 1 is expected.
The definition and the properties of the low-energy LL are valid for any lattice, in any dimension, and are not restricted to nearest-neighbor coupling. We simulated thoroughly many different "canonical" lattices, for which details are provided in Appendix C. The construction of the high-energy DLL, however, used explicitly in our case the chiral symmetry of the honeycomb lattice and hence the central symmetry of its DOS. We will show in the last part of the paper that when this property cannot be used, there remains a general procedure which consists in "demodulating" the Schrödinger equation around a local maximum of the dispersion relation of the Laplacian. This procedure can be applied to any lattice, leading to a DLL that relies on the specific band structure of the given lattice.

IV. HIGH-ENERGY LANDSCAPE FOR NONSYMMETRICAL DOS
To illustrate this, we now focus on the hexagonal lattice [ Fig. 3(a)] whose DOS is not symmetric [Fig. 3(c)]. In Fig. 3(d), we display the states corresponding to the four highest eigenvalues. Similar to what was observed for the honeycomb lattice [ Fig. 1(g)], we first note that the high-energy eigenstates are spatially oscillating with a period equal to the nearest-neighbor distance. Without on-site energy disorder, the highest eigenvalues are located at the vertices K of the first Brillouin zone (BZ) [Figs. 3(b) and 3(c)]. The correspond- To remove the rapidly oscillating part of the highenergy eigenstates, we define the envelope φ of an eigenmode ψ whose wave vector is close to k max by ψ n = e i k max · r n φ n , where r n is the position of site n (see Appendix E for the derivation of the landscape around a local maximum of the dispersion relation). By injecting φ in Eq. (1), we obtain a "demodulated" equation that reads for k max = 4π where, in addition, the band structure has been inverted by changing the signs of both the couplings and the on-site energies, and where the notation n + a denotes the site reached from site n by a translation of a vector a [see Fig. 3 We then compute the landscape u associated to this Hamiltonian, i.e.,Ĥ 4π 3ax u = 1, and obtain a complex confining potential 1/u whose absolute value is plotted in Fig. 3(e). The comparison between the deepest wells of 1/|u | ≡ (1/|u n |) n∈[ [1,N]] , see Fig. 3(d), and the locations of the aforementioned localized high-energy states clearly shows here again a direct match between the two sets (see Figs. 4 and 5 in Appendix B). Moreover, the proportionality between the minima of 1/|u | in the basins and the actual energies still holds (see Fig. 10 in Appendix B). We have performed extensive simulations on eight types of lattices (three 1D and five 2D) for various disorder strengths, all exhibiting Pearson coefficients larger than 0.96 and even larger than 0.98 in most cases.

V. CONCLUSION
Born a decade ago, the localization landscape theory has proven its remarkable efficiency to bring in a more accessible form the information contains in a Hamiltonian [28]. In this work, we have extended its scope to discrete systems described by a tight-binding Hamiltonian, with a focus on 2D lattices. The low-energy part of the spectrum is described by a discrete extension of the effective confining potential defined for continuous systems. It bears the same efficiency as its continuous counterpart in predicting the localization regions and the corresponding energies and hence the density of states [29]. More challenging is the construction of the dual confining potential that acts on the upper part of the spectrum. When the lattice is equipped with chiral symmetry, like the honeycomb lattice, the high-energy theory is directly deduced from the low one. When this symmetry is not present, we have proposed a general procedure to build the dual localization landscape. Our method is efficient, robust, and very general but not yet completely universal. It has yet to be extended to situations like the one encountered with the kagome lattice. In this case, the DOS is not symmetric and the high-energy states lie on a flat band: the definition of k max remains a challenge. Finally, it is also worth noting that interesting properties appear in the center of the band for many of the new 2D materials. Although this energy range does not fall directly into the frame of our approach, we were able to pinpoint the localization regions in this case thanks to an approach inspired by the L 2 -landscape method [14]. This last one is able to give information about the localization subregions, but predictions about the eigenenergies, IDOS, etc. are still missing despite its efficient numerical implementation [30]. Predicting the eigenenergies and the IDOS will be the topic of future investigations.

APPENDIX A: ENERGY SHIFT
For the simulations presented in the paper, we considered the Hamiltonian given by Eq. (1) in the main body of the paper, but written in a more concise way: where V n is the on-site potential at site n and t is the coupling constant between neighboring sites. Additionally, V n = W ν n , where ν n is an i.i.d. random variable with uniform law in [−0.5, 0.5], with W being therefore the disorder strength for a given lattice and V shift the energy shift that avoids negative eigenvalues. Finally V shift = − min(E ) + W/2, where E is the TABLE I. Smallest and largest eigenvalues of the tight-binding Hamiltonian without on-site potential and with V shift = 0. The couplings are all t = 1, except for the 1D dimer chain and the 1D chain with second neighbor coupling cases where they are explicitly written.
The considered values are shown in Table I.
To compute the dual landscape, we use V shift = max(E ) + W/2. Note that for the 1D chain with second neighbors coupling the energy is where the corresponding k max are the solutions of Then, values shown in Table I correspond to the particular case t 2 = t 1 / √ 8.

APPENDIX B: ACCURACY OF THE PREDICTIONS
In this section, we quantify the quality of the localization prediction of the localized states in a situation of small (Fig. 4) and large disorder (Fig. 5). This is done by calculating the distances between the position of the maximum of an eigenstate and the positions of all the wells of the effective potential. The wells are ordered by their depths and we then find the rank of the well corresponding to the minimum distance and report it in the 2D histograms of Figs Similarly, we quantify energy predictions for all the methods presented in the main text. In Fig. 6, we plot the Pearson coefficient of the linear regression as a function of the number of minima taken to compute the slope for the honeycomb lattice and in Fig. 7 the same quantity for the hexagonal lattice. We see that, whatever the disorder, the maximum correlation is obtained when we chose the 3% states with the lowest (or highest) eigenvalues. In Fig. 8, we plot a 2D histogram showing the distribution of the eigenstates in the eigenvalues-IPR plane. We clearly see the more localized states for higher disorder and the asymmetry of the density of states for the hexagonal lattices.  Figure 9 shows that, in every case considered in the paper, the eigenvalues and the minima of the effective potential are highly correlated (Pearson coefficient very close to 1).
Finally, Fig. 10 is similar to main text Fig. 2(d), except that the calculations are done using the explicit demodulation. This explains why there are more cases in Fig. 10

APPENDIX C: EXHAUSTIVE STUDY OF VARIOUS 1D AND 2D LATTICES AND SIMULATIONS DETAILS
In Table II, we list the different lattices studied in this study with their properties. The values of k max shown here for the 1D chain with second neighbor coupling correspond to the case t 2 = t 1 / √ 8. For each lattice, we computed 100 different configurations; see Table III for details. In the case of symmetric DOS, we calculated the landscape prediction FIG. 10. Proportionality factor between V shift − E and the minima of 1/|u | for the different 1D (squares) and 2D (circles) lattices studied. The linear fits are done on states with the 3% highest energy. The dashed horizontal lines show the limit expected using the approximate form Eq. (3). Each symbol corresponds to a disorder strength. The results are plotted as a function of the mean IPR calculated over the 3% states found for a given disorder strength.
both using the symmetry of the DOS and using the explicit demodulation. In Figs. 2(c) and 2(d) of the main text, 3% of the eigenstates correspond to 64 states, and because we consider 100 configurations, the slopes are calculated using 6400 points. Figure 11 shows the counting function for the hexagonal and honeycomb lattices with a strong on-site disorder (W = 240 and W = 480 for the honeycomb and hexagonal lattices, respectively). Blue lines correspond to the solution of Eq. (1), while orange lines are the minima of both confining and dual confining potentials.  Let us assume that k max is a wave vector at which the dispersion relation E ( k) of the Hamiltonian without potential (i.e., minus the Laplacian) exhibits a local maximum. Then one can write locally the dispersion relation as

APPENDIX D: IDOS ESTIMATE
where A is a definite positive 2 × 2 matrix whose eigenvalues are the inverse of the effective masses in both directions. One can then write any eigenfunction ψ of the full Hamiltonian (i.e., with potential V = W ν) as where φ is an envelope function satisfying the following equation: with E being the energy of ψ. The local maximum of the dispersion relation E ( k) at k max implies that m∈ n r m e i k max · r m = 0, where the sum is taken over all interacting neighbors of one site of the lattice assumed to be located at r = 0. One can thus define a new HamiltonianĤ The energy of a plane wave φ = exp( j k · r) for this Hamilto-nianĤ k max is therefore One can therefore apply the localization landscape formalism to this Hamiltonian shifted by a quantity W/2.