Theoretical Studies of Electron Transfer and Electron Localization at

J. Phys. Chem. B 1997, 101, 2459-2475

2459

Theoretical Studies of Electron Transfer and Electron Localization at the Semiconductor-Liquid Interface B. B. Smith* and A. J. Nozik National Renewable Energy Laboratory, 1617 Cole BouleVard, Golden, Colorado 80401 ReceiVed: September 30, 1996; In Final Form: January 7, 1997X

We present new models and simulations of electron transfer (ET) at the semiconductor-liquid interface (SLI). The simulations are of a “first principles” molecular dynamics type and therefore incorporate electronic structure calculations. An In(H2O)62+/3+ redox species, water, and InP semiconductor system are focused on. We discuss the problem of electron localization at this interface, especially as it relates to ET. The study allows the mechanism of the ET process to be analyzed. Rate constant calculations are performed with the dynamics of the entire system incorporated. We apparently present the first calculation of electron transfer coupling matrix elements for the SLI.

1. Introduction

2. An Analytical Model for Semiconductor-Liquid ET

The semiconductor-liquid interface (SLI) is an extremely challenging electron transfer (ET) system both experimentally and theoretically. The experiments are difficult to control due to problems in preparing clean, stoichiometric, and defect-free semiconductor surfaces. Contradictory results for nominally identical systems are prevalent. Theoretically, simple models like those possible with metal systems are lacking. For nearsurface processes semiconductors tend to require sophisticated treatments. Calculations are often very sensitive to surface structure, but predicting surface reconstruction geometry is extraordinarily difficult. Interfacing a solvent with this problematic solid system by forming an SLI makes matters even more difficult. Experimentally, SLI energetics and microscopic structure are in general inadequately understood, since it is often impossible to use many of the techniques of experimental surface science for the SLI. Despite the difficulties, one can develop microscopic models for the SLI. Recently, in ref 1, we developed a SLI ET simulation that accounted microscopically for the full electronic structure of the semiconductor, redox species, and solvent. The determination of the electronic coupling of the redox molecule with the semiconductor was inherent in the simulation. Surface states, chemisorption, surface defects, and other near-surface phenomena are addressable. In ref 1 we also generalized an Anderson Hamiltonian based electron transfer model (ABETT) to the SLI. The ABETT we generalized in ref 1 was one advanced several years ago by one of us (B.B.S.) and J. T. Hynes for electron transfer (ET) at metal-liquid interfaces;2 the latter model2 built on foundations laid by Schmickler,3 Sebastian,4 and others. All these models are generalizations of the well-known Anderson Hamiltonian approach to chemisorption.5,6 Here we further enhance and analyze both the SLI model and simulation and report new results and insights into the mechanism of SLI ET. At the end of the paper we discuss some key experiments to do to investigate the theoretical issues we explore here. The ABETT model serves to introduce the simulation and to illustrate other points below. We discuss it next.

In our ABETT model1 we assume the redox species has a single s-orbital to (from) which electron transfer occurs from (to) the semiconductor. The electronic occupation of this orbital monitors the ET process as in refs 1-4. As eq 16 will show, the occupation of this orbital can be related to the ground state electronic energy of the system and also to the solvent and semiconductor configuration which, in turn, determines the nuclear energy. These electronic and nuclear energy components will determine the total energy of the system analogously to ref 2. Free energy curves similar to those of conventional ET theory7,8 can be constructed from these electronic and nuclear contributions, and the ET rate constant can then be calculated using these free energy curves.2 The formalism is best suited for electronically adiabatic transfer.1,2 The total Hamiltonian (H) for the model system just discussed can be divided into parts

X

Abstract published in AdVance ACS Abstracts, March 15, 1997.

S1089-5647(96)03003-9 CCC: $14.00

H ) He + Hsolv ) Heredox + Hesemi + Heint + Hsolv

(1)

where He is the electronic part of the Hamiltonian and Hsolv is the portion having to do with the solvent. As implied by eq 1, He consists of redox species and semiconductor terms, as well as a term describing the interaction between the redox species and semiconductor. We assume in our ABETT model that Hsolv is of the wellknown quadratic form2 2

HS ) PS /(2M0) + 1/2kS(∆E)2

(2)

where kS is the system force constant2 defined by eq 3 and ∆E is the difference in the energy of interaction of the redox species (oxidized and reduced forms) with the solvent and semiconductor at fixed solvent and semiconductor nuclear configuration. ∆E is directly analogous to the solvent coordinate used widely in modern homogeneous ET theory.9,10 By employing eq 2 we are assuming symmetric free energy curves (of equal and parabolic curvature) and thus2,11

kS ) kBT〈(∆E)2〉-1 ) (2Er)-1

(3)

where kB is Boltzmann’s constant and Er is the reorganization energy. The system mass (M0) can be expressed as2,11 © 1997 American Chemical Society

2460 J. Phys. Chem. B, Vol. 101, No. 14, 1997

kBT/M0 ) 〈(d∆E/dt)2〉

Smith and Nozik

where t is time. For the electronic Hamiltonian (He in eq 1) we will use a spin-dependent Anderson Hamiltonian. We write this Hamiltonian as

He ) naσ(∆E + σ) +

∑k

{knkσ + (Vakc†aσckσ + Vkac†kσcaσ)} (5)

where σ denotes the electron spin, “a” denotes the redox species state, and “k” refers to the electron states of the semiconductor. σ is the effective redox species energy level of spin σ, and k is the energy of the electrons in the semiconductor. More specifically,

σ ) a + U〈na-σ〉

(6)

a is the energy of the redox species orbital in vacuum, which is essentially equal to the orbital ionization energy. U is the Coulomb integral magnitude used to model the on-site Coulomb repulsion (represented by the second term of eq 6) of the redox species level. The niσ are the occupation number operators, and the c†iσ and ciσ the creation and annihilation operators. The electrode-redox species coupling is governed by the matrix elements Vak and Vka. Note the model Hamiltonian (eq 5) differs from the original Anderson-Newns Hamiltonian5,6 in that the effective energy eff of the redox species (the adatom of the Anderson-Newns model) is not σ alone, but rather:

eff ) ∆E + σ

(7)

This is the conventional ABETT prescription, and it models how the effective energy of the redox species orbital is affected by the solvent fluctuations.1-4 ∆E’s value changes with the polarization of the solvent (and semiconductor fluctuations). When ∆E attains its equilibrium value for the reduced species, then eff is the effective energy of the reduced species orbital equilibrated to the solvent and semiconductor. When ∆E fluctuates to the oxidized value we obtain the eff for the oxidized species. The breakdown of eq 5 into the three He terms of eq 1 is obvious. One can write the (projected) local density of states of the redox species as1-6

F(,∆E) ) RLDOS )

|〈mσ|a〉|2δ(-mσ) ∑ m

(8)

For later purposes we call F(,∆E) the redox species local density of states (RLDOS) (in our simulation the summation for F is over orbitals, and all the many redox species orbitals must also be summed1). The integral of F over energy () at fixed ∆E gives the average occupation, or the ground state occupation number, 〈n(∆E)〉, of the single orbital of our ABETT model’s redox species at that ∆E:

〈n(∆E)〉 )

∫f() F(,∆E) d

(9)

f() is the Fermi-Dirac distribution function. Actually, for the model we are using, 〈n(∆E)〉 must in general be calculated as1

〈n(∆E)〉 )

∫f() F(,∆E) d + 〈naσ(∆E)〉L

〈naσ(∆E)〉L ) [1 - dΛ(L)/d]-1

(4)

(10)

so as to include localized redox species states, where 〈naσ(∆E)〉L is the average occupation number of localized redox species states. As follows directly from the full Anderson model5,6

(11)

and the localized energies L are the roots of the equation:

L - σ - Λ(L) ) 0

(12)

The 〈n(∆E)〉 of eq 10, which can be rewritten as 〈naσ(∆E)〉, must be solved for self-consistently when spin is included so that 〈naσ(∆E)〉 ) 〈na-σ(∆E)〉 (cf. eqs 10 and 6).1 The average occupation of the orbital 〈n(∆E)〉 is, of course, directly related to the formal charge of the redox species and thus indicates whether the redox molecule is oxidized, reduced, or is in some intermediate charge state. The Anderson Hamiltonian gives F(,∆E) as5,6

F(,∆E) ) π-1[∆()/(( - ∆E - σ - Λ())2 + ∆2())] (13) where

∑k |Vak|2δ(-k)

(14)

∫

(15)

∆() ) π and

Λ() ) π-1P d′ [∆(′)/( - ′)]

∆() is directly related to the energy width a discrete state acquires on its interaction with a continuum of states (in this case, the discrete state of the redox species interacts with the continuum of states in the solid and, potentially, the solvent), Λ() is directly related to the energy shift a discrete state acquires with such an interaction. The broadening (∆) and shift (Λ) are manifested in the 3D simulations below, although the shape and behavior of F are typically much more complicated than in the simple one-dimensional ABETT model we introduce below. We discussed this topic in detail in ref 1. Using the Hellman-Feynman theorem, the ground state electronic energy for the system is2,4

E0(∆E) )

∫〈n(∆E)〉 d∆E

(16)

In refs 1 and 2 we showed that the total free energy of our model system can be written as

F(∆E) ) 1/2kS(∆E)2 + E0(∆E)

(17)

where the first term stems from eq 2 and E0(∆E) is eq 16. Using eq 17 we can construct free energy curves that are of the typical form assumed in traditional electronically adiabatic ET theory. However, unlike in traditional theories, the electronic coupling that determines the barrier shape and height is calculated directly from a microscopic formulation which takes detailed account of the interaction of the redox species with the continuum of states in the solid and gives E0(∆E). These dependencies can have very pronounced influences on ET rates and characteristics in ways not included in traditional theory.2-4 One of our main tasks in extending the ABETT model of ref 2 to the semiconductor-liquid interface was to evaluate for this new system the quantities ∆() and Λ(). In ref 1 we did this by using a one-dimensional model (note the model as given in eqs 1-17 is general and can address 3D systems). What this accomplished was to have an explicit analytical model, directly relatable to molecular orbital theory, for the electronic coupling and the electronic structure of the semiconductor and redox species. We considered a one-dimensional chain of atoms coupled, within the tight-binding approximation, only to their

Electron Transfer at the Semiconductor-Liquid Interface

J. Phys. Chem. B, Vol. 101, No. 14, 1997 2461

Figure 1. Schematic diagram of the 1D model described in the text (in this figure Hii is abbreviated as Hi, and Hij as Hj).

Figure 2. Simple band structure predicted by the 1D model.

nearest neighbor by a hopping matrix element Hij (this is the same matrix element which in the tight-binding hydrogen molecule ion model describes the frequency with which an electron hops back and forth between the H atoms, ω ∼ Hij/h; Hij has the same fundamental physical interpretation here). We considered only one other matrix element, the on-iste matrix element Hii. Hii is essentially the atomic orbital energy (Coulomb integral) and is often modeled as the isolated orbital ionization energy.12 An end atom of the chain is the model single-orbital redox species and it has its own Hii and Hij, which we call Hii′ and Hij′. Hii′ is directly related to σ of eq 7. The electronic structure of the solvent is left out in this model. The 1D chain of atoms is shown schematically in Figure 1, and the reader will recognize the direct relationship of the model to simple Huckel molecular orbital theory. To mimic onedimension the s-p bonding found in compound semiconductors, we assume the sign alternates for the Hii and Hij in the fashion shown in Figure 1. This model has been applied to compound semiconductors.13,14 Note two simple bands are produced with this model (see Figure 2) with band widths and a band gap that are functions of Hii and Hij (see Appendix A). This simple 1D compound semiconductor model yields rather complicated formulas for ∆() and Λ() given in Appendix A. Appendix A shows that ∆() and Λ() are functions of Hij′, Hii′, Hii, and Hij. Thus ∆() and Λ() depend on the electronic structure of the semiconductor and redox species. Using these models we find (see Appendix A) that F(,∆E) also contains the same functionality, as well as solvent dependence (via ∆E), i.e.,

F(,∆E) ) f(,Hij′,Hii,Hij,a≈Hii′,∆E)

(18)

The important point for this paper is that the electron transfer, as mapped by 〈n(∆E)〉, depends on the semiconductor’s and redox species’ electronic structure, and the key influence of the solvent is that its polarization shifts eff via eq 7 and eff changes 〈n(∆E)〉 via eqs 10 and 13. The assumptions embodied in eqs 1-18 and the physical reasons for the variation of the redox species’ effective energy via eq 7 were discussed in ref 1. But it is worthwhile to briefly reiterate. As the dipolar solvent undergoes its thermal fluctuations, the electric field the redox species experiences changes. This change in field changes the energy of the redox species orbital. If the redox orbital is in quantum communication with the solid (if wave function overlap is sufficient and of the proper symmetry), electron density can potentially be exchanged. If

the energy of the redox species orbital is above (below) that of the semiconductor conduction (valence) band then electron density can potentially be exchanged with the solid (this argument is modified for surface states in the gap in obvious ways). The fluctuation in the solvent and semiconductor configuration must be such that energy is conserved in the ET process (in the absence of phonon-assisted hopping). As electron density is transferred to the redox species, 〈n(∆E)〉 changes and therefore the charge on the redox species changes, which in turn affects the average polarization of the solvent ∆E. The electron density transferred to the redox species screens the nuclei in the redox species and thus also effects the redox species orbital energy. The charge on the redox species, its electronic occupation, the orbital energy, and the solvent fluctuations are thus all mutually coupled. Equations 1-18 represent a simplistic model of this complex interrelationship which is modeled in great detail in the simulation below. Using eqs 10, 13, 16, and 17 and equations from Appendix A, we can calculate the free energy surface F(∆E) which will yield the rate constant k. The calculation of k is possible using transition state theory (TST) as follows. In TST, the rate constant is the equilibrium averaged one-way flux across the barrier top:2

kTST ) 〈j+〉R ) 〈(PS/MS)θ(PS)δ[δ∆E‡]〉R

(19)

The subscript R means that the average is normalized with respect to the reactant state and θ is the Heaviside step function. Combining this with a consideration of the extrema of the full Hamiltonian, one obtains:2

kTST ) (ωR/2π) exp(-∆F‡/kBT); ωR ) (M0-1(∂2F/∂∆E2)∆ER)1/2 (20) The adiabatic activation free energy is obviously:

∆F‡ ) F(∆E‡) - F(∆ER)

(21)

∆ER is the value of ∆E corresponding to the reactant well minimum (the average ∆E for the reactant), and ∆E‡ is the corresponding value for the barrier top. The reactant well frequency ωR is (from eqs 3, 17, and 20):2

ωR ) M0-1/2[E0′′(∆ER) + (1/(2Er))]1/2

(22)

E0′′ denotes the second derivative of E0 with respect to ∆E. We show results in Figure 3 for 〈n(∆E)〉, calculated using eqs 10, 13, A4, and A5, for the 1D model compound semiconductor and redox species, using Hii ) 1.9, Hij ) 0.7 eV, and H′ij listed in the caption. If we use a typically assumed value for the reorganization energy (Er) of 1 eV, we obtain for F(∆E) representative results like that shown in Figure 4 (where we used eqs 10, 13, 17, A4, and A5). The activation energy ∆F‡ for the forward and reverse reactions can be calculated numerically from curves like those in Figure 4 by utilizing eq 21. To account for the displacement of the reactant and product minima (which is in part due to the difference in solvation energies of reactant and product, see below), we simply vary a + ∆E as a parameter (this difference is calculated in detail in the simulation below). By employing simulation results, Er can be calculated using eq 3 and M0 by using eq 4, or Er and M0 can be estimated in other ways. Then the k can be calculated from eq 20 and compared to k from the simulation (see below). The total, as opposed to local, k calculation is described by eq 27 below, and the local k from eq 20 can be used in eq 27. In


Smith and Nozik 3. The Simulation (Background)

Figure 3. 〈n(∆E)〉 as calculated by the 1D model for Hii ) 1.9 eV, Hij ) 0.7 eV, and U ) 0. All curves are actually plotted against ∆E + a, the choice of the zero is obviously arbitrary, and this convention allows easy extension of the results to arbitrary a values. (a) Results for Hij′ ) 0.8, 1.0, and 2.0 eV. The curves can be identified by the fact that they approach 1 faster with ∆E as H′ decreases. (b) The solid curve is the same as the Hij′ ) 2 eV curve of (a). The dashed curve shows the localized contributions to 〈n(∆E)〉, via eq 11.

Figure 4. Free energy curves calculated using the 1D model. The results of Figure 3 (a) were used to construct these curves. Curves were shifted to center on zero. (a) The barrier height (measured from the well minima to the center peak) decreases monotonically as Hij′ increases. Results for Hij′ ) 0.8, 1, and 2 eV are shown. (b) Two nonequilibrium curves (reaction free energy * 0) obtained by varying the vertical shift (the two thin-lined curves) and compared to an equilibrium curve (heavy-lined curve). All curves calculated using Hij′ ) 2 eV.

ref 1 we showed that while our ABETT model qualitatively matched our full simulation results, our comparison of the electronic density of states of the two methods1 clearly shows the ABETT’s simplicity.

Having described an analytical model of the system which sets the general physical background, we now describe the simulations and their results. Recently we presented a simulation of SLI ET.1 There we coupled electronic structure calculations to molecular dynamics simulations in two ways to examine the semiconductor-liquid ET system. The same general methods were used here, except the second technique was augmented. These two techniques are as follows: (1) In the first method we simply calculated the electronic structure of the configurations of the system generated by conventional molecular dynamics (MD). Our conventional MD used internuclear potentials which were a fixed function of internuclear distance, as is typical. Details of this technique are described in ref 1. Here we primarily used potentials that were approximate fits to the interatomic potentials calculated using the PM3 Hamiltonian.15 This simulation method was mainly used as a simple check of the second method which employed PM3 in a different way. We also compared results from conventional MD simulations using potential functions constructed from combined ab initio, semiempirical, and empirical techniques. (2) Fixed interatomic potential functions, as used in conventional MD, are not adequate to accurately describe certain important aspects of the ET system. The electronic structure should, ideally, be calculated self-consistently during the simulation by allowing the electronic wave functions to adjust to the changing nuclear configurations. This simulation technique is typically called “first principles” molecular dynamics. One of the most rigorous self-consistent or first principles MD methods is the Car-Parinello scheme.16 We use a “first principles” technique here which we call molecular orbital molecular dynamics (MO-MD) since we use semiempirical electronic structure techniques incorporating a molecular orbital picture. Like other such methods, MO-MD results in potentials that are much more complicated than a simple fixed function of internuclear distance. Moreover, MO-MD calculates the actual transfer of electron density. Conventional MD does not do this at all; to do so requires an electronic structure calculation like that incorporated in MO-MD. The MO-MD electronic structure calculations also allow bond formation and bond dissolution to be computed so that chemisorption can be treated. Let us use the simple ABETT model above to explain some key features of the MO-MD technique. If we have a sophisticated enough Hamiltonian we can perform a simulation. How would the ABETT Hamiltonian (eq 1) have to be modified to perform an actual simulation? First, in a realistic atomic orbital based SLI model many orbitals for the redox species must be considered, and the same is true for the semiconductor and solvent. We would also have to add the dependence of the Hii and Hij on interatomic geometry for each atomic orbital of the basis set. The repulsive internuclear potential must also be added. In short, the Hamiltonian must actually be able to calculate quantities like ∆E from nuclear and electronic contributions, rather than simply parameterizing the Hamiltonian with such quantities, as in eq 2. Once such a Hamiltonian is constructed, the simulation can be performed as follows. Initial coordinates, for a full 3D system, and momenta must be assigned. Then an electronic structure calculation must be done which calculates the system energy and thus forces via

forcei ) -∂E/∂xi

(23)

where E is the system energy and xi is an arbitrary spatial atomic

Electron Transfer at the Semiconductor-Liquid Interface coordinate. The system can then be propagated a short time, dt, with these calculated forces. Then using the new configuration, another electronic structure calculation is done and so on. This is the same technique used in all first principles MD. The description just given captures the essential features of MO-MD,1 regardless of which semiempirical electronic structure method we might use. In ref 1 we primarily used the extended Huckel method supplemented by including nuclear repulsion and self-consistent techniques.1 Here we used the PM3 Hamiltonian, from the MNDO family of Hamiltonians,17 to calculate the electronic structure. The PM3 method was compared with results of both EHT and ab initio results (the latter only on limited configurations for details of electronic structure and forces). It has previously been found that even the simple EHT method can give reasonable chemisorption geometries if appropriate precautions are taken, and EHT has been used to generate the water-surface potential for conVentional metalliquid system MD.22 We reported EHT-based electronic density of states (DOS) in ref 1, and showed they acceptably matched established results for InP and water. MNDO related methods have in general had better success than EHT for well-known reasons. Namely, all MNDO methods are self-consistent field methods which include electrostatic repulsion and exchange. Conventional EHT does not include these factors, nor is it selfconsistent, and these flaws can only be corrected by rather ad hoc procedures.1 While MNDO typically yields results superior to EHT, it is also true that MNDO remains a research tool regarding solid state systems. It is well-known that parametric adjustments are necessary to bring MNDO into agreement with certain features of the electronic DOS and other quantities for certain systems. In all cases, before settling on specific PM3 parameterizations, we checked the results with our ab initio calculations in the same fashion reported in ref 1. We found little adjustment necessary from the conventional MOPAC 93 values, and we intentionally kept adjustments to a minimum; the standard MOPAC 93 parameters satisfactorily reproduce results reported here. To avoid extreme computational demands, we did not use a full scale electronic structure calculation. This method is discussed in detail in ref 1. Here we note that in our MO-MD simulations the interatomic potentials only for two layers of the solid (more layers were used in test runs) and a varying number of waters (typically less than 20, up to 40 in test runs) nearest the redox species and surface were calculated via electronic structure methods during the simulation. To save CPU time, the rest of the interatomic potentials were modeled by otherwise fixed functions that approximately matched the results from PM3 calculations (e.g., the point charge water model of ref 1 adequately fits PM3 results and was employed). Thus, part of the solvent experienced a changing simulated InP surface potential, and the charge transfer and bond energetics could be monitored (this bounded “first principles” technique is similar to the subspace Hamiltonian technique18). However, the solvent modeled with the fixed interatomic potentials also experienced the changing “first principles” potential to a limited degree since we included its interaction with the changing charge distributions on atoms addressed by “first principles”. In turn, the “first principles” nucleii experienced the internuclear potentials of the “fixed potential” nucleii. To further reduce computational demands, the system was propagated with conventional MD (with PM3 potential fits) in instances where it could be ascertained that electronic effects were not pronounced. Ewald sums were used to mimic extended system electrostatics in conventional ways1,2,10,19 and Born von Karman boundary conditions20 (transversely to the surface) were used for the wave

J. Phys. Chem. B, Vol. 101, No. 14, 1997 2463 functions. Two walls were used. The front wall was the semiconductor slab, and the back wall was an uncharged dielectric continuum. The dimensions of the simulation box were about 15 Å × 15 Å (corresponding to a slab face of about 10 semiconductor atoms on a side) and about 25 Å in length. A total of 182 waters and two to eight layers of semiconductor atoms were used (a few large electronic structure calculations were done to check the size dependence of the results). Dangling bonds were passivated with hydrogen atoms for the 100 surface simulation. We will focus here almost exclusively on these 100 simulations, although we have used other faces and even tried reconstructive geometries. In one set of simulations the semiconductor atoms were fixed in their equilibrium positions, and in a second set they were allowed to move. The methods employed for the dynamics of the semiconductor atoms is discussed in Appendix B. To preface discussions in section 4, note that if semiconductor atoms are partially charged (either intrinsically or from other effects) the dynamical fluctuations of the lattice are roughly analogous to, and have a similar effect as, those of the solvent polarization fluctuations. The simulations that we will focus on here used water as the solvent, In(H2O)62+/3+ as the redox species, and the 100 face of InP. Hexaquo-coordinated indium is a nice redox species for theoretical studies because it is possible to exclude d-orbitals from consideration and it has orbital symmetries with the In of the semiconductor. Computer time can thus be saved. Other hexaquo-coordinated redox species, including transition metals, can be expected to yield similar theoretical results, and indeed we have found similarities with our previous Fe(H2O)62+/3+ studies.1 Since In forms so many hydrolysis products, it is not at all a straightforward species to deal with experimentally.21 We regard our studies of the In2+/3+ system as an illustrative theoretical system only. While in principle it would be possible to measure the In2+/3+ ET rate constant, the reaction would in most instances be an intermediate one in a series leading from In3+ to In+ or to In0. (Aquated In2+ is highly reactive22 and is often in practice formed by reaction of In3+ with hydrated electrons.) Further, In2+ may in some systems exist not as an aquo complex but as a complex composed in part or wholly of supporting electrolyte. It is important to recognize that the energetics of this interface, often the object of focus in experimental work, are determined directly from the simulation from the microscopic description of the atoms involved. In our simulation, if the redox species’ molecular orbitals participating in the transfer of electron density are energetically located in the gap and there are no surface states, no electron density is exchanged. One way in which exchange can occur is the following. Redox species oxidation, for example, can occur if the participating ET orbitals of the reduced redox species partner are elevated above the Fermi level and the conduction band edge by solvent and/or semiconductor fluctuations and if we are simultaneously able to reach the barrier peak of the adiabatic potential energy surface so that energy is conserved. In the case of phonon-assisted electron transfer (see discussions below), the latter condition can be relaxed. Reduction can occur in an analogous way. Of course, both reduction and oxidation can occur in each band. Depending on many factors (whether the semiconductor is p- or n-type, whether or not the system is light-excited, the size of the band gap, and other energetics of the system), the electron can transfer to or from one or all of the following: surface states, conduction band, valence band, and surface states within the bands. The simulation showed that our indium system energetics is such that participation by both the conduction and valence band is


Smith and Nozik

possible if we assume intrinsic InP 100 in dark conditions (without “photoexcited” states). But in general we find that the hydrogen-passivated termination promotes ET domination by the conduction band. Since the hydrogen-passivated InP 100 ET is more clear cut, we used that surface in all simulations reported below, except where noted. Energetic issues and comparisons to real world measurements are discussed in detail in Appendix G. 4. Simulation Results One of the quantities of interest from the simulation is the ET rate constant. To calculate it from the simulation, we will first assume that conventional transition state theory (TST), as it is applied to ET, can be employed. The standard expression, which stems from eq 19, for the local TST ET rate constant (by local we mean the rate constant calculated for an ion held at a fixed perpendicular distance from the electrode and above a particular site) is

k ) (ωRsim/2π) exp(-∆F‡sim/kBT)

(24)

Note that this expression is the same as eq 20, but we rewrite it with “sim” subscripts to emphasize that now the quantities ωR and ∆F‡ will be calculated from the full simulation. Once these quantities are calculated, k follows directly from eq 24. However, eq 24 is only applicable to electronically adiabatic reactions, i.e., those in which the ET can be envisioned as occurring on a single potential energy surface. This is, of course, the same limit to which our simulation technique is best suited. Some of the fastest rates of electron transfer typically occur in this adiabatic limit, and these are of interest in our experimental research efforts.23 The adiabatic limit is also of theoretical interest since the study of near surface phenomena (where adiabaticity is expected to occur) is virtually nonexistent for the SLI. We thus opt to explore this limit here. To establish when the electronically adiabatic assumption is appropriate for our system, we report our calculations of the electronic coupling matrix element (Γ) next. We show values for the electronic coupling matrix element Γ for our hexaquo-coordinated In ion and hydrogen-passivated 100 InP crystal face system in Figure 5. Figure 5 plots Γ as a function of center to center distance between the redox species In atom and a certain semiconductor surface atom (see figure caption). Note all curves fall off with distance roughly exponentially. Orientation of the redox species would appear to have a rather large effect on Γ, and some of the curves in Figure 5 are averages over configuration (see caption). We find that the effect of an applied field on Γ is rather small (see Figure 5), even for fields as large as 106 V/cm. However, these calculations do not account for the fact that the surface may well show bias dependent reconstruction which should significantly affect Γ (cf. later remarks). Further, applied fields can have a significant effect on other components of the rate constant, such as the free energetics. Also, our calculations are only on ground state systems. Γ is difficult to calculate for solid-liquid systems. In Appendix C we describe our methods for this system and our confidence in the results. No previous calculations of Γ appear to exist for the SLI. We want to try to assess, by using Figure 5, under what conditions electronic adiabaticity is a good assumption. The usual textbook two-state model24 used in conventional ET theory has a splitting equal to 2Γ between the lower and upper energy surface (see Figure 6). It is frequently assumed that if this splitting is sufficiently greater than kT, the upper energy surface is practically inaccessible, and the system can be thought of as

Figure 5. Plot of the electronic coupling matrix element (Γ) between the In(H2O)62+/3+ redox species and hydrogen-passivated 100 InP as a function of the center to center distance between the In atom of the redox species and a certain semiconductor surface atom, as detailed below. (Curve 1) Plots of Γ as a function of distance between a central P-terminated 100 InP and a hexaquo-coordinated In with its waters in Th symmetry and with the peak of the water octahedron pointing perpendicular to the surface above the P. (Curve 2) Same as curve 1 except the face of the water octahedron is parallel to the InP. (Curve 3) Same as curve 1 and 2 except averages over three redox species’ configurations from the MD are used. (Curve 4) Same as curves 1-3 except InP was In-terminated. (Curve 5) Same as curve 1 except a 106 V/cm field was applied, directed out of surface plane.

Figure 6. Schematic diagram of the textbook two-state model typically applied to electron transfer reactions.

moving on the lower energy surface only. Since the latter assumption defines adiabatic ET theory, large Γ suggests adiabaticity (a similar way to analyze adiabaticity is to compare electronic time scales {τelectronic = h/(2π‚2Γ)} versus nuclear ones). Another way to assess ET adiabaticity is to use the full Landau Zener (LZ) theory25 and see if the LZ transmission coefficient is near unity (the adiabatic limit) for a given Γ. A more exacting analysis, not possible with our technique here, would ask whether the coupling terms between nuclear kinetic energy operators were sufficiently small. We find from Figure 5 that at distances greater than about 7 Å 2Γ approaches kT and it is unlikely, using either of the first two assessment methods, that ET could any longer be adiabatic. Our calculations suggest that for adiabaticity in our particular system the redox species must be so close to the electrode that there is no intervening solvent, or at most one water. To maintain electronic adiabaticity without extremes of surface interaction, we situated the In of the redox species 5.3 Å from the electrode surface atoms (again, for reasons mentioned above we wanted to examine the adiabatic ET limit). From our calculations it is apparent that the ion can indeed reach these distances thermally. In fact, the calculations show that given sufficient time the species can fully chemisorb. We computed ∆F‡ by calculating the free energy curves as a function of ∆E. To do this we used a modification of a method that has become conventional in ET theory, namely, umbrella sampling employing the equation26,27



F(∆E) ) -kBT ln(P(∆E)) ) -kBT ×

∫

∫

ln{ dΩ exp(-βV)δ(∆E0-∆E)/ dΩ exp(-βV)} (25) where P is essentially the probability of finding configurations with a given numerical value of ∆E equal to ∆E0, V is the full system potential energy, and the relation -kBT ln(P(∆E)) is directly related to the well-known thermodynamic expression for free energy in terms of partition functions.28 The second half of eq 25 is the conventional expression for calculating P, from classical MD, by employing umbrella sampling.26,27 Ω represents all the system coordinates. The conventional method needs modification when applied to first principles MD, and we discuss this in Appendix D. Figure 7 shows the calculated free energy curve for the case where the semiconductor atoms were held fixed and the center of the indium redox ion was centered over and 5.3 Å from a P-terminating semiconductor atom. Note the adiabatic “hump” in Figure 7 (the barrier top) stems directly from the electronic structure calculation. In absence of electronic coupling of the two states, we would of course obtain two parabolic shaped diabatic free energy curves corresponding to the reduced and oxidized species free energy. As in homogenous theory, one of the key players in the vertical displacement of our minima (the difference between reactant and product average free energy) is the difference in the solvation free energies of the redox species. It is in the electronic coupling region of the curves (the region of the barrier peak) that the homogeneous analogy breaks down significantly because the redox species is electronically coupled to a completely different system, a semiconductor, not another redox species. The latter difference is studied in some detail below. The forward and reverse reaction activation energy for use in eq 24 is simply the difference between the reactant/product well minimum and the barrier peak. As mentioned, we did two kinds of simulations, one in which the semiconductor atoms were fixed and a second in which they vibrated. From Figure 7 it can be seen that for the fixed semiconductor atoms case ∆F‡ for the forward reaction was about 0.3 eV and the reverse ∆F‡ about 0.6 eV. For the two vibration conditions, the ∆F‡ deviated by about 20% at most (see Figure 7). Obviously the reorganization energy can also be calculated from the curves, directly by its definition.7,8 The reorganization energy is about 1.7 eV for the forward reaction and 1.5 eV for the reverse (roughly the same for both vibration conditions), as can be inferred from Figure 7. The curves are not perfectly quadratic in ∆E nor are they symmetric; conventional theory assumes both. ωR in eq 24 can be obtained in conventional ways. The easiest method is to employ the second derivative of F(∆E) weighted by the system mass M0:

ωR ) (M0-1(∂2F/∂∆E2)∆ER)1/2

(26)

Analogous methods have been used profusely in homogenous theory.2,11,27,29 The system mass M0 can be calculated using eq 4.2,11,27,29 For the fixed semiconductor atom conditions, ωR was about 620 ps-1 (for the forward reaction); the corresponding vibrating case ωR was about 60% lower. The results reported above show the effect of introducing lattice vibrations was not spectacularly large for ∆F‡ and ωR, and thus by eq 24 the kTST were not dramatically different. The largest kTST difference was for the forward reaction: an order of magnitude difference (lower for the vibrating system). The 20% difference in the ∆F‡ accounts for most of this difference, and 20% is perhaps comparable to the error in our vibrating

Figure 7. Free energy surfaces calculated from the MO-MD. The solid curve is for the fixed semiconductor atoms case. The dashed curve is the curve calculated with the dynamics of the semiconductor atoms included.

∆F‡ estimate, since the lattice vibrations were not modeled wholly from first principles (see Appendix C). If the surface was assumed to be reconstructed or had oxygen or other adsorbates, the effects would be different and potentially larger. In the semiconductor-liquid case the rate constant is conventionally assumed to be “bimolecular” in the concentration of redox molecules and electrons.30 If we follow these notions and wish to calculate the total electron flux through the surface to thus obtain the observed (total) rate constant, we have to form an equation like the following,1 similar to that first proposed by Levich7

ratei ) ki [C]0i [Ce]0i )

∫∫∫kL (z) Ci(z) Ce (,τ) dz dτ d i

i

(27)

for both i ) forward and i ) reverse ET reactions, where ki is the total ET rate constant, kLi is the local rate constant discussed above, z is the redox species’ perpendicular distance from the electrode, dτ represents integration over the whole volume of the semiconductor, [Ce]0i is typically taken to be the carrier (electron or hole) concentration near the semiconductor surface, [C]0i is typically taken to be the bulk redox species concentration, Ci(z) and Cei(,τ) are the corresponding calculated spatially dependent concentrations, and is energy (Cei(,τ) is actually a density of states per unit volume per unit energy). While this bimolecular assumption is theoretically unsettling,1 we are interested in reporting what our simulation predicts for the conventionally calculated rate constant. Due to the computational cost of solving eq 27 with our methods, we opt for finding an upper bound to this conventional rate constant by looking at the adiabatic transfer close to the electrode modeled above. To make this calculation we assume that only density in the first two layers of semiconductor atoms is important, as our simulations suggest for small redox species separation distances. We further only take into account flux from the redox species at 5.3 Å. In this case we get a rate constant of about 10-15 cm4/s for the cathodic current in the conduction band. We calculated this cathodic rate constant from the anodic one as follows. We first calculated the anodic rate constant by integrating over the unoccupied electronic DOS near the conduction band edge to give Cei(τ)/[Ce]0i (cf. eq 27). Conventionally, at equilibrium the cathodic conduction band current/ rate is equal to the anodic conduction band current (which is assumed to be independent of applied potential/band-bending)


Smith and Nozik

times exp(-eη/kBT); this last term accounts for the change in surface carrier concentration with band-bending, where η is the so-called SLI “overpotential”.30,31 Thus, if this conventional approach is taken, our cathodic rate constant follows from the anodic one. This calculation assumed fixed semiconductor atoms, and also assumed the redox species concentration had no z dependence (Ci(z) ) [C]0i). Details of this calculation are discussed in Appendix G, especially concerning how accurately the simulation results match various experimental measurements. Effect of Solvent and Semiconductor Dynamics. We now examine the impact of dynamics on the ET. We find that dynamical effects for our semiconductor-liquid ET system stem from three basic sources (1) the vibrations of the semiconductor lattice, (2) the changing conformations and vibrations of the redox molecule, and (3) solvent dynamics. Points 1 and 2 seem to have received little attention regarding the SLI; point 3 has been previously studied for metal-liquid systems32,33 but certainly not to the degree it has been for homogeneous systems. Studies addressing near surface solid-liquid electronic effects (like ours) are still extremely rare, and these effects may significantly influence all three effects. The physics behind all three dynamical effects will be discussed below, but for now we are interested in their collective effect on the rate constant. It is well-known that when the dynamics of a system is considered, the transition state theory rate constant is often no longer valid. Transition state theory assumes that every time the potential energy barrier maximum is reached, the system moves directly to products. However, when dynamical effects are considered, it turns out that recrossings of the barrier can occur, i.e., upon reaching the barrier maximum the system can move back toward the reactant state or reform the reactant(s). Thus, within this framework, kTST is an upper bound for the rate constant. In the simplest expressions that account for dynamical impacts, the effect on kTST is expressed as

k ) κkTST

(28)

where κ is a transmission coefficient which accounts for recrossing and is formally defined as the fraction of systems that reach the activated complex (transition state) which directly form products. Kramers was among the first to describe recrossing events as arising from frictional effects.34 Many have followed with a deluge of more elaborate stochastic frictional formulas.35 The effects of solvent dynamics is often conceived of as a solvent friction (on the reaction coordinate(s)). Affects of the electronic structure on dynamics is often modeled as an analogous electronic friction. To obtain κ from the simulation, we do not have to appeal to any model, like a frictional model. We will calculate κ directly in the fashion that it has been done in other simulations10,27 as follows. We first constrain the system to be in the transition state and let it come to equilibrium in this state. Then for random time steps in this equilibrium we save the momenta (p) and coordinates (r). We then run unconstrained trajectories of the system using our sampled p and r to start the trajectories. We run these trajectories backward in time to reach the reactant state and then forward in time from the transition state. Using this method we can watch the system evolve from the reactant state to the transition state and then see if the system ultimately produces products or moves back into the reactant well. We can thus see what fraction of those systems that reach the transition state produce products and then calculate κ by its definition above. Obviously, there will be a large number of these p and r sets that will correspond to the transition state, and for good statistics we want to sample many of these. Figure

Figure 8. ∆E trajectories calculated as described in text. All curves are for the fixed semiconductor atoms case.

TABLE 1 ion distance from surface (Å)

case

κ

5.3 5.3 7.0

fixed semiconductor atoms moving semiconductor atoms fixed semiconductor atoms

0.43 ( 0.03 0.05 ( 0.01 0.92 ( 0.04

8 shows the main output of this kind of simulation. ∆E is plotted as a function of time, and thus ∆E maps the system trajectory in time. Near ∆E ) 0 the system is in the transition state, when ∆E is roughly -2 eV the system is in the reactant state (oxidized redox species), and when ∆E is roughly +2 eV it is in the product state (reduced). Note that in Figure 8a the recrossing is pronounced: the transition state (∆E ) 0) being recrossed numerous times; while in Figure 8b there is essentially no recrossing, i.e., the TST assumption is met. In Figure 8c the reactant returns to the initial state; product production (ET) does not occur. The average of many such graphs/trajectories for given system conditions gives κ, after some statistical corrections (see Appendix E). A discussion of some issues involving this κ calculation appears in Appendix E. We calculated κ both for the condition of fixed and vibrating semiconductor atoms: (1) In the fixed case the κ was below the transition state value but not markedly so (see Table 1). In this case only the dynamics of the solvent and redox species affect κ (there is no semiconductor nuclear dynamics, and the effect of the dynamic electron semiconductor subsystem appears minimal). There was a clear trend with coupling Γ. As the redox species was constrained to be closer to the interface, and Γ necessarily increased, the transmission coefficient went down,

Electron Transfer at the Semiconductor-Liquid Interface as seen in Table 1. The barrier is broader as the coupling increases, and the solvent dynamics increases the recrossings of the barrier for well-known reasons.2,10,27 In exploratory calculations we found that the calculated κ also depends on the site over which the redox species was held fixed, and this again seems to be primarily due to the effect of changing Γ. (2) In the second case (vibrating semiconductor atoms) κ decreased by about an order of magnitude (see Table 1). While in this case Γ still plays a role in affecting κ since Γ changes as the atoms vibrate, it appears the main physics of the process is better described as stemming from dynamical disorder, as we discuss below. The dynamics study also allowed us to investigate the time scale of electron transfer which, e.g., allows investigation of whether it is possible for ET to occur before hot electrons thermalize. From graphs like those in Figure 8, it is apparent that the average time of ET, as measured by the time to go from reactant to product well, is on the order of picoseconds. For the systems we studied, this time spanned the range of 1-6 ps, with 1 ps values dominating. Mechanism and Localization. Our dynamics studies also permitted us to learn about the mechanism of ET in these systems by monitoring the electron localization. This mechanistic study will allow us to discuss some details of the three main dynamical effects mentioned above. To begin a discussion of mechanism, it is important to note certain features of the simulation. It is easiest to calculate the charge transfer for the oxidation of the redox species because the electron evolves from a spatially well-localized molecular state of the redox species. To calculate the reverse process (reduction), in principle, a highly delocalized “free” electron must be propagated from the semiconductor bulk to the surface, which is obviously a more difficult calculation. We find that, as long as the redox species is separated from the electrode by one solvent diameter or less, typically not only does the transferred electron evolve from a localized redox species state, but it also initially ends up in a spatially welllocalized site on the semiconductor surface. In this case the transferred electron density is typically localized near a single semiconductor atom close to the redox species. But with increasing separation between the surface and redox species, the transferred electron did spread its density over the slab. Figure 9 illustrates the kind of data we used to analyze localization phenomenon for our simulation. In Figure 9 one notes that the majority of transferred electron density for the seven-layer 100 InP slab used for this particular simulation is localized on the first two layers (the separation distance was intermediate; see caption). In still other cases the transfer was more delocalized over the slab. In any event, when we speak below of “localized” transfer in our simulations, it means transferred electron density localized essentially near one atom. All this behavior is in accord with investigations by Wolynes regarding the metal-liquid interface,36 which showed transfer was more localized the smaller the separation distance. Now the equations of motion for the system, be they classical or quantum mechanical, are symmetric under time reversal. This, of course, is the celebrated law of microscopic reversibility. Thus, with some possible reservations discussed below, we are led to conclude that the exchange of electrons in either direction (oxidation or reduction) can be from a spatially localized state on the semiconductor surface to a spatially localized state on the redox species. Of course, in the case of oxidation eventually this density must be transported through the bulk semiconductor. This transport is directly affected by the applied field, and eventually the initially highly localized density can be trans-


Figure 9. Atomic charges (Mulliken method) from the MO-MD simulation. The three-dimensional results are displayed in a quasi onedimensional format. Only representative atoms are displayed. The dashed lines connect open circle atomic charges before the oxidative ET, as indicated. The solid lines connect solid circle atomic charges immediately after a unit electron enters the semiconductor (later it is transported in). As can be seen, this particular initial transfer was not highly localized (with respect to a single atom), but it did localize on essentially the first layer. The In redox atom was 1 Å further out in the water than the 5.3 Å of the text. The surface was also In-terminated instead of P. Atomic charges are not observables, and they can vary significantly with the basis set. However, we found that contour plots of electron density showed essentially the same information (in a less transparent way).

Figure 10. Plot of the occupation (n) of the redox species calculated via the Mulliken method. The occupation goes from about 0.5 near the transition state to near 1.0 in the oxidized state over about a 1 ps interval (as we calculated it for the figure, n ) Mulliken redox species charge - reduced species charge). The run was made starting in a prepared transition state like the curves in Figure 8. The solid line is the occupation calculated with solvent absent, and the dashed line is with solvent present (see text). In this example the similarity between the two curves starts to disappear near 1 ps.

ported to the back of our slab by applied fields of the correct sign, as we also found with our EHT simulations.1 Let us look at our ET between localized sites and see what conclusions it brings about the ET mechanism (we later return to the question of delocalized transfer). In this case we find that the electronic occupancy of the redox molecule is sometimes significantly correlated with its configuration. We say this because we sometimes get the same trend in redox species occupancy if the solvent is removed and we let the redox species execute the same motion it did with solvent present. We found this both for simulations in equilibrium and near the transition state, although the fluctuations in the equilibrium state occupancy are of course much less dramatic. We show a graph of this phenomenon in Figure 10 where the occupancy of the redox species is plotted as a function of time. The general trend can be such that the occupancy of the orbital tends toward +2 (+3) at configurations closer to the In(H2O)62+ (In(H2O)6+3) equilibrium configuration.


Smith and Nozik

Thus, analogous to homogenous theory, we find that a redox species configurational change can result in a change in the electronic occupancy of the molecule. But the configurational dependence is only very roughly analogous. A simple analytical model of the configurational dependence of the occupation for our system can be developed by extending the Anderson model above, and it brings some insight, as we now describe. In the ABETT model above (eqs 1-22), ∆E includes the displacement of the redox species ligands. Some of the physics of this displacement related to charge transfer is not apparent from the Hamiltonian (eq 5). Let us look at this physics by taking a simple model of the redox species in which only a single vibrational mode of the molecule is accounted for. For the time being we will neglect the solvent. To explicitly account for the redox species we write a modified form of eq 5 often employed in other contexts37,38 which reads

H ) naσσ(Q) +

∑k {knkσ + (Vakc†aσckσ + Vkac†kσcaσ)}

(29)

where Q is the displacement of some harmonic mode of the redox species, which is assumed to be most important (one could also imagine the redox species is a diatom). In Appendix F we show that it can be found that the exchange of electron density caused by the change in bond length, due to vibration, can be attributed to an electron-phonon coupling, where the phonon is that of the redox species vibration. If we were to go on with the analysis we would find the shifts and broadening in the electronic states. But the key point is madesthe occupation of the orbital and the redox species configuration (via its vibration) are coupled. We find that such a conclusion is also had from a full scale calculation of a gas phase adatom, within our PM3 method. Within the model just described one can also imagine putting in the solvent modes. If explicit redox species modes can be identified that contribute to the ET, then one could imagine mapping the electronic occupation of the orbital as a function of Q of the redox species and modes Qi of the rest of the solvent. Such a mapping is extremely difficult, as we now begin to discuss. While the qualitative expectations of the simple model (discussed in Appendix F) are sometimes met for this system in the simulation, we also sometimes find that when the solvent is removed (in simulations like those in Figure 10) the occupancy as a function of time does not correlate at all well with the redox species. In other words, the solvent is sometimes the dominant player in changing the orbital occupancy of the molecule. For example, there can be such large changes in the field acting on the redox species stemming from the solvent that the redox molecule can be in a perfect configuration for ET but the field can prevent it. Also, most simulations did not have supporting electrolyte included. When ion pairing between the redox species and supporting electrolyte is included (by introducing a single counterion into the simulation box), we find that the noncorrelated behavior just spoken of is especially pronounced and the supporting electrolyte field can significantly influence the ET. In general these counterion effects occur in the fashion one might imagine. A supporting electrolyte ion situated in such a way that its field pulls electron density on the ion often produces reduction in instances of redox and even solvent configuration in which it would otherwise probably not occur. But this effect can be quite complicated. Of course, the supporting electrolyte also influences the curvature, shape, and barrier height of the free energy curves. If we inquire further about the mechanism, still keeping the semiconductor atoms fixed, then the obvious question that arises is, does the redox species’ motion primarily cause transfer of

electron density, are the solvent and counterions more involved, or is it the result of all players? What we have just said suggests the latter, and indeed we find in general this is the case. As we said in section 1 and as our simulation results have shown: the charge on the redox species, its electronic occupation, the orbital energy, and the solvent fluctuations are all mutually coupled. While one might hope to separate out a few reaction coordinates comprised of redox species bond length, or normal modes of vibration, or particular solvent polarizations, etc. to describe the ET, it is extremely difficult to do this in practice. In section 5 we have a final summary regarding the “causality” we have been talking about in this system. We reach the expected answer, namely, that we are dealing with a stochastic system. But first let us add a new player. The dynamics picture is even more complicated when vibrating semiconductor atoms are added: there is now a new player. Again we discuss this topic for the case of close redox species distances where localized ET occurs (see description above). If the distances between the atoms is disordered due to thermal vibrations, then at any instant of time all sites will be inequivalent to some degree. The main effect of the difference in on-site energies is that the possibility of the electron being exchanged with the semiconductor is in general decreased because it is more difficult for resonant transfer to occur. Resonant transfer requires both the semiconductor and solvent to fluctuate so that the energy is the same regardless of whether the unit electron density is on the semiconductor or redox species. In this picture the chance of ET becomes proportional to these two roughly independent fluctuations. Reference 39 addresses directly analogous phenomena in a one-dimensional tight binding solid model where the same general picture appears. While there is a possibility for the ET to be facilitated with favorable changes in hopping matrix element related terms, we find that the effect of the disorder seems typically to decrease ET relative to the fixed lattice case. For example, recall κ decreased by about a factor of 10 for this dynamical semiconductor case (see Table 1). This same picture seems to hold for larger slabs (seven layers). Apparently this is because the dynamical disorder at the surface, in the spatial region where surface phonons occur, is the most important since at small separation distances the electrons are initially transported and localized there. But the bulk disorder deserves more attention in our simulations (see Appendix C). With sufficient disorder all sites can become localized. (In a tight binding model, for example, if the on-site matrix elements alone are sufficiently different on each site, then we have the famous conventional Anderson localization.40,41) Then electron transport can only be of the hopping form, from one localization site to the next as a result of thermal activation, which can be phonon-assisted. In hopping models the electron transport rate constant is of the form41

k ) V exp(-2ar) exp(-∆/kT)

(30)

where a is the localization length, r is the distance between localized sites, ∆ is the hopping activation energy, and V is a normalized frequency depending on the phonon density of states and electron-phonon coupling. This k has a form similar to nonadiabatic ET formulas;25 note the exp(-2ar) term describes the overlap of two wave functions,41 analogous to Γ. It is thus of particular interest to see if hopping behavior might be pertinent to the SLI rate constants. It is well-known that localization is very dependent on dimensionality, and that 1D models are a particularly “pathological” case. Thus we opt to use the full simulation and



calculate the participation ratio of disordered media theory to determine the hopping potentialities of our system. The participation ratio for electron localization is (a similar quantity holds for phonon localization):41

P∝[

∑i |ψi|2]2/[∑i |ψi|4]

(31)

The participation ratio (eq 31), where ψi is the wave function/ orbital of interest on sites i, measures the degree of electron localization (P can be generalized for a sum over many orbitals over real space); an electron state is more localized the closer eq 31 is to zero. We find that indeed localized sites or traps for electrons can appear on the surface due to surface disorder. However, at least in the system we examined, the traps are not sufficiently deep to require true phonon assistance for escape nor are the traps that are also present in the disordered solvent. Only if the hopping length is small will the cluster adequately model hopping. A significantly larger cluster might show different characteristics. At any rate, our findings by no means preclude the possibility of hopping transport in other systems, especially those composed of amorphous semiconductors where we fully expect it. The role of semiconductor phonons in SLI ET is essentially unexplored. Further study and development of the simulation is needed to assess well the role of lattice vibrations in the ET. In particular a good handle on how the surface actually reconstructs in the presence of solvent and under bias is needed. Perhaps a definitive answer to the problem must await this development since the role of surface phonons appears so important. 5. Discussion Summary of Mechanism. We said we would give a final summary of causality in this system. Perhaps the easiest way to explain what we find is as follows. If we imagine the whole system in equilibrium with the reactant redox molecule, for transfer to be effected a fluctuation is necessary to tip the scales and cause the occupation to be closer to the product state. As Figure 8 shows, this fluctuation to the transition state, as expected, happens on a relatively fast time scale; it is not something that slowly builds over time. It is also, obviously, a very rare event since there is a substantial barrier to be surmounted. This escape from a metastable state, is of course the subject of a vast literature.35 Once the transition state is reached dynamics operate again to either trap the system and move it to product or send it back to reactants. In our particular system, the occupation of the orbital must be stabilized by fluctuations of the solvent and semiconductor. The redox molecule itself (in rough accord with our simple model above, eq 29) must also form either a stable trap that deepens with time or, in the reverse process, a trap that dissipates. As Figure 8 shows, the process is different depending on the “initial conditions” of the transition state. While our simulation has identified the major players of the ET system, characterized the interactions, and allowed rate constant expressions to be formulated and the simple ABETT model captures the physical essence of the ET, there is still an air of mystery about the process. That air of mystery is common to all stochastic systems. The underlying principle of such systems is that the random forces stemming from the many body dynamics present the key inputs in getting the system away from states of local stability.35 The random forces are rather weak, but there are infrequent fluctuations of the kind we have just mentioned that induce excursions over the barrier to another state of local stability. The complex many body interplay, the

random forces, is reduced in stochastic models to a heat bath or noise. In our case the biggest source of many body noise is the solvent and perhaps the semiconductor dynamics. Simulations of the baths like ours enables them to be modeled realistically rather than by phenomenological guesses, and one can probe key questions about them. The generic bath becomes a detailed physical entity, reflecting the real world as accurately as possible, and simple descriptions can be checked. The bath can thus, with continued study, be described in more “deterministic” terms. Finally, it is apparent from our simulations that the diffusive motion of the redox species itself may play a fundamental role in the ET. Constraining the redox species to be a fixed distance from the electrode (as we do to simplify and decrease computational effort) may introduce an artificial bias into the calculations. Other Models for the SLI. Since our system is of the stochastic type, one can hope to develop a model for it by using the generalized Langevin equation based formalisms that have been extensively developed of late.2,10,35,42 It is possible in principle to model many of the dynamical effects we talked about with the frequency-dependent frictions of this approach. We can begin, for example, with the Grote-Hynes formalism where the ET rate constant is modeled by

k ) κGH kTST ) κGH (ωR/2π) exp(-∆F‡/kBT)

(32)

where κGH is the dynamical transmission coefficient of GroteHynes theory (compare eq 28). Adapting Grote-Hynes theory2,42 to the semiconductor-liquid system we can write

κGH ) ωB(s + M0-1{ζsolv(s) + ζsemi(s) + ζelec(s)})-1 (33) where ωB is the barrier frequency. One models each effect by a separate frequency (s) dependent friction. The ζsolv(s) is the well-known solvent friction of the generalized Langevin equation approach to ET.2,42 ζsemi(s) accounts for friction stemming from the dynamical fluctuations of the semiconductor lattice as we discussed in ref 1. We also include an electronic friction term (ζelec(s)) which stems from the electronic dynamics of the semiconductor, solvent, and redox molecule. The key input to this model, the frictions, can be obtained from the simulation. Then the results can be compared to the k, calculated directly from the simulation, to gauge the success of the model. The activation energy and barrier frequency can also be put in from the simulation, and many of the effects we mentioned will influence these as well. Polaron Model of ET. It may seem to the reader that a polaron model might be able to capture some of the physics of the processes discussed above. This type of model was actually embraced by the Russian school of ET43 and its basic framework lies, often unacknowledged, at the heart of conventional ET theory. That polaron-type model concerns the liquid side of the interface. What we are suggesting here is a polaron model for the solid; the polaron is an extensively studied solid state concept. We propose that the transferring electron in the solid may be trapped in phonon modes at the surface, for both directions of ET. Further, such a polaron model can feed directly into a hopping model like eq 30. But to describe anything like what we talked about in sections 3 and 4, we must address the surface polaron problem in the semiconductor not the usual bulk polaron problem.44 We find that there does not appear to be a bound polaron state near the semiconductor surface within the conventional dielectric continuum polaron formalism,45 unless very large electron-phonon coupling is assumed. These conclusions can be reached by substituting

2470 J. Phys. Chem. B, Vol. 101, No. 14, 1997 static and high-frequency dielectric functions appropriate to semiconductors in conventional models.45,46 We note that the use of nonlocal dielectric functions47 might alter the conclusions. But if a local dielectric continuum polaron model is used, these trapping states apparently appear so close to the surface that the continuum model itself is questionable, in accord with the conclusions of ref 45 for a different system. Our final conclusion is that models like those we have been working on, that explicitly account for the electronic structure, are the most promising for answering questions about the disordered media/ polaron mechanism of electron transfer. We stress, however, that for some systems simple continuum polaron models might be quite insightful. Model for Delocalized ET. Above we said we would return to the question of delocalized SLI transfer. If a highly delocalized “free electron” is propagating in the bulk semiconductor, it is at least conceivable that at some point it could suddenly collapse to a localized state in the redox species (or vice versa). It is known that free electrons can be nonradiatively trapped by local phonon modes into localized sates, at least in bulk systems. The process is often conceived of as a four-step one:48 (1) the development of a quantum or thermal fluctuation, with the distortion caused by it increasing as time progresses; (2) the emergence of a local level for an electron and the trapping of the electron there; (3) penetration of the barrier to trapping by the composite system comprised of the electron and the lattice deformationshere the deformation would be of the redox species-solvent-semiconductor complex; (4) relaxation of the nonequilibrium state. If a similar process can occur for collapse into a redox species state, then not only may there be highly delocalized transfer (delocalized to local and vice versa), but there may not be time in this sudden collapse to even conceive of an activated complex/transition state. It seems far more likely that, due to the numerous effects of symmetry breaking by the interface, the free electron state would first collapse inside the semiconductor near the surface into a localized state and then subsequently transfer to the redox species. In the latter case a transition state seems feasible. In any event, our studies show that near surface transfer is from localized to the localized state (in the sense we used the word localized above) for both directions of transfer due to microscopic reversibility. Promising avenues of attack for detailed investigation of highly delocalized SLI ET is to use the nonradiative trapping theoretical techniques that have been developed48 and methods that address scattering from impurities. Surface State Modification. The solvent and the redox species itself can create, annihilate, or modify surface states. The role of surface states in semiconductors is enjoying a resurgence of interest, and many new ideas are bearing fruit. We will report in detail on the effects of surface states on the ET in another paper. For now we note that with hydrogen passivation of the semiconductor the Fermi level is not pinned by our water solvent. The redox species can produce surface states, but they are often inextricably intertwined with the redox species orbital occupation; in other words, the states are associated both with the redox species and semiconductor, and these states are similar to the localized kind given in eq 11 above. Final Remarks. It should be recalled that we regard the particular SLI system studied here as theoretically illustrative only (see above). Only a systematic study of many systems will allow the general accuracy of MO-MD simulations to be assessed. We will attempt to focus our future simulations on very diverse systems and compare theory and experiment as rigorously as possible.

Smith and Nozik We believe our results here bring a number of insights into the mechanism of ET and localization at the SLI. Our results suggest a number of new SLI theoretical avenues to be investigated. All this has been detailed above. But our models and simulations also suggest a number of key experiments to help SLI understanding. We discuss a few of these below. 1. As we briefly discussed in section 4, ion-pairing of a counterion with the redox species produces large effects in our simulations. Supporting electrolyte has shown a number of interesting experimental effects at metal-liquid interfaces.21 Systematic study of the effect of counterion pairing, and counterion effects in general, on SLI ET seems an almost open field. 2. For more accurate simulations STM studies of the semiconductor surface of the SLI system are absolutely crucial given the extreme difficulty of accurately predicting reconstruction. One paramount goal of such research should be to study how the semiconductor surface reconstructs in the presence of different solvents and as a function of applied bias. Such studies already exist for the metal-liquid interface,49 and it is wellknown that even metal-liquid interfaces display bias dependent reconstruction. Semiconductors are in general much more prone to reconstruction, and the effects should be even greater. 3. Surface-resolved vibrational studies of the redox species or solvent allow rigorous comparison of theory and experiment, since surface modes are often the result of very local structure and can be well-correlated to calculations. 4. If one uses a disordered semiconductor in the SLI system, the ET should primarily be from localized states of the semiconductor to localized redox species states. The differences in behavior between this system and an ordered semiconductor might be very useful for examining hopping transport effects. Temperature dependent studies would be especially welcome. Acknowledgment. This work was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Office of Energy Research, Division of Chemical Sciences. Appendix A In our 1D semiconductor model1 we assumed the sign alternates for Hi and Hj in the fashion shown in Figure 1. (In this appendix, Hii is abbreviated as Hi and Hij as Hj.) This mimics in one-dimension the s-p bonding found in compound semiconductors. Two simple bands are produced with this model (see Figure 2) with bandwidths (see Figure 2)

EB ) (Hi2 + 4Hj2)1/2 - |Hi|

(A1)

and with band gap (see Figure 2)

Egap ) 2|Hi|

(A2)

The energy equation (the band equation which produced Figure 2) is

(k) ) ((Hi2 + 4Hj2 sin2 k)1/2

(A3)

The 1D semiconductor model we employ is that proposed by ref 14. This 1D model yields the following formulas for the Anderson model quantities14

∆() ) and

Hj′2| + Hi| (2 - Hi2)1/2 2Hj2(4Hj2 - 2 + Hi2)1/2

(A4)



Λ()) Hj′2( + Hi){[((2 - Hi2)/(2 - Hi2 - 4Hj2))1/2 -

The method we employed to calculate Γ for the SLI was to view the semiconductor and redox species as a supermolecule in which charge is transferred from one localized region of the supermolecule to the other. In other words, our basic outlook was that the transfer was “intramolecular”. As stated in the text, if the redox species is close to the semiconductor, the charge transfer is quite localized, so this view should be reasonable. Even when one solvent molecule exists between the semiconductor and redox species the transfer remains localized. To calculate the matrix element at the closest distances (i.e., 5 Å range) we did not use intervening water. The Figure 5 caption describes the site over which the redox specie’s In atom was located, whether In or P termination was used, and other conditions of the different Γ calculations. The cluster size for the semiconductor was 15 atoms (2 layers) for most Γ calculations, adding a third layer did not appreciably affect Γ. The redox species In(H2O)62+/3+ and four other waters were present on the surface. Exploratory calculations were made with more waters. The effect of about 10 additional waters seemed relatively small. The actual method we used is as follows. The electronic coupling matrix element (Γ) is by convention defined as 1/2 of the energy separation at the adiabatic seam between upper and lower adiabatic surfaces which gives25,51

1]/2Hj2} {for 2 > Hi2 + 4Hj2} ) -Hj′2( + Hi)/2Hj2 {for Hi2 < 2 < Hi2 + 4Hj2} ) Hj′2( + Hi) {[((Hi2 - 2)/(Hi2 + 4Hj2 - 2))1/2 1]/2Hj2} {for 2 < Hi2} Equations A4 and A5 when substituted in eq 13 give the projected density of states (F) on the redox species. Appendix B Lattice dynamics is an extensively explored field. In classical simulations it was long ago determined that pairwise interactions are often not sufficient to describe the dynamics.20,50 Threebody rather than two-body interatomic interactions must, at least, be considered for many semiconductors. Further, it is difficult to accurately simulate the surface phonon modes in a classical simulation. We use here a combination of classical and first principles MD to simulate the lattice vibrations. We took an InP 100 hydrogen-passivated semiconductor slab eight layers thick and analyzed the interatomic potentials for all atomic layers including the surface. In all 100 crystal face simulations of the text, except where noted, we used the potentials so determined to model all but the first two layers (which were simulated via “first principles” MD). Thus, in these simulations only the first two surface layers of atoms moved in accord with first principles forces. An important question is how to model the interaction of the “bulk” atoms (deeper than two layers inside the semiconductor) with the solvent and redox species. Conventional electrostatic image formulas can account for static effects. But the nuclear charges of the compound semiconductor atoms are screened by the electronic subsystem, and the screening changes with vibration (as charge exchange occurs between atoms), i.e., the screening is dynamic and is unlikely to be welladdressed by simple dynamic image force formulas. However, in full exploratory first principles calculations we seem to find that the surface layers have the most important interaction with the solvent and redox species. Thus to first order it may be sufficient to include only these first two-layer interactions via first principles. In the general case, it is too computationally expensive to do otherwise with our present methods. Appendix C The calculation of the electronic coupling matrix element is difficult for any system. The matrix element can often only feasibly be calculated at an ab initio level for tens of atoms. But Γ can be sensitive to redox species configuration, and thus Γ should be averaged over a configurational ensemble. Further, ideally, a large number of solvent molecules (and in our case semiconductor atoms) should be included in the calculation. Such calculations rapidly become prohibitive. While there is ample reason for concern, at least for homogenous systems it has been found that even methods as simple as EHT can give reasonable agreement with ab initio methods, and the addition of solvent sometimes shows small effects. Regardless, fundamental questions about the level of rigor and system size needed for accurate estimates linger. Only a few calculations for Γ exist for the metal-liquid interface. We have found no calculations of this quantity for the SLI. To reduce uncertainty, we opted to calculate at the ab initio level.

Γ ) |HAB - SAB(HAA + HBB)/2|/(1 - SAB2)

(C1)

where Hij ) 〈ψi|H|ψj〉 and i and j are equal to A or B. SAB ) 〈ψA|ψB〉. H is the total electronic Hamiltonian, and ψA is the diabatic state for the reactant system and ψB for the product. In our case the reactant state has a unit electron located in the semiconductor; this unit electron is located on the redox molecule in the product state. It has been shown that it is possible to calculate Γ using output from UHF LCAO-MO electronic structure calculations.51 One starts with two nonorthogonal diabatic UHF wave functions. Then the overlap matrix

D ) B†SA

(C2)

is calculated. S is the atomic orbital overlap matrix, and A and B are the eigenvectors of the diabatic states mentioned above. These three matrices are conventional output in LCAOMO programs. A corresponding orbital transformation is done using D ) Udˆ V† where U, dˆ , and V† are related as follows:

∫

dˆ ) bˆ †aˆ dτ ) U†DV

(C3)

That is, the constraint is made that the overlap is diagonal (dˆ is diagonal), where the two sets of spin orbitals of the nonorthogonal diabatic UHF wave function we designate as a and b. The a and b are linearly transformed so that

aˆ ) aV

(C4)

bˆ ) bU

(C5)

Then using the matrices A ˆ ) AV, B ˆ ) BU, and T (Tii ) ˆ TB ˆ †) can be ∏j*i dˆ jj), a generalized density matrix (P ) A calculated. Then this density matrix, U, and V† can be used to compute HAB, and U and V† are used to calculate SAB.51 The calculation of HAA and HBB is straightforward. Using all of these quantities in eq C.1 yields the electronic coupling matrix element Γ. This basic method has been implemented in HONDO,51,52 and we employed it with some modifications. The entire procedure must be carried out for both R and β spin

2472 J. Phys. Chem. B, Vol. 101, No. 14, 1997 occupied molecular orbitals. To create the initial localized states, Boys’ localization method or others can be used. Feeding into the algorithm converged results from calculations in which point charges are used to localize density also works. Ideally, the basis set dependence of the results should be investigated and also compared with Γ calculated via other methods (we primarily used STO-nG, with n ) 3-6). The current computational expense prohibits this for us. The addition of many more layers of semiconductor is a necessity for examining the dependence of the matrix element on cluster size. The latter is also currently computationally prohibitive. More waters should also be added to the calculation. The addition of semiconductor atoms should be more important than the latter. Pseudospectral methods may be well-adapted for addressing much larger systems.53 Appendix D In a classical MD ET simulation, in effect, the charge on the redox species (q) is varied in increments so that it spans the range q ) qred + z, where qred is the reduced redox species charge and z varies between 0 and 1 (the oxidized charge is qox ) qred + 1). This allows a broad range of configurations and momenta of the solvent to be explored so a sufficient amount of phase space can be sampled. Without this technique only a small part of the free energy curve (see Figure 7) could be constructed (thermally accessible regions) without resorting to ridiculously long computational times. So the dynamics is done with a given fixed effective charge on the redox species. Then at each time step the oxidized and reduced species are substituted into the frozen configuration and the difference in their interaction energy with the system is calculated, i.e., ∆E is calculated. Several features of this conventional free energy curve method (mentioned in the text) need modification for application to first principles MD, as next discussed. First, the ∆E calculation cannot, in principle, be done the same way in first principles MD. When a species with oxidized or reduced charge is substituted in (when the total system charge is made equal to the oxidized or reduced charge), the charge redistributes between the semiconductor and the redox species since the electronic wave functions adjust to the nuclear configuration. So the charge on the substituted ion has the potential to become very different than its oxidized (or reduced) value. One way to force the charge on the redox molecule to become qox or qred is to substitute for two given appropriate orbitals their sum and difference. A second way is to use Boys’ localization method, or other location methods, to fix the charge on the redox species. A third method is to use point charges to localize charge, the point charges can then be removed, and the converged point charge produced wave function can be used to start the calculation. The density will remain localized as in the point charge case if the nuclear configuration is appropriate. In practice we find the third technique is often the easiest. With any of these localization methods, after the selfconsistent solution is reached it must be verified that the coerced q condition remains. While we had to use one or the other of these techniques for some configurations, often we did not have to. That is, sometimes simply assigning a total charge on the system of a given value resulted in a qox or qred charge on the redox species. Another difference arises when performing the umbrella sampling. We first applied the sampling constraint (the umbrella potential) only to the waters that were treated classically and thus the umbrella forces could be calculated easily. (The calculation of the adiabatic umbrella forces is in principle the

Smith and Nozik most difficult part of the first principles sampling method.) Thus the outer waters, in effect, moved in a potential that corresponded to having a qred + z charge on the redox species. The actual charge on the ion computed from the electronic structure calculation was different; but the sampling still fully served its purposesto get the system biased away from its “thermal trap” and sample more of phase space. After an initial equilibration under the constraint just discussed, the umbrella sampling charge was altered until it and the redox species charge calculated from the electronic structure calculation were self-consistent (actually, until the ion electric fields were equated). In this manner configurations were generated that sampled over the free energy surface. We generally found that the waters close to the ion that were treated quantum mechanically did not have atomic charges significantly different from those in bulk water. Thus, as one might guess, it typically proved possible to treat all the waters except those directly around the ion as classical for purposes of the initial sampling to get the system out of the thermal trap. Later in time some of the other waters were treated quantum mechanically again, and the same self-consistent process as above was undertaken. Appendix E It is nontrivial to identify the transition state for our SLI system. But a working definition is that at the transition state the total system energy is the same with either the reactant or product redox species placed in a given solvent and solid configuration. (This definition does correspond to the maximum of the barrier peaks in our free energy curve of Figure 7. But it is not guaranteed that the true saddle point, corresponding to the transition state, can be accurately found by such a definition. Peculiarities of how ∆E is defined can introduce error.) It proved possible to first run classical MD simulations with potentials made to fit the PM3 results and then use the umbrella sampling potential to force the system near the transition state (near ∆E ) 0) (see also Appendix D). Using these configurations and momenta as a starting point, we then instituted the first principles MD and again used the sampling potential to bias this system to the new transition state (which typically proved to be close to the old). Attempts to find the transition state by the usual quantum chemical methods54 appear to give similar results, although our sample size is small since the computational cost is extreme. The fact that the trajectories, started from the transition state configurations we found, decayed both to reactant and product (see Figure 8) is significant evidence that we were sufficiently close to the transition state for the system. After the transition state has been located and constrained equilibrium runs are done in this state, the calculation of κ is possible (see text). The calculation of κ has been extensively explored for homogeneous systems.10,11,27,35 Most microscopic calculations of κ have been done in classical MD simulations. We used the same basic method, except that the constrained dynamics was done in the fashion described here and in Appendix D. In ref 55 it was found that the bias in the sampling statistics that is introduced because the trajectories from the transition state have d∆E/dt * 0 (whereas the initial distribution has d∆E/dt ) 0) can be corrected for by using the formula

κ ) 〈〈(∆E˙ )θ [∆E(t)]〉e′〉eTS/〈〈∆E˙ θ [∆E˙ (t)]〉e′〉eTS (E1) where e′ means the averaging is to be done over the velocity distribution and eTS means a configurational average is to be done at the transition state. We utilized the same correction. Issues in κ calculations surround whether it is dynamics in the



reactant well or near the barrier top that are most controlling. These issues have been extensively discussed in the literature.10,27 There have been κ studies for metal-liquid systems, but these simulations have been classical MD.32,33

above. But the qualitative physics is anticipated to be roughly the same. We will present a more detailed analysis elsewhere.

Appendix F

Here we talk about issues related to the energetics of the simulation and rate constant calculations; we also make comparisons to experiment. In ref 1 we showed many figures of the electronic density of states of 100 InP, water, redox species, and the combined system. We do not have space for that here; we simply note the results were of similar quality to that described in ref 1. The calculated band gap of the 100 InP semiconductor, when our periodic boundary conditions and Ewald methods are used, is about 1.8 eV using the PM3 Hamiltonian. This is quite acceptable (1.4 ( 0.1 eV encompasses the range of experimental bulk InP direct band gap measurements); even some very sophisticated methods match experiment more poorly. Without full account of many body effects, band gap disagreement is expected. It has been proved that molecular dynamics and electronic structure methods, combined and separately, can do a reasonable job of predicting energies of solvation (practitioners, however, must be thoroughly familiar with pitfalls). Many papers and books have addressed this question, and some calculated values, even for dielectric continuum methods, can result in errors as small as 10% relative to experiment.56-60 Some of these methods, however, are highly parametrized. MOPAC 93, the premier MOPAC program, has standardized routines for calculating such quantities. But here we calculated the free energy of solvation from our PM3-based molecular dynamics by calculating the total energy of the system with and without the ions. For the In3+ ion the calculated solvation energy was about -2900 kJ/mol (experiment ) -3989 kJ/ mol57). Experimental results for In2+ are apparently, not unexpectedly, unavailable. Redox potentials can be determined if one computes the difference in the total calculated energy for the reduced and oxidized gas phase species and adds to this the solvation energy difference between the two species. The zero of energy for this definition is the energy of a free electron in vacuum. The next problem comes in trying to compare this calculated result to experiment, since electrodes are used whose absolute potentials are not rigorously known. The calculated redox potential was about 0.0 V (NHE), whereas a typical experimental value61 for In2+/In3+ is -0.49 V (NHE). (To convert the calculated value to NHE we assumed the absolute potential of the NHE was 4.5 eV, the value typically estimated.62) While this difference may seem large, even ab initio results have apparent errors in redox potential magnitudes of a volt or more when correlation and exchange is not wellaccounted for. MNDO methods are corrected semiempirically for exchange. A systematic study using our first principles PM3 molecular dynamics will have to be completed to reliably assess our method versus experiment. PM3 does not accurately predict ionization energies for higher oxidation states, thus cancellation of error may be occurring in our result. To more accurately predict the redox potential, high-level ab initio methods should be used, or better semiempirical corrections should be instituted (especially for correlation and exchange contributions to the gas phase energies). It is of interest to know, especially in conventional SLI ET outlooks, how the energy of the redox species lines up with the semiconductor band edges. It is easiest to state what we find as follows. The CB and VB edge of the semiconductor we take to be the LUMO and HOMO, respectively, that we calculate

When the conformation of a molecule changes, electron density redistributes. Within the context of tight binding theory, this means that the on-site (and hopping) matrix elements Hii change. That is, when the charge redistributes, the screening of the nuclei changes and the orbital ionization energy changes, and in turn Hii ≈ the orbital ionization energy. When the Hii of the redox species changes, the orbital occupation also does via eq 10 because Hii is directly related to a (and a directly affects 〈n(∆E)〉 via eqs 7 and 13). In a real redox species the change in Hii (and Hij) is quite complex.1 But from our argument just given, the conformational changes and the vibration of a molecule are expected to give rise to electron density exchange between the semiconductor and redox species. Our electronic structure calculations show the change in 〈n(∆E)〉 with conformation can be dramatic, as discussed in the text. Now the physics of this 〈n(∆E)〉 change can be stated differently and with new physical content, as follows. If we start with the Hamiltonian given in the text (eq 29) and assume the vibration (or conformational change) is harmonic, then Q in eq 29 of the text can be written:

Q ) [h/(πm*Ω)]1/2(b + b†)

(F1)

m* is the effective mass of the electron, Ω is the vibrational frequency, and b and b† are boson creation and annihilation operators. If one further assumes that the change in Q is very small (note this is not necessarily the case in our conformational change problem) to first order, one has37,38

σ(Q) ) σ(q) + σ′(q)Q

(F2)

where q is the equilibrium bond length in the presence of the solid. The other variables in eq F2 are described in the text. It is easy to show by substituting eqs F1 and F2 into eq 29 of the text that an electron-phonon coupling constant results with the value:

C ) (2 h/m* Ω)1/2(σ + U′naσ + Un′aσ) ) δσ

(F3)

If one assumes in eq 13 of the text that ∆ and Λ vary slowly with , then one can write (cf. refs 37 and 38):

δna )

∫F() ) ∫[∂F()/∂σ]δσ ≈ Cf(F)

(F4)

Here δna is the change in the number of electrons in the redox species orbital during one period of vibration. (For a metal, the last equality of eq F4 is [C F(F)];38 for our 1D semiconductor model, the dependence is much more complex, but C, with approximations, is still the prefactor as indicated in eq F4.) Thus the occupation of the redox orbital changes by δna with a vibration, as the bond length changes from equilibrium to its vibrational maximum. The new physical content is that this orbital occupation change can be viewed as arising from electron-phonon coupling, where the phonon is the local one of the redox species. In our case there are substantial conformational changes as well as vibration. It is quite questionable that such effects can be treated well in the first order scheme

Appendix G


Smith and Nozik

from the cluster with our supplemental Ewald sums and the imposed periodic boundary conditions. The energy level of the redox species, in the sense it is used in the ABETT model, can be though of as the average value of the eigenvalues to which the electron is transferred (we actually must deal with all redox species orbitals as we discuss elsewhere1). The absolute energy of the redox species eigenvalues is affected by the solvent (and semiconductor) fluctuations and these are calculated explicitly by the electronic structure method. For our system we find that the electron is transferred to and from a roughly atomic-like 5s-orbital of the In, although there is even some apparent exchange with the water ligands. The average energy of this orbital for the In(H2O)63+ lies about 0.4 eV above the conduction band (CB) edge (LUMO) for the hydrogen-passivated InP surface and about 0.8 eV below the CB edge for In(H2O)2+ at the 5.3 Å separation distance. (Note that these orbital energies change with redox species separation distance from the electrode.) The calculated LUMO and HOMO lie at about 1.3 V (NHE) and -0.5 V (NHE). All this is in very rough accord with energetics measured from experiment since the band edges for InP are placed at about +0.7 and -0.7 V versus NHE,30 and the calculated and measured E0’s are as reported above. But this could very well be fortuitous, and anyway the conduction band edges of InP on this scale is at best only known to within 0.2-0.4 eV. Also SLI experiments are sometimes not particularly well-controlled, and it is quite questionable that values which stem mainly from polycrystalline results (not a single crystal 100 face) are even applicable to our simulation. Further, there is almost certainly reconstruction of the real surface, an oxide layer, contaminants on the surface, pH dependencies, etc.; all of which can have a significant impact on the energetics. Let us review how we calculated the rate constant of section 4. We used the conventional model of SLI ET which is equivalent to the diode equation of semiconductor physics30

i ) i0 exp(-qη/kT) - i0

(G1)

where i0 is the exchange current density, η is the semiconductor overvoltage (assumed to fall only in the semiconductor space charge layer), and q is the fundamental unit of charge. By conventional thought this expression can be replaced with

i ) (qnsurface cox kox) - (qnunocccredkred)

(G2)

where q is the fundamental unit of charge, nsurface is the concentration of electrons at the surface, nunocc is the unoccupied surface concentration of states, and cox,red is the concentration of oxidized and reduced redox species, respectively. kred and kox are the second order rate constants which we calculate from eq 27 of the text by using the first order rate constant calculated with eq 24 of the paper and the two concentrations. Equation 27 is basically equivalent to an equation first proposed by Levich.7 We determined nunocc using our calculated density of states, as discussed in the text. We then calculated kred using eq 27 of the text, where we assumed cox ) cred and employed the other simplifications discussed in the text. At equilibrium i ) 0, and by eq G1, G2, and eq 27 of the text we can then can find exp(-η/kT) since we can determine the first order rate constant for the reduction from eq 24 of the text using our simulation (see text). We get a η of about 0.3 V by this equilibrium calculation. About the only thing that can be said about this result is that it is not unreasonable; it is not 10.0 V for example. It can also be said the field in this case is too low to significantly affect the electronic coupling matrix element (see text). It also

does not affect the fact that at close redox species-surface separation the transferred electron density is initially near one atom (see text for effects of the field on how the location of the transferred electron density evolves with time). Only a systematic study of predicted η and k, for systems in which the equilibrium η is experimentally established, could determine the success of the method. The only real merits of this method are that it follows conventional thought, and serves to show how our calculated first order rate constant can fit into such a framework to give a second order rate constant. In subsequent papers we will detail shortcomings of conventional theory. In this paper we are most interested in trying to model fluctuations in energy realistically, and in knowing relative energies within the model. We are not particularly concerned with trying to reproduce experiment. As we say in the text, ultimately we regard the system of the paper as illustrative only and chose it mainly because of its theoretical tractability. The mechanism of the transfer is of most interest here. In our opinion mechanistic studies must come before accurate rate constant formulations. It seems to us there is no rigorous mechanistic model in the conventional SLI literature. References and Notes (1) Smith, B. B.; Nozik, A. J. Chem. Phys. 1996, 205, 47. (2) Smith, B. B.; Hynes, J. T. J. Chem. Phys. 1993, 99, 6517. (3) Schmickler, W. J. Electroanal. Chem. 1986, 204, 31. (4) Sebastian, K. L. J. Chem. Phys. 1989, 90, 5056. (5) Muscat, J. P.; Newns, D. M. Prog. in Surf. Sci. 1978, 9, 1. (6) Newns, D. M. Phys. ReV. 1969, 178, 1123. (7) Levich, V. G. In AdVances in Electrochemistry and Electrochemical Engineering; Delahay, P., Ed.; Interscience: New York, 1966; Vol. 4, pp 249 (8) Marcus, R. A. In Annu. ReV. Phys. Chem. 1964, 15, 155. (9) Warshel, A. J. Chem. Phys. 1982, 86, 2218. (10) Zichi, D. A.; Ciccotti, G.; Hynes, J. T.; Ferrario, M. J. Phys. Chem. 1989, 93, 6261. (11) Carter, E. A.; Hynes, J. T. J. Phys. Chem. 1989, 93, 2184. (12) McGlynn, S. P. Applied Quantum Chemistry; Holt, Rinehart, and Winston: New York, 1972. (13) Davidson, S. G.; Levine, J. D. 1970, 25, 1. (14) Davidson, S. G.; Huang, Y. S. Solid State Commun. 1974, 15, 863. (15) Stewart, J. J. P. J. Comput. Chem. 1991, 12, 320. (16) Laasonen, K.; Sprik, M.; Parrinello, M.; Car, R. J. Chem. Phys. 1993, 99, 9080. (17) Levine, I. N. Quantum Chemistry; Simon and Schuster: Princeton, NJ, 1991. (18) Menon, M.; Allen, R. Prog. Surf. Sci. 1987, 25, 317. (19) Curtiss, L. A.; Halley, J. W.; Hautman, J.; Hung, N. C.; Nagy, Z.; Rhee, Y. J.; Yonco, R. M. J. Electrochem. Soc. 1991, 138, 2032. (20) Desjonqueres, M. C.; Spanjaard, D. Concepts in Surface Physics; Springer-Verlag: Heidelberg, Germany, 1993. (21) Bard, A. J. Encyclopedia of Electrochemistry of the Elements; Dekker: New York, 1973; Vol. VII. (22) Cotton, F. A.; Wilkinson, G. AdVanced Inorganic Chemistry; John Wiley and Sons: New York, 1988. (23) Rosenwaks, Y.; Thacker, B. R.; Ellingson, R. J.; Burr, K. C.; Tang, C. L.; Nozik, A. J. J. Phys. Chem. 1994, 98, 2739. (24) Cohen-Tannoudji, C.; Diu, B.; Laloe, F. Quantum Mechanics; John Wiley and Sons: Paris, 1977. (25) Newton, M. D. Chem. ReV. 1991, 91, 767. (26) Kuharski, R.; Bader, J.; Chandler, D.; Sprik, M.; Klein, M.; Impey, R. J. Chem. Phys. 1988, 89, 3248. (27) Smith, B. B.; Staib, A.; Hynes, J. T. Chem. Phys. 1993, 176, 521. (28) Hill, T. Statistical Mechanics; Dover: New York, 1956. (29) Zichi, D.; Ciccotti, G.; Hynes, J. T.; Ferrario, M. J. Phys. Chem. 1989, 93, 6261. (30) Morrison, S. R. Electrochemistry of Semiconductors and Oxidized Metal Electrodes; Plenum Press: New York, 1980. (31) Memming, R. Electroanal. Chem. 1988, 11, 1. (32) Rose, D. A.; Benjamin, I. J. Chem. Phys. 1994, 100, 3545. (33) Straus, J. B.; Calhoun, A.; Voth, G. A. J. Chem. Phys. 1995, 102, 529. (34) Kramers, H. A. Physica 1940, 7, 284. (35) Hanggi, P.; Talkner, P.; Borkovec, M. ReV. Mod. Phys. 1990, 62, 251. (36) Morgan, J. D.; Wolynes, P. G. J. Phys. Chem. 1987, 91, 874. (37) Persson, B. N. J.; Persson, M. Solid State Commun. 1980, 36, 175.

Electron Transfer at the Semiconductor-Liquid Interface (38) Ueba, H. J. Chem. Phys. 1982, 77, 3759. (39) Vekhter, B. G.; Ratner, M. A. J. Chem. Phys. 1994, 101, 9710. (40) Adler, D.; Schwartz, B.; Steele, M. Physical Properties of Amorphous Materials; Plenum: New York, 1985. (41) Elliott, S. R. Physics of Amorphous Materials; Longman: London, 1984. (42) Grote, R. F.; Hynes, J. T. J. Chem. Phys. 1980, 73, 2715. (43) Dogonadze, R. R. Reactions of Molecules at Electrodes; Wiley: New York, 1971. (44) Kuper, C. G., Whitfield, G. D., Eds. Polarons and Excitons; Plenum Press: New York, 1963. (45) Elmahboubi, A.; Lepine, Y. Surf. Sci. 1993, 280, 231. (46) Elmahboubi, A.; Lepine, Y. Surf. Sci. 1994, 303, 409. (47) Smith, B. B.; Koval, C. A. J. Electroanal. Chem. 1990, 277, 43. (48) Quantum Tunnelling in Condensed Media; Elsevier: B. V., 1992. (49) Furtak, T. E. Surf. Sci. 1994, 299/300, 945. (50) Abraham, F.; Batra, I. Surf. Sci. 1985, 163, L752. (51) Farazdel, A.; Dupuis, M.; Clementi, E.; Aviram, A. J. Am. Chem. Soc. 1990, 112, 4206.

J. Phys. Chem. B, Vol. 101, No. 14, 1997 2475 (52) Dupuis, M.; Marquez, A. HONDO 95.3 from CHEM-Station; IBM Corporation: Kingston, N.Y., 1995. (53) Zhang, L. Y.; Friesner, R. A. preprint, 1996. (54) Baker, J. J. Comput. Chem. 1986, 7, 385. (55) Carter, E. A.; Ciccotti, G.; Hynes, J. T.; Kapral, R. Chem. Phys. Lett. 1989, 156, 472. (56) Nagono, Y.; Sakiyama, M.; Fujiwara, T.; Kondo, Y. J. Phys. Chem. 1988, 92, 5823. (57) Marcus, Y. Ion SolVation; John Wiley and Sons Ltd: New York, 1985. (58) Cramer, C. J.; Truhlar, D. G. Science 1992, 256, 213. (59) Zauhar, R. J.; Morgan, R. S. J. Comput. Chem. 1990, 11, 603. (60) Miertus, S.; Scrocco, E.; Tomasi, J. J. Chem. Phys. 1981, 55, 117. (61) CRC Handbook of Chemistry and Physics, 74th ed.; CRC Press: Cleveland, Ohio, 1993-94. (62) Pearson, R. J. J. Am. Chem. Soc. 1986, 108, 6109. (63) Lohmann, F. Z. Naturforsch 1967, 22a, 843. (64) Lowrey, A. H.; Cramer, C. J.; Urban, J. J.; Famini, G. R. Comput. Chem. 1995, 19, 209.

Theoretical Studies of Electron Transfer and Electron Localization at

Recommend Documents