Correlation Effects in Scanning Tunneling Microscopy Images of

Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Correlation effects in STM images of molecules revealed by quantum Monte Carlo Matteo Barborini, Sandro Sorella, Massimo Rontani, and Stefano Corni J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00710 • Publication Date (Web): 06 Oct 2016 Downloaded from http://pubs.acs.org on October 8, 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 34

Journal of Chemical Theory and Computation

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 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

we visualize the electronic correlation embedded in the simulated STM images, highlighting the many-body features that might be observed.

1

Introduction

Scanning tunneling microscopy (STM) 1–5 and spectroscopy (STS) are able to gain respectively topological and spectroscopic information regarding the local density of states of extended surfaces 6–10 and nano-objects, such as quantum dots 11–22 and single molecules 23–33 . To probe the intrinsic molecular states, it is essential that the molecules are decoupled from the substrate, either through physisorption 27,34,35 or by depositing an insulating overlayer (for example NaCl) between the conductive substrate and the molecule 24–26,28 . With proper energy resolution and ideal shape 36 of the STM tip, the probability of the electron addition (or removal) is ruled by the square modulus of the wave function of the electron (or hole) quasiparticle added to (or removed from) the interacting system 12,13,37–41 . The quasiparticle wave function (QPWF), also known as Dyson orbital 39,42,43 , experiences the whole many-body correlation arising from the interaction of the added (removed) electron with the other electrons already present in the molecule. For this reason, the QPWF is also used to interpret the photoelectron spectroscopy (PE) data 42 , although the perturbations and time scales associated with the electron removal are quite different from those of STM. Whereas the QPWF square moduli provided by STS experiments are usually interpreted in terms of molecular orbitals, as obtained for example from density functional theory (DFT) calculations 25,34,44–48 , several physical phenomena could be explained only taking into account the full effect of electronic correlations 11–13,38,41,49 . This is the case e.g. of the Kondo ground state, which exhibits the collective behaviour of electrons including those in the conductive substrate and STM tip 49–52 . In this paper we focus on a distinct, recently predicted few-body effect inherent to the isolated molecule – a sub-nanometer distortion of the QPWF with respect to its uncorrelated shape induced by strong electron correlation 39 .

2

ACS Paragon Plus Environment

Page 2 of 34

Page 3 of 34

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


In a recent work 40 by some of us the QPWF has been built from state-of-the-art coupled cluster and configuration interaction calculations with single and double excitations (CCSD and CISD, respectively), which partially take into account the effects of electronic correlation. By comparing the correlated STM images based on these QPWFs with those obtained from unrestricted Hartree-Fock (UHF) orbitals, it was shown that certain metallorganic molecules exhibit significant distorsions of the single-particle orbitals induced by few-body correlations. These orbital alterations are robust against molecule-substrate interaction (NaCl) and large enough to be possibly measured by comparing the images of two similar metallic complexes (Cu- and Zn-saloph), only one of the twos being sensitive to correlation. Unfortunately, the CCSD and CISD methods suffer from a notoriously Ne6 scaling of the computational cost as a function of the number of electrons, Ne . Besides, these methods might not be accurate enough to properly treat correlation, as suggested by the frequent requirement of the more expensive CCSD(T) method to reliably predict molecular properties. To overcome the limitations of truncated multideterminantal techniques, in this paper we compute the QPWFs using quantum Monte Carlo (QMC) methods 53–55 , which are stochastic algorithms to evaluate averages of physical observables over correlated many-body wave functions. The first advantage of QMC lies in its computational implementation: The algorithms scale at most as the fourth power of Ne and are highly parallelizable without loss of efficiency, overcoming the large computational prefactor and making the description of sizeable systems feasible 56 . The second advantage is that, to represent electronic states, it is possible to use quite complex variational wave functions with both a multideterminantal fermionic part and a Jastrow factor. The latter is a bosonic term, namely symmetric under particle permutations, which is able to describe many-body correlation by coupling the electronic variables 57 . By using advanced wave function 58–60 and structural 61–63 optimization procedures it is nowadays possible to study through QMC several electronic 64–68 and structural properties of molecular compounds 69–71 within chemical accuracy. In the following, by means of efficient and accurate QMC methods, we define a procedure

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 ACS Paragon Plus Environment

Page 4 of 34

Page 5 of 34

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


paper is the reformulation of the QPWF in a form suitable for evaluation via QMC, which is illustrated in section 3.1.

2

Quasiparticle wave function

Low-temperature STS allows to measure the spectral density of states, N (r, ω). For a holelike excitation the spectral density is: 2 X ˆ N (r, ω) = h ¯ hNe − 1, ν|ψ(r)|Ne , 0i δ(¯hω − E0 (Ne ) + Eν (Ne − 1)),

(1)

ν

where hNe − 1, ν| is the νth excited state of the system with Ne − 1 electrons, which includes the molecule and possibly the substrate, |Ne , 0i is the reference ground state, Eν (Ne − 1) ˆ and E0 (Ne ) are the energies respectively associated with these two states, and ψ(r) is the annihilation operator that destroys an electron in the position r. Throughout the paper the vector r includes the spin index, unless explicitly stated otherwise. If the difference between the energy of the ground state, E0 (Ne − 1), and that of the first excited state, E1 (Ne − 1), of the molecule is larger then the energy resolution of the STM tip for a small transport 2 ˆ window centered at the Fermi energy, then eq 1 reduces to N (r, ω)dω = hNe − 1|ψ(r)|Ne i ,

where only the ground states hNe − 1| and |Ne i are left.

As discussed in refs 39,40, and previously proposed by Rontani and Molinari for quasitwo-dimensional quantum dots 37 , the spectral density is thus proportional to the square modulus of the quasiparticle wave function (QPWF) that for a hole-like excitation is defined as: ˆ ϕhQPWF (r) = hNe − 1|ψ(r)|N e i.

(2)

Together with this hole-QPWF (hQPWF), we also define the electron quasiparticle wave function (eQPWF) as the probability amplitude of finding an electron, added to the molec-

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

Page 6 of 34

ular ground state by the field operator ψˆ+ (r), in the position r: ϕeQPWF (r) = hNe + 1|ψˆ+ (r)|Ne i.

(3)

STS experiments aim at linking the measured current flowing through the microscope tip to the square modulus of the QPWF using different setups. For example, in the so called ‘constant current’ operational mode the probing tip scans both the planar location and apparent height of the molecule while keeping the tunneling current constant through a feedback loop. Therefore, the isosurfaces of the QPWF square modulus signify the apparent height of the molecule as seen by the tip, i.e., the vertical distance from the reference plane at which the probability weight assumes a certain fixed value 40 .

3

Quantum Monte Carlo compliant formulation of QPWFs

Quantum Monte Carlo (QMC) methods 53,54 are stochastic algorithms developed to evaluate the mean values of physical observables over complex trial wave functions, ΨT (¯r, {¯ α}), that ˆ e , which depend are approximate guesses of the eigenfunctions of the electronic Hamiltonian H on a set of variational parameters {¯ α}. The idea is to calculate the observables as mean values of local quantities, estimated for a set of uncorrelated electronic configurations ¯r, which are extracted according to a probability proportional to the modulus square of ΨT (¯r, {¯ α}). In order to establish a QMC procedure to compute the QPWFs defined in eqs 2 and 3, we have to evaluate the following integral expression of the eQPWF (see the Supporting Information for the derivation):

ϕeQPWF (r) =

s

Ne + 1 NNe NNe +1

Z

d¯r(Ne ) ΨNe +1 (¯r(Ne ) , r)ΨNe (¯r(Ne ) ),

6


(4)

Page 7 of 34

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


where NNe =

R

d¯r(Ne ) |ΨNe (¯r(Ne ) )|2 and NNe +1 =

R

d¯r(Ne +1) |ΨNe +1 (¯r(Ne +1) )|2 are the two

normalization factors of the many-body wave functions, and ¯r(Ne +1) = {r1 , . . . , rNe , rNe +1 } and ¯r(Ne ) = {r1 , . . . , rNe } are the electronic configuration vectors containing the electronic coordinates. Unfortunately, these normalization factors cannot be directly calculated using QMC. Here we propose a workaround based on the ‘correlated sampling technique’ 73,74 . To proceed we rewrite eq 4 as the integral of a local function of the electronic configurations, multiplied by a density probability related to the square modulus of one of the two wave 1/2

functions. This can be done by multiplying and dividing the integrand by ΨNe (¯r(Ne ) )/NNe so that it becomes the product of two terms,

ϕeQPWF (r) = Q

Z

d¯r(Ne )

ΨNe +1 (¯r(Ne ) , r) ΠNe (¯r(Ne ) ), ΨNe (¯r(Ne ) )

(5)

which are respectively the ratio between the two wave functions, and the probability density ΠNe (¯r(Ne ) ) = |ΨNe (¯r(Ne ) )|2 /NNe associated with the ΨNe (¯rNe ) many-body wave function. The constant Q that appears in the equation is the ratio between the two normalization factors: Q=

s

(Ne + 1)

N Ne . NNe +1

(6)

In this way the integral 5 is evaluated as the mean value of the ratio ΨNe +1 (¯r(Ne ) , r)/ΨNe (¯r(Ne ) ) calculated for each electronic configuration, sampled according to the probability ΠNe (¯r(Ne ) ) through a random walk:

ϕeQPWF (r) = Q

ΨNe +1 (¯r(Ne ) , r) ΨNe (¯r(Ne ) )

.

(7)

Π Ne

The wave function ΨNe that appears in the ratio is calculated on the same ¯r(Ne ) electronic configuration generated during the stochastic sampling, whereas the ΨNe +1 is calculated adding to the Ne electrons another electron at the position r, where we want to evaluate the QPWF. In this way, we estimate the QPWF using the QMC sampling of the |ΨNe |2 in 7



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 8 of 34

the whole configurational space with a relatively low error bar. Note that the still unknown normalization factors of the QMC wave functions are hidden in the prefactor Q, which is dealt with in the next section.

3.1

The normalization prefactor Q

To estimate the normalization prefactor Q appearing in eq 7 we combine the independent samplings of both ground states with Ne and Ne + 1 electrons, as explained in this section. Consider the square modulus of the eQPWF, Ne +1 NNe NNe +1

|ϕeQPWF (r)|2 =

R

×

R

d¯r(Ne ) ΨNe +1 (¯r(Ne ) , r)ΨNe (¯r(Ne ) ) ×

(8)

d¯r(Ne +1) ΨNe (¯r(Ne ) )ΨNe +1 (¯r(Ne ) , rNe +1 )δ(r − rNe +1 ),

which may be written as the product of two integrals to be computed independently. The first integral I1 over Ne electronic positions is evaluated with the same procedure described above, i.e., by sampling the probability associated with the square modulus of ΨNe (¯r(Ne ) ). The other integral I2 is dealt with by multiplying and dividing by ΨNe +1 (¯r(Ne +1) ) and including the NNe +1 normalization factor, so it is replaced by the integral I2 =

Z

d¯r(Ne +1)

ΨNe (¯r(Ne ) )δ(r − rNe +1 ) ΠNe +1 (¯r(Ne +1) ), ΨNe +1 (¯r(Ne ) , rNe +1 )

(9)

which is computed stochastically by sampling a certain number of electronic configurations ¯r(Ne +1) according to the probability density ΠNe +1 (¯rNe +1 ) =

I2 =

|ΨNe +1 (¯ r(Ne +1) )|2 : NNe +1

✘ rN✘ ΨNe (r1 , . . . , rNe , ✘ e +1 )δ(r − rNe +1 ) ΨNe +1 (¯r(Ne +1) )

.

(10)

ΠNe +1

The ratio between the two wave functions is evaluated each time the Ne + 1 electron of the ¯r(Ne +1) configuration vector is in r; specifically, ΨNe +1 is estimated using all the electrons in the configuration whereas ΨNe is computed on the ¯r(Ne ) configuration obtained by removing

8


Page 9 of 34

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


the rNe +1 electron once it is located in r. Furthermore, this expression is obviously generalized considering that each one of the Ne + 1 electrons may be located in r: 1 Ne + 1

*N +1 + e X ΨNe (r1 , . . . , ri , . . . , rNe +1 )δ(r − ri ) ΨNe +1 (¯r(Ne +1) ) i=1

.

(11)

ΠNe +1

Here the ratio between the two wave functions is evaluated each time the ith electron of the ¯r(Ne +1) configuration satisfies the condition δ(r − ri ). In conclusion, the square modulus of the eQPWF is the product

|ϕeQPWF (r)|2 =

ΨNe +1 (¯r(Ne ) , r) ΨNe (¯r(Ne ) )

Π Ne

*N +1 + e X ΨNe (r1 , . . . , ri , . . . , rNe +1 )δ(r − ri ) ΨNe +1 (¯r(Ne +1) ) i=1

ΠNe +1

(12)

between two terms, each one being separately evaluated through a distinct QMC integration. Comparing eq 12 with eq 7 we obtain that the prefactor Q is a constant function of r: Q2 = Q2 (r) =

DP

ri ,...,rNe +1 )δ(r−ri ) Ne +1 ΨNe (r1 ,...,✚ i=1 ΨNe +1 (¯ r(Ne +1) )

D

ΨNe +1 (¯ r(Ne ) ,r) r(Ne ) ) ΨNe (¯

E

E

ΠNe +1

=

N (r) . D(r)

(13)

Π Ne

Here some comments are in order. Firstly, N (r) and D(r) are functions which only differ in their amplitudes (by a factor Q2 ) and have the same nodes. For this reason in those values of r in which N (r) = D(r) ≡ 0 the ratio in eq 13 is indeterminate. Secondly, to evaluate the numerator N (r), one would need to sample configurations satisfying exactly the δ(r − ri ) condition, which occurs with almost vanishing probability. In order to overcome these difficulties there are two possible ways to proceed in the evaluation of Q. In the first approach, since the prefactor Q must not depend on r and N (r) = Q2 D(r), R R also the V N (r)dr = Q2 V D(r)dr equivalence holds, leaving us with the possibility of

evaluating the constant as:

R N (r)dr . Q = RV D(r)dr V 2

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

Page 10 of 34

By integrating N (r) and D(r) over the same arbitrary volume V we obtain a meaningful statistical expression where both quantities represent random variables with finite variance. R † This is formally equivalent to considering a quasiparticle operator Z † = ψr dr that creates r∈V

a particle in the volume V. Note however that only the variable r is limited to the region

V, while the other electronic variables are not limited in such region. Consequently, within this scheme we introduce a discretization error in evaluating the denominator D(r) in eq 13 due to the finite mesh of the grid inside the volume V. In the second approach, we calculate the denominator only in one point r (for which D(r) 6= 0) and the numerator’s mean value N (r) in the volume V around r. Unfortunately, as V → 0 the number of stochastically extracted configurations that actually have an electron in V decreases, enlarging the stochastic fluctuations of the numerator and consequently of Q. In order to compare the two approaches we will apply both of them to lithium atom in Section 5.

4

Jastrow antisymmetrized geminal power

The trial wave functions used in quantum Monte Carlo calculations are written as products of a fermionic part and a Jastrow factor 53,55 , which is a bosonic term. In this work we use the Jastrow antisymmetrized geminal power 75 (JAGP) wave function already applied to study the structural and electronic properties of various molecular systems 56,75–77 . The particularity of this wave function is that its fermionic part, the antisymmetrized geminal power 78 (AGP), although rather compact is a constraint multideterminantal expansion that includes various molecular excitations, as described in refs 67,70,76. For this reason, it has been shown that even at the variational Monte Carlo level the JAGP is able to correctly describe various molecular properties recovering a high level of electronic correlation 62,66–68,70,71 . For molecular systems of Ne electrons in a spin singlet state, i.e. Ne /2 = Ne↑ = Ne↓ , the AGP is

10


Page 11 of 34

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


written as the antisymmetrized product, 

ΨAGP (¯ x) = Aˆ 

Ne /2

Y i=1



ΦG x↑i ; x↓i  ,

(14)

of two-electron geminal pairing functions in a spin singlet state (Aˆ is the antisymmetrization operator): 1 ΦG (xi ; xj ) = φG (ri , rj ) √ |↑ii |↓ij − |↑ij |↓ii . 2

(15)

The spatial part φG (ri , rj ) of the pairing functions,

φG (ri , rj ) =

M X X

λµa νb ψµa (ri ) ψνb (rj ) ,

(16)

a,b=1 µ,ν

is a linear combination of products of two atomic orbitals multiplied by the λµa νb coefficients, that can be seen as weights which modulate the superposition of two atomic orbitals centred on the ath and bth atoms and in which the indices µ and ν represent the quantum numbers (n, l, lz ) of the orbitals. Here M is the number of atoms in the molecular compound. In order to describe spin polarized systems, Ne↑ > Ne↓ , it is necessary to generalize eq 14 introducing S = Ne↑ −Ne↓ single electron wave functions in the antisymmetrized product 79,80 . In this case, we obtain the general expression 

ΨGAGP (¯ x) = Aˆ 

↓

Ne Y i=1

ΦG x↑i ; x↓i

S Y s=1



Φs x↑N ↓ +s  ,

(17)

s = 1, . . . , S

(18)

with the single electron wave functions

Φs (xi ) =

M X X a=1

µ

λsµa ψµa (ri ) | ↑i

being essentially independent “molecular orbitals” built as linear combinations of the atomic orbitals in the basis set. 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

Page 12 of 34

The bosonic term that appears in the JAGP, i.e. the Jastrow factor 57,62,76 , does not depend directly on the spin coordinates and is written as the product of three terms, J = J1 J2 J3/4 , the first of which is the one body Jastrow factor: PM,N ¯ ¯r = e a,i e {−(2Za )3/4 ξ((2Za )1/4 ria )+Ξa (ri )} , J1 R,

(19)

where ria = |ri − Ra | is the distance between the ith electron and the ath nucleus, Za is the atomic charge. This factor includes both the homogeneous interaction between the electron and the nuclei through the function ξ (r) = B2 1 − e−r/B , and the non homogeneous term P built as a linear combination of atomic orbitals Ξa (ri ) = µa gµa χµa (ri ). The second term

is the purely homogeneous two body factor

J2 (¯r) = exp

(

Ne X

)

ξ (rij ) ,

i

Correlation Effects in Scanning Tunneling Microscopy Images of

Recommend Documents