An electron gas calculation of the nuclear magnetic shielding tensor of 3 H2

We calculate the interaction nuclear magnetic shielding tensor for the 3Σu + state of H2, using the electron gas theory that we have presented previously. We confirm the finite direct contribution of the exchange energy, and our general results agree qualitatively with a simple Heitler-London calculation. Comparison with experimentally derived results for 129Xe suggest that an electron gas calculation on that system would be valuable.


Introduction
Recently we presented the essential requirement for electron gas calculation of the magnetic shielding tensor of nuclei [1], a gauge-invariant approximation of the ground state energy density of a locally homogeneous electron gas in the presence of a constant magnetic field and a nuclear magnetic dipole. We pointed out that the theory would probably have greatest utility in the calculation of the part of the magnetic shielding due to interatomic interactions in systems of closed-shell components at intermediate distances from each other, where electron exchange and overlap are important, but long-distance induction and dispersion forces are not. In this paper we apply the electron gas theory to the prototype for interacting closedshell systems, the non-bonding triplet state of H2. This system is an electron 'desert' rather than gas, so we do not expect the results so much to be quantitatively correct as to, if qualitatively so, encourage the application of the theory to electron-rich heavier systems where the electron gas approximation holds more nearly true. It turns out our results are qualitatively correct, to the extent that can be determined, and thus the electron gas approach to nuclear magnetic shielding shows promise in its first test.
We begin with a very brief review of the origin and nature of our approach. Gordon and Kim [2], in their calculation of interatomic forces, introduced the idea of writing the interaction energy of a system of closed-shell components as a functional of the electron density by approximating at each point in space the energy density by that of a uniform interacting electron gas of the same density. This electron gas approximation is equivalent to the assumption that locally the kinetic and potential energy operators commute [1]. This is most nearly true in regions of large and slowly varying potential, that is in regions far from atomic cores where electron density is nevertheless high. The second step in the Gordon-Kim approach is to approximate the true density of the system that appears in this energy density functional by the sum of the densities of the non-interacting components. Heller et al. [3] have pointed out that the Hohenberg-Kohn theorem [4] implies that the interaction energy density thus calculated differs from the actual interaction energy density 0026-8976[91 $3.00 9 1991 Taylor & Francis Ltd to second order in the non-additivity. The interaction energy density can be integrated directly to obtain the interaction energy.
We should like to make several points about this approach. First of all, the theory still requires the ground state densities of the non-interacting components to be found accurately by some other method. This is considerably easier than calculating the density of the whole system, however, because of the high symmetry of the components, which makes the trade-off advantageous. Note that an electron gas approach (e.g. ordinary Thomas-Fermi theory) would in general give rather poor results for the isolated component de~sities, since near atomic cores the electron gas assumption is rather badly violated. Indeed, the success of the Gordon-Kim method rests to some extent on the subtraction, at the energy density stage, of the rather incorrect isolated component contributions to the energy. What is left in the energy density integrand is then the contribution from regions of electron overlap, where the electron gas approximation holds best.
The second point is that the success of the method relies on the non-additivity being relatively small, which is why it is not useful for systems where strong bonding interactions exists, or where very close approach of the component is important.
Thirdly, the assumption of additivity means that long-range induction and dispersion forces are not included in this approach, which is why the asymptotic behaviour of the system is not described well. This is pointed out by Harris [5], who also showed how induction and part of the dispersion forced could be included in an electron gas approach, albeit at the price of more work. We have not included long-range forces here.
Finally, it is a little paradoxical that the approach works at all: the assumption of additivity is inconsistent with the Hohenberg-Kohn and Hellmann-Feynman theorems, and yet it is precisely that inconsistently that gives non-trivial answers, since strict additivity would give an interaction energy of strictly zero. Harris and Heller [6] have discussed this point more fully.
That a Gordon-Kim type of approach would work when magnetic fields are present is assured by the Rajagopal-Colloway extension of the Hohenberg-Kohn theorem [7], which proves that the ground state energy of the inhomogeneous interacting electron gas in the presence of a magnetic field is a unique functional of the electron density and current. Of course, it remains to find that functional! Harris and Cina [8] constructed a kinetic energy density functional, for the case of a slowly varying magnetic field, which was a functional of the density and current, in the hopes of using additive densities and currents to write the energy density. However, the functional was miltivalued in places and could not be so used. Vignale and Rasolt have also recently constructed an exchange correlation energy density as a functional of the density and current [9], but its calculational usefulness has not yet been determined. Cina and Harris [10] also found an energy density for the case of a slowly varying magnetic field as a functional of the density and explicitly of the magnetic field, and used it to calculate the diamagnetic susceptibility of the triplet state of H2, but they found that the lowest-order direct exchange energy contribution inconveniently diverged, and could not be included in the energy density. Stephen [11], in considering the case of a purely homogeneous electron gas many years ago, noted this problem, and found that screening the exchange resolved it. Our previous work is the construction, in a manner similar to that of Cina and Harris, of a complete energy density for the locally homogeneous electron gas, in the presence of a specific kind of magnetic field: that in an ordinary nuclear magnetic resonance (NMR) experiment, a super-position of a constant field and that due to a nuclear magnetic dipole. We found that by not treating this field as slowly varying (locally constant in the sense the potential is treated), it is not necessary to screen the exchange to obtain a finite lowest-order direct exchange contribution to the energy. We believe this to result from the fact that the magnetic field is 'spatially limited' in the sense that the lowest-order term in the square of the field does not have amplitude over all space.
Our energy density functional is correct to first order in the nuclear dipole moment and in the strength of the constant field. If we find the ground state densities, to similar order in the fields, of closed-shell components, we can use the Gordon-Kim approach to calculate the interaction energy of the interacting components as a function of the fields and the arrangement of the components. The coefficient of the term in the interaction energy bilinear in the two fields is then the interaction nuclear magnetic shielding tensor as a function of the nuclear position vectors.

The spin-polarized energy density functional
In our previous work we required the system to be closed-shell, i.e. to have total electron spin S = 0, so that if we assumed that the ground state wavefunction could be written as a product of one-particle orbitals then a single Slater determinant of these orbitals would suffice. All physical properties could then be calculated using the physical space part of the orbitals, with multiplicative factors of two in the right places to account for the double occupancy of each physical space orbital. The triplet state of H2 is not, of course, a closed-shell system, but it is completely spin-polarized, meaning the azimuthal component of spin is at its maximum, that is m s = S. For spin-polarized as for closed-shell systems, just one Slater determinant is needed, and all physical properties can be found from the physical space orbitals. However, each single-particle orbital is singly instead of doubly occupied, because the spins of the electrons are parallel. For a locally uniform electron gas this means that at a particular point, with local potential V(r) and local Fermi energy ~F(r), the total spin-polarized electron gas density p~ (r) = 89 where p(r) is the total electron density when each physical space orbital is doubly occupied. To convert the results of our previous work to those appropriate to a spin-polarized electron gas, we need to insert factors of two in appropriate places. For example, from the ordinary field-free Thomas-Fermi-Dirac energy functional [12] V TED [p(r)], which gives the local potential from the local density, we can find VT TFD, the local potential for a spin-polarized electron gas, as a function of the spin-polarized density: Carrying through the procedure logically consistent with this on the kinetic and exchange energy densities given in (37) where the variable a is defined as The functions D 1 , D 2, Te I and Te 2, as given in [1], are where j0 (x) is the spherical Bessel function of order zero and si (x) is the sine integral function. The functions D l , D 2, Te I and Te 2 are finite everywhere. The functions Xe ~ and Xe 2 we defined in [1]. We cannot write them in closed form, but they have been tabulated [13] as functions of their single variable. The electron gas potential energy density is just that from ordinary mean field theory: where Vn is the nuclear Coulomb potential. The first term comes from the interaction of the electrons with the nuclei, where the second is the classical electron-electron repulsion energy. No modification of the form of this functional is necessary for spin-polarized systems. Note that this is just the electronic potential energy; to arrive at the total potential energy we should add the nuclear repulsion. The latter does not, however, contribute to the magnetic shielding tensor.

Densities of the isolated hydrogen atoms
What we have said in the previous section is still completely general. Now we must specialize to the triplet H 2 system that we are considering and determine the input to the energy density functional, which is the densities of the isolated hydrogen atoms in the presence of the two magnetic fields. This is the place where each application of the electron gas theory requires serious work; however, once the work is done for any particular set of components, the magnetic shielding tensor of a system consisting of any non-bonding arrangement of those same components can be calculated with one three-dimensional integral over the energy density functional. For this symmetric system the nuclear magnetic shielding is identical for each proton, so we arbitrarily designate the atom with the nuclear magnetic dipole under consideration as Ha. The other we call H b. We write the densities, in the presence of the dipole p and constant magnetic field B as the field-free ground state density p0 plus some small correction due to the fields: where R is the internuclear separation. The p0b are the usual hydrogen atom ground state densities, in atomic units: We calculate the additions pl,b(/~ , B) using ordinary time-independent perturbation theory, in the manner of Dalgarno and Lewis [l 4]. The energy density of the interacting system is determined by using p~ + Pb for the density p(r) in the energy density functionals above. To find the interaction energy density, we then need to subtract the energy densities of the isolated components. We find the isolated energy density for Ha by inserting p~ into the energy density functionals above. For Hb, however, the isolated energy density is found by inserting pO alone into only the field-free part of the functionals above, because when H b is infinitely far from H a no contribution to its energy can come from terms involving/~, and only the field-free parts of Pb, t[p(r)] and x[p(r)], contain no/l.
We now describe the calculation of the P~.b(/~, B). The effects of the dipole on Ha are a special case (with internuclear separation going to zero) of the effects of Hb, SO we begin with the latter. We locate the atom at the origin. A magnetic dipole at R and constant magnetic field lead to the following extra terms in the Hamiltonian: mc r 3 q-~mc B" 1 + 2mc~ r3 We have dropped terms higher than first order in ~u and B, and ignored the spin-magnetic-field interaction. The notation IR indicates angular momentum about R, that is IR = (r -R) • p. The first two terms of (10) are the paramagnetic terms: to evaluate the change that they make to the wavefunction in the general case, to first order in p and B, we should need to do second-order perturbation theory, involving in the traditional approach a sum over all the excited states. However, because the ground state density must be real, the paramagnetic terms cannot contribute in first order, and, because the ground state of the hydrogen atom has zero angular momentum, they contribute nothing in second order. We need only deal then with the final, diamagnetic term in (10). This requires only first-order perturbation theory, since the diamagnetic term already contains/~ and B explicitly. Therefore the perturbation theory equation that we must solve for ~b~, the first-order change in wavefunction, is

{ I[p • (r-R)]'(B • r)}
(~b ~ --E~ = Eb ~ 2C 2 r3 ~k ~ ( 1 1) in atomic units, where O ~ is the lowest eigenfunction of the field-free Hamiltonian ~o and E ~ is its eigenvalue. The first-order energy shift Eb' can easily be calculated as el = e 2" +N+N +R g3 The term multiplying c-2/, 9 Bwe call E 1 , it tends to 89 as R ~ 0. The term multiplying c 21~zB z (if we let R lie along the z axis) we call E~, and it of course tends to zero as R does. Expanding both sides of (1 l) in spherical harmonics and equating coefficients gives us a second-order inhomogeneous differential equation, in just the radial coordinate r, for each coefficient functionft(r) multiplying the spherical harmonic of order l in the perturbation wavefunction expansion. These we solve by numerical integration [13]. We orthogonalize the l = 0 coefficient to the field-free wavefunction specifically, and, by that and the orthogonality of the spherical harmonics, ~b~ is made orthogonal to ~b ~ In practice we need include only the first few spherical harmonics and coefficient functions in ~b 1, because for higher l the coefficient functions become small and the spherical harmonics oscillate quickly, cancelling each other out. Our experience is that the first nine are more than sufficient. (As a side comment, we should point out that the retention of more and more spherical harmonics in the wavefunction changes it hardly at all, with one interesting exception: very near to and centred on the dipole, the spherical harmonics and their coefficient functions add coherently in a smaller and smaller region, building what for the complete solution is probably a strong singularity. However, this singularity is such that spherically symmetric integration over it with any non-singular function gives zero. We have noted before this property of Dalgarno-Lewis-type perturbation wavefunctions for singular perturbations [15], and the problems of integration over them has been discussed in a more general context by Pitzer et al. [16].) To calculate the first-order change in density for Ha is considerably easier than for Hb. The dipole is now located at the origin, where the atom is, so the perturbation equation that must be solved is Note that the gauge origin of the constant magnetic field has been changed from that in (1 l). Although this gauge transformation can affect the perturbed wavefunction Ola, it cannot affect the resulting perturbed density pal. The first-order energy shift is found to be 89 9 B, by direct calculation or by letting R ~ 0 in (12). We follow the same logic as in the determination of Ob above--the expansion in spherical harmonics, finding coefficient functions, and so on. However, the expansion of the perturbation stops with l = 2, and only three coefficient functions are non-zero, so the process is considerably abbreviated.

Results
Because of the symmetry of the H2 system, the interaction nuclear magnetic shielding tensor a 1 (R) has just two independent components. If the z is the interwhich are identical.
nuclear axis then these are a~, which we label a I , and alx~ or ayy, We define a~ = 1 1 0.1 ~(axx + yy), in accordance with the usual notation for diatomic molecules; the anisotropy Aa ~ is then conventionally defined as a I -a~. The scalar interaction magnetic shielding a ~ is one-third of the trace of the full interaction magnetic shielding tensor. Note that this is all that survives rotational averaging, and hence is what is determined in ordinary gas phase experiments.
Figure l(a) shows a l, alll and Aa I for the triplet H2 system, as functions of internuclear separation R, calculated with the electron gas theory. It should be borne in mind that the full magnetic shielding tensor a(R) would be a 1 (R) plus the shielding tensors of the two isolated hydrogen atoms. The latter are just 89 -2 times the unit tensor for H~, and zero for Hb. Figure l  interaction magnetic shielding from the kinetic, exchange and potential energies. Note the significant contribution of the exchange interaction energy. Figure 2 shows the portions of the exchange energy contribution where the fields come in directly and indirectly, the 'direct' portion comes from integration of that part of the exchange energy density x[r, F, B, p(r)] that contains the field explicitly, while the 'indirect' portion contains the fields only by virtue of the field dependence of the non-interacting densities pa(r, F, B) and pb(r, R, F, B). As mentioned in section 1, Cina and Harris [10] showed that the direct exchange contribution to the electron gas energy diverges if the magnetic field is treated as slowly varying, and cannot therefore be included in an electron gas calculation. From figure 2 it can be seen that leaving out the direct exchange energy contribution here would be a significant error. It should be recalled that no screening of the exchange is included in this theory. The finite direct exchange energy contribution is the first interesting result of this calculation. The second is the existence of a minimum in the deshielding effect of the interaction at intermediate distances near 'contact' (R = 2a0). There is no rigorous calculation of the triplet shielding tensor available for comparison with the electron gas results, but there is one other ab initio calculation at these short distances, by Marshall and Pople [17] using the Heitler-London molecular wavefunction. The atomic basis orbitals Xs that they use are simply hydrogen orbitals multiplied by a phase factor: in atomic units, where ~b H is the hydrogen orbital and R, is the location of nucleus s. These atomic orbitals are correct to first order in B. This follows from the observation that the hydrogen orbitals are unchanged to first order in B if the gauge origin is the location of the nucleus. Moving the gauge origin to -R,, which is equivalent to keeping it at zero and moving the nucleus to R~, is a gauge transformation, and modifies the wavefunction by the phase factor given above. Of course, the molecular wavefunction constructed from these orbitals is not correct to order B, because of the inexact formulation of the former from the latter by the Heitler-London prescription.
It becomes more accurate at large separations because the Heitler-London wavefunction becomes exact asymptotically.
Marshall and Pople used their wavefunction to calculate the magnetic field induced at the nucleus with the magnetic dipole by the external field; the coefficient of the term linear in the external field is then the magnetic shielding. (It is not readily apparent in Marshall and Pople's paper that this is what they are doing, but careful consideration of their calculation shows its formal and exact equivalence to use of the Biot-Savart law to find the induced magnetic field from currents in the wavefunction induced by the external magnetic field.) The electron gas approach, by contrast, uses a density that has both y and B in it to find the interaction energy, whereupon the coefficient of the term biIinear in the fields is the interaction magnetic shielding tensor. This 'energy' approach would naturally give the same answer as the 'induced-field' approach used by Marshall and Pople, were the wavefunctions exact. If approximate wavefunctions were used that satisfied the Hellmann-Feynman theorem [18], even trivially, for the two fields, then the two approaches would also give the same answers.
The difference between the two approaches is not therefore necessarily an indication of the accuracy of either, but we cautiously take it to be so, since the methods of calculation differ significantly.
A comparison of the scalar magnetic shielding calculated by Marshall and Pople with the electron gas results is shown in figure 3, along with, for comparison, the results of the perturbative calculation of Rummens [19], made to order 1/R ~~ which is good at large distances. The agreement between the electron gas and Marshall and Pople calculations is qualitative, but certainly not quantitative. (It is rather encouraging that there is even qualitative agreement, given the doubtfulness of the electron gas assumption for a two electron system.) Although it is not easily seen in the figure, the    Jameson's [20] empirically derived interaction chemical shift as a function of , function shielding is rather more negative at intermediate distances than the Marshall and Pople shielding, and l~h~ latter does not display a minimum. We suggest, however, that this minimum is a real feature, missing in the Marshall and Pople calculations, although its position and depth for the triplet H 2 system we do not expect to get right with electron gas theory. We base this first on the physical intuition that the loss of spherical symmetry, resulting from the interaction, that reduces diamagnetic shielding begins to reverse as the interacting atoms get very close. There is also indirect experimental evidence that this minimum exists, at least in heavier systems. Jameson [20] has extracted from experimental data obtained [21] in her laboratory an empirical estimate of the interaction chemical shift, as a function of internuclear separation, for pairs of 129Xe atoms, a real interacting closed-shell system. Her results are shown in figure 4. Note that, in parallel to what we find in the triplet H2 system, there exists a minimum in the interaction deshielding near 'contact' (the equilibrium separation Ro). Figure 4 also shows the function of the xenon interatomic potential Jameson found that best reproduces the experimental data. While there are no compelling reasons for the interaction chemical shift to depend in a simple way upon the interatomic potential, it is interesting that, of the class of functions that do so, the one that works best is very similar to our interaction magnetic shielding as shown in figure 1. We feel that it would be worthwhile to perform an electron gas calculation for the 129Xe system, where the plethora of electrons lends the electron gas assumption considerably more accuracy in the region of importance. With a good xenon interatomic potential [22], the temperature and density dependence of the magnetic shielding of 129Xe could then be directly calculated and compared with Jameson's experimental results. Our preliminary work on triplet H 2 suggests there is a good chance this comparison will be favourable, since the minimum in the interaction magnetic shielding function, at least the most obvious feature demanded of it for agreement with experiment, which is missing in previous theoretical formulations, is here found even in the electron desert limit of the electron gas. We wish to make three further comments. First of all, it should be borne in mind that the electron gas results are expected to be best at intermediate distances. At large separations, where overlap is small (here where R is greater than about 4), the main contributions to the energy are from induction and dispersion forces, of which the electron gas theory takes no notice, since it uses the additive density approximation. At small distances, the non-additivity is large; furthermore, the potential changes rapidly and the Thomas-Fermi condition no longer holds (or more precisely is even less accurate) in the region of overlap, whence most of the interaction energy comes. We estimate this region as being R less than about 1. The short-distance region, we could argue, is less important for ensembles of atoms, since close approaches are unlikely because of the strongly repulsive nature of the potential at short distances. As for the asymptotic region, as Harris has shown, it is in fact possible to include induction forces in the electron gas theory--at the price of a fair amount more work in the calculation of the constituent component densities. It may also be possible to make some kind of interpolation between the electron gas results at intermediate distances and the results of perturbative calculations, which are reliable at asymptotic distances.
Our second comment is that the NMR spectrum of spin-polarized atomic hydrogen at low densities has been recorded by Johnson et al. [23]. Since the interatomic potential of triplet H2 is known, the density dependence of the nuclear magnetic shielding for spin-polarized hydrogen could be calculated directly from our results, perhaps using the perturbative results for asymptotic distances, and compared directly with experiment. (Johnson and colleagues have not actually published the density dependence of the magnetic shielding, but their published results make clear that it could easily be measured with their technique, if they have not already done so.) Finally, of some interest are the results of our electron gas theory applied to direct calculation of the complete (not interaction) nuclear magnetic shielding of just the hydrogen atom itself, because they are curiously good considering that a single electron in a Coulomb potential is very far from an electron gas, and in fact the Thomas-Fermi condition is satisfied nowhere.
This latter should exactly cancel with the exchange energy, since for just one electron We add these potential energies to the results of integrating the kinetic and exchange energy density functionals above using the hydrogen atom unperturbed density. The The electron gas approach gives a kinetic energy that is too small (the exact value is 89 but this is compensated for by the exchange energy also being too small. Thus the remarkable accuracy of the total energy seems somewhat fortuitous. The exact nuclear shielding tensor for the hydrogen atoms is 89 -2 times the unit tensor. The electron gas results are, in units of c -2, Here the electron gas theory overestimates the exchange energy, but this is compensated for by an overestimation of the kinetic energy, leading to a final result that is surprisingly good.