3268
J . Phys. Chem. 1988, 92, 3268-3273
Bipolarons in Metal-Metal Halide Solutions E. S. Fois, A. Selloni, M. Parrinello,* and R. Car International School for Advanced Studies, Strada Costiera 1 1 , 3401 4 Trieste, Italy (Received: August 27, 1987; In Final Form: December 22, 1987)
A novel molecular dynamics method is used to follow the adiabatic dynamics of two electrons solvated in molten KC1. The electrons are treated quantum mechanically within the local spin density approximation. A coupled set of Newtonian and time-dependent Schrodinger-like equations is used to describe the evolution of the ions and of the Kohn-Sham orbitals. We find that parallel spin electrons repel each other and form separate F-center-like states. Antiparallel spin electrons, instead, attract each other and coalesce into a single bipolaronic complex. The electrons sit mostly in an ionic cavity which is surrounded by cations. The diffusion of the bipolaron, while bound, occurs on an ionic time scale. However, dissociation processes occur during which the electrons can acquire a high mobility leading on average to a large electronic diffusion.
I. Introduction The physics of extremely dilute metal-metal halides solutions has recently been thoroughly in~estigated,'-~ and it has been found without ambiguity that in the extremely dilute limit the valence electrons released by the excess metal atoms localize in states that are the liquid-state counterpart of the F centers in solids. In such states the electron surrounds itself mostly by cations, replacing an anion in the charge-ordered structure of the fluid. At higher concentrations different F centers can interact, and experimentsbs indicate that spin pairing is operative and that various forms of aggregates can occur. This has been confirmed by numerical path integral (PI) simulations9 that have demonstrated a tendency for pairs of electrons to form spin-paired bielectronic complexes. By analogy with solid-state concepts we call these states bipolarons.I0 In fact two antiparallel spins tend to localize in the same spatial position, since the gain in electron-ion interactions is larger than the loss in Coulombic energy that is brought about by the electronic repulsion. PI calculations are however limited in that they have great difficulty in treating the triplet state and cannot give direct information on the dynamics of the electrons, which is an important facet of the physics of these systems. In particular, questions such as the nature of diffusion processes, or the excitation spectra, and the recombination dynamics can barely be answered by PI calculations. Although analytical continuation of the PI calculations11J2can provide a possible way around the latter problem, it is not clear how the result is sensitive to numerical uncertainties in the input data. Besides, it eventually yields as output only thermally averaged correlation functions. From this a detailed model for the dynamics of the electrons has to be deduced in a rather indirect way. Here instead we propose and apply to the bipolaron case an alternative approach that, although approximate, is capable of providing a detailed and physically very ~
~
(1) Warren, Jr., W. W. In The Metallic and Non-Metallic States of
Matter: An Important Facet of the Chemistry and Physics of the Condensed Stare; Edwards, P. P., Rao, C. N., Eds.; Taylor and Francis: London, 1985. (2) Warren, Jr., W. W.; Campbell, B. F.; Brennert, G F. Phys. Reu. Lett 1987, 58, 941. (3) Pfeiffer, E.; Garbade, K.; Freyland, W. J. Electrochem. SOC.,in press;
Proceedings of the Vth International Symposium on Molten Salts, Honolulu,
HI,1987, in press. (4) Parrinello, M.; Rahman, A. J . Chem. Phys. 1984, 80, 860. ( 5 ) Selloni, A.; Camevali, P.; Car, R.; Parnnello, M. Phys. Rev. Lett. 1987, 59, 823. (6) Nicoloso, N.; Freyland, W. Z . Phys. Chem. (Munich) 1983, 135, 39. (7) Nicoloso, N.; Freyland, W. J. Phys. Chem. 1983, 87, 1997 (8) Freyland, W.; Garbade, K.; Heyer, H.; Pfeiffer, E. J . Phys. Chem. 1984,88, 3745. (9) Parrinello, M.; Rahman, A., to be published. (10) Alternatively, these states can be considered the liquid-state analogue
of the F' center in ionic solids (see, e.&: Stoneham, A. M. Theory ofDefects in Solids; Clarendon: Oxford, 1975). (11) Thirumalai, D.; Berne, B. J. J . Chem. Phys. 1983, 79, 5029. (12) Nichols 111, A. L.; Chandler, D., to be published.
0022-3654/88/2092-3268$01.50/0
vivid picture of the bipolaron behavior. In our approach we confine ourselves to the Born-Oppenheimer (BO) approximation. Namely, we shall assume that electrons always remain in the ground state pertaining to the instantaneous ionic configuration. A further approximation will be that of treating the electronic quantum many-body problem in the local spin density (LSD) appro xi ma ti or^.'^*'^ In so doing, we will loose information on the full many-body wave functions and focus only on the number and magnetic momentum electron densities. The validity of the combined BO and LSD approximations will be proved in the following for the bipolaron in molten salts. We should just add here that the present approach is of general use and can be applied to other physical systems, and within the above-mentioned approximations it does not have a negative sign problem as an exact simulation of a many-fermion system would. For instance, here we will treat the two-particle singlet and triplet states with the same ease, in spite of the fact that these two states have different symmetry relative to particle coordinates exchange. The singlet state is in fact even, while the triplet is odd. This leads in the latter case to nonpositive defined probabilities and therefore to severe difficulties in the numerical evaluation of the relevant expectation values, within the PI approach. Our way of solving the LSD equations will be rather different from the conventional approach in that in the course of the simulation we will not explicitly solve the self-consistent Schriidinger-like equations, which are the Euler equations of the minimum problem associated with LSD theory, but rather we will introduce a dynamics for the combined time evolution of electrons and ions. This dynamics will lead to Newtonian equations for the ions which are treated classically and to time-dependent self-consistent Schriidinger equations for the LSD single-particle orbitals. Under appropriate conditions this dynamics coincides with the BO dynamics for the combined quantum classical system. We will show that these conditions are to a high degree of accuracy satisfied here, and therefore our approach provides a viable and computationally convenient way of finding the BO-LSD trajectory. 11. Theoretical Framework
In this section we derive and discuss the equations of motion that are used in our study. Fundamental to our approach is the Born-Oppenheimer or adiabatic approximation which allows decoupling the electronic and ionic degrees of freedom under the assumption that the typical ionic kinetics energies, of the order of kT, are much smaller than the typical electronic excitation energies A. In this situation it becomes possible to use a classical (13) See, for instance: Theory of the Inhomogeneous Electron Gas; Lundqvist, S.. March, N. H., Eds.; Plenum: New York, 1983. and references therein. (14) Gunnarson, 0.;Lundqvist, B. I. Phys. Reo. B: Solid State 1976, 13, 4'71.
0 1988 American Chemical Society
The Journal of Physical Chemistry, Vol. 92, No. 11, 1988 3269
Bipolarons in Metal-Metal Halide Solutions approximation for the ions and thus describe them in terms of their positions (RJ. The equations describing such an adiabatic motion are
H((RJ) 'Po((RJ) = EO((R1) 'Po((RJ)
(2)
where MI is the mass of the Ith ion, VI, the interionic potential, and Eo((RJ)the electronic ground-state energy. This is obtained by solving the Schrodinger equation ( 2 ) with the electronic Hamiltonian H((R1) = T +
VeI
+ Ve,
(3)
where T = -( 1 /2)xiV; is the kinetic energy operator, VeI(r,(R)) = ~ p e 1 ( r i , ( Rthe J ) electron-ion potential, and V, the electronelectron interaction. The set of equations (1) and (2) simply states that the electronic system is always in the ground state corresponding to the instantaneous ionic configuration (RJand, correspondingly, the ions move on the BO potential surface Eo((R)). Using the Hellman-Feynman theorem, one can express eq 1 in a form which is more convenient for numerical computation. This is
MI-d2Rl = --aVII((RJ) - Jdr a u e I ( r m no(r) dt2 aRI a RI
(4)
where no(r) is the electronic ground-state density. Sometimes it can also be convenient to replace eq 2 with the time-dependent SchrMinger equation -i
a*(O = H(t) *(t), at
*(to)
= %((R(to)t)
(5)
where to is a given initial time and H(r) = H((R(t))). The set of equations (4)and (5) has been extensively discussed in the literature.Is For adiabatic motion eq 5 is equivalent to eq 2, since the time variation of H ( t ) is then so small with respect to A-l that no electronic transition is induced, and thus * ( t ) coincides with @,((R(t))) except for a phase factor. One can remark however that if for some ionic configuration the electronic excitation energy A becomes of the order of kT, then *(t) can deviate substantially from a,((R(t))), thus indicating the occurrence of an electronic transition. This does not mean however that the set of equations (4) and (5) is capable of describing nonadiabatic motion, since in this case the assumption of a single, well-defined classical trajectory for the ions (eq 4) breaks down. Thus, the set of equations ( 2 ) , (4) or (4), (5) has in general the same range of validity, both being restricted to a description of adiabatic motion. It is mostly a matter of computational convenience to use either eq 2 or eq 5 to construct dynamically the BO potential surface for the ionic classical motion. So far no approximation has been introduced in our treatment of the many-electron system; Le., H in eq 3 is the full many-body Hamiltonian and, correspondingly, a0or \k is the exact groundstate wave function for the system. However, such an exact treatment is impossible in practice, since handling and storing a many-body wave function poses enormom computational problems. An accurate and practically very convenient solution of the quantum many-body problem can however be obtained by using the local spin density (LSD) approximation to spin density functional (SDF) theory.14 Using this approach, one can introduce a set of single-particle wave functions $Ji,"(r) (here a = t, 1 is the spin-label) and define in terms of them a density n(r) = Cld12+ 1
Cld12 i
nt(r)
+ nr(r)
(6)
and local spin polarization (15) See, for instance: Tully, J. C. In Dynamics of Molecular Collisions, Part B Miller, W. H., Ed.; Plenum: New York, 1976, and references therein. More recent references can be found, for instance, in: Kirson, 2.;et al. SurJ Sei. 1984, 137, 527.
The calculation of the ground-state energy Eo((RJ)of the electronic system corresponding to the ionic configuration (R) amounts to a minimization of the functional 1 E = - -2Ei.aJ d r 4*i,aV24i,a+ Jdr n(r) ue1(r,(RJ)+ l S d r dr' n(r) n(r') + Jdr 2 Ir - r'l
cxc(n(r),{(r)) (8)
with respect to the (+i,a(r)).Here txc(n(r),{(r)) is the exchangecorrelation energy per particle. For the set (diJr)J which minimizes the functional (8), n(r) and {(r) coincide with the true ground-state density b ( r ) and spin polarization lo(r) of the system, respectively. The problem of minimizing E with respect to the (+i,n(r))is equivalent to solving the Kohn-Sham (KS) equations13
where uH and uXc," are the Hartree and spin-dependent exchange-correlation potentials, respectively, and the 4i,aare the KS orbitals. Within SDF-LSD theory, eq 2 and 4 describing the adiabatic dynamics of the system thus become equivalent to the set of eq 4 and 9. The problem, although enormously simplified, is still rather formidable. Indeed, for each instantaneous ionic configuration (RJone must solve eq 9 (one for each spin direction) self-consistently, a procedure which is quite time-consuming. By analogy with eq 5, we now replace eq 9 with the corresponding time-dependent equations
subject to the initial condition $i,a(to)
= 4'i,a((R(tO)l)
(lob)
where (+oi,a((R(to)))] is the set of KS orbitals which minimize (8) for the ionic configuration at to. We remark that h, in (loa) depends on t not only through the ionic coordinates (R(t)Jbut also through the time-varying electronic charge density n(t) and spin polarization {(t). However, if the coordinates (R(t))vary slowly on the time scale A-l, then so do n(t) and {(t),which thus coincide with the ground-state density no((R(t)))and spin polarization Ib((R(t)]) for the instantaneous ionic configuration. Hence, eq 10 is equivalent to eq 9 for adiabatic motion. It should be noted that even under adiabatic conditions it is not necessary for the ($i,a(t)) derived from eq 10a to coincide with the instantaneous KS eigenstates (q5°i,a((R(t)))).In fact, any unitary transformation within the subspace of the doi," preserves the electronic charge and spin density and leads to identical forces on the ions and identical electronic energies. If such a rotation does occur, the KS orbitals can be simply obtained by diagonalizing the small matrices hyJ = ($i,alhal$,,a)where i and j run over the occupied can then states N". The corresponding unitary transformation This be used to obtain the KS orbitals doi,,as doi,, = procedure is obviously much less time-consuming than that of diagonalizing the ha matrix, whose size M is as large as the size of the basis set used to represent the +Oi,,, and in general M >> N".
eJ xg,vJ$j,a.
111. Model Calculation
In this section we discuss some technical aspects of our calculations for two solvated electrons in KCl. We solve the coupled eq 4 and 10 relevant to the BO condition of our system as follows. For the classical equations (4) we use the standard Verlet algorithm.I6 For the time-dependent KS equations (10) we use instead the Trotter formula (16) Verlet, L. Phys. Rev. 1967, 159, 98
3270 The Journal of Physical Chemistry, Vol. 92, No. 11, 1988
Fois et al.
$i,a(t+At) = exp(-ihaAt)$i,a(t) LZ
exp(-iV,At/2)
exp(-iTAt) e~p(-iV,At/2)$~,,(t)
+ ...
(1 1)
+ +
where V, = ueI uH uXc,,. We perform the matrix multiplications in (1 1) using fast Fourier transforms (FFT) to switch from real space, where exp(-iV,Aht/2) is diagonal, to reciprocal space, where exp(-iTAf) is diagona1.l’ One advantage of eq 11 over other schemes is that the time evolution of the $ is automatically unitary. Therefore, if the wave functions are initially orthonormal, they will maintain this property during all subsequent time evolutions. This remark is particularly relevant in the case of parallel spin electrons where orthonormality is a necessary consequence of the Pauli principle. This makes the simulation of triplet state just as easy as that of the singlet state. The time step Ar to be used in the integration procedure (1 1 ) should be on the time scale of the electronic motion, Le., rather short. On the other hand, the difference in time scale between the electronic and ionic motions allows us to use a somewhat larger time step A T to integrate eq 4. We typically use Af = 1.2 X s and AT = 20Af. With these values the total energy is conserved during a run of 2 X IO5At. within 1 X In our system periodic boundary conditions are imposed, each unit box (a cube of length 26 au) containing 32 K+ and 30 C1ions. We represent each electronic wave function on a discrete mesh of 163 points in the cubic cell. For the ion-ion interaction potential we take the same Born-Meyer potentials used in ref 4 and treat the long-range part by the Ewald summation method using 163reciprocal lattice vectors, Le., all the G vectors compatible with the choice of the mesh in direct space. The overall convergence was verified for a few typical ionic configurations by using 323 grid points and reciprocal lattice vectors. Our electron-ion potentials are a smoothed version of those used in ref 4, Le. 1 u* = f- erf ( r / R * )
r
where the core radius is R+ = 3.7 au for K+ and R- = 2.2 au for C1-. Such smoothing is required by our use of a discrete mesh to represent the $,,, and does not introduce any significant quantitative efect. This was explicitly checked in the case of a single solvated e l e c t r ~ nwhere ,~ we studied the effect of varying R+and/or the mesh size. For the exchange-correlation part of the electron-electron interaction u,,,, we use an interpolation*8 of the results of ref 19. An important feature of the calculations is the preparation of the initial state. We proceed as follows. First we prepare a well-equilibrated liquid of 32 K+ and 30 C1- in a neutralizing uniform background a t a temperature of about 1300 K and a density of 1.5 g/cm3. At a given time to we remove the uniform background and calculate the ground state of the two-electron system corresponding to the ionic configuration (R(to)l and the chosen value of the total spin. ( S = 0 and 1 for the singlet and triplet states, respectively.) We do this by minimizing the energy functional (8) with a steepest descent method.*O The set of functions +Oi,, ((R(to))) thus obtained and the stored ionic configuration ((R(to)))are used as initial conditions for the run, which proceeds using eq 4 and loa. This procedure was followed in order to avoid as much as possible bias from a particular choice of initial conditions. Given the fact that eq 4 and 10a lead to a correct dynamics only if the system is adiabatic, it is necessary to check the motion adiabaticity. This is done by projecting the current electronic state (qi,(t))onto the ground state ($Oi,,((R(t)))] obtained by minimizing the functional (8) for the fixed ionic configuration (R(t)J. For a strictly adiabatic motion this projection is unity. In our simulations this condition is verified except for rather infrequent events during which for a relatively short time deviations (17) Feit, M. D.; Fleck, J. A.; Steiger, A. J. Comput. Phys. 1982,47,412. (18) Perdew, J. P., Zunger, A. Phys. Reo. B: Condens. Matter 1981, 23, 5048. (19) Ceperley, D. M.; Alder, B. J. Phys. Reu. Lett. 1980, 45, 566. (20) Car, R.; Parrinello, M.; Andreoni, W. In Microcluster; Sugano, S., Nishina, Y.,Ohnishi, S., Eds.; Springer-Verlag: West Berlin, 1987; Springer Ser. Mater. Sci. No. 4, p 134.
Figure 1. Contour plots of the projected electronic density for the S = 0 (left panels) and the S = 1 (right panels) states at time to = 0 (upper panels) and tF = lo5 au = 2.4 ps (lower panels). The dots are the projections of the ionic positions. TABLE I: Averaged Kinetic (EK), Potential ( E P ) ,Hartree (EH), Exchange-Correlation ( E x c ) and , Total ( E ) Electronic Energies for the S = 0 and S = 1 States‘ S=O
E, Ep E” E,, E
0.085 -0.296 0.122 -0.193 -0.283
(0.005) (0.015) (0.009)
(0.007) (0.010)
S = l 0.106 -0.216 0.074 -0.192 -0.227
(0.011) (0.042) (0,018) (0.014) (0.031)
one electron 0.052 -0.128 0.041 -0.094 -0.129
(0.008) (0.027) (0,011) (0.009) (0.020)
’In the third column the same quantities for a single electron also treated with BO-LSD are reported. The zero of our energy scale is such that E P = EH = 0.0 for uniformly spread electrons. Numbers in parentheses are estimated errors. For the S = 0 case, these were obtained by analyzing independent sections of a very long (- 106At) run. For the shorter (-3 X 105Ar)S = 1 and one-electron runs, errors were taken equal to the rms fluctuationsof the corresponding quantities. All values are computed excluding the initial section of the trajectory, typically of length 105At. Atomic units are used throughout (1 au = 27.21 eV). up to 10-20% are observed. These events are related to pseudo level crossings. Similar effects were observed in the case of a single solvated e l e c t r ~ nand , ~ as in that case, they are more likely to occur (at least the more dramatic ones) in the initial part of the run. However, in contrast to what observed in ref 5 , in the present calculations the electrons are seen to return quickly to the BO state after each event, so that here no special action is taken to deal with such situations. IV. Singlet vs Triplet State In this section we present the results of two parallel calculations for electrons with paired and unpaired spins and focus in particular on the relative energies of their states. In our LSD approach Pauli’s exclusion principle requires the two single-particle wave functions to be orthogonal in the triplet state. Instead, no restriction must be imposed for the singlet state. This largely accounts for the differences in the time evolution of the two states. The qualitative features of our results are summarized in Figure 1 where we show contour plots of p(x,y) = Sdz n(x,y,z) for the S = 0 (left panels) and S = 1 (right panels) states. For both cases two plots are given, one corresponding to the initial time to and the other to a subsequent time tF = 2.4 X s, representative of the equilibrium configurations. At time to the initial ionic
The Journal of Physical Chemistry, Vol. 92, No. 1 I , 1988 3271
Bipolarons in Metal-Metal Halide Solutions
,
(a:
‘I
,/’
,,,’
i
2.1.
~
’
I I 2 3 4
1.0
1,-
0.
2.
4.
6.
1
8.
lo.
J r(a.u.)
12.
Figure 2. Radial ion-ion pair correlation functions in the molten salt. configuration is the same for both runs. As is apparent from Figure 1, the singlet and triplet states look very similar at to, both being quite delocalized. The situation is very much different at tF. By this time the two electrons with antiparallel spins have paired, giving rise to a single-center localized state which looks like a somewhat spread-out F-center state, whereas the two electrons with parallel spins have repelled each other and formed two distinct F-center states. In order to make the differences betweeen the singlet and triplet more quantitative, in Table I we give the total energies of the two states together with the corresponding partial contributions of kinetic, potential, Hartree, and exchange-correlation energies. In order to be consistent and possibly cancel some of the errors by taking approximate differences, we have repeated the one-electron calculation of ref 5 using the LSD approximation. The results of this more approximate calculation are also shown in Table I and are very similar to those already obtained. It is apparent that the total energy of the singlet state is lower than that of two noninteracting electrons, whereas the opposite is true for the triplet. Therefore, two antiparallel spin electrons can form a bound state, whereas there is a repulsive interaction for the parallel spin case. The bound singlet state is what we call a bipolaron. Its estimated binding energy is of the order 0.7 eV, in reasonable agreement with recent experimental estimates. It is interesting to contrast the various contributions to the total energy for the singlet and triplet states. The larger spread of the electronic wave functions causes the kinetic energy to be lower in the S = 0 than in the S = 1 state. This energy gain however is completely offset by the larger Coulombic repulsion of the singlet state, since the two electrons are on average much closer. The exchange-correlation energies are not very much different for the two cases, whereas the largest energy gain for the S = 0 state comes from the electron-ion interaction. A direct quantitative comparison between the present and the PI results for the singlet state9 is hindered by the different choice of the zero of the energy scale in the two cases. However, all quantities which do not depend on such choice turn out to be rather similar in the two calculations. In particular, the average kinetic energy is 0.085 and 0.106 au in the present and in the PI calculations, respectively. Quite similar are also the electron-ion radial pair correlation functions (to be discussed below), with an average coordination of -4 in both calculations.
V. Structural Properties of the Bipolaron We now turn to a discussion of some static properties of two solvated electrons with antiparallel spin. As shown in the previous
0
2.
4.
6.
8.
10.
12.
r (a.u.)
Figure 3. (a) Radial pair correlation functions between the solvated electrons and the ions. (b) Same as above between the center of mass of the two electrons and the ions. Curves referring to the scale on the right give the coordination of the center of mass. section, the typical configuration for the two electrons is that of a bound complex-the bipolaron-which looks qualitatively similar to the F center. Our main purpose in this section is to characterize the distribution of ions around the bipolaron. It is useful to characterize first the structure of the ionic liquid. This is done in Figure 2 where we show the radial pair correlation function g(K+
