After some years of controversy, it was recently demonstrated how to obtain the correct long-distance limit [point-dipole approximation (PDA)] of pseudo-contact nuclear magnetic resonance chemical shifts from rigorous first-principles quantum mechanics [Lang et al., J. Phys. Chem. Lett. 11, 8735 (2020)]. This result confirmed the classical Kurland–McGarvey theory. In the present contribution, we elaborate on these results. In particular, we provide a detailed derivation of the PDA both from the Van den Heuvel–Soncini equation for the chemical shielding tensor and from a spin Hamiltonian approximation. Furthermore, we discuss in detail the PDA within the approximate density functional theory and Hartree–Fock theories. In our previous work, we assumed a relatively crude effective nuclear charge approximation for the spin–orbit coupling operator. Here, we overcome this assumption by demonstrating that the derivation is also possible within the fully relativistic Dirac equation and even without the assumption of a specific form for the Hamiltonian. Crucial ingredients for the general derivation are a Hamiltonian that respects gauge invariance, the multipolar gauge, and functional derivatives of the Hamiltonian, where it is possible to identify the first functional derivative with the electron number current density operator. The present work forms an important foundation for future extensions of the Kurland–McGarvey theory beyond the PDA, including induced magnetic quadrupole and higher moments to describe the magnetic hyperfine field.
I. INTRODUCTION
Nuclear magnetic resonance (NMR) observables, such as resonance frequencies and linewidths, are strongly affected by the presence of unpaired electrons. These perturbations encode information about the molecular structure and the electronic structure and, therefore, represent a precious source of information in modern chemistry. As a consequence, the study of paramagnetism in NMR (pNMR) has attracted increasing attention,1 from both experimental and theoretical standpoints.
For closed-shell systems, which have a nondegenerate ground state usually well-separated from any excited states, it has been known for many decades that the chemical shielding tensor can be calculated via the equation originally derived by Ramsey.2,3 Moon and Patchkovskii derived an equation for pNMR chemical shifts for individual Kramers doublets in terms of spin Hamiltonian (SH) parameters,4 an approach that was generalized to arbitrary spin degeneracy by Pennanen and Vaara.5 In both cases, the chemical shielding tensor is defined as the mixed second derivative of the Boltzmann-averaged energy. Van den Heuvel and Soncini calculated the shielding as the derivative of the averaged magnetic hyperfine field6 and obtained expressions that are also valid in the strong spin–orbit coupling limit. They successfully applied their theory to pNMR shifts in lanthanide systems. Shortly after, it was found by the same authors that the chemical shielding tensor can actually be considered as the mixed second derivative of the electronic Helmholtz free energy.7 This resulted in one of the most general equations for chemical shieldings developed so far and showed that earlier treatments based on the Boltzmann-averaged energy instead of the free energy can lead to wrong conclusions. The Van den Heuvel–Soncini equation can be employed practically in any situation, as long as the basic assumption of the different time scales of electronic dynamics and nuclear spin dynamics and the weak-field approximation (i.e., that shieldings are independent of the magnetic field strength) are fulfilled. The same authors also used their theory to obtain expressions in terms of SH parameters for a single spin state subject to zero-field splitting (ZFS).8 Gendron et al. were the first to evaluate the Van den Heuvel–Soncini equation in a sum-over-states (SOS) manner, without recourse to SH parameters.9
An alternative route to the calculation and interpretation of pNMR shifts that has been developed since the 1960s separates isotropic chemical shifts into “contact” shifts, i.e., shifts arising from the Fermi contact (FC) interaction, and “pseudocontact” shifts (PCSs), arising from all other interactions. The isotropic shifts result after rotational averaging of the full chemical shielding tensor and only include contributions to the latter that have a non-zero trace. Starting from the SH approximation, McConnell and Robertson10 used the point-dipole approximation (PDA) to derive an equation for PCSs in terms of the g values of the paramagnetic center and structural parameters. Their work was extended by Jesson.11 Kurland and McGarvey generalized their results under weaker assumptions and showed that the previous equations were special cases of an equation that expresses the PCSs in terms of the susceptibility tensor and structural parameters,12
Here, Δχax and Δχrh are parameters that describe the anisotropy of the susceptibility tensor of the paramagnetic center, θ and ϕ are the polar and azimuthal angle of the direction of the nucleus in the principal axis system (PAS) of the susceptibility tensor, and R is the distance between the nucleus and the paramagnetic center. It should be stressed that compared to the Van den Heuvel–Soncini equation, which is valid for any distance, Eq. (1) is more approximate and only strictly valid in the long-distance limit.
The dependence of PCSs on geometrical parameters, as given in Eq. (1), has been used extensively to extract information about the molecular structure of metalloproteins and of small complexes. This dependence can also provide insight into conformational heterogeneity.13–15 More recently, the improved computational efficiency and quality of relativistic wavefunction-based quantum chemistry (QC) methods16 made it possible to rapidly and reliably calculate the magnetic susceptibility of metal complexes from first principles.17–19 Therefore, the dependence of PCSs on the magnetic susceptibility has been used to gain a deeper understanding of the metal coordination, for instance, in paramagnetic metal complexes,20,21 including single ion magnets,22 and even to refine the structure of the active site of cobalt(II)-bound enzymes.23
Equation (1) for the PCS can also be expressed via
in terms of the pseudocontact contribution to the full chemical shielding tensor,
This is one of the central equations around which the current work revolves. We shall call it the Kurland–McGarvey (KM) equation, although it did not appear in the full tensor form in the original paper. Incidentally, the same equation was already derived by McConnell in his treatment of the NMR chemical shielding in closed-shell molecules arising from distant electrons.24 Equation (3) has a clear physical content that can be intuitively understood as follows: For a nucleus K at position RK that has a sufficiently large distance from the paramagnetic center, the magnetic hyperfine field due to induced currents can be described in the PDA, i.e.,
where ⟨μ⟩ is the average induced magnetic moment. The latter can be expressed in terms of the “single-molecule” magnetic susceptibility (which has SI units of m3) via
Inserting this into Eq. (4) gives
which, compared with the expression ⟨BHF(RK)⟩ = −σT · B for the induced magnetic field in terms of the chemical shielding tensor, immediately leads to Eq. (3).
In 2016,25 it was proposed that the susceptibility tensor χ in Eq. (3) should be replaced by an unsymmetrical tensor χ′, where in the SH approximation, .15 This can lead to significantly different results compared to the original KM equation.1,26,27 The derivation of the equation containing χ′ started from the QC expression for the chemical shielding in terms of SH parameters and assumed that the only relevant contribution to the A tensor is the spin–dipolar one, in particular neglecting the orbital contribution induced by spin–orbit coupling (SOC). Several papers found apparently good agreement of this approach with experimental results.28–30 However, as already noted before,31 the uncertainty in the calculated SH parameters (g and D tensors) can lead to uncertainties that are of the same order of magnitude as the expected difference of the two approaches. Therefore, the better performance of the equation with the unsymmetrical χ′ tensor in some cases could also possibly be ascribed to fortuitous error cancellation. A recent experimental study27 compared the performance of the two versions of the KM equation for an S = 1/2 copper(II) ion in an axial coordination environment, where there is no zero-field splitting (ZFS) and it was possible to experimentally obtain the g tensor under the same conditions as the NMR data. The results were clearly in favor of the classical equation with the symmetrical χ tensor.27 This led to our reinvestigation of the derivation leading to the unsymmetrical χ′ tensor. Based on rigorous first-principles quantum mechanics and without invoking the SH approximation, we found that the KM equation can be directly derived from the Van den Heuvel–Soncini equation.31 In our derivation, we only assumed that the PDA holds and that the Hamiltonian is the effective nuclear charge SOC Hamiltonian, including gauge correction terms. Thus, it became clear that the neglect of contributions beyond the spin–dipolar interaction was not justified. Within the SH approximation, the KM equation is exactly equivalent to the approximation of the A tensor in terms of the g tensor,32
where γ is the gyromagnetic ratio. We found a proof for this long-distance limit as well and could verify it via QC calculations.31 The current work contains a detailed discussion of the proofs for Eqs. (3) and (7), the latter also in the context of the approximate Kohn–Sham density functional theory (DFT) and Hartree–Fock (HF) methods. Figure 1 gives an overview of the different ways to calculate the chemical shielding and the approximations involved. Furthermore, we go beyond the simple effective nuclear charge Hamiltonian of our previous work31 and generalize the results to arbitrary Hamiltonians that fulfill a certain gauge invariance condition, including the fully relativistic Dirac Hamiltonian.
Overview of different approaches for calculating PCSs and the approximations involved. Note that all the equations are given in atomic units. The application of the SH approximation (red solid lines) is straightforward, while the correct application of the PDA (blue dashed lines) has been unclear in the recent past. The top left arrow refers to the derivation of Eq. (3) involving Van Vleck’s equation for the susceptibility from the Van den Heuvel–Soncini equation. The bottom right arrow refers to the derivation of Eq. (7), the long-distance limit of the A tensor. Both derivations are discussed in depth in the present article.
Overview of different approaches for calculating PCSs and the approximations involved. Note that all the equations are given in atomic units. The application of the SH approximation (red solid lines) is straightforward, while the correct application of the PDA (blue dashed lines) has been unclear in the recent past. The top left arrow refers to the derivation of Eq. (3) involving Van Vleck’s equation for the susceptibility from the Van den Heuvel–Soncini equation. The bottom right arrow refers to the derivation of Eq. (7), the long-distance limit of the A tensor. Both derivations are discussed in depth in the present article.
II. DERIVATION OF THE LONG-DISTANCE LIMIT OF THE QUANTUM-CHEMICAL EXPRESSION FOR THE CHEMICAL SHIELDING
In this section, we present in detail our derivation of the KM equation from rigorous first-principles quantum mechanics, which we could only roughly sketch in our previous work31 due to space limitations. The most important non-standard symbols used throughout this article are given in Table I. Figure 2 illustrates the meaning of the different position vectors used in our derivations.
Symbols used throughout this article. Non-boldface versions of symbols representing vectors are used to denote their magnitudes. Furthermore, uppercase letters (such as in ) denote quantities summed over all electrons. The superscript is used for Cartesian indices.
| H | Electronic Hamiltonian (including relativistic effects such as SOC). |
| H(0) | Hamiltonian in the absence of magnetic fields. |
| Hnonrel | Nonrelativistic Hamiltonian [Eq. (8)]. |
| Hnonrel,(0) | Nonrelativistic Hamiltonian in the absence of magnetic fields [Eq. (9)]. |
| Eigenstates of H(0). | |
| Eigenvalues of H(0). | |
| Eigenstates of H(0),nonrel. | |
| Eigenvalues of H(0),nonrel. | |
| ξA(r) | Effective nuclear charge SOC factor [Eq. (11)]. |
| B | Homogeneous magnetic field. |
| MK | Magnetic dipole moment of nucleus K. |
| (∂H/∂μ)(0) | Derivative of the Hamiltonian with respect to some parameter μ, evaluated at B = 0, MK = 0. |
| O | Gauge origin (chosen at the “location” of the paramagnetic center). |
| K | Index of the nucleus whose chemical shielding we are interested in. |
| A | Index of another nucleus that, e.g., creates SOC. |
| R = RO − RK | Vector from the magnetic nucleus under consideration to the gauge origin. |
| i, j | Electron indices. |
| riL = ri − RL | Position of electron i relative to point L. |
| f(r) | The vector-valued function defined in Eq. (33). |
| Spin density matrix in the atomic orbital (AO) basis. | |
| z | Spatial part of the SOC operator [Eq. (88)]. |
| Perturbed orbital coefficients for spin channel σ ∈ {α, β}. | |
| (g|h) | Coulomb repulsion integrals [Eq. (91)]. |
| χ(r) | Scalar field that defines a gauge transformation. |
| (δH/δA(r))(0) | Functional derivative of the Hamiltonian with respect to the magnetic vector potential evaluated at A(r) = 0. |
| H | Electronic Hamiltonian (including relativistic effects such as SOC). |
| H(0) | Hamiltonian in the absence of magnetic fields. |
| Hnonrel | Nonrelativistic Hamiltonian [Eq. (8)]. |
| Hnonrel,(0) | Nonrelativistic Hamiltonian in the absence of magnetic fields [Eq. (9)]. |
| Eigenstates of H(0). | |
| Eigenvalues of H(0). | |
| Eigenstates of H(0),nonrel. | |
| Eigenvalues of H(0),nonrel. | |
| ξA(r) | Effective nuclear charge SOC factor [Eq. (11)]. |
| B | Homogeneous magnetic field. |
| MK | Magnetic dipole moment of nucleus K. |
| (∂H/∂μ)(0) | Derivative of the Hamiltonian with respect to some parameter μ, evaluated at B = 0, MK = 0. |
| O | Gauge origin (chosen at the “location” of the paramagnetic center). |
| K | Index of the nucleus whose chemical shielding we are interested in. |
| A | Index of another nucleus that, e.g., creates SOC. |
| R = RO − RK | Vector from the magnetic nucleus under consideration to the gauge origin. |
| i, j | Electron indices. |
| riL = ri − RL | Position of electron i relative to point L. |
| f(r) | The vector-valued function defined in Eq. (33). |
| Spin density matrix in the atomic orbital (AO) basis. | |
| z | Spatial part of the SOC operator [Eq. (88)]. |
| Perturbed orbital coefficients for spin channel σ ∈ {α, β}. | |
| (g|h) | Coulomb repulsion integrals [Eq. (91)]. |
| χ(r) | Scalar field that defines a gauge transformation. |
| (δH/δA(r))(0) | Functional derivative of the Hamiltonian with respect to the magnetic vector potential evaluated at A(r) = 0. |
Position vectors occurring in our derivations. Dashed lines denote position vectors with respect to some arbitrary coordinate origin 0, while full lines are position vectors relative to other positions in the molecule. Red lines are vectors describing the position of an electron with index i. The vectors are illustrated for the NiSAL–HDPT complex from our previous work.20 For this example, where the paramagnetism arises due to a single Ni ion, the gauge origin O is chosen at the position of this ion. A is another nucleus that, e.g., causes SOC (in this example, a N atom), and K is the nucleus whose chemical shielding we are considering (in this example, a H atom).
Position vectors occurring in our derivations. Dashed lines denote position vectors with respect to some arbitrary coordinate origin 0, while full lines are position vectors relative to other positions in the molecule. Red lines are vectors describing the position of an electron with index i. The vectors are illustrated for the NiSAL–HDPT complex from our previous work.20 For this example, where the paramagnetism arises due to a single Ni ion, the gauge origin O is chosen at the position of this ion. A is another nucleus that, e.g., causes SOC (in this example, a N atom), and K is the nucleus whose chemical shielding we are considering (in this example, a H atom).
A. Quantum-chemical Hamiltonian
Hartree atomic units (ℏ = me = e = 4πϵ0 = 1) are used in the following, where the Bohr magneton has the value μB = 1/2 and the vacuum permeability is given by μ0 = 4πα2, with α being the fine-structure constant. The nonrelativistic QC Hamiltonian in the presence of a magnetic field with vector potential A is given by
Here, Hnonrel,(0) is the nonrelativistic Hamiltonian in the absence of magnetic fields. In order to go beyond “spin-only” magnetism, it is necessary to extend this Hamiltonian and include spin–orbit coupling (SOC). For simplicity, we assume it in this section to be of the effective nuclear charge form,33
If this operator is added to the Hamiltonian, a so-called “gauge correction”34 must be added in the presence of magnetic fields, which guarantees gauge-invariant observables. The gauge correction for the effective nuclear charge SOC operator has the form
Note that this term also arises when making the substitution p → p + A(r) in Eq. (10), i.e., when using the proper expression for kinetic momentum instead of canonical momentum. If the vector potential describes the sum of an external homogeneous magnetic field B and the magnetic field created by a nucleus at position RK with magnetic moment MK,
the complete Hamiltonian can be written as
The last term in Eq. (14) (quadratic in MK) contributes to nuclear spin–spin coupling constants but is not relevant for the current work. Note that we use the notation rL = r − RL to denote a position relative to a point L at position RL. The point O occurring in Eq. (13) is the “gauge origin” of the vector potential for the homogeneous magnetic field.35
The field-free zeroth order Hamiltonian is given by
For the purpose of the present discussion, only SOC is added to the nonrelativistic zeroth order Hamiltonian. However, it is also possible to add the direct electronic spin–spin coupling operator without having to change any of the arguments that are presented in the following. The first derivatives of the Hamiltonian can be interpreted as the operators for the electronic magnetic moment, Mel = −∂H/∂B, and the magnetic hyperfine field created by the electrons at the position of the nucleus, BHF(RK) = −∂H/∂Mk. The derivative with respect to the nuclear magnetic moment can be written as the sum of four terms: The paramagnetic spin–orbit (PSO) term is defined as
the spin–dipolar (SD) term is defined as
the Fermi contact (FC) term is defined as
and the spin–orbit gauge correction term is defined as
Likewise, the derivative with respect to the external magnetic field can be written as the sum of three terms: an orbital part involving the total angular momentum operator with respect to the chosen gauge origin RO,
a spin part involving the total spin operator,
and a gauge correction term given by
Finally, the second derivatives are given by
B. Chemical shielding and susceptibility as second derivatives of the electronic free energy
Van den Heuvel and Soncini derived the following expression for the chemical shielding tensor in a molecule with arbitrary degeneracy:
Here, F = −kBT ln Z is the electronic Helmholtz free energy of the molecule, is the canonical partition function belonging to the zeroth order Hamiltonian, and β = 1/kBT. These expressions assume that the eigenfunctions and eigenvalues of the relativistic zeroth order Hamiltonian are known,
where ν is an additional label for states that are degenerate at zero field. For completeness, we present a simplified derivation of Eq. (25) in Appendix A. The basic assumptions in Eq. (25) are that the timescale of transitions between different electronic states is much shorter than that of the nuclear spin dynamics, meaning that the nuclei sense the average effect of all thermally populated electronic states, and that the shielding is independent of the field strength (weak-field approximation). These assumptions are fulfilled in most cases.
C. Long-distance limit of the operator derivatives
We now investigate the limit of Eq. (25) for the case that the magnetic nucleus is at a large distance from the paramagnetic part of the molecule that gives rise to the chemical shielding. McConnell performed a similar analysis24 based on Ramsey’s expression for nondegenerate ground states.2,3 Choosing the gauge origin used for the vector potential of the external magnetic field in the region of space where the unpaired electrons are located, the vector
is a good measure for the distance between the nucleus and the paramagnetic center. In principle, the gauge origin can be chosen arbitrarily and should not influence any observables because the effective nuclear charge Hamiltonian in Eq. (14) (including the gauge correction term) is well-behaved under gauge transformations. For example, McConnell chose the gauge origin for the external field at RK.24 However, this choice makes the derivation more complicated.
If we only consider the effect of electrons i located at the paramagnetic center, in the long-distance limit, riO can be considered small compared to R and then riK = riO + R ≈ R. Therefore, it is possible to expand the inverse powers of riK occurring in Eqs. (16), (17), (19), and (24) in a Taylor series around riO = 0 to obtain24,31
Inserting these expansions into the expressions for the Hamiltonian derivatives and retaining only the terms with inverse powers of R up to three, one obtains
and
with
Note that Eq. (32) cannot be obtained as a derivative of Eq. (31) because the latter is here evaluated at zero field. However, as we will demonstrate in Sec. IV B, Eq. (31) is also true at non-zero fields. The detailed derivation of Eqs. (31) and (32) is not trivial and is presented in Secs. II D and II E. Note that the FC contribution to the hyperfine field vanishes if there is no overlap between the nucleus and the electronic system that is responsible for the chemical shielding. It can be observed that Eqs. (31) and (32) contain terms scaling as R−2 due to the presence of the first term in f(r), in apparent violation of the PDA, for which the lowest nonvanishing inverse power is expected to be R−3. We will demonstrate below that these two R−2 operators give rise to terms that are of equal size but of opposite sign such that they exactly cancel each other.
D. Derivation of Eq. (31)
This section is quite technical and may be skipped without compromises in conceptual understanding. A basic identity used in this section is the commutator of the nonrelativistic field-free Hamiltonian with an arbitrary local one-electron operator,
which can also be written as
A local operator is one that is “diagonal” in the position representation, which is equivalent to saying that it is “multiplicative.” That is, when acting on a wavefunction in the position representation, the latter is multiplied by another function. A prominent example of a local operator is the potential energy. Another useful expression is the product of two Levi-Civita symbols in terms of a determinant of Kronecker deltas,
From this equation, it follows that a product of two cross products can be written as
Furthermore, an expression for the contraction of two Levi-Civita symbols over a common index can be obtained from Eq. (36),
It is also useful to note that the commutator of many-particle operators that are the sum of one-electron operators can be reduced to the commutator of the corresponding one-particle operators,
where only terms scaling up to R−3 are kept. In the derivatives of the Hamiltonian with respect to the magnetic field that occur in the van Vleck equation for the susceptibility tensor [Eq. (27)], there are no momentum operators except for the one that is part of the orbital angular momentum operator. It is therefore necessary to remove the remaining momentum operators from Eq. (40). By expressing the total momentum operator as
The third term in Eq. (40) can be written using the Levi-Civita symbol as
Antisymmetric combinations of terms such as occur in the orbital angular momentum operator. Hence, using Eq. (38),
Using Eq. (34), the corresponding symmetric combinations are obtained via
Here, we used Eq. (38), and the term involving N vanishes because of R ×R = 0. Equations (42) and (47) can now be inserted into Eq. (40) to obtain
with f(r) defined in Eq. (33). In order to achieve a cancellation of the R−2-scaling terms introduced by the f(r) function in the sum-over-states expression Eq. (25), there should be a commutator with the full effective nuclear charge Hamiltonian H(0) instead of Hnonrel,(0). We therefore anticipate that another contribution to the long-distance limit of will produce the term . In order to know what to look for, we now calculate this commutator. Since the one-electron spin operator commutes with any spatial operators, we leave it out of the expressions for now and only consider one-particle operators [see Eq. (39)]. We can write
Using this, the commutators with the two terms in fl(r) are, using Eq. (38),
The term involving two cross products can be rewritten using Eq. (37) as
When inserted into Eq. (51), this gives
We now turn our attention to the long-distance limit of the gauge correction Eq. (19) (again just the one-particle part with omission of the spin operator),
Here, we used Eq. (29) and only kept terms scaling up to R−3. When combining Eqs. (50), (53), and (54), the following expression is obtained:
Upon multiplication of this expression with and summation over the Cartesian components p and electrons i [using Eq. (39)], we obtain
We see that the long-distance limit for the gauge correction indeed produces the desired commutator with the SOC operator.
Using Eq. (30), the long-distance limit of the spin–dipolar part (retaining terms scaling up to R−3) can straightforwardly be obtained as
E. Treatment of the commutator terms and derivation of Eq. (32)
Given two arbitrary Hermitian operators A and B, it is straightforward to show that (see Appendix B)
which is valid for arbitrary n. A similar equation, valid only for nondegenerate ground states, can be formulated in terms of the polarization propagator and was, for example, used by Geertsen in his proof of the gauge invariance of the second order energy in the presence of magnetic fields.38
We now consider the case and , for which we need to evaluate the commutator on the right-hand side of Eq. (58). Since local operators commute with other local operators and with spin operators, only the orbital part of , which is proportional to the total angular momentum, contributes to this commutator. Using
together with Eq. (38), one obtains
for the first term and
for the second term of the commutator with [defined in Eq. (33)]. Using Eq. (37), the term involving two cross products can be rewritten as
Inserting this into Eq. (61) leads to
Furthermore,
The long-distance expansion of the diamagnetic operator operator Eq. (24) [using Eq. (29) and retaining only terms scaling up to R−3] is given by
F. Long-distance limit of the chemical shielding
In the first term of Eq. (25), the term arising from the commutator in Eq. (31) vanishes since for an arbitrary operator A,
Equation (58) shows that the contributions to the second and third term of Eq. (25) arising from the commutator terms in Eqs. (31) and (32) exactly cancel each other. Together, these results mean that inserting Eqs. (31) and (32) into Eq. (25) for the chemical shielding gives
with defined in Eq. (27). Equation (67) is exactly the KM equation [Eq. (3)] written component-wise. This shows that an exact QC treatment confirms the KM long-distance expression and indicates that the different results predicted by the use of the unsymmetrical χ′ tensor are an artifact that derives from the neglect of the orbital contribution to hyperfine coupling.
III. THE LONG-DISTANCE LIMIT IN THE SPIN HAMILTONIAN FRAMEWORK
It is often a good approximation to assume that all electronic states that make sizable contributions to the paramagnetic shift can be described by a SH of the form
Note that we employ the “IAS” convention for the magnetic hyperfine term in the SH. The A tensor in this convention is related to the one in the also often used “SAI” convention by transposition. In our previous work,31 we argued that within the SH approximation, the validity of the KM equation [Eq. (3)] implies that in the long-distance limit, the complete tensor A can be expressed in terms of the tensor g,31,32
which is Eq. (7) written in atomic units. This expression was previously only known to hold for the SD contribution to the A tensor in terms of the spin contribution to the g tensor.26 In this section, we present a detailed derivation of Eq. (69) at the level of degenerate perturbation theory up to second order (DPT2). It will turn out that it is crucial to use gauge-invariant expressions for the SH parameters, i.e., to properly include the gauge correction terms [Eqs. (19) and (22)].
A. Exact eigenstates
We now consider eigenstates of the nonrelativistic field-free Hamiltonian, which are also eigenfunctions of the total spin operator and its z component,
Using DPT2, the g and A tensors can be written in terms of the nonrelativistic eigenfunctions and eigenvalues as34,39–41
and
Equations (74) and (79) involve a sum over all excited states. Hence, we will call those terms the sum-over-states (SOS) contributions to the g tensor and the A tensor, respectively. The gauge correction to A was only considered in a few ab initio studies so far.42,43 In the latter article, it was called “diamagnetic orbital contribution.” The fact that it is usually neglected in QC calculations indicates that its contribution is often small in standard situations. However, without this contribution, the A tensor is not gauge invariant.
The spatial integrals that are necessary to evaluate Eq. (78) are very similar to the integrals needed for calculating the contribution of SOC to NMR indirect nuclear spin–spin couplings.42 In practice, there is the problem that the integrals diverge for K = A. One, therefore, usually follows the pragmatic approach to leave out the SOC contribution of the nucleus for which the property is calculated.42,44 We expect this to be an excellent approximation for the numerical results in our previous31 and present work since the probe nucleus is far away from the electronic system of interest. Hence, SOC due to this nucleus should be negligible. Note that the necessity for such an ad hoc approach is an artifact of the approximate effective nuclear charge Hamiltonian employed in this section. For the Dirac Hamiltonian discussed below, no such divergence issues arise.
The Fermi contact contribution to A is easily seen to vanish in the long-distance limit, where riK ≠ 0. However, note that there are examples of extended electron spin delocalization, e.g., across π systems. In these cases, our basic assumption of a spatially confined paramagnetic center is not strictly applicable. For the spin–dipolar contribution, it was already known before that26
which can also immediately be obtained by using Eq. (30). This leaves us with the investigation of the long-distance behavior of the SOS contribution and the contribution due to the gauge correction. For an arbitrary Hermitian operator A and an arbitrary Hermitian singlet operator B, i.e., one that does not mix different total spins and magnetic sublevels of nonrelativistic zeroth order states, there is a relation similar to Eq. (58),
This equation is derived in Appendix B. Inserting Eq. (48) into Eq. (79) and applying Eq. (81) lead to
The two “nonphysical” commutator terms scaling like R−2 cancel when adding Eqs. (82) and (83), which concludes the proof of Eq. (69) at the level of DPT2. Note that only the gauge-invariant sum of the gauge correction term and the SOS term fulfills Eq. (69). In contrast to that, the SOS term alone, Eq. (82), has an unphysical R−2 long-distance contribution.
B. Approximate methods: Kohn–Sham density functional theory and Hartree–Fock
In practice, one has no access to exact eigenstates of the Hamiltonian and has to resort to approximate methods instead. In our previous work,31 we observed that for unrestricted Kohn–Sham (UKS) density functional theory (DFT) with the BP86 exchange–correlation functional, Eq. (69) was only fulfilled for very large basis sets. In order to understand this behavior, we now discuss the PDA for UKS and unrestricted HF (UHF). In this section, we use Greek letters μν… for atomic orbitals (AOs), ij… for occupied molecular orbitals (MOs), ab… for virtual/unoccupied MOs, and pq… for general MOs. In UKS and UHF, the expressions for the SOC45,46 and gauge contributions to the g tensor and the A tensor corresponding to Eqs. (73), (74), (78), and (79) are given by
with the spatial part of the SOC operator being, for convenience, defined as
is the spin density matrix, given in terms of orbital coefficients by
The perturbed spin density matrices in Eqs. (85) and (87) are defined with respect to the orbital Zeeman and the PSO part of the Hamiltonian, respectively. Note that we assume the use of a standard (field-independent) AO basis set. In this case, the g tensor depends on the chosen gauge origin. For a single local paramagnetic center, it has been demonstrated that choosing a gauge origin close to the paramagnetic center usually yields good results.47 The gauge origin dependence can be removed by employing gauge-including AOs (GIAOs).48,49 However, this makes the theory significantly more complicated at not much added conceptual insight. For simplicity, we therefore employ field-independent basis sets with the center of electronic charge as the gauge origin in our calculations. The assumption of a field-independent basis set is made throughout the derivations in this section.
The first-order perturbed orbital coefficients can be obtained by solving the coupled-perturbed self-consistent field (CP-SCF) equations, which for purely imaginary perturbations (using the “magnetic Hessian”45) read
where cHF is the amount of HF exchange. Note that the perturbed orbital coefficients are linear in the perturbation operator such that one can consider different additive contributions to the perturbation operator separately. Here and in the following, we use the notation
for Coulomb repulsion integrals. For pure functionals (cHF = 0), the magnetic Hessian is diagonal45 and the CP-SCF equations are completely decoupled. This means that one can directly write down the solutions to Eq. (90),
The field-free Kohn–Sham (or Fock) operator hKS,σ, which in the chosen one-electron AO basis has the orbitals |pσ⟩ and orbital energies as eigenfunctions and eigenvalues, is given by45
The terms are (from left to right) the kinetic energy, the external (nuclear) potential, the Coulomb operators, the exchange operators, and the exchange–correlation potential. For pure functionals, cHF = 0, while for HF, cDF = 0 and cHF = 1. For hybrid DFT functionals, both cHF and cDF are non-zero. The only operators in Eq. (93) that are non-local (i.e., that can contribute to the commutator with a local operator) are the kinetic energy and the exchange operators. This means that for arbitrary local operators f, i.e., operators that multiply a wavefunction with the function f(r), one has in analogy to Eq. (35),
Therefore, in the long-distance limit, one can write the perturbation operator relevant for the A tensor as
which is analogous to Eq. (48). The commutator term in Eq. (95) can be dealt with in a similar way as in Eq. (81) for exact states.
From Eq. (92), one immediately sees that, for pure functionals, such a commutator term gives a contribution to the perturbed orbital coefficients of the form
We will now show that the same expression also holds for HF and hybrid functionals. One can write a one-electron matrix element of the commutator with an exchange operator as
Assuming that it is possible to write
(for which completeness of the one-particle basis is a sufficient but not necessary condition), it follows that
In the second step, pσ were restricted to virtual orbitals since there is an exact cancellation of terms for occupied orbitals. A matrix element of the commutator with the Kohn–Sham operator is given by
Comparing this with the CP-CSF equations, Eq. (90), shows that the contribution of a commutator term to the perturbation operator (∂h/∂μ)(0) leads to a contribution to the first-order orbital coefficients of the form
which is the same result we also obtained for pure functionals [Eq. (96)]. If Eq. (98) is fulfilled, the corresponding contribution to the contraction of the perturbed orbital coefficients with the SOC operator can be written as
In this step, Eq. (98) is necessary not only in the presence of HF exchange but also for pure functionals. In our specific case, with the local operator being defined through the function f(r) given in Eq. (33), the commutator term in Eq. (95) leads to a contribution
Inserting this into Eq. (87) gives
In comparison, inserting Eq. (55) into Eq. (86) for the gauge contribution to the A tensor at the DFT and HF level gives
Adding Eqs. (106) and (107) leads again to a cancellation of the unphysical terms. This shows that Eq. (69) is also fulfilled at the level of DFT and HF, but only if the condition Eq. (98) is fulfilled. This condition will, in general, not be true for commonly used finite basis sets.
To summarize our findings, the derivation of the physically correct long-distance behavior of the A tensor at the level of UKS or UHF requires the validity of Eq. (98). For pure functionals, where Eq. (96) is always true, this condition enters only once in the derivation [for Eq. (103)]. In contrast, the presence of HF exchange also requires the condition in the derivation of Eq. (102). Therefore, one could expect that problems associated with the incompleteness of the one-electron basis will be more pronounced in the presence of HF exchange. In order to test this hypothesis, we present results for a bare proton at different distances from a CO+ radical at the UHF level. The computational details were chosen identical to the ones in our previous paper31 except that we used extremescf convergence and lowered the CP-SCF convergence threshold to 10−9 in ORCA.50,51 The results are shown in Fig. 3, where we plot the “relative difference”31 between the exact and approximate APSO/SOC + Agauge,
Here, ‖ ⋅‖ is the Frobenius norm and A(g) denotes the long-distance functional relationship between the A tensor and the g tensor given in Eq. (69). One can observe that the basis set convergence is not slower than for the pure BP86 functional that we used in our previous study.31 Actually, the results of both methods are extremely similar. This demonstrates that, other than expected, the additional approximation that is present when using HF exchange together with incomplete basis sets apparently does not lead to an additional deviation from the correct behavior. In both cases (BP8631 and UHF), one can observe only for the largest decontracted cc-pV6Z basis set the physically expected behavior that the error associated with the PDA tends toward zero for sufficiently large distances. We note that the slower basis set convergence at larger distances might be related to the fact that the commutator terms (which must cancel in order for the PDA to hold) scale as R−2 and therefore increase in importance relative to the “physical” contributions. Therefore, the cancellation must be more accurate at larger distances in order to obtain the same accuracy for the total A tensor, explaining the heavier basis set demands.
Relative difference between the exact sum of PSO/SOC and gauge contributions to the A tensor and a long-distance approximation using the sum of orbital and gauge contributions to the g tensor.
Relative difference between the exact sum of PSO/SOC and gauge contributions to the A tensor and a long-distance approximation using the sum of orbital and gauge contributions to the g tensor.
IV. GENERALIZATION TO OTHER HAMILTONIANS
A. The Dirac Hamiltonian
In this section, we will generalize the results obtained with the approximate effective nuclear charge Hamiltonian (Sec. II) to the fully relativistic Dirac Hamiltonian. We will restrict the discussion to the one-electron Dirac Hamiltonian. The generalization to many electrons (Dirac–Coulomb Hamiltonian) is straightforward. The one-electron Dirac Hamiltonian in the presence of magnetic fields has the form
Note that the speed of light occurring in this equation is equal to the inverse fine-structure constant in the atomic units used throughout this work, c = 1/α. However, we will explicitly leave it in the expressions as is common in the literature. Also note that in this section, σ denotes the vector of Pauli matrices and not the chemical shielding tensor.
In the absence of fields, the Dirac Hamiltonian takes the form
The derivatives of the Dirac Hamiltonian, using the vector potential of Eq. (13), are given by
All second derivatives vanish, i.e., the second term of the Van den Heuvel–Soncini equation [Eq. (25)] is zero in contrast to the nonrelativistic and effective nuclear charge Hamiltonians. We will now investigate whether the long-distance limit of Eq. (112) fulfills the same equation as the corresponding derivative of the effective nuclear charge Hamiltonian, Eq. (31). Inserting Eq. (29) into the off-diagonal term of Eq. (112), one obtains
In the following, we will make use of the expression for the commutator of the field-free Dirac Hamiltonian with a local operator,
Using this equation, we obtain
Equations (115) and (116) are the two contributions to the commutator between hD,(0) and the local operator fl defined in Eq. (33). Furthermore,
[a special case of the equation
which is valid for arbitrary vectors a,b,c,d], after comparison with Eqs. (112) and (113), we obtain
which is identical to the equation one obtains in the case of the effective nuclear charge approximation [Eq. (31)]. Since the commutator of the derivative defined in Eq. (111) with the local operator fl [occurring on the right-hand side of Eq. (58)] is zero, the commutator term in Eq. (120) does not contribute to the chemical shielding tensor. This concludes the proof of the long-distance limit of the chemical shielding tensor using the Dirac Hamiltonian, which again confirms the KM equation.
B. Arbitrary Hamiltonians that respect gauge invariance
The fact that the same long-distance expression for the derivative of the Hamiltonian with respect to the magnetic moment is obtained for both the effective nuclear charge Hamiltonian [Eq. (31)] and the Dirac Hamiltonian [Eq. (120)] suggests that there must be a general principle that is independent of the choice of a specific Hamiltonian. In particular, the same function fl(r) is present in the commutators. Therefore, the equation for this function [Eq. (33)], which is independent of the chosen Hamiltonian, needs to be explained. The observed cancellation of commutator terms for the chemical shielding tensor is reminiscent of the cancellation of terms when applying a gauge transformation to a gauge-invariant property.34 Therefore, in the following, we will focus on Hamiltonians that respect gauge invariance.
Under changes of the gauge of the vector potential, A(r) → A′(r) = A(r) + ∇χ(r), we assume that the Hamiltonian will transform according to52
where χ is a local operator that multiplies the many-electron wavefunction with . Equation (121) is fulfilled for the nonrelativistic, effective nuclear charge, and Dirac Hamiltonians considered in this article. If HΨ = EΨ, it follows that H′e−iχΨ = Ee−iχΨ, i.e., eigenenergies stay unchanged, while the eigenfunctions acquire a (position-dependent) phase factor under the gauge transformation. We note that it is possible that the Hamiltonian and wavefunction are transformed under the change of gauge via a more complicated unitary transformation than that assumed in Eq. (121).53 It should be straightforward to extend the following argument to such cases as well. Using the Baker–Campbell–Hausdorff (BCH) expansion, Eq. (121) can be written as
Alternatively, one can consider the Hamiltonian as a functional of the magnetic vector potential A and write the gauge-transformed Hamiltonian as a Taylor series around the reference vector potential (see Appendix C),
Comparing Eqs. (122) and (123) and equating terms that involve equal powers of the gauge function χ, we obtain
and similar for higher orders. Equation (124) means that a contraction of δH/δA(r) with the gradient of some scalar field can be written in terms of a commutator with the Hamiltonian. This is a quite remarkable result. When considering the reference point A = 0, it means that the form of (δH/δA(r))(0), which describes how the terms of the Hamiltonian linear in A look like, is constrained by the field-free Hamiltonian H(0). The first functional derivative of the Hamiltonian appearing in Eq. (124) can be identified with the particle number current density operator;54 see also Appendix D. As an illustration, the first and second functional derivatives of the Hamiltonians considered in this work are presented in Table II. It is illustrative to use the functional derivative from Table II to show that Eq. (35) that we used earlier is a special case of Eq. (124).
One-electron integrals over the first and second functional derivatives with respect to the vector potential for the Hamiltonians considered in this work. For the nonrelativistic and effective nuclear charge Hamiltonian, the one-electron functions are two-component spinors (spin orbitals), whereas for the Dirac Hamiltonian, they are four-component spinors.
| Hamiltonian | ||
| Hnonrel | ||
| HD | 0 |
| Hamiltonian | ||
| Hnonrel | ||
| HD | 0 |
We will now consider the derivatives of the Hamiltonian with respect to the external magnetic field and the nuclear magnetic moment. Without specifying the Hamiltonian, we cannot write down these derivatives explicitly. However, we do know the dependence of the vector potential on these parameters [Eq. (13)]. Using the chain rule, one can write
Note that we consider the derivatives in Eqs. (126) and (127) at finite field as opposed to before. In order to obtain the long-distance limit, it is convenient to transform the vector potential of the nuclear magnetic moment to the multipolar gauge.55–58 This gauge expresses the magnetic vector potential in terms of the magnetic field. Expanding the fields in Taylor series around the origin RO leads to the vector potential in the multipolar gauge given by
where the gauge function is defined as58
More explicitly, the vector potential in the multipolar gauge is given by58
i.e., it is defined in terms of the magnetic field and its spatial derivatives at the chosen origin.
Equation (128) can be inverted to obtain
If the series expansions in Eq. (129) and Eq. (130) are truncated, approximations to Anuc(r) are obtained from Eq. (131) that are most accurate in the vicinity of the origin RO and become increasingly worse with increasing distance from RO. When doing this, the gauge function should be truncated, for consistency, at an order that is higher by one than that at which the vector potential in the multipolar gauge is truncated. This is because the gradient operator in Eq. (131) reduces the order in rO by one.
The exact dipolar magnetic field created by the nucleus is given by59
The fields Anuc and Bnuc and the gradient of Anuc at the origin RO are given by
Inserting these expressions into Eqs. (129) and (130) and using Eq. (131), one obtains for the vector potential up to order ,
where f(r) is the same function we have already obtained earlier in a different way [Eq. (33)]. We now recognize that this function is a long-distance approximation to the derivative of the gauge function that transforms between the Coulomb and the multipolar gauge [Eq. (129)],
Of course, Eq. (136) could also have been obtained by a direct calculation of the long-distance limit of the nuclear vector potential [using Eq. (29)], but the use of the multipolar gauge [Eqs. (129)–(131)] is more straightforward.
From Eq. (136), one can calculate the derivative of the vector potential with respect to the nuclear magnetic moment,
Taking the derivative of this expression with respect to a component of the homogeneous magnetic field gives
Evaluation of Eqs. (140) and (141) at zero field (B = 0, MK = 0) immediately yields Eqs. (31) and (32). The remainder of the derivation of the long-distance limit of the chemical shielding (cancellation of commutator terms, etc.) proceeds exactly as before, once again confirming the KM equation, Eq. (3). The preceding derivation demonstrates that the long-distance KM equation is valid for arbitrary Hamiltonians that fulfill the gauge invariance property of Eq. (121), which encompasses the cases of the effective nuclear charge and Dirac Hamiltonians that we treated explicitly in Secs. II and IV A.
V. DISCUSSION AND CONCLUSIONS
In the present work, we discussed in detail the rigorous quantum-mechanical proofs of the KM equation, Eq. (3), and the (in the SH approximation) equivalent long-distance expression for the A tensor, Eq. (7).31 We also extended the proof of Eq. (7) to the case of the approximate DFT and HF methods, which showed that the equation is only valid in a sufficiently complete basis. This explains the computational results that we reported previously, where extremely large Gaussian basis sets were needed for the sufficient cancellation of the unphysical terms that scale like R−2 with the distance between the magnetic nucleus and the paramagnetic center.31
The major drawback of our previous QC proof of the KM equation was the reliance on an approximate effective nuclear charge SOC Hamiltonian.31 However, we have shown here that a proof for the fully relativistic Dirac Hamiltonian is possible along very similar lines. This led to the question whether it is possible to derive the KM equation without reference to any specific Hamiltonian. This is, indeed, the case, and we were able to show that it can be derived solely based on the assumption of gauge invariance. This is the most generally valid proof of the expression that exists to date. Our new derivation also provides an interpretation of the function f(r) that arose naturally in the derivation for the effective nuclear charge Hamiltonian: it defines the gauge transformation of the vector potential of the nuclear magnetic moment from the Coulomb to the multipolar gauge. The latter is parameterized by the magnetic field and its derivatives at a chosen expansion point. This leads to the following alternative viewpoint of our results: if the Hamiltonian transforms properly under gauge transformations, it does not matter which gauge is chosen. Therefore, one could use the multipolar gauge from the very beginning instead of Eq. (13) to define the Hamiltonian. Then, the unphysical R−2 terms would never arise and no cancellation between commutator terms would be necessary. Therefore, the KM equation arises very naturally from the use of the multipolar gauge. This choice of gauge also has the advantage that the long-distance expansion is already “built in” in the form of the series expansion [Eq. (130)]. The PDA can simply be derived by truncating this series expansion after the first term.
Our general derivation based on gauge invariance also sheds new light on the DFT and HF results. The commutator terms arise because of a gauge transformation from the multipolar to the Coulomb gauge. The fact that these terms do not exactly cancel is then a simple consequence of the fact that approximate methods using finite basis sets are not gauge invariant. For other magnetic properties, such as magnetizabilities and chemical shieldings of closed-shell molecules, this problem has been effectively solved by introducing GIAOs as basis functions.48,49 However, in the way GIAOs are commonly implemented, they can alleviate the gauge problem only for changes in the gauge origin for a homogeneous magnetic field. In our case, the problems arise because of more general gauge transformations in a non-homogeneous magnetic field due to the magnetic nucleus. Therefore, GIAOs will not be able solve the slow basis set convergence observed in Sec. III B and our previous work.31
The developed theoretical framework based on the multipolar gauge also seems promising for extending the KM theory beyond the PDA. The physical picture behind Eq. (3) is that the external magnetic field induces currents in the paramagnetic center that create a pure dipole field, which is then probed by the magnetic nuclei in the molecule. By truncating the expression for the vector potential in the multipolar gauge at higher order, one could include the effect of the induced quadrupole moment and even higher magnetic moments. This would lead to a more accurate description of the induced magnetic field and, therefore, to more accurate chemical shifts. We assume that the inclusion of higher induced magnetic multipole moments could be especially beneficial for systems where the paramagnetic center cannot easily be associated with a single atom, e.g., in oligonuclear metal complexes. In such molecules, there is no “natural” location in space for the induced dipole moment. However, even in mononuclear metal complexes, contributions beyond the PDA can become significant, in particular for nuclei close to the metal center. In contrast to the direct computational evaluation of Eq. (25),9 a formulation in terms of induced magnetic multipole moments could yield accurate results while retaining the philosophy of the KM equation, meaning that the chemical shieldings of different nuclei are not independent, but are probes for the same quantity: the magnetic field due to currents induced in the paramagnetic center. Such an approach also has the advantage that quantities such as the susceptibility can be calculated using truncated molecular models that possibly do not even include the nuclei for which the PCSs are required.25 This could make it possible to target systems such as proteins that would otherwise be intractable. Furthermore, the direct evaluation of Eq. (25) with the nuclear vector potential in the usual gauge [Eq. (13)]9 has potentially the same demanding basis set requirements for large distances as we observed for the A tensor. In contrast, the susceptibility tensor usually converges quickly when increasing the basis set size.
We conclude with a final note on the topic of gauge invariance. Hamiltonians that are not exactly gauge invariant are widespread. For example, the gauge-invariant expression for the A tensor in the framework of the relativistic second order Douglas–Kroll–Hess theory (DKH2) following from an fπFW transformation turns out to be divergent.43 Therefore, only the expression resulting from the so-called fpFW transformation can be practically used. However, this expression is not gauge invariant. For the spin–orbit mean field (SOMF) approximation to the SOC operator that is widely used in ORCA,60 there has also not been developed a gauge correction to date. However, we expect that it should be possible to formulate a gauge correction for the SOMF operator from the full Breit–Pauli Hamiltonian. One can also probably expect approximate gauge invariance when combining the SOMF operator with the effective nuclear charge gauge correction. The derivations presented in this paper demonstrate that the use of expressions that are not gauge invariant can in the worst case lead to completely unphysical results, in our case the presence of terms that scale as R−2 with the distance between the magnetic nucleus and the paramagnetic center. Therefore, when a Hamiltonian is used that violates gauge invariance, the results should be thoroughly checked.
SUPPLEMENTARY MATERIAL
See the supplementary material for raw data for all calculated contributions to the SH parameters.
ACKNOWLEDGMENTS
L.L. thanks Willem Van den Heuvel for invaluable feedback on the manuscript. The authors acknowledge financial support from the Max Planck Society. E.R., G.P., and C.L. acknowledge financial support from the Fondazione Cassa di Risparmio di Firenze, the Ministero della Salute, through Grant No. GR-2016-02361586, the Italian Ministero dell’Istruzione, dell’Università e della Ricerca through “Progetto Dipartimenti di Eccellenza 2018-2022” to the Department of Chemistry “Ugo Schiff” of the University of Florence, the PRIN 2017A2KEPL project “Rationally designed nanogels embedding paramagnetic ions as MRI probes,” and the European Commission through H2020 FET-Open project HIRES-MULTIDYN (Grant Agreement No. 899683).
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
DATA AVAILABILITY
The data that support the findings of this study are available within the article and its supplementary material.
APPENDIX A: DERIVATION OF THE CHEMICAL SHIELDING TENSOR
The basic assumption often used to arrive at practical equations for the chemical shielding tensor is that transitions between electronic states are much faster than the nuclear spin dynamics. Then, the nuclei sense the average effect of the ensemble of electronic states in thermodynamic equilibrium. The chemical shielding tensor σ in the weak-field approximation is defined via the following relation between the average hyperfine magnetic field (i.e., the average field generated by the orbital and spin magnetic moments of the electrons at the position of the nucleus) and the applied magnetic field:
This can be rewritten as
This equation was already used in the past as a starting point to calculate the chemical shielding tensor.6 The Boltzmann average of the magnetic hyperfine field can be written as
where we used the operator representing the hyperfine field given by . Using Wilcox’s equation for the derivative of an exponential operator,61
one can write
Furthermore, since the zeroth order states are time-even and is time-odd, ,7 which leads to
Using the product rule to calculate the derivative and employing Eqs. (A5) and (A6), Eq. (A2) can be written as
This equation is identical to Eq. (8) of Van den Heuvel and Soncini.7 The remainder of the derivation proceeds exactly as in their case. The traces in Eq. (A7) are evaluated in the eigenbasis of H(0) together with the insertion of a resolution of the identity in the second term. This gives
The integral over ω can now be evaluated as
After swapping the summation indices in one of the terms, one immediately obtains Eq. (25). The connection of our derivation to Van den Heuvel and Soncini’s definition of the chemical shielding tensor as a second derivative of the free energy becomes clear by recognizing that the average hyperfine field as defined in Eq. (A3) can also be written as
an equation that can be straightforwardly derived using Eq. (A4).
APPENDIX B: DERIVATION OF EQS. (58) AND (81)
While Eqs. (58) and (81) look very similar at first sight, they differ in important details, which makes separate derivations necessary: in Eq. (58), there is an additional summation over the degenerate manifold belonging to energy level n, which is not present in Eq. (81). In the latter, the sum is only over those excited states having the same spin as the ground state and only involves the principal component (M = S). Furthermore, the operator B is restricted to be a singlet operator, which is not necessary for the derivation of Eq. (58).
We start with the derivation of Eq. (58). We first let H(0) in the numerator act on the eigenfunctions and cancel the resulting energy differences with the denominator to obtain
Using completeness of the set of eigenfunctions, one can write
After inserting this into Eq. (B1), it can be seen that one term in the first part of the equation cancels another term in the complex conjugate, leading to Eq. (58).
Equation (81) can be derived as follows: By acting with the zeroth order Hamiltonian in the commutator on the zeroth order eigenfunctions, one obtains
The sum over excited-state projectors can be written as
where PSS is the orthogonal projector on the subspace of states with total spin S and magnetic sublevel M = S. Since B is a singlet operator by assumption, it does not change the spin quantum numbers of the bra state on which it acts from the right, i.e.,
This means that
APPENDIX C: FUNCTIONAL DERIVATIVES AND TAYLOR SERIES EXPANSION OF THE HAMILTONIAN
One can write the Hamiltonian in terms of its matrix elements in a complete basis {|I⟩},
The matrix elements are simple numbers, i.e., their functional dependence on the vector potential can be written in terms of a usual (functional) Taylor expansion around some reference potential,
This expansion allows us to define the functional derivatives of the Hamiltonian itself,
etc., which leads to a Taylor series expansion for the Hamiltonian
APPENDIX D: IDENTIFICATION OF THE FUNCTIONAL DERIVATIVE OF THE HAMILTONIAN WITH THE PARTICLE NUMBER CURRENT DENSITY OPERATOR
In the following, we generalize an argument given by Epstein for the nonrelativistic Hamiltonian.62 We assume that the expectation value of a local operator f can be expressed as an integral over the electron density,
For the Hamiltonians considered in this work (nonrelativistic, effective nuclear charge, Dirac), this implies that the electron density for a single electron is given by ρ(r) = ψ†(r)ψ(r).
If both the wavefunction and the operator are time-dependent, the time derivative of the expectation value ⟨f⟩(t) = ⟨Ψ(t)|f(t)|Ψ(t)⟩ is obtained from Eq. (D1) as
One can also express the time derivative of the expectation value via the generalized Ehrenfest theorem,
Here, we used Eq. (124), which is valid if the Hamiltonian respects gauge invariance, and performed integration by parts. A comparison of Eqs. (D2) and (D3) shows that
which is true for arbitrary functions f(r, t). This means that
which is simply the continuity equation if we identify
as the particle number current density. Hence, we can identify the first functional derivative of the Hamiltonian with respect to the magnetic vector potential with the particle number current density operator.
