Decay Rate of Correlated Real-Space Delocalization Measures

Subscriber access provided by UNIV OF NEBRASKA - LINCOLN

Article

Decay rate of correlated real space delocalization measures: insights into chemical bonding and Mott transitions from Hydrogen chains Alfonso Gallo-Bueno, Miroslav Kohout, and Angel Martín Pendás J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00139 • Publication Date (Web): 02 Jun 2016 Downloaded from http://pubs.acs.org on June 7, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Decay rate of correlated real space delocalization measures: insights into chemical bonding and Mott transitions from Hydrogen chains A. Gallo–Bueno,† M. Kohout,‡ and A. Martín Pendás∗,† Departamento de Química Física y Analítica. Facultad de Química. Universidad de Oviedo. 33006 Oviedo. Spain., and Max Planck Institute for Chemical Physics of Solids. Nöthnitzer Strasse 40, 01187 Dresden. Germany. E-mail: [email protected]

Abstract We study in this contribution the spatial decay rate of real space localization and delocalization indices in correlated systems. To that end we examine Hubbard and quantum chemical models of simple cyclic hydrogen chains, showing that all descriptors of delocalization converge quickly towards the infinite chain limits. It is then shown that the localization index may be understood as a generalization of the standard order parameter in Mott insulator transitions, and that the origin of the enigmatic sigmoidal profile of delocalization indices in chemical bond breaking processes lies in the nonlinear mapping between intersite distances and correlation parameters. Although the ∗

To whom correspondence should be addressed University of Oviedo ‡ Max Planck Institute †

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

long-range asymptotic decay of delocalization indices is exponential, we show that as the correlation parameter decreases quantum mechanical interference sets on and a switch to an oscillating pattern, related to core chemical concepts like resonance or mesomerism, appears.

1

Introduction

Electron delocalization, an intimate consequence of quantum mechanical behavior, lies behind most of chemical phenomena. In fact, the commonly held principle that the chemical bond implies overlap could be rephrased as there is no chemical bond without electron delocalization, a more general statement that makes no direct reference to the orbital paradigm. The basic prototypes of covalent and ionic bonding use delocalization channels in order to occur. Neither the former, loosely characterized by electron sharing, nor the latter, associated to charge transfer, are possible if electrons are only allowed to dwell in disconnected regions of the many particle Hilbert space. 1 Even dispersion interactions are subject to this rationale. From the chemist’s perspective, it is in the real space where chemical phenomena occur, although the standard approach to chemical bonding is based on the molecular orbital (MO) description. 2 Actually, during the last few decades, an alternative paradigm, based upon analyzing orbital invariants constructed from reduced density matrices (RDMs), has been developed. These techniques are known as Quantum Chemical Topology (QCT), 3 and among them, the Quantum Theory of Atoms in Molecules (QTAIM) by R. F. W. Bader and coworkers, 4,5 based on the scrutiny of the electron density is best known. QCT has delivered over the years a profusion of information about chemical bonding, including measures of electron delocalization, after the work of Bader and Stephens, 6 who introduced the delocalization index (DI(A, B) or δ A,B ) between two disjoint regions A, B of real space. The QTAIM partitions the space into domains associated to atoms and DIs have usually been obtained over this kind of domains, although this is not necessary. Delocalization indices 2


Page 2 of 32

Page 3 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


have been appropriately generalized to include many centers 7–9 or to be applicable to extended systems. 10 The DI between atomic regions may be shown to be directly related to the standard Wiberg-Mayer 11,12 bond order, thus quantifying basic core chemical concepts. Although they have been computed in many systems, usually under the single determinant Hartree-Fock (HF) or Kohn-Sham (KS) approximations, the exact behavior of DIs is not known in general circumstances. This is particularly true as concerns their response to changes in the A, B distance, since in most cases this implies an appropriate treatment of electron correlation. We have already shown 13 that the profile of DIs along reaction coordinates reveals the nature of chemically relevant interactions, with sigmoidal profiles signaling the formation or breaking of chemical bonds and exponential shapes detecting non-bonding clashes. If read in terms of plain physical quantities, these findings tell us about the different patterns of electron delocalization taking place in both kinds of chemical processes. Another almost unexplored territory lies in how DIs decrease with inter-center distance in (large) molecules or extended systems. This is a relevant issue, for their algebraic or exponential decay rate might be related to the metallic or insulating nature of the system. This relation stems from the seminal work of W. Kohn, 1 who showed that the degree of delocalization of the many particle wave function of a system in its ground state is intimately connected to its response to electric fields and thus to its insulating or metallic character. Kohn’s theory of the insulating state has been reformulated by Resta, revitalizing the field. 14 The link between ground state delocalization and conductivity is particularly interesting in chemistry, for it interweaves very naturally with the standard chemical notions, quantifying the vague but important concept that associates metals to exceptionally delocalized systems. Independently from Resta, several authors 15–17 showed for tight binding (TB) models that the decay rate of the off-diagonal elements of the first order RDM is algebraic for metals, their power law depending on the dimensionality of the system, and exponential for insulators, the exponent increasing with the gap. For single determinant approximations, including KS determinants, DIs are constructed from domain-averaged squared 1RDMs. It

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

has been previously shown 15 that, both in finite molecules and in extended systems, DIs decay algebraically/exponentially for metallic-/insulating-like materials, both in Hückel or TB models and in actual simple systems computed at the HF or KS levels. It was also found that DIs show an oscillatory behavior in metallic-like TB models, vanishing whenever simple selection rules were satisfied. This behavior, related to quantum mechanical interference, persists in HF or KS calculations only in one-dimensional chains, disappearing in 2D and 3D systems. It was also proved that the oscillations may be easily interpreted in terms of the classical Pauling resonance theory or mesomerism, establishing a link between physical behavior and chemical reactivity. To close the loop, we examine in this contribution how the decay rate of DIs behaves under geometrical changes. In extended systems, this is coupled to the possible shift between a metallic and an insulating state, a metal-insulator transition (MIT). In molecules, to the shift from a sigmoidal to an exponential DI shape. Little is known about how DIs behave in these cases, since examining their long-range behavior in a geometrical rearrangement normally implies entering strong correlation regimes, which are difficult to deal with in large molecules. To keep the discussion as simple as possible, the focus will be placed on the 1D cyclic hydrogen chain. This is no doubt one of the simplest systems where a possible metal to insulator transition may occur, induced by the localization of electrons in their original potential wells as the H–H distance increases. This collective change is known as a Mott transition, 18 which is difficult to model with standard theoretical chemistry tools. In a simple chemical language, increasing the H-H distance in the H-chain is equivalent to doing so in the H2 molecule. The simplest qualitatively correct solution in this case involves mixing two Slater determinants or, similarly, performing a complete active space CAS[2,2] calculation. Since the size of the CAS space explodes combinatorially if we consider a CAS[n,n] model when the number of atoms in the chain, n, increases, the problem soon becomes intractable with usual techniques. Standard KS density functional theory (DFT)

4


Page 4 of 32

Page 5 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


calculations are of no help here, although several corrections may be applied to partly correct the problem. 19 Out of the basic quantum chemical methods, difficulties to deal with 1D (or their 2D, 3D analogues) H chains appear from: (i) the lack of size-consistency in truncated CI expansions, as exemplified by configuration interaction with single and double excitations (CISD); (ii) convergence problems as interatomic distances increase, so common in coupled-cluster expansions; 20 (iii) little resemblance of the reference determinant with the final solution in single-reference methods; or (iv) a combination of all of them. In extended systems, standard DFT methods also face most of the above problems. Intense efforts have been devoted in the last years to deal with strongly correlated systems, with varying success. Among the battery of techniques that have been devised, we may cite the density matrix renormalization group (DMRG), 21 together with other tensor and cluster product methods, 22 density matrix embedding theory (DMET), 23,24 variational two-electron reduced density matrix methods, 25 constrained-pairing mean-field theory (CPMFT) 26 and other symmetry breaking and restoration procedures, auxiliary field quantum Monte Carlo (AFQMC), 27 the dynamical cluster approximation (DCA), 28 reduced density matrix functional theory (RDMFT), 29 the pair coupled cluster doubles (pCCD), 30 approximations to antisymmetrized products of interacting geminals, like the AP1roG method, 31 etc. Noteworthy are also the everyday more successful experimental quantum simulations of strongly correlated systems using ultracold fermionic atoms in optical lattices. 32,33 Given the importance of strongly correlated systems, an also large number of simplified models habe been proposed to cope with them, usually formulated in second-quantized language. For the problem treated here, the Hubbard hamiltonian 34 is the simplest one, since it allows for controlling the strength of electron correlation through a continuous parameter. When this parameter vanishes, the Hubbard model with nearest neighbour couplings falls onto the Hückel or TB one, so a relatively simple comparison with standard quantum chemical results is also at hand. The Hubbard model has been extensively studied in one, two, three and infinite dimen-

5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

sions, at zero and finite temperatures, in all kind of lattices. 35 In 1968, Lieb and Wu 36 (LW) showed that the 1D Hubbard chain admits an analytical solution and that there is no Mott transition for this model, the ground state remaining a gapped antiferromagnetic singlet at any non-zero value of the correlating parameter. The existence of the LW solution has made the 1D Hubbard chain a cherished system in which to check the power of all of the above-mentioned methods, and many papers examining the performance of DMRG, DMET, DCA, etc, have been published. We refer the reader to a specific monograph on the 1D Hubbard model for futher information. 37 An even more interesting situation exists in 2D or 3D, where no analytical solutions have been found. In 2D, for instance, the state of the art has been internationally assessed, and the Simons Collaboration on the many-electron problem has recently published their conclusions after comparing the results of a large number of methods. 38 Recent results 39 conclude that there is indeed no MIT at low T for U > 0 in the half-filled 2D Hubbard model. An important topic when dealing with phase transitions is the choice of adequate order parameters that may sense their onset. Interestingly, a possible order parameter for the MIT within simple lattice models is the so-called double occupancy, D, that measures the probability of finding two electrons at a given node in the lattice (of necessarily opposite spin, for only one spinless state is available per site). Together with spin(less) two-point pair correlation functions, or with the decay of the off-diagonal terms of the 1RDM, these quantities have been repeatedly used to detect long-range order. 26,40–42 At a first order MIT, D drops to exactly zero, remaining so in the insulating phase. As it will be shown, D is directly related to the real space (de)localization measures used in QCT, so another unexpected relation between seemingly unrelated concepts appears that deserves being considered. The aim of this work is to explore and compare the behavior of DIs in both the Hubbard model at different correlating parameters and in simple correlated descriptions of hydrogen chains. In the process, the origin of the sigmoidal shape that has been found for DIs with

6


Page 6 of 32

Page 7 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


interatomic distance in homolytic bond breaking processes, 13 will be proved to lie in the nonlinear relation between the Hubbard correlating parameter and the interatomic distance. It will also be verified that the oscillatory distance decay of the DI that characterize 1D metals in single determinant descriptions vanish as correlation increases and that its onset may be understood as a transition from a tail-overlap-dominated regime to one in which quantum mechanical interference becomes prominent. The paper is organized as follows. First the models and descriptors used will be presented, together with some of the basic knowledge related to them in the context of the present work. Then, some results will be shown, first convincingly proving that the combined use of quantum chemical and Hubbard results allows to consider the rest of the data with reliability. Finally, the decay rate of DIs will be examined.

2 2.1

Delocalization measures and the 1D Hubbard chain Real space Localization and Delocalization indices

Delocalization measures in real space, after the works by Bader and Stephens 6 in 1974 and Fradera, Austen, and Bader 43 in 1999 are interpreted as a way to quantify the number of pairs of shared (i.e. delocalized) electrons between two real space domains of an N electron system. Two descriptors are introduced: the localization index, LI or λ, and the delocalization index, DI or δ. Starting from the spinless diagonal first and second order RDMs (i.e. the density and P ˆ + (r)Ψ ˆ σ (r)|Ψi the pair density, normalized to N and N (N −1), respectively), ρ(r) = σ hΨ|Ψ σ P ˆ + (r1 )Ψ+′ (r2 )Ψ ˆ σ′ (r2 )Ψ ˆ σ (r1 )|Ψi, the LI and the DI are obtained and ρ2 (r1 , r2 ) = σ,σ′ hΨ|Ψ σ σ

by domain-averaging the exchange-correlation density, ρxc (r1 , r2 ) = ρ(r1 )ρ(r2 ) − ρ2 (r1 , r2 ): A

LI(A) = λ =

Z

dr1 A

7

Z

dr2 ρxc (r1 , r2 ). A


(1)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Similarly, DI(A, B) = δ

A,B

=2

Z

dr1 A

Z

Page 8 of 32

dr2 ρxc (r1 , r2 ).

(2)

B

Since ρxc is also the second order cumulant density, it is size extensive,

R

dr1 ρxc (r1 , r) =

ρ(r), so when an exhaustive partition of space into disjoint regions is performed λA + P (1/2) B6=A δ A,B = nA , where nA is the average electron population in domain A, nA = R drρ(r). A As already noticed, the DI is a real space extension of the Wiberg-Mayer (WM) bond or-

der (BO). 11,12 DIs and WM bond orders can be considered equivalent if we identify averaging over a real space domain with Mulliken condensation. 15,44 DIs have been successfully used to quantify covalent bond orders for several decades, showing how a key chemical concept may be derived from orbital invariant quantities. The interacting quantum atoms (IQA) approach 45,46 allows us to show that DIs are related to the covalent component of interatomic interaction energies. The cumulant part of the two-particle reduced density can also be obtained from the second order cumulant density matrix, which norm has been used as an entanglement measure. 47,48 Interestingly, LIs and DIs also admit a statistical interpretation. Using the density and P P pair density operators, ρˆ(r) = i δ(r − ri ), ρˆ2 (r1 , r2 ) = i6=j δ(r − ri )δ(r − rj ), we can R construct domain restricted particle (and pairs of particles) operators, 49 n ˆ A = A dr ρˆ(r), 1R R nA d dr dr2 ρˆ2 (r1 , r2 ). With this, = 1 2 A 2 A λA = nA − var(nA ) = nA − h(nA − hnA i)2 i = nA − (hn2A i − n2A ),

(3)

since nA = hˆ nA i. Similarly, δ A,B = −2cov(nA , nB ) = −2 (hnA nB i − hnA ihnB i) .

(4)

Notice that only in the case that the electron population in a region does not fluctuate at

8


Page 9 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


all will λA be equal to nA . Taking into account that bond orders are naïvely used to identify bonded atoms, the link between DIs and WM BOs, and thus the connection between the fluctuation of the electron population, introduces new insightful interpretations in chemical bonding theory. For instance, whenever the DI between two regions is non-zero, i.e. whenever they are bonded, the population of one responds to a fluctuation in the population of the other and vice versa. In other words, two regions are bonded together if there is a delocalization channel between them. In those circumstances, any perturbation in the electron population of one of them will be followed by a change in the population of the other. How standard bonding concepts change in cases of strong correlation deserves definitely further study. The interpretation of chemical bonding in terms of the fluctuation of electron populations is being actively developed in terms of the concept of electron number distribution functions (EDFs). 50–52 Very succintly, once a partition of space, e.g. into m atomic regions, has been chosen, an EDF is the statistical distribution function corresponding to partitioning the N electrons of a system into these m domains. Its basic quantity is p(n1 , n2 , . . . , nm ), the probability of finding exact integer number of electrons n1 , n2 , . . . , nm in regions 1, 2, . . . , m. Algorithms for computing these p’s from one- and multi-determinant wave functions have been developed. 52 Apart from the 1- and 2RDMs, two-center (or two-site) variances and covariances may also be computed from EDFs, i.e. from the probabilities p(nA , nB ) of finding given numbers of electrons in the A and B regions. Moreover, in the same way as the (co)variance in the case of two regions senses two-center bonding, the i−th cumulant moments of the EDF are related to multi-center bonds. 9 It becomes clear now how the standard metal-insulator transition (MIT) order parameter, the probability of double occupancy of a site, D = p(2), is related to real space descriptors. D is nothing but one of the components of a one-center EDF. On the one hand, D only vanishes for a Hubbard system when the variance at a site is zero. This is clear, but stems naturally from the average site population, hni = 2p(2) + p(1) = 1, and the sum rule 9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

p(0) + p(1) + p(2) = 1. If p(2) = 0, then p(1) = 1 and p(0) = 0. Since the average site population in the half-filled Hubbard chain is 1, using Eq. 3 we see that, in terms of order parameters, var = 1 − λ, plays the same role as D, vanishing when D = 0. Actually, the variance is more general and it only vanishes if the electrons are completely localized. It is also possible to resolve EDFs in spin components 51 and to define p(↑↓), etc. It is also immediate to show that for the ground state of H-chains the symmetry p(2) = p(0) guarantees that var = 1 − λ = 2p(2) = 2D. With this, many of the results found in the physics literature about the Hubbard model can be directly read on the chemical scale. Let us stress that we need the two-particle 2RDM to get population variances. A final comparative note is also due. DIs may also be understood as real space domain-averages of the cumulant part of the pair correlation function, which is widely used in the physics community. 42 Fig. 1 shows the evolution of the DI with the internuclear distance in the singlet and triplet states of the H2 molecule, computed at the qualitatively correct CAS[2,2] level. This wellknown behavior has been published many times 13 and is reproduced here for convenience. In a homodiatomic there is only a possible exhaustive partition of the space into two equivalent atoms that corresponds to the one in the QTAIM: the surface separating both regions is the plane that bisects the internuclear axis. In H2 , the full EDF has three components, p(2, 0) = p(0, 2), p(1, 1), corresponding to the probabilities that the two electrons are found in one, the other, or each of the two atoms, respectively. At the HF level, the two electrons are statistically independent, so p(2, 0) = 1/4, and p(1, 1) = 1/2. 53 This leads to DI= 1, i.e. to a full electron pair being shared between the atoms, to a bond order equal to one, and to the ideal covalent bond model. The inability of HF to dissociate the molecule correctly makes the HF DI remain unchanged with the internuclear distance a. Only the introduction of electron correlation restores the expected behavior, and already at the CAS level (and beyond), the DI correctly decays to zero as the internuclear distance increases. In EDF terms, this is done at the expense of increasing p(1, 1), which tends to 1 in the dissociation limit. This decreases the covariance, so at the

10


Page 10 of 32

Page 11 of 32

0.9 0.8 0.7 0.6 δH,H

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.5 0.4 0.3 0.2 0.1 0.0 1.0

2.0

3.0

4.0 5.0 a/bohr

6.0

7.0

8.0

Figure 1: DI(H,H) versus the internuclear distance, a, for the X1 Σ+ g (sigmoidal curve in 3 + red) and b Σu (blue) states of the H2 molecule computed at the CAS[2,2]//6-311++G(d,p) level. correlated level, the DI at equilibrium is about 0.85, close but not equal to the ideal single bond unit value. Notice that in a diatomic molecule, var(nA ) + cov(nA , nB ) = 0, so δ H,H = −2cov(H, H) = 2var(H) = 4p(2) = 4D is the direct analogue of the extended order parameter D. As previously commented, the sigmoidal shape of the ground state DI is characteristic of (homolytic) bond breaking processes. When a bond breaks into open-shell fragments, many examples 13 show that its associated DI decreases slowly at first, then accelerates its drop, and after an inflection points it slows again its falling rate. In all the cases, the DI at the inflection point turns out to be extremely close to half its maximum value, where the bond is "half-broken". This chemically intuitive analogy may be modelled in terms of EDFs. For two-electron links, a transfer parameter t ∈ [0, 1] measuring the breaking process monitors the DI. t = 1/2 implies a half-maximum DI. In our case, the inflection is found at a ≈ 2.90 a.u. Similarly, the exponential behavior of the triplet is typical of non-bonded contacts, like the one here

11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

presented, or of the dissociation towards closed-shell molecules like He2 . The origin of the ground state sigmoid is basically unknown. However, since the appearance of sigmoidal shapes is usually associated to cooperative effects, the relation between the DI and a phase transition order parameter is intriguing.

2.2

The 1D Hubbard chain

One of the simplest models proposed to gain insight into the difficult strong correlation regime in the thermodynamic limit was proposed by Hubbard 34 and has been extensively studied so far. The literature is vast, since almost any new methodology is tested and validated with it. A part of it can be found in a monograph, 37 and in a number of contributions containing extensive authoritative recent information. 38,39,47,54 In its basic formulation the Hubbard hamiltonian models a one-dimensional chain composed of identical single energy level sites (that can thus host two opposite spin electrons at most), subject to nearest neighbor interactions. In the half-filled band case, each site accommodates an electron on average, corresponding to an N -sites, N -electron problem. Usually, periodic boundary conditions (PBC) are imposed and N is made to tend to infinity. Two electrons are supposed to interact strongly only when located at the same site via a positive interaction energy U > 0, the so-called on-site Coulomb repulsion. They are also allowed to delocalize between neighboring sites through a kinetic-like parameter t, known as the hopping energy, which is envisioned to be determined by the effective overlap of their site states. Although no exact link exists between the t, U Hubbard parameters and those of a conventional quantum chemical Hamiltonian, there have been many attempts that try to fit or approximate them. Particularly interesting in this sense is the work of Spalek and coworkers 55 who, through variationally determined renormalized single-particle wave functions, provide a chemically intuitive way to establish that connection.

12


Page 12 of 32

Page 13 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


With this, the Hubbard hamiltonian is now written in second quantized form as

H = −t

X

+ (c+ iσ cjσ + cjσ ciσ ) + U

X

ni↑ ni↓ .

(5)

i

hi,ji,σ

In the above expression, the c+ iσ and cjσ are fermionic operators that create or annihilate a σ spin electron at site i, respectively. They are subjected to standard anticommutator + + relations {ciσ , c+ jσ ′ } = δi,j δσ,σ ′ , {ciσ , cjσ ′ } = {ciσ , cjσ ′ } = 0. The hi, ji sum runs over first

neighbors only, with each term describing the hopping of an electron from site j to site i. niσ = c+ iσ ciσ is the spin-σ number operator at site i. Notice that the anticommutation relation guarantees that the maximum occupation of a site is 2. Finally, the second term adds a Coulombic repulsion energy (a positive energy U ) to each site that is doubly occupied. t and U play opposing roles, so the model is conveniently described by the U/t dimensionless correlating parameter. Small U/t values favor hopping, thus delocalized solutions. As the correlating parameter increases, so does the penalty for electron hopping. At a MIT, the system ceases to delocalize and each electron sits permanently at the same site. It is important to note on passing that the MIT parameter D = p(2) equals hni↑ ni↓ i, so D = de/dU , e being the energy per site. This means that if D = 0 in the Mott insulating state, the energy ceases to change with U . From the quantum chemical point of view, it is interesting to examine the two-site Hubbard system, which is an analogue of the H2 molecule in a minimal basis. This textbook model is immediately solvable in the 6D space formed by the 42 states constructed by arranging two electrons in two nodes, which reduces to dimension 4 in the Ms = 0 spin sector. √ Using r = U/(4t) as correlating parameter and defining κ = 1 + r2 − r, the two lowest

states are a singlet, at energy E = −2tκ, and a triplet at energy E = 0. The probability of any electron configuration, thus the EDF in this model, may be read directly from the square of the coefficients of the eigenvectors. The singlet is a mixture of Heitler-London-(HL) like covalent and ionic forms, while the triplet is fixed at its HL-like state. Notice that, since no

13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

√ explicit overlap appears, the triplet, for instance, is 1/ 2(↑↓ + ↓↑). With this, it is clear that D = p(2) vanishes for the triplet and that without overlap but with the antisymmetry restrictions, no interaction may occur, so E = 0 constantly. From the singlet eigenvector, one has p(2) = 1/2 (κ/(1 + κ2 )), which decays smoothly towards zero from its 1/4 starting value at r = 0 as r increases. The statistics of the electron distribution for the singlet and triplet states follows the same trend as in a minimal basis full configuration interaction (FCI) under Mulliken condensation, and has the same r = 0, ∞ limits. The triplet is non-bonding since p(2) = 0 always, and the singlet changes from the non-interacting p(2) = 1/4 that we already examined at r = 0 towards p(2) = 0 as r increases. In this sense, the Hubbard model correctly captures the basic population correlations of the H2 molecule. Notice also that the bond breaking, chemically parameterized by an increasing H-H distance, is projected as a process in which r grows. In passing, the map between U = 0 and the TB approximation is met when the Hückel α and β parameters are set to 0 and −2t, respectively. Since most of the concepts used in practical chemistry 2 originate from independent particle descriptions, the U = 0 limit makes the Hubbard model an adequate starting point towards constructing improved tools in chemical bonding theory in cases of strong correlation. Surprisingly, Lieb and Wu (LW) showed in 1968 36 that the Hubbard infinite 1D chain has analytic solution and that the half-filled ground state is a gapped antiferromagnetic singlet at any U 6= 0 value. The gap only vanishes at U = 0. There is no MIT transition, and the double occupancy D ceases to behave as an order parameter, reaching zero only at its r → ∞ limit. No analytic solutions are known in 2D or 3D, but simulations 35,55,56 show that in these cases MITs develop at well defined r values. It has also been shown 57 that, although not obvious from the Bethe Ansatz, all densitydensity correlations, i.e. the hni nj i expectation values for any two sites i, j, decay exponentially with intersite distance, as expected for an insulator. However, since in the U = 0 limit the Hubbard solution falls into a metallic state with algebraic decay and oscillations in its DIs, it remains interesting to examine the behavior of this model and compare it with

14


Page 14 of 32

Page 15 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


standard quantum chemical results.

3

Results

In this Section some results concerning the decay of DIs in Hubbard and real hydrogen chains are considered. To that end, Hubbard cyclic chains for lattices of increasing size with 4m + 2, m = 1, 2, . . . sites have been exactly solved using PBC, examining their evolution with the correlating parameter r. The open (no-PBC) m = 0 two-site chain has also been included. 1D chains show very clear bond length alternation (precursors of Peierls distortions, as already shown 15 ), so to simplify as much as possible, only aromatic-like cycles are considered. Hubbard calculations have been performed using Lanczos diagonalization within the SPINPACK code. 58 As it will be shown, DIs (or density-density correlations) get quickly saturated towards the LW limit. The same H-cycles at different H-H distances a have been solved in regular polygon geometries with standard quantum chemical tools. Results at the HF, unrestricted HF (UHF), valence CAS, singles and doubles configuration interaction (CISD), and FCI levels using the GAMESS 59 code will be shown. Focusing on qualitative rather than on quantitative reasoning, simple 6-311+G, or even STO-3G basis sets in the FCI calculations, have been used. Delocalization measures have been obtained from QTAIM space partitionings, using the PROMOLDEN 60 code. QTAIM indices are known to be rather basis set independent. 61 Succinctly, LIs and DIs are obtained after expanding ρxc in the orbital basis. 62,63 With this, as Eqs. 1 and 2 show, our measures are reconstructed from R overlap integrals restricted to real space atomic domains, A dr φi (r)φj (r). These are called

atomic overlap matrices (AOM), and are numerically obtained. They are available for many

computer codes, 64 including PROMOLDEN. DIs or LIs are then obtained from the AOMs either with PROMOLDEN or equivalently with the EDF 52 program.

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3.1

Saturation rate

The existence of the LW analytic solution allows to study the speed at which the quantum chemical or finite Hubbard chain results approach the infinite limit. There has been previous interest in the literature about the evolution of LIs or DIs in H-chains as regards their possible behavior as order parameters in MITs. For instance, Baranov and Kohout 10,65 showed that symmetry broken (unrestricted KS, UKS) DFT LIs experience an abrupt transition from low to high LI as the lattice parameter of the infinite H-chain increases. Since many models, like the Gutzwiller Ansatz, 66 show a MIT in the 1D H chain, these results might simply signal that UKS displays a MIT transition, but do not help much in showing the proper behavior of the LI as r changes. Fig. 2 shows how size extensive valence CAS DIs, as well as finite Hubbard DIs and LI, behave with the lattice parameter a and 4r = U/t, respectively, in cyclic chains with 4m + 2 nodes (the linear m = 0 case is also included in the upper diagram). Do not confuse the correlation strength r with the lattice parameter a. The r(a) relation has been obtained by numerically matching the exact diagonalization δ results obtained for the m = 1, 2, 3 finite Hubbard chains to the CAS δ’s. The saturation is very quick, at any a or r value, but only after comparing with the m = ∞ LW limit does it become clear that this rapid convergence is not an artifact of finite size models. Notice that it is possible to compare with analytic results at U/t = 0. In this case, δ 1,2 = 0.444, 0.419, 0.413, 0.409 for m = 1, . . . , 4, respectively, and the n → ∞ value tends to 4/π 2 ≈ 0.405. All the CAS results show an inflection point at about a ≈ 3.6 bohr that is absent in the DI versus U/t Hubbard plots, much as shown in the H2 case. To understand the origin of the sigmoid from the CAS and Hubbard data a r(a) map is generated. This has been done before with different models, 55 and as the middle plot in Fig. 2 reveals, the mapping is non-linear and quickly convergent. It appears that r grows slowly at small internuclear distances and faster as the distance increases. Excellent fits can be obtained in our limited a window to exponential forms like r = αexp(−β/a). The r(a) non-linearity explains the sigmoidal shape of DIs as follows. A DI will show 16


Page 16 of 32

Page 17 of 32

0.90 m=0 m=1 m=2

0.80 0.70

δ1,2

0.60 0.50 0.40 0.30 0.20 0.10 0.00 1

2

3

4

5 6 a/bohr

0.50

7

8

9

6

m=1 m=2 m=3

0.40

4

a(U/t)

δ(U/t)

δ

a

1,2

0.30

10

0.20 2 0.10 0.00

0 0

2

4

6

8

10

U/t

0.50 m=1 m=2 m=3 m=∞ LW

0.40 0.30 λ

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.20 0.10 0.00 0

2

4

6

8

10

U/t

Figure 2: Saturation of LIs and nearest neighbour DIs with size for 4m + 2 sites cyclic chains. The upper diagram shows the evolution of DIs in valence CAS//6-311+G H-chains and the middle and bottom pictures those of DIs and LIs for Hubbard chains, respectively. The middle plot also shows a mapping of the a, U/t relation, that should be read on the right scale, obtained from matching the CAS and Hubbard results. The bottom diagram contains 17bohr. also the m → ∞ LW LIs. All distances are in ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 32

a zero second derivative with respect to a due to the opposite signs of both terms in the following expression, d2 δ = da2

dr da

2

d2 δ dδ d2 r + ≡ ((r′ )2 δ ′′ + δ ′ r′′ ), dr2 dr da2

(6)

so that the inflection point condition will be met when δ ′′ /δ ′ = −r′′ /(r′ )2 . From the presented data, this occurs at about U/t ≈ 4 or r ≈ 1, when the nearest neighbor DI is very close to half its maximum value (at r = 0). This is an important insight, since it was previously shown 13,67 that important chemical processes occur close to points where δ ≈ δ max /2. The saturation data show that the convergence of actual QTAIM LIs or DIs with n is similar to that shown by Hubbard models towards the LW solution. We are thus in a position to state that the evolution of either the LI or the nearest neighbor DI in the infinite H chain will be sigmoidal with a, due to the non-linear r, a map smoothly decaying with r, and that the discontinuity in the LI shown by UKS data 65 points to a fictitious MIT due to symmetry breaking.

3.2

Decay rate of DI(1,j)

The inter-site decay rate of DIs is a much less explored subject with possibly relevant repercussions on chemical thinking. Let us label the nodes of the cyclic chains from 1 to n and recall from a previous contribution 15 that in the Hückel, TB, or Hubbard at U = 0 models, δ i,j shows an oscillatory pattern, vanishing whenever i + j is even. This is a result of exact interference cancellation, that may be successfully interpreted in terms of Pauling resonance structures. In fact, 15 many important chemical effects, like mesomerism in conjugated systems, i.e. the alternating propagation along the change of substituent effects, rely on this oscillatory propagation of electron delocalization. In the n → ∞ limit, the envelopes of all non-vanishing DIs decay with an inverse square law, δ 1,j → 4/ (π 2 (j − 1)2 ). This signals the closure of the gap and the transition to a metallic extended system. It was also shown in

18


Page 19 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


that work that, although the long-range decay of the DIs for finite (gapped) chains might be exponential, their behavior for intermediate j values converges quite rapidly with n towards the inverse square law. Exact vanishing of the odd j DIs is only observed if no direct through-space overlap between the nodes is allowed, as in the Hückel/TB cases. Real HF/KS calculations show that the oscillations persist in real systems when the geometries (a values) lie in a range where mean-field calculations are accurate. However, exact cancellation conditions are now not met, and odd j DIs do not vanish, although their values are considerably smaller than their neighboring even j partners. This is, in our opinion, sufficient to warrant mesomeric effects: non-monotonous decay of influence (delocalization) along the chain. Interestingly, as it has been already reported, 15 exact cancellations do also occur in 2D and 3D Hückel/TB models, but seem to disappear in periodic single-determinant KS computations. A possible rationalization lies in the larger number of short-range overlaps in 2D and 3D cases, which destructively work against interference and oscillations. The effect of electron correlation on this picture is not well known. On the one hand, in finite systems, correlation does not influence the presence of a gap, so the long-range decay of δ should be exponential. On the other, the exact, true thermodynamic limit of the H chain is not available, but the LW solution also shows long-range exponential decay at any r. This does not preclude interesting chemistry at mid-range distances. By naïvely interpolating these ideas, a possible scenario appears in which, at large a or r values DIs do not show any anomaly, decaying exponentially such that δ 1,i > δ 1,j if i < j. For small a or r parameters, however, a transition towards an oscillatory pattern, converging on the j-parity rule, should be found. Fig. 3 shows the decay of δ 1,j with the intersite H-H distance for a H14 cyclic chain with a = 1.84 bohr for several levels of theory. Each type of result has also been fitted to a polynomial decay law to show the shift in tendencies. As expected, once the zero j values have been removed from the Hückel results, we get the slowest decreasing model, with an

19



1

δ1, j = (j−1)−fδ1,2 0.1

δ1, j

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 32

0.01

0.001

0.0001

Hückel f = 1.43 HF f = 2.88 UHF f = 2.88 CAS f = 4.35 CISD f = 4.09 FullCI f = 4.20

2

4

8

j Figure 3: Evolution of δ 1,j with j, equivalent to the intersite distance for the H14 cycle at different levels of theory, with a = 1.84 bohr. Notice the double logarithmic scale. A linear fit of each set of data, with explicit description of the f exponent is also shown, so that δ 1,j ≈ (j − 1)−f δ 1,2 . Be aware that the fit is done with both high and low values. Odd j δ’s vanish at the Hückel level, so only even j’s are included in that case. Some unreliable data due to loss of accuracy in small magnitude numerical integrations at j = 8 have not been included. The UHF data are indistinguishable from the HF ones. All calculations performed with 6-311+G basis sets except the full CI one, done with the STO-3G basis set.

20


Page 21 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


exponent f ≈ 1.43. This tends to exactly 2, the theoretical value in the n → ∞ limit. As it can be clearly seen, all the other theory levels predict oscillations compatible with mesomeric or medium-range delocalization effects, much as in a metallic-like system, but with all odd j DIs not vanishing now. Among the realistic models, the mean-field HF solution is the one with slowest decay, and the main effect of including electron correlation is a general decrease of all DIs. Interestingly, the situation is quite similar to that already examined in the H2 molecule, where the DI converges toward the exact value in a damped oscillatory pattern as the amount of correlation considered increases. As exemplified here, we have found that in this regime UHF calculations reproduce, essentially, the HF landscape. Inclusion of only static correlation at the CAS level has a large impact on DIs. As f tells us, the CAS data decay fastest of all. Addition of dynamical correlation at the CISD level dampens the CAS decrease and consideration of all possible excitations at the FCI level (even with a minimal basis) again reacts back slightly. In any case, at this internuclear distance where the HF approximation is still reasonable, introduction of correlation seems to alter quantitatively, but not qualitatively, the oscillatory decay. Fig. 4 shows two examples (with 6 and 14 sites) of the effect of the increase of r on the oscillatory pattern of DIs in the Hubbard chains used here. The behavior exemplified is general. At sufficiently large r, the canonical decay expected for an insulator is recovered and δ 1,j decreases on increasing j. Decreasing r takes us to a point where a first crossing is met. This signals the onset of quantum mechanical interference in electron delocalization and once this barrier is crossed, the δ 1,j order is no longer canonical. Successive crossings end up with the Hückel alternation, with vanishing odd j DIs. Interestingly enough, the first crossing between DIs on decreasing r occurs again, independently of m, very close to r = 1 or to U/t = 4. As expected, the Hubbard chains show a very interesting transition. Whatever the number of sites, the DIs go from an asymptotic regime with large r, decaying smoothly with intersite distance, to a regime at small r, where interference effects give way to a non-monotonous

21



0.45

j=2 j=3 j=4

0.40 0.35

δ

1,j

0.30 0.25 0.20 0.15 0.10 0.05 0.00 0

2

4

6

8

10

U/t

0.06

j=2 3 4 5 6 7 8

0.05

1,j

0.04 δ

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 32

0.03 0.02 0.01 0.00 0

2

4

6

8

10

U/t Figure 4: Decay of δ 1,j for 4m + 2 cyclic Hubbard lattices with m = 1, 3 versus U/t.

22


Page 23 of 32

oscillatory decay. Since these oscillations are one of the fingerprints of mesomeric effects, we firmly think that electron correlation plays against it. This adds to other studies showing that the classical rules used in light element chemistry may not hold at all with heavy element analogues as electron correlation strengthens its role. 68

1

U/t=0.25 U/t=2 U/t=4 U/t=6 U/t=8

0.1

1,j

0.01

δ

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.001

0.0001

1e-05

1e-06 2

3

4 j

5

6

7

8

Figure 5: Evolution of the decay rate of δ 1,j for the 14 sites cyclic Hubbard chain with U/t. Notice the doubly logarithmic scale. The transition from the oscillatory to the monotonous decay regime in the 14 sites chain is also shown in Fig. 5. Not surprisingly, r ≈ 1 (or U/t ≈ 4) is found to be the approximate boundary between the two situations. At U/t = 8 delocalization is hampered so much that the asymptotic infinite n limit has been fully entered, which is found to be exponential. 57 How this image readjusts in real H chains is the next focus. Fig. 6 shows the quantum chemical results on the H10 chains, at several levels of theory and at varying a values. Most of the relevant points regarding the behavior at small H-H distances have already been commented. The parity rule (odd-even j alternation) is obeyed at all theory levels, with damped CAS, CISD, and FCI alternation. 23



1 j=4 0.1

j=2

δ1, j

0.01

j=3

0.001 f = 1.38 f = 2.41 f = 3.33 f = 2.80 f = 3.45

0.0001 1e−05 1e−06 2

4

8

16

1

δ1, j = (j−1)

−f 1,2

δ

0.1

δ1, j

0.01 0.001 Hückel f = 1.38 HF f = 2.52 CAS f = 4.12 CISD f = 2.62 FullCI f = 4.69

0.0001 1e−05 1e−06 2

4

8

16

8

16

1 0.1 0.01

δ1, j

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 32

0.001 f = 1.38 f = 2.89 f = 7.54 f = 2.72 f = 9.09

0.0001 1e−05 1e−06 2

4

j Figure 6: Decay of δ 1,j for the H10 ring at three H-H distances: 2.7 (top), 3.5 (middle) and 24 5.5 (bottom) a.u., and various ACS levels of theory, are only labeled in the middle plot. Paragon Plus which Environment Notice that the fit includes all the data. All calculations performed with 6-311+G basis sets except the full CI ones, done with the STO-3G basis set .

Page 25 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


As a increases, the mean-field, single-determinant HF model becomes unable to localize the electrons in their atomic domains (or sites), and the DI continues oscillating. At a = 3.5 bohr, close to r = 1 according to the r, a mapping, the parity alternation of DIs no longer holds for both the CAS and the FCI descriptions, but continues at the CISD level. At larger distances all signs of oscillatory behavior disappear for the CAS and FCI, and continue at the HF and CISD levels. A first obvious conclusion is immediate. In order to qualitatively follow the evolution of electron delocalization, size consistency is a must and the CISD method soon ceases to behave properly. A second important point is that, once the consistency issue is assumed, quantum chemical and Hubbard data are surprisingly compatible, pointing towards the generality of our results. Although density-density correlations are known to decay exponentially in the long range regime at any r 6= 0 for the Hubbard infinite chain, 57 this by no means excludes interesting phenomena at short- or mid-range. We think to have convincingly shown that this is the case and that delocalization of electrons in real chains will suffer a transition from oscillatory to monotonous decay as correlation is increased.

4

Conclusions, perspectives

We have examined in this work the spatial decay rate of real space delocalization measures as electron correlation increases. After showing that delocalization indices (DIs) decay according to inverse power laws in metallic-like cases and exponentially in insulating-like moieties 15 for single-determinant and Hückel models, we go one step further, examining the validity of those insights when electron correlation is switched on. The consideration of electron correlation as a variable opens the possibility of examining transitions between insulating and metallic phases upon simple geometrical rearrangements. These Mott insulating transitions (MIT) play an important role in contemporary condensed matter physics. Correlation also poses important problems for standard quantum chemical

25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

procedures, since MITs are collective phenomena (involving many particles) in the strong correlation limit (that needs from considering a combinatorial number of Slater determinants). It is therefore important to combine standard calculations with models of strong correlation to ascertain the reliability of the results. By using the standard Hubbard hamiltonian at half-filling, the amount of correlation can be controlled by a simple parameter r. The systems analyzed here are restricted to simple finite and infinite H rings, examined under several standard quantum chemical theoretical levels and the Hubbard approximation. The latter admits an analytic solution in the infinite 1D chain, as found by Lieb and Wu 36 who showed that no MIT exists and that the system is an insulator for any r 6= 0. A first insight is that the standard order parameter used to detect MITs, the double occupancy probability of a site, that vanishes at the MIT, is directly related to the localization index (LI). In fact, the latter is a generalization that might be used in general systems. LIs and DIs show a quasi-universal sigmoidal behavior with intersite distance in bond breaking processes 13 . It is shown that the origin of this sigmoid lies in the non-linear relation between inter-site distance and the r correlating parameter, and that its inflection point is found at r ≈ 1. It has also become clear that DIs and LIs converge quickly with chain size. We can confidently state that the LI in the H chain will decay smoothly to zero as r or a increases and that the jump shown in a previous work 65 is due to symmetry breaking. Although our conclusions are based on real space indices computed with modest basis sets, and DIs depend weakly on the latter, further analyses should be done to check the role of residual basis set errors. Although DIs must decay exponentially at long range for the infinite Hubbard chain, light has been cast on the behavior of the metallic-like oscillations shown in mean-field models, that survive at short and midrange when r is small, to further disappear as r increases. At about r = 1 they vanish and a clear exponential behavior sets in. Due to the link between DIs and core chemical concepts like bond order, bond conjugation, mesomerism, etc, these results

26


Page 26 of 32

Page 27 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


imply a possible wide impact in chemistry. Conjugation and mesomerism, for instance, will not survive in strongly correlated analogues of alternant hydrocarbons. Extension to further dimensions and/or systems should be done. As soon as the 1D chain is abandoned in favor of 2- or 3D analogues, true MIT will appear. After the analyses set forth, it comes out that the variance of the electron population at a given site (VAR, related to the LI) may be used as a general order parameter in real simulations. At a MIT, the VAR will show a jump, but will not vanish exactly, as it happens in the Hubbard model (where, for instance, VAR is null in the triplet state of H2 ). VAR will remain non-zero when computed quantum chemically, due to exponentially decreasing, but non-vanishing overlap of same spin electrons onto neighboring atomic domains. Finally, the results shown here have implications about the profile of DIs in chemical processes. Inflection points in those cases where sigmoidal shapes appear mark the switch from short to long range delocalization. At smaller a, r values the regime is bonded, covalentlike, while at larger ones it turns to the non-bonded, exponential overlap region. Given the a, r map, increased correlation will suppress covalency. This has been put forward repeatedly, but stems very neatly from our results.

Acknowledgement AGB thanks the FICYT for a Ph. D. grant that has made this research possible. The authors also thank the Spanish MINECO, grant CTQ2012-31174, for financial support.

References (1) Kohn, W. Phys. Rev. 1964, 133, A171. (2) Gimarc, B. M. Molecular structure and bonding. The qualitative molecular orbital approach; Academic Press: New York, 1979. (3) Popelier, P. L. A.; Brèmond, E. A. Int. J. Quant. Chem. 2009, 109, 2542. 27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(4) Bader, R. F. W. Atoms in Molecules; Oxford University Press: Oxford, 1990. (5) Matta, C. F., Boyd„ Eds. The Quantum Theory of Atoms in Molecules: From Solid State to DNA and Drug Design; Wiley-VC, 2007. (6) Bader, R. F. W.; Stephens, M. E. Chem. Phys. Lett. 1974, 26, 445. (7) Mundim, K. C.; Giambiagi, M.; de Giambiagi, M. S. J. Phys. Chem. 1994, 98, 6118– 6119. (8) Giambiagi, M.; de Giambiagi, M. S.; Mundim, K. C. Struct. Chem. 1990, 1, 423. (9) Francisco, E.; Martín Pendás, A.; García-Revilla, M.; Álvarez Boto, R. Comput. Theor. Chem. 2013, 1003, 71. (10) Baranov, A. I.; Kohout, M. J. Comput. Chem. 2011, 32, 2064. (11) Wiberg, K. B. Tetrahedron 1968, 24, 1083. (12) Mayer, I. Chem. Phys. Lett. 1983, 270, 97. (13) García-Revilla, M.; Popelier, P. L. A.; Francisco, E.; Martín Pendás, A. J. Chem. Theory Comput. 2011, 7, 1704. (14) Resta, R. Phys. Rev. Lett. 1998, 80, 1800. (15) Gallo-Bueno, A.; Martín Pendás, A. Phys. Chem. Chem. Phys. 2016, 18, 11772. (16) Taraskin, S. N.; Fry, P. A.; Zhang, X.; Drabold, D. A.; Elliot, S. R. Phys. Rev. B 2002, 66, 233101. (17) Taraskin, S. N.; Drabold, D. A.; Elliot, S. R. Phys. Rev. Lett. 2002, 88, 196405. (18) Mott, N. F. Rev. Mod. Phys. 1968, 40, 677. (19) Ruzsinszky, A.; Perdew, J. P.; Csonka, G. I.; Scuseria, G. E.; VYdrov, O. A. Phys. Rev. A 2008, 77, 060502R. 28


Page 28 of 32

Page 29 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(20) Bartlett, R. J.; Musial, M. Rev. Mod. Phys. 2007, 79, 291. (21) White, S. R.; Martin, R. L. J. Chem. Phys. 1999, 110, 4127. (22) Marti, K. H.; Bauer, B.; Reiher, M.; Verstraete, F. New J. Phys. 2010, 12, 103008. (23) Knizia, G.; Chan, G. Phys. Rev. Lett. 2012, 109, 186404. (24) Knizia, G.; Chan, G. J. Chem. Theory Comput. 2013, 9, 1428. (25) Mazziotti, D. Phys. Rev. Lett. 2008, 101, 253002. (26) Tsuchimochi, T.; Scuseria, G. E. J. Chem. Phys. 2009, 131, 121102. (27) Shi, H.; zhang, S. Phys. Rev. B 2013, 88, 125132. (28) Maier, T.; Jarrel, M.; Pruschke, T.; Hettler, M. Rev. Mod. Phys. 2005, 77, 1027. (29) Lathiotakis, N. N.; Helbig, N.; Gross, E. Phys. Rev. B 2007, 75, 195120. (30) Henderson, T. M.; Bulik, I. W.; Stein, T.; Scuseria, G. E. J. Chem. Phys. 2014, 141, 224104. (31) Limacher, P. A.; Ayers, P. W.; Johnson, P. A.; De Baerdemacker, S.; Van Neck, D.; Bultinck, P. J. Chem. Theory Comput. 2013, 9, 1394. (32) Cirac, J. I.; Zoller, P. Nat. Phys. 2012, 8, 264. (33) Cocchi, R.; Miller, L. A.; Drewes, J. H.; Koschorreck, M.; Pertot, D.; Brennecke, F.; Köhl, M. Phys. Rev. Lett. 2016, 116, 175301. (34) Hubbard, J. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 1963, 276, 238. (35) Baeriswyl, D., Campbell, D. K., Carmelo, J. M. P., Louis, E., Eds. The Hubbard Model: Its Physics and Mathematical Physics; Springer US, 1995. 29



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(36) Lieb, E. H.; Wu, F. Y. Phys. Rev. Lett. 1968, 20, 1445. (37) Essler, F. H. L.; Frahm, H.; Göhmann, F.; Klümper, A.; Korepin, V. The OneDimensional Hubbard Model ; Cambridge, Cambridge., 2005. (38) LeBlanc, J. P. F.; et al., Phys. Rev. X 2015, 5, 041041. (39) Schäfer, T.; et al., Phys. Rev. B 2015, 91, 12509. (40) Sinitskyi, A. V.; Greenman, L.; Mazziotti, D. A. J. Chem. Phys. 2010, 91, 014104. (41) Chen, Y.-H.; Tao, H.-S.; D-X., Y.; Liu, W.-M. Phys. Rev. Lett. 2012, 108, 246402. (42) Rubin, N. C.; Mazziotti, D. A. J. Phys. Chem. C 2015, 119, 14706. (43) Fradera, X.; Austen, M. A.; Bader, R. F. W. J. Phys. Chem. A 1999, 103, 304–314. (44) Matito, E.; Poater, J.; Solà, M.; Duran, M.; Salvador, P. J. Phys. Chem. A 2005, 109, 9904. (45) Martín Pendás, A.; Blanco, M. A.; Francisco, E. J. Comput. Chem. 2007, 28, 161. (46) Francisco, E.; Martín Pendás, A.; Blanco, M. A. J. Chem. Theory Comput. 2006, 2, 90. (47) Juhász, T.; Mazziotti, D. A. J. Chem. Phys. 2006, 125, 174105. (48) Skolnik, J. T.; Mazziotti, D. A. Phys. Rev. A. 2013, 88, 032517. (49) Ziesche, P. In Many-electron densities and reduced density matrices; Cioslowski, J., Ed.; Kluwer Academic Publishers, 2000; Chapter 3, p 33. (50) Francisco, E.; Martín Pendás, A.; Blanco, M. A. J. Chem. Phys. 2007, 126, 094102. (51) Martín Pendás, A.; Francisco, E.; Blanco, M. A. J. Chem. Phys. 2007, 127, 144103. (52) Francisco, E.; Martín Pendás, A. Comput. Phys. Commun. 2014, 185, 2663 – 2682. 30


Page 30 of 32

Page 31 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(53) Martín Pendás, A.; Francisco, E.; Blanco, M. A. Phys. Chem. Chem. Phys. 2007, 9, 1087. (54) Rubin, N. C.; Mazziotti, D. A. Theor. Chem. Acc. 2014, 133, 1492. (55) Kurzyk, J.; Wojcik, W.; Spalek, J. Eur. Phys. J. B 2008, 66, 385. (56) Georges, A.; Kotliar, G.; Krauth, W.; Rozenberg, M. J. Rev. Mod. Phys. 1996, 68, 13. (57) Essler, F. H. L.; Frahm, H. Phys. Rev. B 1999, 60, 8540. (58) Schulenburg, J. program package SPINPACK. 2001. (59) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S. J.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. J. Comput. Chem. 1993, 14, 1347. (60) Martín Pendás, A.; Francisco, E. A QTAIM/IQA code (Available from the authors upon request). (61) Jablonski, M.; Palusiak, M. J. Phys. Chem. A. 2010, 114, 12498. (62) Martín Pendás, A.; Blanco, M. A.; Francisco, E. J. Chem. Phys. 2004, 120, 4581. (63) Martín Pendás, A.; Francisco, E.; Blanco, M. A. J. Comput. Chem. 2005, 26, 344. (64) Keith, T. A. The AIMAll program. 2015; The code is avalaible at http://aim. tkgristmill.com. (65) Baranov, A. I.; Kohout, M. Acta Crystallogr. A 2011, 67, C115. (66) Guztwiller, M. C. Phys. Rev. Lett. 1963, 10, 159. (67) Martín Pendás, A.; Francisco, E.; Blanco, M. A. Faraday Discuss. 2007, 135, 423. (68) Frenking, G., Shaik, S., Eds. The Chemical Bond: Fundamental Aspects of Chemical Bonding; Wiley-VC, 2014. 31


Correlation hampers resonance effects

Journal of Chemical Theory and Computation Page 32 of 32

1

U/t=0.25

~1,j Bond order

1,j

δ

Os Low U Mesocillatory /t meri decay sm s ets in

U/t=4

0.1

U/t grows. Correlation increases

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

0.01

U/t=8 H1

0.001

0.0001

H14 1e-05

1e-06 2

Hj

Large U/t Exponential regime No Mesomeric effects

ACS Paragon Plus Environment 3 4 5

j

6

7

8

Decay Rate of Correlated Real-Space Delocalization Measures

Recommend Documents