Decay rate of real space delocalization measures : a comparison between analytical and test systems

Please note that technical editing may introduce minor changes to the text and/or graphics, which may alter content. The journal’s standard Terms & Conditions and the Ethical guidelines still apply. In no event shall the Royal Society of Chemistry be held responsible for any errors or omissions in this Accepted Manuscript or any consequences arising from the use of any information it contains. Accepted Manuscript


Introduction
Real space theories of the chemical bond 1,2 have provided a physically sound alternative to the molecular orbital (MO) paradigm over the last two decades, 3 incorporating orbital invariant descriptors endowed with chemical meaning to the chemical bonding toolbox.Among several proposals, commonly gathered together under the Quantum Chemical Topology (QCT) umbrella, the Quantum Theory of Atoms in Molecules (QTAIM) proposed by Bader and coworkers stands out by its own. 1 This theory is commonly known as the topological analysis of the electron density, ρ, since the real space is partitioned using the topology induced by ρ and many of its insights are obtained from the density itself or from its succesive derivatives (gradient, laplacian. . .).However, the QTAIM can also take profit of other reduced density matrices (RDMs).
A key distinguishing feature that separates the QTAIM from other techniques is its energetic face.All the standard components of the Coulomb Hamiltonian may be examined over QTAIM real space domains, providing a unique route to bind the physicist and the chemist points of view.In order to do so, the non-diagonal first order RDM (1RDM) as well as the diagonal 2RDM are needed.The first is a standard ingredient of the orthodox QTAIM, and the second was shown to provide a measure of electron delocalization in a seminal paper by Bader and Stephens 4 .The 2RDM was finally added to the energyrelated descriptors of the QTAIM in the interacting quantum atoms (IQA) scheme. 5partamento de Química Física y Analítica.Facultad de Química.Universidad de Oviedo.33006-Oviedo.Spain.Fax: 34 985 103125; Tel: 34 985 103037; E.mail: ampendas@uniovi.esAs the set of systems for which the QTAIM was applied increased (including molecules, clusters, and solids), it soon became clear that ρ and/or ∇ 2 ρ contain a wealth of information about chemical bonding, and that a simple classification of systems into shared-shell or closed-shell types, in quite good agreement with the standard covalent/non-covalent (including ionic) bonding models was possible.This knowledge is now mainstream, reaching general chemistry textbooks, and its success encouraged researchers to look for features in ρ peculiar to metallic systems, an enterprise with little initial success.The density of conducting materials seemed not different from that of standard covalent ones.A proposal that nonnuclear maxima (or non-nuclear attractors, NNAs), known in the Li 2 molecule since 1956, 6 and found in other lithium clusters 7,8 might signal metallic behavior was received with hope.Its chemical image matched well with the qualitative idea that conducting electrons were transferred to interstitial positions in crystalline lattices.However, when reliable calculations of the topology of the electron density in solid alkali metals were available, this initial hope vanished. 9Only Li displayed NNAs at equilibrium geometries.Since then, experimental densities have demonstrated that NNAs may appear in Beryllium, 10 but are absent in the metals of largest conductivity, like Cu, Al, or Ag.Whatever the origin of NNAs, 11 they do not signal conductivity.Other attempts have shown that the density of metallic systems is characterized by its interstitial flatness, 12 again in agreement with conventional chemical wisdom, but no salient feature of the density determining conductivity has ever been found.
From the purely theoretical side, the naïve difference between electrical conductors and insulators lies in their exci-tation spectra, 13 i.e. in the nature of their excited states, far from the real space realm.However, a seminal paper by Kohn in 1964 14 showed that insulators and conductors also differ essentially in the organization of electrons in their ground state.In the former, the wave function is composed of manybody building blocks localized in disconnected regions of the many-particle configuration space.Long forgotten, Kohn's theory was reformulated by Resta in 1998, 15 in what today is known as the modern theory of polarization, offering a new view deeply linked to Berry phases. 16Resta has shown that the finiteness or divergence of Kohn's localization tensor (LT) is the key to conductivity, and that for one-determinant descriptions, the LT is closely related to Boys theory of localization, 17 very familiar to the quantum chemical audience.Application of the LT (or total position spread tensor, TPS) to chemical problems has been pioneered by Evangelisti and coworkers. 18An important point for what follows here is that Resta's formulation lies in real space.Thus, conductivity may not leave scars in the plain density, but should lead to a recognizable imprint if we examine other RDMs in real space.
Notwithstanding the role that the TPS should play in this important problem, in this contribution we will focus on how the standard real space measures of (de)localization may be related to metallic behavior in molecular systems.To that end, we will start recalling some known results based on Kohn's nearsightedness principle. 19They show that the decay behavior of the 1RDM, ρ(r r r; r r r ′ ), determines the locality of all relevant observables.We will then relate the 1RDM to the real space delocalization index (DI or δ) defined within the QTAIM by Bader and Stephens, 4 δ A,B .Armed with this, we will examine the analytical decay rate of δ for Hückel and tight binding (TB) models of metals and insulators, comparing with simple calculations in toy systems.Our results will show that, as expected, the decay rate of delocalization measures differs in insulating-or metallic-like systems, being exponential in the former and algebraic in the latter.Other interesting links, like that between the well known oscillations of δ in conjugated molecules, clearly related to resonance and chemical behavior, and Friedel oscillations in metals will be put forward.
The paper is organized as follows.Section 2 will be devoted to the relation between 1RDMs and delocalization indices.Then we will briefly present the models and computational details we have used.Sections 4 and 5 will examine finite and extended analytical models, while we will comment on single-determinant results on simple systems in Section 6.We will end with some prospects and conclusions.

Decay rate of 1RDMs and delocalization indices
As commented, the decay rate of the density matrix, a fundamental issue after Kohn's insights on nearsightedness 19 has been widely studied in the physical literature, usually under a one-electron picture within density functional theory (DFT), or within a tight binding (TB) Hamiltonian approximation.For instance, Goedecker 20 showed that assuming the electronic structure of an isotropic 3D metal to be dominated by its free-electron band structure, the 1RDM ρ(r r r; r r r ′ ) = ρ(|r r r − r r r ′ |) = ρ(s s s) decays algebraically at zero temperature, where k F is the Fermi vector modulus, related to the valence electron density k 3 F /(3π 2 ) = N el /V .As it can be seen, ρ(s) oscillates on decaying like s −2 , with zeros at sk F ≈ 4.49, 7.73, 10.90, etc.For reasonable valence density values, it can readily be found that these zeros are close to lattice vectors.As we will see, this oscillatory behavior, which is closely related to the well known Friedel oscillations of metals, 13 has close relatives in finite molecules.
Taraskin and coworkers 21 have refined these results for 1D to 3D TB metals in simple linear, square, or cubic cells, showing that ρ decays as s −(d+1)/2 , d being the dimensionality of the system.These authors 22 have also shown that, for two bands TB models of insulating lattices the 1RDM falls exponentially with s, where the inverse decay length λ depends on the gap, ∆, scaling linearly with it as ∆ → 0. ρ(s) turns out to be anisotropic, showing its slowest decay along the (1,1) or (1,1,1) diagonals in 2D or 3D, respectively.Effective λ values have been shown to lie between 1-5.Once these results have been presented, we turn to delocalization measures in real space.Several indicators have been proposed over the years, among which the electron localization function (ELF) of Becke and Edgecombe, 23 very popular in theoretical chemistry after the work of Savin and Silvi, 24 is probably best known.Other possibilities like the electron localizability indicator (ELI) introduced by Kohout, 25 valid for correlated descriptions, also exist.All of these are local descriptors, bearing no decay information, and do not serve our purposes.
Fortunately, the DI introduced by Bader and Stephens 4 within the QTAIM describes how many pairs of electrons are shared (thus delocalized) between two finite regions A and B in real space: Page 2 of 10 Physical Chemistry Chemical Physics

Physical Chemistry Chemical Physics Accepted Manuscript
In this expression, the integrand is the standard exchangecorrelation density, ρ xc (r r r 1 , r r r 2 ) = ρ(r r r 1 )ρ(r r r 2 ) − ρ 2 (r r r 1 , r r r 2 ).δ A,B is a scalar parameter between any two regions (notice that if A = B we usually talk about a localization index) that due to the extensivity of ρ xc , adds to the total electron population (N), 1/2 ∑ A,B δ A,B = N.
The DI may be understood as a simple generalization of the Wiberg-Mayer (WM) bond order, 26,27 to which it reduces upon identification of real space averaging with atomcentered, Mulliken-like basis set condensation.Actually, for single determinant expansions with one-electron spinorbitals φ i , where S A i j = A dx x x φ * i (x x x)φ j (x x x) is the domain restricted overlap integral between spinorbitals i, j.To compare this expression to the WM bond order, which is usually written as with sums running over primitive functions χ µ centered on the A or B nuclei, and P denoting the density matrix written in terms of primitives, that can immediately be recast as δ A,B if the Mulliken condensation is made.Although any orthogonalization procedure will destroy the original adscription of primitives to centers, it is very often the case that orthogonalization tails are not dominant, and that one can still formally assign the new orthogonal primitives to nuclei.Although this lies at the core of many of the problems of Mulliken or Löwdin population analyses, it plays no role in the following.
The chemist bond order is a measure of delocalization that the DI simply puts into a proper physical context.We will use the above condensation procedure soon in what follows.DIs have been widely used, providing a number of interesting insights.Particularly important in this context is the general finding that electron correlation tends to decrease the covalent bond order well below the standard integer numbers used by chemists.
It is also important to recognize that, being a domain condensation of the exchange-correlation density, the DI reflects the two-domain statistics of electron populations.For instance Here, n A , n B are the domain electron counts, so that n A is the average electron population in region A. Delocalization, as sensed by the DI, is a measure of the fluctuation of electron populations.In chemical terms, two regions display a non-vanishing mutual bond order (if we like, they are bonded) when fluctuations in the electron population of one of them are sensed in the other, and vice versa.
At this point we also notice that most of the known results about the decay rates that we have commented above are based on TB Hamiltonians or effective one-electron formulations within DFT.In such cases, which we can assimilate to one-determinant expansions in a theoretical chemistry context, the exchange-correlation density reduces to its Fock-Dirac expression, This means that the decay rate of DIs with A − B distance, see Eq. 3, should allow us to distinguish between metalliclike and insulating-like behavior in not only extended but also finite systems.
It is our purpose to show with the help of Hückel and TB model hamiltonians that the above insights hold indeed for molecules and solids.DIs should fall algebraically in metalliclike systems, possibly showing Friedel-like oscillations, and exponentially in insulating-like molecules, with decay lengths depending on the gap.

Models, computational details
We will restrict to the simplest possible cases that can be solved both analytically and modeled via single determinant, Hartree-Fock (HF) or Kohn-Sham (KS) DFT, expansions.To simplify as much as possible, we will consider homoatomic A n linear chains of growing size with one electron per node to model metallic-like cases, and heteroatomic (AB) n ones with also one valence electron sites to understand insulating-like behavior.We will obtain DIs from analytical Hückel solutions with Mulliken condensation, and compare them to HF results in H and LiH chains obtained with the GAMESS 28 code using 6-311G(p) and 6-311+G(p) basis sets, respectively.In these cases, DIs for QTAIM topological partitions have been computed through our PROMOLDEN 29 program.
We have also obtained TB solutions for linear, square, and simple cubic one-electron per site extended lattices, and compared the decay rate of their DIs to that obtained from hydrogen lattices computed through the all-electron, full-potential linearized augmented plane wave (FP-LAPW) code ELK. 30 QTAIM DIs from Elk solutions were obtained through the DGrid code. 31e will start considering our finite analytical models.Then we will generalize to 1D-3D extended systems, and finally we will compare results with HF and KS real data.

The H ückel homoatomic chain
The Hückel homoatomic chain is an excellent semi-empirical model not only of hydrogen chains, but also of the π skeleton of alternate conjugated hydrocarbons, where a p, instead of an s function is placed at each node.We will freely switch between the s-H chain and the p-alternate hydrocarbon interpretations in what follows.
Let us label the n nodes of the chain with Latin indices, and build each one-electron function φ µ = ∑ i c i µ χ i , where χ i denotes each node's primitive and the orbital index µ runs from 1 to n.We can both consider open-ended or closed chain conformations.Both admit well-known analytical solutions, so to simplify, we will stay with the open-ended, linear chain case.This is characterized by a Hamiltonian matrix H H H = αI I I + βT T T , where α, β are the standard Hückel Coulomb and resonance parameters, respectively, and T T T i j is a Toeplitz tridiagonal adjacency matrix, with elements equal to 1 whenever |i − j| = 1 and equal to zero otherwise.Toeplitz systems are easily diagonalized by discrete Fourier transforms. 32To simplify further, let us assume that n is even.Then, the eigenvalues of H H H and its associated spinorbital coefficients are Similar solutions may be obtained for a closed chain, now by solving a circulant matrix problem.Using Mulliken's condensation, δ i, j = 2(∑ µ c i µ c j µ ) 2 (change the prefactor from 2 to 4 if the sum runs over occupied orbitals).Notice that the DI is built up from trigonometric functions, so a clearly oscillating behavior is expected.Fig. 1 shows all the DIs for both open and closed n = 6 chains, which may be seen as models of hexatriene and benzene, respectively.DIs obtained with this simple prescription have already been reported in several model systems. 33here are several points to be commented.First, notice that the DI between nodes separated by an even number of edges is exactly zero, which is valid for any value of n and for cyclic or open chains.A chemically appealing connection between electron delocalization via DIs and mesomerism thus appears.
The resonance link is very clear when the covariance interpretation of the DI is taken into account.For instance, it is straightforward to check that on building the standard Pauling resonance structures of the hexatriene analogue, if the charge of node (atom) 1 (at one edge) is altered, then only those charges of atoms 2, 4, or 6 will also be found altered in the possible resonance schemes.This means that only the 1even populations will display non-zero covariance, thus nonvanishing DIs.This interpretation may be generalized to other dimensions.
The ortho (or 1,2) DI in Hückel's benzene is 4/9, so adding the classical σ bond order would add to a total C-C bond order of 1.44, different from the naïve value 1.5.HF or DFT C-C DIs in benzene have been calculated many times, giving values clustered around 1.4.The 1,4 (para) DI, or PDI is quite large in benzene (although smaller than in the open chain), and has been successfully related to aromaticity in real calculations. 34econdly, DIs in the open-ended chain show the expected bond order alternation of alternant hydrocarbons, with an oscillatory pattern of partial double (if the σ component is added) bonds, in good agreement with chemical wisdom.If the open chain is taken as a model for H n , DIs predict the Peierls distortion (dimerization) of the hydrogen chain. 13If, on the contrary, it is understood as an alternant hydrocarbon model, then DIs predict bond length alternation.Finally, this very simple example shows that DIs decay slowly in chains: the 1,6 value is as large as 0.0908.
Let us examine now the infinite chain (n → ∞) limit.It is easy to show that This analytical expression has several interesting readings.
For instance, the open chain does not lead to bond equalization at its ends.The 1,2 and 2,3 DIs tend to 0.721 and 0.259, respectively.Equalization is however achieved far from the edges, i.e. when n ≈ n/2, where the DI tends to 4/π 2 ≈ 0.405.As our main objective is regarded, if we take s = | j − i|, δ(s) → 16/(π 2 s 2 ) (when non-zero), decaying algebraically as s −2 , in agreement with Taraskin and coworkers. 21 Similar expansions can be performed with the mid-chain first neighbor DI, so the DI may be used as an indicator of the gap for large chains, as Fig. 3 shows.

The H ückel AB heteroatomic chain
A model for a chain insulator may be easily constructed by interpenetrating two different α homoatomic lattices.Since all the physics is contained in ∆α = α − α ′ , we can arbitrarily set one of them (e.g.α ′ ) to zero.This is a model for the valence electrons in LiH, for instance.Let us construct a chain with n (n even) sites, and order them such that the first n/2 are A (α) and the second n/2 B (with α ′ = 0).Then the Hückel matrix is square-blocked, where T T T i j is again a Toeplitz tridiagonal matrix.Splitting eigenvectors into A and B components, the eigensystem is easily solved with a simple generalization of the Coulson-Rushbrook 35 theorem.The set of eigenvalues is where the plus/minus sign differentiates the occupied/virtual solutions.Similarly, with c i µ as in the homoatomic chain and τ = α/ε µ .It is not difficult to show that this is a gapped system.In the infinite n limit, using a 1/n = γ expansion similar to that used before, ∆(γ) = α + 2π 2 γ 2 /α + O(γ 3 ), and the gap approaches α on growing chain size, the faster the larger the α value.Fig. 4 shows the evolution of DI(1, j) for the n = 10 heteroatomic chain with three different α values.Notice how the results collapse on the homoatomic ones if α = 0, and how the metallic-like oscillations get damped for small values of α to completely disappear as this parameter increases.This is a very interesting insight.It is also pretty clear that heteroatomic DIs decay much faster than homoatomic ones.All DIs converge extremely fast to the n → ∞ limit.For instance, with α = 2, DI(1,2) attains the limiting 0.368 value at n = 8, with just 4 AB units.
Fig. 4 also shows the onset of the exponential decay of DIs for even pretty small n values.Our δ 1, j = δ(s) falls off exponentially with exponent λ approximately equal to 1.5, 1.8 for α = 3, 4 respectively.Our finite chain results support the proportionality between the gap and λ in the small gap limit.In this case, ∆ ≈ α ≈ 2λ.The faster the decay rate, the larger the gap.This is a valuable insight in molecular calculations.

Periodic analytical models
The calculation of DIs from TB models in one to three dimensions has been pioneered by R. Ponec, who first presented a simple calculation, 36 later extended and reformulated. 37His second paper actually provides tight binding results under the Mulliken condensation approximation discussed previously.DIs from DFT calculations over QTAIM or ELI real space log δ 1,j j Fig. 4 Top: δ 1, j for the n = 10 AB heteroatomic chain at α = 0 (black), 1 (green), and 2 (blue).Oscillations rapidly disappear as α grows from 0. Bottom: Evolution of the logarithm of δ 1, j for a n = 20 chain with α = 3 (solid red with crosses) and 4 (dashed blue with dots).
domains are available since the work of Kohout and coworkers. 25,38However, all these authors have been more interested in first or second-neighbor DI values than in the decay rate of the indices.We will focus here on this second aspect for homoatomic lattices, referring the reader to the above-mentioned papers for further details.Imposing periodic boundary conditions (PBC) on a lattice with one primitive function χ per site allows us to use Bloch's theorem 13 to immediately write the one-electron Bloch functions as where R R R runs over real space lattice vectors, N is the total number of sites, and k k k runs over the first Brillouin zone (BZ).Under a nearest neighbors TB (or Hückel) hamiltonian, the above Bloch ansatz leads to the one-electron eigen-values: where R R R n only covers nearest neighbors.
In a 1D lattice with lattice parameter a, where −π/a ≤ k < π/a, we will have ε k = α + 2βcos(ka), which may be compared with our previous finite 1D results.To obtain the DI between two lattice sites, let us center our reference frame at one of them (the 0 site, with R R R = 0 0 0).The band (orbital coefficient) of this site is independent of k, c 0 k = 1/ √ N.That of site located r lattice parameters away, c r k = 1/ √ Ne ikra .Using then the Mulliken condensation scheme, and integrating over the BZ, dk e ikra 2 = 4 sin 2 (πr/2)/(r 2 π 2 ).
(17) Of course, this result is equivalent to our previous infinite n limit, and shows that only if r is odd the DI does not vanish.PBC also leads to bond equalization.All nearest neighbor DIs are equal to 4/π 2 , independently of the lattice parameter a.This is of course a flaw of the TB hamiltonian.It is also relevant to comment on the on-site localization index (half the diagonal δ 0,0 value), which turns out to be equal to 1/2, showing that half of the electron population is localized, half delocalized over the full lattice.Notice also that the sum rule . Integration over the BZ in a 2D square lattice of lattice parameter a is again trivial, since the Fermi surface is a perfect square.Taking an arbitrary site of the lattice as origin O, we will label any other site with Cartesian coordinates (ra, sa) with the (r, s) integer pair, which reduces to Our results show a more complex landscape than that provided by Taraskin and coworkers, who would describe an inverse third power decay.Here we show that the decay rate depends on the direction, following an inverse fourth power law envelope along the (1,0) direction.The localization index of each site is again equal to 1/2.The oscillatory pattern found in 1D is seen here to propagate in 2D.From a given site, the network of nodes with non-zero DIs resembles a check board.This behavior is clearly related to the ability of this lattice to be decomposed into two 45 • rotated interpenetrating alternate sublattices, like in alternate hydrocarbons.It seems that DIs between elements of the same sublattice vanish.Again, this may be understood in terms of charge fluctuation (covariance) if the allowed Pauling resonance structures are examined.Several interesting investigations regarding this should be undertaken.On the one hand, it would be interesting to check the behavior of frustrated lattices.On the other, it would be of great interest to study also the chemical consequences (as noticed with mesomerism in 1D) of these patterns.
For the time being, it is relevant to shift to the crystalline coordinates of the sublattices.This can be done by using two new orthogonal coordinates p = r + s, q = r − s.With this, with p, q both odd, and δ = 0 otherwise.We thus see that the square lattice behaves as a Cartesian product of two 1D networks.With this expression, it is straightforward to show that the sum rule adding to the site population of 1 electron is also fulfilled.Nodes along the (1,1) diagonal belong to the same sublattice (r + s is even, or q = 0).The decay along p =constant-odd or q = constant-odd diagonals follows an inverse square power law, and if particular relations between p and q are satisfied along a nodes sequence, intermediate power laws also appear.We have found it difficult to obtain an analytical angularly averaged decay rate.
The non-trivial shape of the Fermi surface in the 3D case precludes an analytical integration over the BZ.Anyway, if we label nodes on the simple cubic lattice by the trio (r, s,t), then which may be reduced to simple numerical quadratures.The symmetry properties of the above expression allow us to assure that δ 0,rst is only non-zero when r + s + t is odd, and we can again consider the lattice as composed of two interpenetrating sublattices such that δ only communicates nodes belonging to different sublattices.We have no analytical decay rates, but clear numerical evidence points towards faster, likely inverse sixth power, decay speed.As an example, δ 0,100 , the nearest neighbors DI, is equal to 0.112 (to be compared to 0.405 and 0.164 in 1D and 2D, respectively).Summarizing, extended TB models show algebraic decay of DIs in gapless homoatomic systems coupled to a very interesting interference cancellation that leads to wild oscillations that may be traced back to Friedel behavior, from the physical point of view, or to resonance and mesomerism, from the chemical one.

Single determinant (HF, KS) results
We will discuss here how the analytical models compare to actual one-determinant (or pseudo-determinant, in the case of DFT) calculations in hydrogen and lithium hydride toy systems.We have chosen interatomic separations for which these methods are known to provide reasonable answers, and a QTAIM space partitioning.We leave the true role of electron correlation aside, that we expect to consider soon elsewhere.All calculations have been performed at fixed geometries.As we will see, the exact interference cancellation behind zero DIs disappears as we allow for the primitive functions to overlap.However, many of the insights developed from the Hückel or TB models remain valid.The exponent of the power law decay, f , is to be compared with the Hückel or TB result, f = 2.This inverse square law is also represented in blue.The lattice parameter, or nearest neighbor H-H distance, is set to 1.84 bohr the theoretical limiting equilibrium parameter for a HF cyclic chain as n grows.
Fig. 5 shows that actual calculations in 1D chains display deep oscillations, and that δ 1,2i+1 values are non-zero, but certainly much smaller than δ 1,2i ones.Both odd and even envelopes evolve algebraically, with exponents larger than 2, but close to it.Several other points may be highlighted.For instance, f decreases approaching 2, as we move from openended finite to cyclic finite, and finally to PBC infinite chains.This is to be expected, since open finite chains differ considerably from the stringent approximations of the Hückel or TB models.We have also found that results in finite chains converge very quickly with size, as in the models, and that our computed values are quite close to those provided by the latter.For instance, DIs δ 1,2(4) computed in the infinite chain are 0.44, 0.04, to be compared with the TB results, 0.39, 0.04, respectively.
Changing the lattice parameter does only introduce quantitative changes in the picture.For instance, at a = 2.5 bohr, probably out of the confidence window where KS-DFT is reliable for this system (see Fig. 6), the PBC chain f value equals 2.21, slightly closer to 2, the value that it should attain at dissociation values of a. Changing the dimensionality alters qualitatively the analytical results.We will show only PBC calculations in square and simple cubic 2D, and 3D H lattices, both computed at a = 2.5 bohr.This is the lattice parameter used by Baranov and Kohout 38 (BK) in a seminal study of first neighbors DIs in solids.We use it here so that the reader may compare our values with those obtained by BK. Results at a = 1.84 bohr do not differ qualitatively from those shown here.Fig. 6 depicts that the decay is algebraic in the three cases, with f values roughly increasing in 2 units as we change dimension.What is noticeable is that oscillations disappear in 2D and 3D, while they widely persist in the TB models.We think that this is due to the increase in the number of neighboring overlaps that contribute to cancelling the destructive interference that lies behind the oscillations found in TB.In 1D, each site's primitive overlaps with 2 nearest neighbors, while in 2D and 3D this number increases to 4 and 6, respectively, or even more if we consider second neighbors.Be it as it may, our results clearly support an algebraic decay of DIs in gapless systems, with f values increasing steadily on going from 1D to 3D.
We now turn to insulating materials.This time we will present HF finite calculations in a (LiH) 9 1D chain and a 9 × 9 LiH square 2D foil, both with fixed Li-H distances equal to 3.0 bohr.In order to avoid as much as possible termination effects, we have built in each case models in which a Li or a H atom is placed at the center of the chain (or foil).Fig. 7 (top) shows the QTAIM DIs.As expected, their decay with distance is extremely fast.So fast, indeed, that we have not been able to obtain numerically reliable values beyond fourth neighbors.Notice that the existence of an AB lattice introduces three types of DIs: Li-Li, Li-H, and H-H.This is the origin of the kinks in the plot.For instance, the second neighbors ( j = 3) DI between Li atoms is almost one order of magnitude smaller than that between H atoms.With such a quick decay we have not enough data to support an exact exponential decay, but we can rule out a slow polynomial one.Notice also that the exponent seems to increase with dimensionality.This also explains why we have not added 3D data to the figure: numerical issues render even fourth neighbors unreliable for them.Decay rates do also depend on the direction, as expected, but numerical problems again preclude us from extracting precise conclusions.Numerical issues are much less important if instead of QTAIM basins we use Eq. 7 together with, for instance, Löwdin's orthogonalization.Fig. 7 (bottom) shows the exponential-like decay of these Löwdin DIs in a (LiH) 17 linear chain.DIs decrease 8 orders of magnitude on traversing the chain.Although the chemical interpretation of these Löwdin indices is prone to severe criticisms, they serve well our purpose of showing the evolution of decay rates.Overall, analytical and real models support an exponential, non-algebraic decay of DIs in gapped systems.Further work is necessary to establish trends.

Conclusions and prospects
The search for real space descriptors that could discriminate metallic from insulating materials has been a recurrent quest in chemical bonding theory in the last decades.After the reformulation of Kohn's theory of the insulating state by Resta, 15 it is now known that electrical conductivity does not leave recognizable scars in the electron density itself.However, the modern theory of polarization, as the reformulation is known, points towards a possible link between electron delocalization measures and the insulating or conducting nature of a material.Previous knowledge in the physical literature had noticed that the decay rate of the 1RDM changes from algebraic to exponential when going from metals to insulators in tight binding models.Here we have shown that since the delocalization index of real space theories of the chemical bond is dominated by the square of the 1RDM, it must follow the same behavior, so a link between bond orders (that is what the DI measures in isolated molecules) and conductivity appears.
To that end, we have solved several Hückel finite and TB extended models.As it turns out, in fairly small molecular chains the shift from polynomial to exponential decay is evident when a gap is forced in the system, and the larger the gap, the faster the decay.Our results show that in metallic-like systems (because the gap does only close when we go from finite to infinite systems) the DI decays algebraically, wildly oscillating due to quantum mechanical interference cancellations that annihilate the DI between nodes belonging to the same alternating sublattice.These oscillations have been shown to be intimately linked to Pauling resonant structures and chemical mesomerism, well known in alternant hydrocarbons.The dimensionality of the system simply changes the decay exponent, larger as we go from 1D to 3D.
Real computations within the Hartree-Fock or Kohn-Sham single (pseudo)determinant schemes together with a space partitioning according to the quantum theory of atoms in molecules, show that DIs decay algebraically in metallic-like molecular systems, and very like exponentially, or at least ex-tremely fast, in insulating-like ones.Oscillations persist after the inclusion of overlap in 1D chains, although the vanishing DIs in the analytical models are now small, but nonzero.Overlap seems to block interference in larger dimensions, making oscillations to disappear (or at least dampen substantially).
Examining the decay of DIs in real molecules and extended materials may provide very interesting clues to their conducting behavior.Since DIs may be computed between any pair of atoms, their decay may be followed along particular directions, making it possible to detect facile or non-facile conductivity channels.This may provide relevant information in material science and molecular electronics.
The impact of electron correlation on these results remains to be ascertained.We expect it to be small in simple systems at geometries close to equilibrium, but it should be important, yielding interesting insights into metal-insulator transitions, when the single (pseudo)determinant approximation ceases to be useful.Work along this direction is under way.

Fig. 2
Fig. 2 Hückel DI(1,i) for the homoatomic linear chain as its size (n sites) grows.Different lines correspond to different length chains.At the end of the chain j = n.

Fig. 3
Fig.3Evolution of the HOMO-LUMO gap against the midchain DI for the Hückel homoatomic chain.The length of the chain increases as the gap approaches zero (and the DI reaches its 4/π 2 limit).All data in a.u.

Fig. 5
Fig.5Logarithmic plot of δ 1, j against the H-H distance in the H 1D infinite (black), cyclic finite with 42 atoms (green), and open-ended finite with 28 atoms (red) hydrogen chain.Linear fittings are superimposed to clearly observe the algebraic decay rate.The exponent of the power law decay, f , is to be compared with the Hückel or TB result, f = 2.This inverse square law is also represented in blue.The lattice parameter, or nearest neighbor H-H distance, is set to 1.84 bohr the theoretical limiting equilibrium parameter for a HF cyclic chain as n grows.

Fig. 7
Fig.7Top: Semilogarithmic plot of the QTAIM δ 1, j against the node index in a LiH 1D chain (with 9 LiH units, solid line), and a LiH 2D square foil (9 × 9, dotted line), both with a = 3.0 bohr.The red (for Li) and black (for H) colors distinguish which atom is placed at the center of the model and labelled as node 1.Only the (1, 0) direction is shown in the 2D case.Bottom: Semilogarithmic plot of Löwdin δ 1, j against the node index j in a LiH 1D chain of 17 LiH units, using also a = 3.0 bohr.DI(Li,H) is shown with crosses, and DI(Li,Li) with squares.Minimum square lines are also shown just to aid the naked eye.