Quantum Chemical Calculation of Molecules in Magnetic Field

For most practical applications, NMR experiments are the situations where molecules are subjected to low field strengths. A detailed review on ab initio calculation of shielding constants and chemical shifts is presented in (13). We only summarize the conceptual framework, which is crucial to the quantum chemical calculations performed by commercially available softwares.
We start by writing the NMR Spin Hamiltonian

$$\hat{H}=-\sum_{K}\gamma_{K}\hbar\vec{B}^{T}(1-\hat{\sigma_{K}})\hat{I}_{K}+\frac{1}{2}\sum_{K\neq L}\gamma_{K}\gamma_{L}\hat{I}^{T}\hbar^{2}(\hat{D}_{KL}+\hat{K}_{KL})\cdot\hat{I}_{L}.$$

(1)

Here, $\gamma_{K}$ is the gyromagnetic ratio, $\hat{I}_{K}$ is the nuclear spin operator, $\hat{\sigma_{K}}$ is the magnetic shielding tensor (due to electron density), $\hat{D}_{KL}$ and $\hat{K}_{KL}$ are dipolar and indirect coupling constants. The notation of (13) has been preserved for the reader’s convenience. The isotropic Hamiltonian can be written as the rotational average of Eq.1,

$$\hat{H}=-\sum_{K}\gamma_{K}\hbar(1-\sigma_{K})B\hat{I}_{K_{z}}+\frac{1}{2}\sum_{K\neq L}\gamma_{K}\gamma_{L}K_{KL}\hbar^{2}\hat{I}_{K}\cdot\hat{I}_{L}.$$

(2)

where $\sigma_{K}=\frac{1}{3}tr(\hat{\sigma_{K}}$) and $K_{KL}=\frac{1}{3}tr(\hat{K}_{KL})$. The second term in the Hamiltonian can be compactly written in terms of the spin-spin coupling tensor J${}_{\textbf{{KL}}}$, which is defined as

$$\textbf{J}_{\textbf{KL}}=\frac{\hbar\gamma_{K}\gamma_{L}\textbf{K}_{\textbf{KL}}}{2\pi}.$$

(3)

We then note that the shielding and indirect spin-spin coupling is four orders weaker than Zeeman and Dipolar Interactions and hence could be treated as first order perturbation. To find the NMR parameters from the Spin Hamiltonian, the energy (E($\vec{B},\vec{M}$)) is expanded in terms of magnetic field and magnetic dipole moment $\vec{M}$. This gives us the following relations, which are solved by variational perturbation theory.

$$\bm{\sigma}_{\textbf{K}}=\vec{B}^{T}\frac{d^{2}E(\vec{B},\vec{M})}{d\vec{B}d\vec{M}_{K}}\vec{M}_{K},\hskip 14.22636pt\hat{K}_{KL}+\hat{D}_{KL}=\hat{M}_{L}^{T}\frac{d^{2}E(\vec{B},\vec{M})}{d\hat{M}_{L}d\hat{M}_{K}}\hat{M}_{K}$$

(4)

It is a common practice to write the magnetic field as a curl of a vector potential, or a gauge potential. For a net magnetic field, including the contribution from the individual dipole moments and the external field, one can associate a total vector potential A${}_{\textbf{tot}}$,

$$\textbf{A}_{\textbf{tot}}=\textbf{A}_{\textbf{O}}(\vec{r}_{i})+\sum_{k}\textbf{A}_{k}(\vec{r}_{i})\,\hskip 14.22636pt\vec{B}_{tot}=\nabla\times\textbf{A}_{\textbf{tot}}$$

(5)

where the second term in summation is the contribution from individual magnetic dipole moment, and $\textbf{A}_{\textbf{O}}(\vec{r}_{i})$ is the vector potential associated with the external magnetic field defined as

$$\textbf{A}_{\textbf{O}}(\vec{r})=\frac{1}{2}\vec{B}\times\vec{r}_{iO}$$

(6)

The subscript of A denotes gauge origin and the notation $\vec{r}_{iO}=\vec{r}_{i}-\vec{r}_{O}$ is adopted. Although the NMR spectra depends largely on the spin Hamiltonian in Eq.1, electronic effects cannot be ignored. For an electron, it interacts with the vector potential due to the external magnetic field as well as the intrinsic field due to the nuclear magnetic moment. Thus, in general, the canonical momentum operator for an electron in magnetic field could be written as

$$\bm{\pi}=\textbf{p}+q\textbf{A}$$

(7)

Hence, in the Hamiltonian, terms linear in A and quadratic in A are obtained. By inserting the definition of A, we get

$$\textbf{A}\cdot\textbf{p}=\frac{1}{2}(\vec{B}\times\vec{L}_{O}),\hskip 11.38092pt\frac{1}{2}\textbf{A}^{2}=(\frac{1}{2}\vec{B}_{ext}\times\vec{r}_{iO})^{2}=\vec{B}\textbf{P}^{\zeta}\vec{B}$$

(8)

here, $\textbf{P}^{\zeta}=\frac{1}{8}(\vec{r}^{T}_{iO}\vec{r}_{iO}\textbf{1}-\vec{r}_{iO}\vec{r}^{T}_{iO}$). The symbol $\zeta$ represents the type of nuclear-nuclear or nuclear-electronic interaction. The magnetic field due to the nuclear spin moment, which interacts with the electron spin is given as

$$\vec{B}(\vec{r})=\frac{-g_{K}\mu_{N}}{c^{2}}[\frac{\hat{I}_{K}}{|\vec{r}-\vec{R}_{K}|^{3}}-\frac{3(\vec{r}-\vec{R}_{K})((\vec{r}-\vec{R}_{K})\cdot\hat{I}_{K}}{|\vec{r}-\vec{R}_{K}|^{5}}-\frac{8\pi}{3}\hat{I}_{K}\delta_{\vec{r},\vec{R}_{K}}]$$

(9)

where, $g_{K}$ is the nuclear g-factor and $\mu_{N}$ is the nuclear magnetic moment This field interacts with the electron through $g_{e}\mu_{B}\bm{s}\cdot\vec{B}$ and gives rise to the following set of nuclear-electronic (ne) and nuclear-nuclear(nn) interactions:

$$\hat{H}_{ne}^{SD}=g_{e}g_{K}\mu_{B}\mu_{N}\bm{s}\cdot\frac{(\vec{r}^{T}_{iO}\vec{r}_{iO}\textbf{1}-\vec{r}_{iO}\vec{r}^{T}_{iO})}{\vec{r}_{iO}^{5}}\hskip 14.22636pt\text{Spin-Dipolar (SD)}$$

(10)

$$\hat{H}_{ne}^{FC}=\bm{s}\cdot\hat{I}_{K}\delta_{\vec{r},\vec{R}}\frac{8\pi g_{K}g_{e}\mu_{N}\mu_{B}}{3c^{2}}\hskip 14.22636pt\text{Fermi-Contact(FC) operator}$$

(11)

$$\hat{H}_{ne}^{PSO}=(\frac{g_{K}\mu_{N}}{c^{2}})\frac{\hat{I}_{K}\cdot\vec{L}_{O}}{|\vec{r}_{iO}|^{3}}\hskip 14.22636pt\text{Paramaganetic Spin Operator (PSO)}.$$

(12)

In the presence of an external magnetic field, the diamagnetic coupling (linear in $\vec{B}_{ext}$) of the external magnetic field, with an associated gauge centered at G, is given by

$$\hat{H}^{DS}_{nn}=\frac{g_{K}\mu_{N}}{2c^{2}}\vec{B}_{ext}\frac{\vec{r}_{iG}^{T}r_{iO}\textbf{1}-\vec{r}_{iO}r_{iG}^{T}}{|r_{iO}|^{3}}\hskip 14.22636pt\text{Nuclear Diamagnetic Shielding operator(DS)}$$

(13)

One can notice that the interaction terms, which scale as $c^{-2}$ could be treated as weak perturbations and thus we obtain the shielding tensor by correcting Diamagnetic Shielding (DS) operator to the second order

$$\bm{\sigma}_{K}=\langle\phi_{HF}|\hat{H}_{ne}^{DS}|\phi_{HF}\rangle-2\sum_{i\neq 0}\frac{\langle\phi_{HF}|\hat{L}_{iO}|\phi_{i}\rangle-\langle\phi_{HF}|\hat{H}^{PSO}_{ne}|\phi_{HF}\rangle}{E_{HF}-E_{i}}$$

(14)

The above expression is due to Ramsey(33). Here the magnetic field strength is a few Teslas, but when the molecules are in strong magnetic fields, the perturbation cannot be treated as “weak”. This problem is discussed in SectionIII.

From Eq.14, the shielding tensor depends on $\hat{H}_{ne}^{DS}$ and $\vec{L}_{iO}$, which depends on the choice of gauge origin. We notice that the shielding tensor (and chemical shift) should not depend on the coordinates of the molecule in space, that is, it should be gauge independent. This can be assured when the calculation is performed with a complete basis set. However, for practical reason (computational cost and power), we often use a optimized basis set (determined through convergence tests). The first term in Eq.14 is just the expectation value of and can be determined with a good precision. But in the second term, the effect of magnetic field may generate a function outside our defined basis set and hence, the accuracy of the second term is compromised. This can be illustrated by considering a basis set consisting of Gaussian Type Oribtals (GTO), with certain Gaussian centers. Further, let the gauge origin span complete set of Gaussian centers. In a complete basis, irrespective of the gauge origin, a function generated upon a gauge dependent operation necessarily lies in the basis. But in an optimized basis, the gauge origin might lie outside the set of Gaussian centers and hence a function outside the basis set could be generated. This introduces the gauge origin problem. Choice of gauge is just a tool for convenient mathematical construction, but the observables are gauge invariant. We can explore the prescription for the gauge origin to be the nuclear position. From elementary gauge transformations, the vector potential transforms as

$$\textbf{A}_{\textbf{O'}}(\vec{r})=\textbf{A}_{O}(\vec{r})+\nabla\cdot f(\vec{r}),\hskip 14.22636pt\nabla\times\nabla\cdot f=0$$

(15)

where $f(\vec{r})$ is a scalar function and O’ is the position vector of the transformed origin. The conditions in Eg.15 ensure gauge invariance. The wavefunction, in case of such a transformation, picks up a phase $e^{-if\cdot r}$, where $f(\vec{r})$ has the following definition.

$$f(\vec{r})=\frac{1}{2}\vec{B}\times(\textbf{O'-O})$$

(16)

So applying this to a one electron Slater orbital $\chi_{lm}$, we obtain an overall wavefunction, which is known as the Gauge Included Atomic Orbital (GIAO).

$$\omega(\textbf{A}_{\textbf{O}})=e^{-\iota\frac{1}{2}\vec{B}\times(\textbf{N}-\textbf{O})}\chi_{lm}$$

(17)

We shall encounter this scheme of writing atomic orbitals in future contexts, however the notation might differ slightly from Eq.17.

The chemical shift tensor is proportional to the shielding constant and can be defined as

$$\bm{\sigma^{\prime}}(\vec{r})=-\vec{B}_{in}\vec{B}^{T}$$

(18)

where, $\vec{B}_{in}$ is the induced magnetic field due to the ring current density $\textbf{j}(\vec{r}$), defined as

$$\vec{B}_{in}=\frac{1}{c}\int d^{3}r^{\prime}\textbf{j}(\vec{r^{\prime}})\times\frac{\vec{r}-\vec{r^{\prime}}}{|\vec{r}-\vec{r^{\prime}}|^{3}}$$

(19)

The current density is the crucial component of these integrals and is treated in two fashions-the GIPAW method (Gauge Included Projector Augmented Wave) or the CTOCD method (Continuous Transformation of Current Density). The GIPAW method, explained by Pickard and Mauri (32) is used in the popular software packages CASTEP(3; 18) and Wien2K(23; 24). The CTOCD-DZ (DZ:Diamagnetic contribution set to zero)(21; 25), described by Ligabue et al(29), is implemented in the software package DALTON(1). It is beyond the limitations of this review, to condense these theories and to discuss their implementation. Nonetheless, we leave the discussion by remarking that these methods are primarily different in the choice of bases they use; the GIPAW method is suitable for plane wave basis set and therefore advised for the solid state, while the CTOCD-DZ uses Gaussian type basis and is more suitable for molecules. The reader is strongly advised to refer to the cited articles above for a detailed explanation of these methods.

In the previous section, we treated the effect of external magnetic field as a weak perturbation. However, in the cases where the magnetic field strengths appear in the magnitudes of 10^4-7T, the problem could no longer be solved using perturbation theory.
We start our discussion with the simplest atom; the hydrogen atom. Kappes and Schmelcher (20) presented basis optimization of London type wavefunctions. For this, we begin by choosing out basis of the kind,

$$\Psi(\textbf{r},\alpha,\textbf{R},\textbf{C})=e^{-\iota\textbf{A}\cdot\textbf{r}}(x-R_{x})^{n_{x}}(y-R_{y})^{n_{y}}(z-R_{z})^{n_{z}}e^{-\alpha(\textbf{r}-\textbf{R})^{2}}$$

(20)

where, C is the position vector chosen as the gauge origin (otherwise spanned by position vector R), and A(C) is the vector potential at position C. The complex phase multiplied at the beginning promises gauge invariance. In our calcuations, $\alpha$ enters as a 3$\times$3 symmetric matrix (since r-R is a three component vector), which we shall use as a variational parameter for basis set optimization.
The external magnetic field could be chosen to be acting in the z-direction, which naturally conserves the z-component of the angular momentum. It is therefore, convenient to work with cylindrical coordinates. This modifies our basis set in the following fashion,

$$\Phi_{k,l}^{m_{\pi}}=\rho^{m+2k}\cdot z^{l}\cdot e^{(-\alpha\rho^{2}-\beta z^{2})}e^{\iota m\phi}$$

(21)

where, $\alpha$ and $\beta$ are the matrix elements, $\alpha_{xx}$=$\alpha_{yy}$= $\alpha$ and $\alpha_{zz}$=$\beta$. The Pauli Hamiltonian matrix (see below) could then be constructed using the basis functions in Eq.21 and then diagonalized to obtain eigenenergies.These eigenenergies are in terms of a given set of variational parameters $\alpha_{i},\beta_{i}$. The ground state energy could be minimized with respect to these parameters iteratively. It is not of our interest to discuss optimization algorithm. A more useful result, which the authors reported, is the dependence of $\alpha_{i},\beta_{i}$ on the magnetic field.

$$\hat{H}_{Pauli}=\hat{H}_{el}-\frac{e}{m}\hat{L}\cdot\vec{B}+\frac{e^{2}}{8m}|\vec{B}\times\textbf{r}|^{2}$$

(22)

Here, $H_{el}$ is the Hamiltonian of the unperturbed molecule.
It was observed by the authors that for low field strengths, $\alpha_{i}=\beta_{i}$, but they start to deviate from each other in the high field regime. This attributed to the different forms of orbital representation in the directions parallel and perpendicular to the field. In the parallel direction, the orbital dynamics is captured by Slater type representation $e^{-c|{z}|}$, while the in-plane orbital dynamics (perpendicular to the direction of field) could be described by the Landau orbitals¹¹1Landau Orbitals are functions that describe Landau levels in a symmetric gauge. It can be represented as $\Psi(x,y)=f(z)e^{c|z|^{2}}$, where $z=x+\iota y$. This is a useful result to choose the type of orbitals we shall use as basis functions for our calculations. The authors report results on dissociation energies and bond lengths as a function of magnetic field strength (upto 1 a.u.), but we postpone the discussion for later sections. The London orbitals were further popularised by Tellgren et al (38) where they described non-perturbative formalism for quantum chemical calculations in magnetic field.
In their article, they describe the London orbitals in the same way as in Eq.20, but with a slightly different interpretation. The phase factor we saw in Eq.20 appeared to addressed the well known gauge origin problem. But in this formalism, the phase is treated as a plane wave with the wavevector having the same form as that of A. Having said so much, we can attempt to formulate a two-electron problem, following the McMurchie-Davidson scheme(30). Let us consider two single electron functions, described by Eq.20, $\phi(\textbf{r,C})$ and $\phi(\textbf{r,D)}$, where C and D are position vectors of the gauge origin in each case. We can define a new “center” or gauge origin for the overlap of these two orbital functions as

$$\textbf{P=}\frac{c\textbf{C}+d\textbf{D}}{c+d}$$

(23)

where, c and d are Gaussian coefficients of $\phi(\textbf{r,C})$ and $\phi(\textbf{r,D)}$ respectively, and $p=c+d$. The composite system of two function is described by the simple product rule

$$\phi(\textbf{r,C})\cdot\phi(\textbf{r,D})=\Omega^{k_{1}}(\vec{r}_{P})=e^{-\iota k_{1}\cdot r}\Omega(\vec{r}_{P})$$

(24)

where $k_{1}=k_{C}-k_{D}$. Next, a magnetic field dependent spherical London orbital (known as the London Hermite-Gaussian wavefunction) is introduced, which describes the composite system. We obtain our energy expectation values in terms of this function.

$$\Lambda_{tuv}^{k_{1}}(\vec{r}_{P})=(\frac{\partial}{\partial P_{x}})^{t}(\frac{\partial}{\partial P_{y}})^{u}(\frac{\partial}{\partial P_{z}})^{v}e^{-p}\vec{r}_{P}^{2}e^{-\iota k_{1}\cdot\vec{r}_{P}}$$

(25)

We notice that in Eq.25, the function is secured from gauge origin problem as $k_{1}=\frac{1}{2}\vec{B}\times(C-D$) depends only on relative separation of Gaussian centers C and D. We can write Eq.25 more compactly in the following form

$$\Lambda_{tuv}^{k_{1}}(\vec{r}_{P})=\Lambda_{tuv}(\vec{r}_{P})e^{-\iota k_{1}\cdot\vec{r}_{P}}$$

(26)

Since the representations in Eq.24 and Eq.25 describe the same system, their equivalence could be established by expanding $\Omega^{k_{1}}$ as shown below

$$\Omega^{k_{1}}=\sum_{tuv}E_{t}^{n_{x}^{c}n_{x}^{d}}E_{u}^{n_{y}^{c}n_{y}^{d}}E_{v}^{n_{z}^{c}n_{z}^{d}}\Lambda_{tuv}^{k_{1}}(\vec{r}_{P})$$

(27)

where the expansion coefficients are field independent. It is then straightforward to show, from Eq.26 that

$$S^{k}_{1}=\sum_{tuv}E_{tuv}^{n_{i}^{c}n_{i}^{d}}\Lambda_{tuv}(k_{1})$$

(28)

where $\Lambda_{tuv}(k_{1})$ is the Fourier transform of $\Lambda_{tuv}^{k_{1}}(\vec{r}_{P})$.

To describe the details of this formalism is beyond the scope of this review and hence we shall summarize the use of this aforementioned formalism by writing the Coulomb repulsion term in the following manner, for the composite system.

$$(\phi_{C}\phi_{D}|\phi_{F}\phi_{G})=\sum_{tuv}\sum_{\tau\nu\mu}E_{tuv}^{n_{i}^{c}n_{i}^{d}}E_{\tau\nu\mu}^{n_{i}^{f}n_{i}^{g}}(\Lambda_{tuv}^{k_{1}}(\vec{r}_{P})|\Lambda_{\tau\nu\mu}^{k_{1}}(\vec{r}_{P}))$$

(29)

The overlap of the London Hermite-Gaussian functions, in the summation, is computed as

$$(\Lambda_{tuv}^{k_{1}}(\vec{r}_{P})|\Lambda_{\tau\nu\mu}^{k_{1}}(\vec{r}_{P})=\frac{2\pi^{5/2}}{pq\sqrt{p+q}}exp[\frac{k_{1}^{2}}{4p}]exp[\frac{k_{2}^{2}}{4q}]\times R_{tuv,\tau\nu\phi}^{0}$$

(30)

where, $R_{tuv,\tau\nu\phi}^{0}$ is the zeroth order (n=0) function described by

$$R_{tuv,\tau\nu\phi}^{n}=\frac{\partial^{t+u+v}}{\partial P_{x}^{t}\partial P_{y}^{u}\partial P_{z}^{v}}\frac{\partial^{\tau+\nu+\phi}}{\partial Q_{x}^{\tau}\partial Q_{y}^{\nu}\partial Q_{z}^{\phi}}\times e^{-\iota\vec{k}_{1}\cdot\textbf{P}}e^{-\iota k_{2}\cdot\textbf{Q}}(-2\alpha)^{n}F_{n}(\alpha(\textbf{P'-Q'}))$$

(31)

where, P’=P- $\frac{i}{2p}\vec{k}_{1}$, Q’=Q- $\frac{i}{2q}\vec{k}_{2}$, $\alpha=\frac{pq}{p+q}$ and $F_{n}$ is the nth-order Boys Function. The authors calculate the required $R^{0}_{tuv,\tau\nu\phi}$ using downward recursion (set of equations described in Eq.16 of Tellgren et al(40)), where we start with $n^{th}$ order function (n$>$ 0) in the recursion and recover the zeroth order function. It is important to note that the effect of magnetic field only enters our calculation through the wavevectors $\vec{k}_{1}$ and $\vec{k}_{2}$. Hence, the interaction terms, like in Eq.29, are a product of recursively solved molecular integrals, that include the effect of magnetic field. That is to say, this method does not treat the effect of magnetic field as a perturbation and therefore is a non-perturbative method. One can appreciate the fact that the theory presented above is suitable for any magnitude of the magnetic field. We rejoice over the fact that such non-perturbative treatment is appropriate to calculate molecular properties in strong magnetic fields.
This theory has been implemented at the Hartree-Fock level of computation at in the software package LONDON. The authors report their results on benzene, the HF molecule and a few diamagnetic systems, but we shall postpone our discussion on real molecules in the next section.
In our discussion so far, we have not considered the structural consequences of subjecting molecules to strong magnetic fields. This problem is addressed by the authors of Tellgren et al(39), where they lay down the theory for computing molecular gradient that enables us to obtain the geometrically optimized structure in high magnetic fields.

We again start with the London type orbitals (they term it as Gaussian Type Orbitals-Plane Wave Orbitals, or GTO-PW orbitals), but now the wavefunctions are expressed in terms of a contracted basis

$$\phi_{\alpha}(\vec{r})=\sum_{\textit{P}}\chi_{p}(\vec{r})U_{\alpha}^{p}$$

(32)

where P is the set of all primitive atomic orbitals and $\phi_{\alpha}$ is the set of real basis functions transformed, from the primitive basis, by the matrix $U_{\alpha}^{p}$. Another linear operator $\hat{L}$ (not to be confused with angular momentum operator) is introduced, which generates a set of primitive functions, which may not be contained in P. This operation is defined as follows.

$$\hat{L}\chi_{p}(\vec{r})=\sum_{\tilde{L}}v_{q}(\vec{r})L_{p}^{q}.$$

(33)

The linear operator $\hat{L}$ acts on the real basis functions to yield the double summation

$$\hat{L}\phi_{\alpha}(\vec{r})=\sum_{\textit{P}}\sum_{\tilde{L}}v_{q}(\vec{r)}L_{p}^{q}U_{\alpha}^{p}$$

(34)

The matrix elements for this linear operation on basis functions, such that $\hat{L}\phi=\psi$, is given by

$$L_{\alpha\beta}=\langle\phi_{\alpha}|\hat{L}|\phi_{\beta}\rangle=\langle\phi_{\alpha}|\psi_{\beta}\rangle$$

(35)

Now the authors consider $\hat{L_{1}}$ and $\hat{L_{2}}$ to be two linear operators mapping the space of primitive functions to space containing linear combination of those. One can represent an overlap integral/matrix element for an operator $\hat{\Omega}$ in the terms of non-primitive basis functions $\psi$ as

$$\Omega_{\alpha\beta}=\langle\psi_{\alpha}|\hat{\Omega}|\psi_{\beta}\rangle$$

(36)

The benefit of using non-primitive basis circumvents the memory issues in computing matrices of larger dimensions that we otherwise get while dealing with primitive basis. According to the authors, developing such a formalism aids encoding the quantum mechanics in object oriented programming languages (here, C++).
The authors treat the linear operator, mentioned above, as the partial derivative with respect to the nuclear coordinate C. They solve molecular gradients for the one dimensional case and employ it within RHF (Restricted Hartree Fock) theory to study Helium clusters. The molecular gradient for the one dimensional case is

$$\frac{\partial E}{\partial C_{x}}=\frac{\partial D^{ji}}{\partial C_{x}}F_{ij}+\sum_{i}(D^{ji}F_{\frac{\partial_{i}}{\partial C_{x}},j}+c.c.)-\langle\chi_{i}|\sum_{N}\frac{\partial}{\partial C_{x}}\frac{Z_{N}}{|\vec{r}-\vec{R}_{N}|}|\chi_{j}\rangle D^{ji}$$

(37)

where $D$ is the density matrix, whereas the elements $F_{ji}$ constitute the Fock matrix. The authors studied He₃, He₄ and He₆ clusters, observed in the atmosphere of white dwarfs. To directly relate the shape of the clusters with the magnetic field strength is complicated, nonetheless, the cluster dimensions (measured as bond distances) reduce with increased magnetic field (see Fig.1).

What causes such a contraction? Well, direct answers to the case of Helium clusters may not be available, but this is the appropriate juncture to discuss a distinguished bonding mechanism, known as the paramaganetic bonding, which we shall encounter in most of our discussions throughout this review.

Lange et al (22) first proposed the paramagnetic bonding for the H₂ molecule, which they attribute to stabilization of anti-bonding orbitals perpendicular to the external magnetic field. In their calculations, authors perform FCI (Full Configuration Interaction) calculation, using an aug-cc-pVTZ uncontracted basis set, where the wavefunction is written as the anti-symmetric product of the spinors $\phi_{p}$ and the coefficients $C_{n}$ are obtained from Rayleigh-Ritz variational theory.

$$|\Psi\rangle=\sum_{n}^{N}C_{n}det|\phi_{p1}\phi_{p2}\phi_{p3}...\phi_{pN}|$$

(38)

This calculation is performed for GIAOs and therefore the set of equations discussed earlier are useful here as well.
In case of a diatomic species, the authors distinguish two kinds of bonds- parallel to the applied field and perpendicular to the applied field. For the singlet case of the H₂ molecule, the parallel bonding is favored due to a greater overlap. The authors confirm this by plotting polar plots of potential energy surfaces in magnetic field (see Fig.1 of Ref.(22)) .

The triplet, on the other hand, behaves conversely; electronic energy is lowered with the increased field strength, with preferred perpendicular orientation. To understand this stabilization mechanism, when the 1$\sigma_{g}$ molecular orbital is written as two 1s Gaussian orbitals, we realise that when the internuclear distance tends to zero, the MOs of the H₂ molecules could virtually be treated as atomic orbitals of He atom.

$$1\sigma_{g/u}=[2\pm 2e^{-a/2(1+B_{x}^{2}+B_{y}^{2})R^{2}}]^{-1/2}\times(1s_{A}\pm 1s_{B})$$

(39)

Therefore, the authors prescribe that the 1$\sigma_{g}$ could be viewed as the 1s orbital of He, while the anti-bonding orbitals could be treated as $2p_{0}$ and $2p_{-1}$, with the former preferring parallel orientation and the latter prefers perpendicular orientation. It is found that the $2p_{-1}$ is stabilized than the $2p_{0}$, in strong magnetic field due to orbital Zeeman coupling. This hypothesis is verified by computing the bond energies for both the parallel and perpendicular orientations. Their plots clearly show that the 1$\sigma_{u\bot}^{*}$ is stabilized as compared to 1$\sigma_{u\|}^{*}$ (Fig.2).

In this subsection, our aim is to discuss the examples where Hartree-Fock methods were employed, which significantly added to the conceptual framework of response of molecules in strong magnetic fields.
We go back in 1999 to discuss the results of Ivanov et al (15) , where they study the ground state of the carbon atom in strong magnetic fields. They perform calculations using the 2D Unrestricted Hartree-Fock method. They start by writing a simple expression for single particle energy of carbon atom in magnetic field as

$$\epsilon_{B\mu}=(m_{\mu}+|m_{\mu}|+2S_{z\mu}+1)\cdot B/2-\epsilon_{\mu}$$

(40)

where, $m_{\mu}$ is the magnetic quantum number, $S_{z\mu}$ is the z-projection of the spin $S_{\mu}$ and $\epsilon_{\mu}$ ia the one-electron energy.

From their calculations, it is easy to identify configurations in the extremes-the field free and $B\rightarrow\infty$ case. The ground state (GS) configuration for the field-free case is $1s^{2}2s^{2}2p_{0}^{1}2p_{-1}^{1}$ and it retains it form till a maximum field strength of 0.1862 a.u. For the latter case, the authors assume the fully polarized configurations, where each magnetic state (identified by the magnetic quantum number) contains one electron. Such a fully spin polarized configuration is given by $1s2p_{-1}3d_{-2}4f_{-3}5g_{-4}6h_{-5}$. While it is easy to identify these two ground states, it is difficult to identify the GS for intermediate field strengths. Hence, it was found that there exist multiple (5) ground states for the intermediate field strengths, which occur for the spin states $S_{z}=-1,-2,-3$. These observations are tabulated in Table1.

The same formalism is extended, in their article in 2000(16) , to all the atoms from hydrogen to neon. Their results lay an important foundation for the results we shall discuss next in the context of these ions. In this article, they investigate the magnetic fields at which the transition from the fully spin polarized configuration happens to the field-free ground state. An important conclusion drawn from this study is that, for Z$\geq$6, there could be two competing fully spin polarized configurations, namely $1s2s2p_{0}$… and $1s2p_{0}3d_{-1}...$. For lighter atoms, there exists a global fully spin polarized configuration, represented by$1s2p_{-1}3d_{-2}4f_{-3}5g_{4}6h_{-5}$.

We now present a more modern Hartree-Fock based method, presented by Boblest et al (2) , which couples with Diffuse Quantum Monte Carlo (DQMC) to compute the energy spectrum of molecules in intense magnetic fields. In this method, authors write single particle spinors $\Psi^{i}$ in terms of Landau orbitals $\Phi_{i}$

$$\Psi^{i}(\rho_{i},\phi_{i},z_{i})=\sum_{n=0}^{N}P_{n}^{i}(z_{i})\Phi_{nmi}(\rho_{i},\phi_{i})$$

(41)

where, $P_{n}^{i}$ is expanded in terms of B-splines as follows

$$P_{n}^{i}(z_{i})=\sum_{\nu}\alpha_{n\nu}^{i}B_{\nu}^{i}(z_{i})$$

(42)

Here the B-spline coefficients $\alpha_{\nu}$ enter as variational parameter and minimization of the energy functional generates a set of equations known as the Hartee-Fock-Roothaan equations, which are solved iteratively. The energies obtained here are then improved by the DQMC, for which the importance sampled imaginary time Schrödinger Equation is solved. Here, the trial wavefunction is constructed simply by the forming a product of all spins up ($\psi_{\alpha}$) and all spins down ($\psi_{\beta}$) states, modulated by the Pade-Jastrow factor J(8).

$$\Psi_{T}(\vec{R})=J\psi_{\alpha}\cdot\psi_{\beta}$$

(43)

The Pade-Jastrow factor is essentially optimized in the Monte-Carlo calculations. The energy is estimated by computing the local energy as $H\Psi_{T}/\Psi_{T}=E_{L}$.
The authors find that the previously reported configurations for the ground and fully polarized state, for 2 to 4 electron systems is consistent with their results. They also report variation of the transition magnetic field (magnetic field at which configuration transition occurs) with the atomic number Z. In such cases, as the Z number increases, the nuclear interactions overwhelm electron-electron interactions. This leads to a saturation of the transition magnetic field at about 0.17058 a.u. Here, we introduce the following notation to denote the atomic configuration of atoms in different field strength regimes. The configuration of Fully Spin Polarized state $1s2p_{-1}3d_{-2}4f_{-3}5g_{-4}6h_{-5}$ is denoted by $|0\rangle$. Then the state $1s^{2}2p_{-1}3d_{-2}4f_{-3}5g_{-4}$ is compactly written as $|1s^{2}...\rangle$, that is the ket contains the orbital that is different from the $|0\rangle$ configuration. For He-like, Li-like and Be-like ions, the variation in the magnetic field strengths for the $|1s^{\alpha}...\rangle\rightarrow|0\rangle$ transition, with the atomic number is presented in Fig.3.

The authors adopt a similar procedure for the Boron and Carbon like ions. For $C^{+}$, they determine that $|1s2p_{0}...\rangle$ is the ground state for field strengths lower than 0.0788 a.u. They also predict a possible transition path, for different magnetic field regimes, for 6-electron systems as $|1s2p_{0}...\rangle\rightarrow|1s...\rangle\rightarrow|2p_{0}...\rangle\rightarrow|0\rangle$. Their results for the neutral atom differ from those reported by Ivanov et al; they also find a $|2p_{0}3d_{-1}...\rangle$ as a possible ground state. These discrepancy needs to be resolved by studying the fully spin polarised states in Neon by other methods including correlation or exchange interactions.

Following the aforementioned formalism of 2D HF, Thirumalai et al (41) in 2014, solved the problem of the carbon in a slightly different fashion. The solution to single electronic hydrogenic system is obtained, without considering any exchange or direct interactions. The wavefunctions obtained are useful for two reasons; first, they work as initial estimates for the (coupled) Hartree-Fock problem and second, they are used to determine solutions to the elliptic partial differential equations (as described in (41) , which give us direct and exchange potentials. The coupled Hartree-Fock problem is solved iteratively to extract eigenstates and eigenenergies.
They observe that there are twelve possible fully spin polarized configurations for the carbon atom, of which only two are reported by Ivanov et al, in their article in 2000. The authors claim that the key success of this method is its low computation cost at non-compromised accuracy. Hence, this method is also viable than FCI, which might end up being too expensive for a 6-electron system. However, the authors point out certain limitations of this method. Here, the relativistic effects are neglected, which might prove important for such high magnetic fields. Secondly, their equations treat the mass of nucleus as infinite in comparison with that of the electron. But in the case of ultrahigh magnetic fields, nuclear effects would also become relevant and influence of nuclear effects should be investigated. Lastly, as most Hartree-Fock methods do, this one also ignored electron-electron correlation effects.

It is well known that closed shell systems are expected to exhibit diamagnetic (induced magnetic field opposes external one) ground state in the zero field caase. However, works of Hegstrom and Lipscomb(11; 12; 10) suggested the existence of molecules like BH (and isoelectronic compounds, due to Fowler and Steiner(5)) exhibit paramagnetism in the ground state, in the field-free case. On the contrary to the behaviour of conventional closed shell molecules, paramagnetic closed shell systems undergo a transition to a diamagnetic excited state at high magnetic fields.
Tellgren et al(38) performed Hartree-Fock based calculations for such molecules in strong magnetic fields. An important output of their work is a simple two level model, which is used to rationalize and fit the calculated energy spectrum. The part of the total Hamiltonian, that couples with the external magnetic field, which is assumed to be perpendicular to the bond, is given as

$$\hat{H}_{B}=\vec{B}\cdot\hat{L}_{\alpha}-\chi_{\alpha\alpha}\vec{B}^{2}=\sum_{i}\vec{B}\cdot\hat{l}_{i}-\chi_{\alpha\alpha}\vec{B}^{2}$$

(44)

where $\hat{L_{\alpha}}$ is the total angular momentum operator and $\hat{l_{i}}$ is the single particle angular momentum operator. To demonstrate this model, the authors consider a closed shell system of 6 electrons such as BH or the CH⁺ ion. The symmetric ground state is therefore written as

$$\Psi_{0}=|\phi_{1}\bar{\phi_{1}}\phi_{2}\bar{\phi}_{2}\phi_{H}\bar{\phi}_{H}|$$

(45)

where the notation is such that wavefunctions indexed with ’H’ denote the HOMO and those with ’L’ denote the LUMO. The doubly degenerate excited states are defined as

$$\Psi_{1x/y}=\frac{1}{\sqrt{2}}|\phi_{1}\bar{\phi}_{1}\phi_{2}\bar{\phi}_{2}(\phi_{H}\bar{\phi}_{Lx/y}+\phi_{Lx/y}\bar{\phi}_{H}|.$$

(46)

Here, let us suppose the $2p_{z}$ orbital is occupied and hence the LUMO is comprised of two $2p_{x}$ and $2p_{y}$, centered on the carbon or boron atom. The Hamiltonian matrix elements between ground and first excited states can then be formulated as

$$\langle\Psi_{0}|\hat{H}|\Psi_{1x/y}\rangle=\langle\phi_{H}|\hat{l}_{x/y}\cdot\vec{B}|\phi_{Lx/y}\rangle$$

(47)

We write $\phi_{H}$ as a linear combination of atomic orbitals, which contribute to hybridization (1s, 2s and $2p_{z}$), and make the substitutions $2p_{y}=p_{0},2p_{x}=\frac{-i}{\sqrt{2}}(p_{1}-p_{-1})$ and $2p_{z}=\frac{1}{\sqrt{2}}(p_{1}+p_{-1})$, where $p_{n}$ are the eigenstates of the angular momentum operator $\hat{l}_{y}$. These substitutions simplify the above equation to yield the matrix elements,

$$\langle\Psi_{0}|\hat{H}|\Psi_{1x}\rangle=\iota\mu B,\hskip 11.38092pt\langle\Psi_{0}|\hat{H}|\Psi_{1y}\rangle=0.$$

(48)

where $\mu$ is the magnetic dipole moment. A similar approach is followed for an N-membered carbon ring with one $p_{z}$ orbital at each site. The candidates considered here by the authors are $C_{4}H_{4},C_{6}H_{6}$ and $C_{12}H_{12}$. These conjugated systems possess non-degenerate HOMO and LUMO, which correspond to rotationally allowed eigenstates (eigenstates of $\hat{l}_{z}$). The HOMO and LUMO are written as a linear superposition of the eigenstates of the single particle angular momentum operator.

$$\Psi_{HOMO}=\frac{1}{\sqrt{2}}(|+n\rangle+|-n\rangle),\hskip 11.38092pt\Psi_{LUMO}=\frac{-\iota}{\sqrt{2}}(|+n\rangle-|-n\rangle)$$

(49)

Upon calculating the Hamiltonian Matrix elements, a two level Hamiltonian can be constructed as follows,

$$\hat{H}=\left(\begin{matrix}-\Delta-\chi_{0}B^{2}&&\iota\mu B\\ -\iota\mu B&&\Delta-\chi_{1}B^{2}\end{matrix}\right)$$

(50)

and the eigenvalues obtained upon diagonalizing the Hamiltonian yield energies of the ground and excited state as described in Eq.51.

$$E_{0/1}=\frac{-1}{2}(\chi_{0}+\chi_{1})B^{2}\mp\frac{1}{2}\sqrt{(2\delta+(\chi_{0}-\chi_{1})B^{2})^{2}+4\mu^{2}B^{2}}$$

(51)

where $\chi_{0/1}$ are expectation values of the ground and excited state magnetizability tensors, $\Delta$ is the energy difference between the two states in the field-free case. It is then possible to write conditions for the magnetic phase transitions. From the definition of magnetizibility, the magnetic phase transition occurs when $\frac{\partial^{2}E_{0/1}}{\partial B^{2}}=2\chi_{0/1}\pm\frac{\mu^{2}}{4}$ changes its sign. Similarly, the critical magnetic field for magnetic phase transition could be obtained from the condition $\frac{\partial E_{0/1}}{\partial B}=0$. The critical magnetic field for such a phase transition (magnetic field at which paramagnetic to diamagnetic transition occurs) is given as

$$B_{c}=\pm\frac{\sqrt{-2\chi_{0}\chi_{1}[\mu^{2}+\Delta t]+|\mu(\chi_{0}+\chi_{1})|\sqrt{\chi_{0}\chi_{1}(\mu^{2}+2\Delta t}}}{\chi_{0}\chi_{1}t^{2}}$$

(52)

where $t=\chi_{0}-\chi_{1}$. Further, imposing the condition for the existence of a critical magnetic field, we get a simple condition

$$\frac{\mu^{2}}{2\Delta}>|\chi_{0}|.$$

(53)

The authors also conclude from this model that the paramagnetic ground state and diamagnetic excited stated only cross when B$\rightarrow$ $\infty$, that is to say that such a crossing is only possible at ultrahigh magnetic fields, not the ones realizable on Earth. In addition, a limitation of this approach is that it does not take into account the energy changes due change in geometry of molecule at high fields.
To demonstrate the accuracy of the two level model, the authors compare exact magnetizability (perpendicular field), obtained from linear response function, with that obtained from a eighth order polynomial fit and the two level model for the four molecules, as tabulated below (Table 2). We can clearly see that the two level model is more accurate fit than the eighth order polynomial fit for all the four molecules.

In most cases discussed above, the electron-electron correlation was ignored. In this section, methods like the coupled cluster theory and density functional theory are discussed in the context of molecules in high magnetic field.
Stopkowicz, et al(36) proposed the Coupled Cluster theory for calculating correlation energy in molecules subjected to intense magnetic fields. The authors write the “coupled cluster” wavefunction as

$$|\Psi_{cc}\rangle=e^{\hat{T}}|\Phi_{0}\rangle$$

(54)

where, $\hat{T}$ is the cluster operator defined as

$$\hat{T}=\hat{T}_{1}+\hat{T}_{2}+...=\sum_{n=1}^{N}(\frac{1}{n!})^{2}\sum_{i,j,k...,a,b,c...}t_{ij}^{ab}\hat{a}_{a}^{\dagger}\hat{a}_{b}\hat{a}_{i}^{\dagger}\hat{a}_{j}$$

(55)

where, a,b and i,j are indices for the virtual and occupied orbitals respectively. $|\Phi_{0}\rangle$ is chosen as Hartree-Fock wavefunction. Therefore, the correlation energy is computed as the Hartree-Fock expectation value.

$$E_{corr}^{cc}=\langle\Psi_{HF}|e^{-\hat{T}}\hat{H}_{N}e^{\hat{T}}|\Psi_{HF}\rangle$$

(56)

That is to say that the total Hamiltonian is the sum of Hartree-Fock energy ($E_{HF}$) and the coupled cluster Hamiltonian ($H_{N}$), which describes correlation effects.

This model is then applied to the simplest case where correlation effects could be studied; the Helium atom. In this case, the diamagnetic to paramagnetic crossover is observed at approximately $0.8a.u.$ of magnetic field. The authors show that the energy of the diamagnetic singlet state increases continuously with magnetic field, while that of the triplet paramagnetic state decreases. This is linked by the authors to the spatial extension of orbitals , which destablilize the singlet. They present the argument by comparing the expectation value of the diamagnetic contribution (Eq.57), in the correlated regime, to that of the Hartree-Fock result.

$$E_{dia}=\frac{1}{8}B_{z}^{2}\langle\Psi|\hat{r}^{2}|\Psi\rangle$$

(57)

This quantity, for HF level of theory, is lower as compared to the CCSD (Coupled Cluster Singles and Doubles) theory. Qualitatively, it means that there is more correlation in the singlet state, where we have 2 electrons in the 1s orbital, whereas the correlation is lowered upon promoting one electron to the 2p orbital, in case of the triplet. In addition, the correlation in the 1s would cause expansion of the orbital, reducing the screening for the 2p orbital. This causes spatial contraction in the 2p orbital. However, a general trend is not observed; for instance the singlet and triplet states show greater spatial extension in case of Hartree-Fock while the singlet of Ne shows greater spatial extension in case of coupled cluster theory calculations. The authors also report the case of LiH molecule. The paramagnetic bonding in LiH is interpreted in two ways; by noting the destabilization of the ${}^{1}\Sigma_{\|}$ and stabilization of ${}^{1}\Sigma_{\bot}$ with increasing magnetic field. Secondly, the binding energy profile for LiH in Fig.4 shows continuous increase in the binding energy for magntic field strength greater than 0.2 a.u. Reconciliation of these two observations confirm paramgnetic bonding in LiH.

Around the same time, in 2015, Furness et al(7) devised a way to determnine magnetic response of molecules using density functional theory, known as the Current Density Functional Theory (CDFT). The authors introduce a formalism that uses the meta-GGA (Meta-Generalized Gradient Approximation) exchange-correlation. Meta-GGA (MGGA) functionals have proven to be very successful in predicting effects of correlation, and therefore could potentially be exploited to address the problems posed by HF based methods. However, traditional MGGA functionals depend on canonical kinetic energy operator, in which case, gauge invariance does not hold. In order to circumvent this issue, the authors replace the canonical kinetic energy by its generalized form including the paramagnetic current density.
The mathematical formulation of CDFT is compactly written as in Eq.58

$$[\frac{\hat{p}^{2}}{2m}+\frac{1}{2m}{\hat{p},\hat{A_{s}}}+V_{ext}+\frac{1}{2}\hat{A}_{ext}^{2}+V_{C}+V_{x}c+\hat{s}\cdot(\hat{\nabla}\times\hat{A}_{s})]\psi_{i}=\varepsilon_{i}\psi_{i}$$

(58)

Here, $A_{s}$ is the sum of external vector potential ($A_{ext}$ and the one due to exchange-correlation interaction ($A_{xc}$. It is important to distinguish between the latter and $V_{xc}$. They are defined, as a function of electron density ($\rho$) and current density ($j$) in the following way:

$$V_{xc}=\frac{\delta E_{xc}(\rho,j)}{\delta\rho(\vec{r})},\hskip 14.22636ptA_{xc}=\frac{\delta E_{xc}(\rho,j)}{\delta j(\vec{r})}.$$

(59)

The aforementioned Kohn-Sham system is solved iteratively using the uncontracted aug-cc-p-CVTZ basis set(citation needed). The authors use the cTPSS, cTPSS(h), B98 and the KT3(Keal-Tozer-3) functionals, available with the software package LONDON. Success of these MGGA functionals is decided by comparison with the methods established before (see Table3).

In case of the H₂ molecule, the KT3 functional is found to yield reasonably accurate number for the dissociation energy and bond length, but an absurd barrier in potential energy is observed. The barrier, in case of B98, is even bigger. Hence, these functionals, although, give reasonable magnitudes cannot be considered reliable to extract the physics of $H_{2}$. In the case of He₂, the B98 functional underbinds, but since it is empirically parametrized, the authors claim that its optimization could be enhanced.

Similar trends are observed in the case of NeHe and NeNe, and are summarised in Fig5. The authors report on the paramagnetic bonding using different functionals and levels of theory, as a measure to benchmark the performance of meta-GGA functionals. This is carried out for $H_{2}$, $He_{2}$, HeNe cluster and NeNe dimer. The paramagnetic bonding is explored by computing the variation in bond lengths. The authors also correctly predict perpendicular paramagnetic bonding (by examination of bond lengths and electron density plots).

Among the more recent calculation methods, Equation of Motion CCSD (EOM-CCSD) was introduced by Hampe et al(9), in 2017. In this method, the authors assume a general excited state wavefunction generated by the action of the action of a general excitation operator $\hat{R}$, which assumes the same form as that of the cluster operator, defined in Eq.55.

$$\hat{R}=\hat{R}_{1}+\hat{R}_{2}+...=\sum_{n=1}^{N}(\frac{1}{n!})^{2}\sum_{i,j,k...,a,b,c...}r_{ij}^{ab}\hat{a}_{a}^{\dagger}\hat{a}_{b}\hat{a}_{i}^{\dagger}\hat{a}_{j}$$

(60)

The energy is also computed in the same way,

$$e^{-\hat{T}}(\hat{H}-E_{cc})e^{\hat{T}}\hat{R}|\phi_{0}\rangle=\Delta E_{exc}\hat{R}|\phi_{0}\rangle$$

(61)

Here, $\Delta E_{exc}$ is the difference between energy of the excited and reference coupled cluster state.
The authors discuss two kinds of excitations in their article; electronic and spin-flip. In case of the former, the excitation occurs in such a way that the total spin of the system is conserved. However, in case of spin-flip excitation, the total spin of the system may not be conserved. Hence, for an overall spin-less system, a spin-flip excitation can give rise to an overall spin 1.

The authors consider C, $H_{2}$ and $CH^{+}$ ion in magnetic fields upto 1a.u. For carbon, they consider $1s^{2}2s2p_{-1}2p_{0}2p_{+1}$ as the reference state, while it is found that the excited states have the configuration $1s^{2}2s^{2}2p_{-1}2p_{0}(^{3}\Pi_{g})$, $1s^{2}2s2p_{-1}^{2}2p_{0}2(^{3}\Delta_{u})$, $1s^{2}2s2p_{-1}2p_{0}3d_{-2}(^{5}\Phi_{g})$. The authors compute the magnetic field at which the state crossovers happen and compare with Ivanov et al(15). Their results are tabulated in Table4.

From the above observations, it is clear that electron-electron correlation also influences crossover magnetic field strengths. However, in the case of $H_{2}$ molecule, the authors find no significant difference between the FCI and EOM-CCSD calculations. For the $CH^{+}$ ion, a similar result is obtained where ${}^{1}\Delta$ is the lowest lying state in high magnetic field, parallel to the molecule.
Therefore, not only does this method of calculation consider correlation effects to a reasonable accuracy, but it is also ovrecomes the limitations of a FCI calculations, of applicability to large systems.