Electron-gas theory of the chemical shift

We construct a gauge-invariant approximation of the energy density of an interacting electron gas in the presence of a non-uniform magnetic field. The interaction is approximated by a slowly varying local potential. The magnetic field is a superposition of a constant field B and that due to a magnetic dipole μ, thus making the energy density suitable for direct calculation of chemical-shift tensors of, for example, interacting closed-shell systems. Unlike a similar theory for the magnetic susceptibility, the field-dependent exchange energy does not diverge.


Introduction
We present the basic requirement for an electron-gas calculation of the chemicalshift tensor: a gauge-invariant approximation to the ground-state energy density of an interacting electron gas in the presence of a constant magnetic field and a magnetic dipole.The energy density is a functional of the local electron density and the magnetic field strengths.Unlike the case of an electron gas in a uniform magnetic field, the direct electron-exchange contribution to the energy density does not diverge.
Probably the most important venue for the use of this energy density is that of systems of closed-shell components interacting in an intermediate range where electron exchange and overlap are important, but strong bonding interactions are not.For these systems an electron-gas approach to calculating that portion of the responses due to the interaction has been introduced by Kim and Gordon [1] in their successful calculation of the forces between the components.They wrote the interaction energy of the system as a functional of the density by approximating at each point in space the energy density by that of a uniform interacting electron gas of the same density.They then approximated the density of the interacting system by the sum of the densities of the non-interacting components, which they obtained from Hartree-Fock calculations.Heller et al. [2] have pointed out that the Hohenberg-Kohn theorem 1-3] implies that the interaction energy thus calculated differs from the actual interaction energy to second order in the non-additivity of the densities.For interacting closed-shell systems the non-additivity is generally small [1].
That a similar theory would work when a magnetic field is present is assured by the Rajagopal-Calloway [4] extension of the Hohenberg-Kohn theorem, which proves that the ground-state energy of the inhomogeneous interacting electron gas is a unique functional of the density and current.The key ingredient is the groundstate energy density of an interacting electron-gas in a magnetic field.It would have to be gauge-invariant, since even the simplest interacting closed-shell system--pairs of closed-shell atoms--has both a para-and diamagnetic response.The natural approach would be to write the energy density as a functional of the density and current, then assume additive densities and currents of the components.These quantities are presumably considerably easier to calculate than the corresponding quantities for the full system.If the components are closed-shell atoms, for example, only the diamagnetic contribution need be found.Harris and Cina ['5] tried this path several years ago.They constructed a gauge-invariant kinetic-energy density that was a functional of the density and current.It turned out, however, that this energy density was not useful for actual computation because it was multivalued in places.(Vignale and Rasolt [6] have recently calculated an exchange-correlation energy density of an interacting electron gas as a functional of the density and current.) The energy density of the interacting electron gas in a magnetic field can also be written as a functional of the density and explicitly of the magnetic field.Cina and Harris [7] did this for the slowly varying (locally constant) magnetic field, and, assuming additivity of the densities of the components, made the first electron-gas calculation of the diamagnetic susceptibility of the triplet hydrogen state.They had to leave out the lowest-order direct-exchange contribution, however, because it diverged.It is a long-known fact that this divergence occurs for an electron gas in a uniform magnetic field [8].This is a result of the long range of the Coulomb force and can be resolved by screening the exchange.It is not obvious, however, that without screening the lowest-order exchange contribution to the energy need diverge for a non-uniform, spatially limited field.(By 'spatially limited' we mean that the lowest-order term in the square of the field does not have amplitude over all space.) An important example of a non-uniform, spatially limited magnetic field is that in a nuclear magnetic resonance experiement: a superposition of a constant magnetic field and that due to a magnetic dipole.With the energy-density functional of an interacting electron gas in such a field, one could calculate the interaction contribution to the chemical-shift tensors of interacting closed-shell systems.This is the functional we have calculated and present here.It is gauge-invariant and computationally useful--i.e,mathematically well behaved.It also contains a finite exchange contribution, proving that in at least one case a finite exchange contribution to the energy can be obtained without screening.
In the next section we present the details of our calculation, but in brief we proceed as follows.First we calculate by time-dependent perturbation theory the propagator for an electron in the non-uniform magnetic field and a constant external potential.A time contour integral converts this to the density matrix of an inhomogeneous electron gas in the same magnetic field and a slowly varying local potential.From this density matrix we obtain the kinetic-and exchange-energy densities as functionals of the local potential.From the density matrix we also determine the relationship between the local potential and the density, which we then use to eliminate in favour of the density the local potential in the kinetic-and exchange-energy densities.Together with the usual potential-energy density, these two then give the complete energy-density functional.Because of the method of construction, which works directly with the fields, the energy density is manifestly gauge-invariant.
Following the section with our calculations, we remark on the finiteness of the exchange contribution and on how the energy-density functional may be used to calculate a chemical-shift tensor.

Calculation of the energy-density functional
We construct the energy-density functional at the level of the Thomas-Fermi-Dirac approximation.We begin with the Hartree-Fock approximation of the ground-state energy of a multi-electron closed-shell system: where e is the absolute magnitude of the electron charge, p(r) is the ground-state electron density, V=xt(r) and A(r) are the external scalar and vector potentials, and p(rlr') is the single-particle density matrix for electrons of the same spin.Were it constructed from the Hartree-Fock spin orbitals, this density matrix would be given by N p(rl r') = ~ ~b*(r)q~,(r'), where the ~b i are the Hartree-Fock spin orbitals of energy Ei.We write the energy as an integral overall space of kinetic, potential-and exchange-energy densities: where we then have, from (1), the following definitions: 1 I e 12 t[r, It, B, p(r)] -2 lim ~m p + -A(r) p(r'l r),  (4 a)

c)
This is a theory for closed-shell systems, so spin enters the picture only in the leading factor of 2, for double occupancy of each orbital, in the kinetic-and exchange-energy densities.We want to get each of these quantities as a functional of the fields It and B and of the electron density p(r).The potential-energy density is already in that form, so we need not consider it further.For the kinetic-and exchange-energy density we need the density matrix.It turns out that density matrix can be obtained from the single-electron propagator through the following time contour integral 19] : where ~ is the Hartree-Fock Hamiltonian, and ev is the Fermi energy, which for neutral systems we may set equal to zero.
We now set about calculating the propagator.First we make a local approximation to the non-local Hartree-Fock potential, so the single-electron Hamiltonian may be written as where f + ~(r).( 7) p(r) V~ is a local approximation to the exchange potential.It will turn out that we need never consider its exact form or its relation to the Kohn-Sham exchange potential 1-10].We next assume that V(r) varies slowly enough in space that we may ignore the fact that in general it does not commute with the kinetic-energy operator J -= (2m)-lip + (e/c)A(r)]2.The propagator may then be written as We evaluate V arbitrarily at the endpoint r instead of at the midpoint 89 + r') because we are ignoring variations of V over short distances, and we use the propagator only for r very near to r'.It is important to emphasize that A(r) is not treated as slowly varying even through V(r) is.This is to some extent inconsistent, since a non-slowly-varying vector potential might be expected to produce a nonslowly-varying density p(r) and hence local potential V(r).We have not fully explored the consequences, if any, of this inconsistency.We shall assume here that they are unimportant.Finally, we restrict ourselves to a chemical-shift calculation; that is, our magnetic field consists of a constant magnetic field B and that due to a magnetic dipole IX.The elements of the chemical-shift tensor are obtained from the energy as That is, the tensor is the coefficient of the I IX[LB I term in an expansion of the energy in powers of the magnetic field strength.Therefore we need keep terms in our calculation of the energy density only up to bilinear order in [ IX I and [ m I.This will be an important point later.Choosing the Coulomb gauge and fixing the magnetic dipole at the origin, the vector potential is given by The kinetic-energy operator can then be written as p2 e e IX x r e 2 (B x r) 3-"~2--mm+~-mcP'(Bxr)+--P'---r'T--+--(IXxr)'--+mc 2mc2 ra ..., (11) where, as promised above, we have dropped terms of order I It [ 2 or I B [ 2 or higher.We denote the first two terms of J as ~--o, and the third and fourth terms as 3-u and 9"-~B respectively.Ordinary time-dependent perturbation theory then gives e -itg-/~ = e-it~~ (12) where The reason that we do not choose the free (no magnetic field) kinetic-energy operator p2/2m as our zeroth-order Hamiltonian is that, in order to include all terms bilinear in I ltl and I B 1, we should have to do second-order perturbation theory on it.We cannot perform all of the required integrals.Luckily, this turns out to be unnecessary.
Because the two terms of 9-0 commute, we may approximate the exponential:

(rtl B 9 II r't') = B 9 (r x r') t -t----;
We may now substitute this back into (8) and use the resulting form of the propagator in (5) to obtain the density matrix.The time contour integral over the free-particle propagator delivers the well known zeroth-order density matrix [11] fcdt where jl is the spherical Bessel function of order 1, and we define the parameter ~ by The t' integral in (20), over the product of free-particle propagators, has been performed by Feynman and Hibbs [12]: The density matrix that we obtain by doing these time integrals in ( 5) is We now calculate the kinetic-and exchange-energy densities.The kinetic-energy density (4 a) is found by applying the kinetic-energy operator given by (11) to the density matrix:

The time integral over this result gives fc dt e-iW/~ fi'dt' (rtl,t')(~t'lr'O) _
Unfortunately we have been unable to perform this entire integral in closed form.However, it turns out that it can be written in the following form: where Xe 1 and Xe 2, given in Appendix B in ( 66) and (67), while not elementary functions or combinations thereof, may be tabulated for their single parameter.
Figure 1 shows Xe ~ and Xe 2 for a range of the parameter.
The kinetic-and exchange-energy densities given in ( 27) and ( 30) are functionals of the local potential V(t).We need them as functionals of the electron density p(t), so we must substitute into them for the local potential its expression in terms of the electron density (and magnetic fields), F[p(t), r, It, B].Remember, however, that we are calculating the energy density only to first order in I It l and I B I. Therefore, for parts of ( 27) and ( 30 The functions Xe 1 and Xe 2 appearing in (38), as given by ( 66) and (67).

+
where When combined with the potential-energy density (4 b) due to the mean-field and external potentials, (37) and (38) give the total energy density of the inhomogeneous electron gas.Integrating the energy density gives the energy, and one can simply (37) The functions D 1 and D 2 are finite for all positive x.We do not have to contemplate inverting (33) to obtain V[p(r), t, It, B] exactly, because we need the latter only to first order in I P l and I B I. We can do that in the following manner.In each term on the right-hand side of (33) where It and B already appear explicitly, we substitute for a the zeroth-order (Thomas-Fermi) result (32).The resulting equation then contains just one V(r) on the right-hand side, and so can be easily inverted to give V[p(r), r, It, B] to first order in I It l and [ B I. After a final expansion in terms of the magnetic field strengths, we obtain for ~ to first order in I It I and I B I the following: where a = [3~2p(r)] 1/3r (36) We now use (35) to write the kinetic-and exchange-energy densities in their final forms as functionals of the electron density: read off from that the chemical-shift tensor, which will be the coefficient of the term bilinear in I P l and I B I.

Remarks
We wish to emphasize that, in distinct contrast with the situation for a constant magnetic field, the exchange energy density given by (38) will give a finite exchange contribution to total energy for any (reasonable) bound system.Since the functions D 1, D 2, Xe 1 and Xe 2 are everywhere finite, for large r the right-hand side of (38) goes to zero as fast as p(r) itself does.To sketch this, we show in figure 2 the exchange energy density as a function of r for the simple case of p(r) = rr-le-z,.
The energy density that we have constructed can, as we mentioned, be used directly in the calculation of chemical-shift tensors of interacting closed-shell systems.To do that, one would first calculate the density of the isolated components of the system in the presence of the fields.If the components are simply closed-shell atoms--the situation we have in mind--this is enormously easier than calculating the full magnetic response of the system, because the atoms have only a diamagnetic response.The density in the energy-density functional we have constructed can then be approximated by the sum of the densities of the isolated components.The error in the energy resulting from assuming additivity turns out to be second order in the non-additivity, as mentioned in section 1.To isolate the interaction part of the energy, one would subtract off from the energy density the energy density of the isolated components.Integrating the energy density would then give the interaction energy of the system, and the coefficient of the term of first order in Ipl and I B I is the part of the chemical-shift tensor due to the interaction.We are presently carry- ing out a calculation of the chemical-shift tensor for the worst candidate among interacting-closed shell systems for an electron-gas calculation, the triplet state of H 2 [13].A somewhat more complicated calculation done for interacting 129Xe atoms could be compared with the experimental observations of Jameson et al. [14].
(Cynthia Jameson [15] has determined one-third of the trace of the chemical-shift tensor for 129Xe at low densities, where pair interactions dominate, as a function of pair separation.Calculating this function would be a good use of the electron-gas theory of the chemical shift.) We wish to make two further comments.First, for this theory--unlike a theory in which both V(r) and A(r) vary slowly--it should be possible to correct the deficient treatment of the variation of V(r) by making a gradient expansion in V(r).Secondly, one could produce an exchange potential for a slowly varying local potential by varying with respect to p(r) the integral of the exchange-energy density given by (38).This might be a starting point for a somewhat more general theory, such as a consistent Thomas-Fermi-Dirac theory with these magnetic fields.

Appendix A
The density matrix Here we sketch the final calculation of the density matrix from (25).

+lR-89189 89 ~ ~2
The variable x is the cosine of the angle between i and R. The tilde above the various quantities indicates that ~ has been set equal to one.Also, in terms of the new coordinates, we note that =IR+ 89189 = ?.
and Xe 2 are shown in figure 1 for a likely range of the The two functions Xe a parameter ~R.

exp h 1 h 2mc (
note that I is ordinary angular momentum).Once again, we have dropped terms of higher order than we need.If we define the ket I rt)---exp ~m l r)(15) and the corresponding bra then the kinetic term in the propagator to the appropriate order in [ It[ and [ B [ can be written as it e(rle-i'r/~l r') = (rtl t'0) --(rtlB 9 II r'0) h 2mc if/ Itxr +-~ dt' mc-e (r(t-t')lP "-7-Ir'(-t')) i e 2 [ pxr + ~ 2m2c 2 (t --t')(rtl(B 9 I)p" ~ [ r'(t')) (B x r).--~ [r'(-t')).(16)Note that with our definitions, (rt [ r'0) is the ordinary free-particle propagator: (rtlr'0)=(rlexp( it p2 h 2--mm/Ir') \ (17) We may now simplify (16) by the use of two relations.The first is simply the position representation of the operator It 9 I : it.I~ -ih]/t I dq~.(18)where (p, would be the azimuthal angular coordinate if It were to lie along the z axis.The second relation is the following, which can easily be verified by direct calculation: m (rt I r't').

f+ 2mc 2
PjPm jX(~I r--~l + ~lr--r'l) Bi/llgijk~'klm h E 4x 3 Ir -rl Ir -r'l where d Jk is the Levi-Civita symbol for the completely antisymmetric tensor.The final integral over ~ can be done in closed form by changing to prolate spheroidal coordinates and making use of an expansion of 1/I r -r l in these coordinates.The algebra is complicated, however.The explicit calculations are sketched in Appendix A. The resulting density matrix is given there by (62).

r
:;...:.,.~,~ .................The exchange-energy density given by (38) as a function of r for the simple case of 1 2r p(r) = ~-e-.The dashed line is that part of the exchange energy density multiplying p 9 B, the dotted line is that multiplying (p x r) 9 (r x B)/r z, and the solid line is that with no magnetic fields in it.