Quantizing relativistic one- and two-electron operators
In chemistry, the energy ranges (typically a few eV) are not high enough that we should consider pair creation and annihilation processes (typically around 1000 MeV). Therefore, we should be able to construct an easier, pragmatic theory rather than a full QED treatment. In this quasi-relativistic (i.e. a `practical relativistic’) theory, we will use non-quantized electromagnetic fields and the Dirac correspondence principle. In this section we will introduce relativistic fermionic field operators and quantize the one- and two-electron operators.
Relativistic fermionic field operators
Starting from the time-independent Dirac equation, it first sight, this one-electron equation doesn’t seem so special, it only seems more complicated because a system of differential equations has to be solved instead of scalar differential equation. We recall that in the case of a free relativistic electron Reiher and Wolf (2015), its energy is doubly degenerate: \[\begin{equation} \require{physics} E(\vb{p}) = \pm \sqrt{m_e^2 c^4 + c^2 \norm{\vb{p}}^2} % \thinspace , \end{equation}\] whose negative-energy solutions (not to be confused by negative energies that we use in standard quantum chemistry - these are not on an absolute scale) may seem unphysical. It is indeed this peculiarity that needs some special consideration.
The existence of the negative-energy continuum means that an electron in a positive-energy state could spontaneously dissolve into the negative-energy continuum with the emission of a photon. Any positive-energy electron would therefore be unstable. Dirac’s solution to this issue was that these negative-energy states should be fully occupied. Let us then build up a second-quantized theory of these relativistic electrons, in a manner that is similar to how QFT would treat the problem.
We start by introducing the second-quantized operators \(\{ \hat{d}^\dagger_{\mathcal{P}} \}\), which are related to the general complete set of one-electron states \(\{ \phi_{\mathcal{P}}(\vb{r}) \}\). Then, we partition this set into \(\{ \phi^+_{\mathcal{P}}(\vb{r}) \}\), which correspond to a positive energy and have related creation operators \(\{ \hat{a}^\dagger_{\mathcal{P}} \}\), and the set of one-electron states with a negative energy \(\{ \phi^-_{\mathcal{P}}(\vb{r}) \}\) whose creation operators are \(\{ \hat{b}^\dagger_{\mathcal{P}} \}\). Having introduced this partitioning, we can write the field operator as:
\[\begin{align}
\hat{\psi}^\dagger(\vb{r}) &= %
\sum_{\mathcal{P}}^\text{all} %
\phi^\dagger_{\mathcal{P}}(\vb{r}) %
\hat{d}^\dagger_{\mathcal{P}} %
\label{eq:general_relativistic_field_operator} \\
&= %
\sum_{\mathcal{P}}^+ %
\phi^{+, \dagger}_{\mathcal{P}}(\vb{r}) %
\hat{a}^\dagger_{\mathcal{P}} %
+ \sum_{\mathcal{P}}^- %
\phi^{-,\dagger}_{\mathcal{P}}(\vb{r}) %
\hat{b}^\dagger_{\mathcal{P}} %
\thinspace .
\end{align}\]
Dirac’s interpretation requires the relativistic vacuum to be fully filled with respect to the negative-energy states, so we can write down the vacuum state as:
\[\begin{equation}
\ket{\text{vac}} = %
\prod_{\mathcal{P}}^- \hat{b}^\dagger_{\mathcal{P}} %
\ket{\emptyset} %
\thinspace ,
\end{equation}\]
which is also referred to as the Dirac sea'. If one of these electrons being in such a negative-energy state would be annihilated, there is ahole’ that behaves as an electron with positive charge: a positron. We will therefore reinterpret the annihilation (removal) of an electron with \(E<0\) as the creation of a positron (i.e. a particle with the same mass as an electron, but opposite charge) with \(E>0\). In our second quantized formalism, we will represent this reinterpretation as follows:
\[\begin{align}
\hat{b}_{\mathcal{P}} \ket{\text{vac}} &\neq 0 \\
&= \hat{c}^\dagger_{\mathcal{P}} \ket{\text{vac}} %
\thinspace ,
\end{align}\]
in which we have introduced the positron-related elementary second-quantization operators \(\{ \hat{c}^\dagger_{\mathcal{P}} \}\). Note that the set of bispinors that these operators refer to is still \(\{ \phi^-_{\mathcal{P}}(\vb{r}) \}\). Since all the previously introduced elementary second quantized operators refer to fermions, the usual anticommutator relations are still fulfilled:
\[\begin{align}
& \comm{ %
\hat{a}_{\mathcal{P}}
}{ %
\hat{a}^\dagger_{\mathcal{Q}} %
}_+ = %
\comm{ %
\hat{c}_{\mathcal{P}}
}{ %
\hat{c}^\dagger_{\mathcal{Q}} %
}_+ = %
\comm{ %
\hat{a}_{\mathcal{P}}
}{ %
\hat{c}^\dagger_{\mathcal{Q}} %
}_+ = %
\delta_{\mathcal{P} \mathcal{Q}} \\
%
& \comm{ %
\hat{a}^\dagger_{\mathcal{P}}
}{ %
\hat{a}^\dagger_{\mathcal{Q}} %
}_+ = %
\comm{ %
\hat{c}^\dagger_{\mathcal{P}}
}{ %
\hat{c}^\dagger_{\mathcal{Q}} %
}_+ = %
\comm{ %
\hat{a}^\dagger_{\mathcal{P}}
}{ %
\hat{c}^\dagger_{\mathcal{Q}} %
}_+ %
= 0 %
\thinspace .
\end{align}\]
We also should note that even though the vacuum is fully filled with respect to the negative-energy electrons, the vacuum is empty with respect to positive-energy positrons. Finally, this means that the field creation operator \(\hat{\psi}^\dagger(\vb{r})\) will be written as
\[\begin{equation}
\hat{\psi}^\dagger(\vb{r}) %
= \sum_{\mathcal{P}}^+ %
\phi^{+,\dagger}_{\mathcal{P}}(\vb{r}) %
\hat{a}^\dagger_{\mathcal{P}} %
+ \sum_{\mathcal{P}}^- %
\phi^{-,\dagger}_{\mathcal{P}}(\vb{r}) %
\hat{c}_{\mathcal{P}} %
\thinspace .
\end{equation}\]
Since we have constructed an appropriate field creation operator, we can start quantizing one- and two-electron operators, just like in the non-relativistic case. We will explore this quantization in the subsequent chapters.
Quantizing the Dirac Hamiltonian
Using the quantization procedure that we have used before, the Dirac Hamiltonian is quantized as follows: \[\begin{align} \hat{h}^{\text{D}} % &= \int \dd{\vb{r}} % \hat{\psi}^\dagger(\vb{r}) % \thinspace h^{c, \text{D}} \thinspace % \hat{\psi}(\vb{r}) \\ &= \sum^+_{\mathcal{P} \mathcal{Q}} % h^{\text{D}}_{\mathcal{P}+, \mathcal{Q}+} % \hat{a}^\dagger_{\mathcal{P}} % \hat{a}_{\mathcal{Q}} % + \sum^-_{\mathcal{P}} \sum^+_{\mathcal{Q}} % h^{\text{D}}_{\mathcal{P}-, \mathcal{Q}+} % \hat{c}_{\mathcal{P}} % \hat{a}_{\mathcal{Q}} \notag \\ & \hspace{12pt} + \sum^+_{\mathcal{P}} \sum^-_{\mathcal{Q}} % h^{\text{D}}_{\mathcal{P}+, \mathcal{Q}-} % \hat{a}^\dagger_{\mathcal{P}} % \hat{c}^\dagger_{\mathcal{Q}} % + \sum^-_{\mathcal{P} \mathcal{Q}} % h^{\text{D}}_{\mathcal{P}-, \mathcal{Q}-} % \hat{c}_{\mathcal{P}} \hat{c}^\dagger_{\mathcal{Q}} % \thinspace , \end{align}\] in which we have introduced the Dirac Hamiltonian matrix elements \[\begin{equation} h^{\text{D}}_{\mathcal{P}+, \mathcal{Q}+} % = \int \dd{\vb{r}} % \phi_{\mathcal{P}}^{+, \dagger}(\vb{r}) % \thinspace h^{c, \text{D}} \thinspace % \phi_{\mathcal{Q}}^{+}(\vb{r}) % \thinspace , \end{equation}\] and analogously for the other types. Note that, similarly to the non-relativistic matrix elements, these are also just scalar numbers.
The quantized Dirac Hamiltonian is quite peculiar, as it contains both electron-positron pair creation and annihilation operators, and in the last term the order of the positron operators is unusual. From these observations follows that that the quantized Dirac Hamiltonian does not even preserve particle number. Furthermore, it has an infinitely negative vacuum expectation value: \[\begin{equation} \ev{ % \hat{h}^{\text{D}} % }{\text{vac}} % = \sum_{\mathcal{P}}^- % h^{\text{D}}_{\mathcal{P}-, \mathcal{P}-} % \thinspace , \end{equation}\] which is something we should really take care of. From coupled cluster theory, we know that a normal ordered operator is nothing but the operator itself minus its vacuum expectation value, whatever the definition of the vacuum is. We then introduce the normal-ordered Dirac Hamiltonian as: \[\begin{align} \normord{ \hat{h}^{\text{D}} } % &= \hat{h}^{\text{D}} % - \ev{ % \hat{h}^{\text{D}} % }{\text{vac}} \\ &= \sum^+_{\mathcal{P} \mathcal{Q}} % h^{\text{D}}_{\mathcal{P}+, \mathcal{Q}+} % \hat{a}^\dagger_{\mathcal{P}} % \hat{a}_{\mathcal{Q}} % + \sum^-_{\mathcal{P}} \sum^+_{\mathcal{Q}} % h^{\text{D}}_{\mathcal{P}-, \mathcal{Q}+} % \hat{c}_{\mathcal{P}} % \hat{a}_{\mathcal{Q}} \notag \\ & \hspace{12pt} + \sum^+_{\mathcal{P}} \sum^-_{\mathcal{Q}} % h^{\text{D}}_{\mathcal{P}+, \mathcal{Q}-} % \hat{a}^\dagger_{\mathcal{P}} % \hat{c}^\dagger_{\mathcal{Q}} % - \sum^-_{\mathcal{P} \mathcal{Q}} % h^{\text{D}}_{\mathcal{P}-, \mathcal{Q}-} % \hat{c}^\dagger_{\mathcal{Q}} % \hat{c}_{\mathcal{P}} % \thinspace , \end{align}\] where we have used that the first three terms were already in normal ordered form, and we have performed the normal ordering on the last term by incorporating the minus sign, consistent with the fermion anticommutator relations. By construction, the vacuum expectation value of this normal-orderd Hamiltonian is zero: \[\begin{equation} \ev{ % \normord{\hat{h}^{\text{D}}} % }{\text{vac}} % = 0 \end{equation}\] and the energy of a positron state \[\begin{equation} \ev{ % \hat{c}_\mathcal{P} % \normord{\hat{h}^{\text{D}}} % \hat{c}^\dagger_\mathcal{P} % }{\text{vac}} % = - h^{\text{D}}_{\mathcal{P}-, \mathcal{P}-} \end{equation}\] is now positive, which is consistent with the fact that the negative-energy electron states have now become positive-energy positron states.
Since the QED-like quantized Dirac Hamiltonian is vacuum-dependent, we must always be vigilant of clearly defining the vacuum state that is used in derivations. Since a general rotation in bispinor space changes the vacuum state, we will follow the floating vacuum approach by Mittleman (1981): the vacuum is always defined by the current set of one-particle states.
Quantizing the two-electron interactions
Following Reiher and Wolf (2015), we find an expression for the electron-electron interaction, correct up to order \(c^{-2}\): \[\begin{equation} V^c(\vb{r}_1, \vb{r}_2) = \frac{q_1 q_2}{r_{12}} - \frac{q_1 q_2}{2} \qty( \frac{\boldsymbol{\alpha}_1 \boldsymbol{\alpha}_2}{r_{12}} + \frac{(\vb{r}_{12} \vdot \boldsymbol{\alpha}_1) (\vb{r}_{12} \vdot \boldsymbol{\alpha}_2)}{r_{12}^3} ) \thinspace , \end{equation}\] which is called the Coulomb-Breit interaction. The naming of this interaction is clear: the first term is the general (instantaneous) Coulomb interaction and second term is called the Breit term, which is in its turn decomposed into a current-current interaction (the Gaunt term) and an extra gauge term.
Let us first digress a little on the Brown-Ravenhall disease, which is also called continuum dissolution. Consider two non-interacting bound-state particles. The combination of these two states is degenerate with an infinite number of combinations of states in which one particle excites to the positive continuum while the other de-excites into the negative continuum, with a negative energy such that the sum of the energies is equal to the sum of the bound-state energies. If any two-electron interaction (Coulomb-only or full Breit) would then be turned on, the system could dissolve into two continua. In order to avoid any problems related to this collapse, we should employ the reinterpretation scheme that was discussed in the previous section (G. Dyall and Faegri 2007). Regardless of which two-electron interactions that are included, we should use the general field operator \(\eqref{eq:general_relativistic_field_operator}\) to quantize the two-electron interaction: \[\begin{align} \hat{g} &= % \frac{1}{2} \int \dd{\vb{r}_1} \int \dd{\vb{r}_2} % \hat{\psi}^{\dagger}(\vb{r}_1) % \hat{\psi}^{\dagger}(\vb{r}_2) % V^c(\vb{r}_1, \vb{r}_2) % \hat{\psi}(\vb{r}_2) % \hat{\psi}(\vb{r}_1) \\ &= \frac{1}{2} % \sum_{ % \mathcal{P} \mathcal{Q} % \mathcal{R} \mathcal{S} % }^\text{all} % g_{\mathcal{P} \mathcal{Q} \mathcal{R} \mathcal{S}} % \hat{d}^\dagger_{\mathcal{P}} % \hat{d}^\dagger_{\mathcal{R}} % \hat{d}_{\mathcal{S}} % \hat{d}_{\mathcal{Q}} % \thinspace , \end{align}\] in which the two-electron matrix elements/integrals are defined as \[\begin{equation} g_{\mathcal{P} \mathcal{Q} \mathcal{R} \mathcal{S}} % = \int \dd{\vb{r}_1} \int \dd{\vb{r}_2} % \phi^\dagger_{\mathcal{P}}(\vb{r}_1) % \phi^\dagger_{\mathcal{R}}(\vb{r}_2) % V^c(\vb{r}_1, \vb{r}_2) % \phi_{\mathcal{S}}(\vb{r}_2) % \phi_{\mathcal{Q}}(\vb{r}_1) % \thinspace . \end{equation}\]
Let us now cover the quantization of the Coulomb, Coulomb-Gaunt and Coulomb-Breit interaction operators. If we only allow for the instantaneous Coulomb interaction to take place, we should use the familiar Coulomb operator, which is given in four-component and two-component form as \[\begin{align} V^c_{\text{C}}(\vb{r}_1, \vb{r}_2) % &= \frac{q_1 q_2}{r_{12}} \thinspace \vb{I}_4 \\ &= \frac{q_1 q_2}{r_{12}} \begin{pmatrix} \vb{I}_2 & 0 \\ 0 & \vb{I}_2 \end{pmatrix} \end{align}\] If we furthermore want to include the instantaneous magnetic interaction operator due to the moving electrons, i.e. the interaction operator due to the unretarded vector potential, we should add the Gaunt operator. In four- and two-component form, it can be written as: \[\begin{align} G_0^c(\vb{r}_1,\vb{r}_2) % &= - \frac{q_1 q_2}{r_{12}} % \boldsymbol{\alpha}_1 \boldsymbol{\alpha}_2 \\ &= - \frac{q_1 q_2}{r_{12}} \begin{pmatrix} 0 % & \boldsymbol{\sigma}_1 \boldsymbol{\sigma}_2 \\ \boldsymbol{\sigma}_1 \boldsymbol{\sigma}_2 & 0 \end{pmatrix} \thinspace . \end{align}\] The product \(\boldsymbol{\alpha}_1 \boldsymbol{\alpha}_2\) that appears in the Coulomb-Gaunt operator has to be seen as a scalar product of matrix vectors: \[\begin{equation} \boldsymbol{\alpha}_1 \boldsymbol{\alpha}_2 % = \alpha_{1,x} \alpha_{2,x} % + \alpha_{1,y} \alpha_{2,y} % + \alpha_{1,z} \alpha_{2,z} \end{equation}\] that act on the corresponding orbitals \(\phi_{\mathcal{P}}(\vb{r}_1)\) and \(\phi_{\mathcal{Q}}(\vb{r}_2)\) separately. We should emphasize that the Coulomb-Gaunt operator is the true unretarded interaction operator between the two electrons.
If, on top of the unretarded interaction Coulomb-Gaunt operator, we want to add a first-order retardation effect, we should use the Breit operator. In four- and two-component form, it reads \[\begin{align} B^c_0(\vb{r}_1,\vb{r}_2) % &= - \frac{q_1 q_2}{2} \qty( % \frac{\boldsymbol{\alpha}_1 % \boldsymbol{\alpha}_2}{r_{12}} % + \frac{ % (\vb{r}_{12} \vdot \boldsymbol{\alpha}_1) % (\vb{r}_{12} \vdot \boldsymbol{\alpha}_2) % }{r_{12}^3} % ) \\ &= - \frac{q_1 q_2}{2} % \begin{pmatrix} 0 % & \dfrac{ % \boldsymbol{\sigma}_1 % \boldsymbol{\sigma}_2 % }{r_{12}} % + \dfrac{ % (\vb{r}_{12} \vdot \boldsymbol{\sigma}_1) % (\vb{r}_{12} \vdot \boldsymbol{\sigma}_2) % }{r_{12}^3} \\ \dfrac{ % \boldsymbol{\sigma}_1 % \boldsymbol{\sigma}_2 % }{r_{12}} % + \dfrac{ % (\vb{r}_{12} \vdot \boldsymbol{\sigma}_1) % (\vb{r}_{12} \vdot \boldsymbol{\sigma}_2) % }{r_{12}^3} % & 0 \end{pmatrix} \end{align}\] and since it includes the Gaunt operator, it includes both instantaneous and retarded interactions. The true retardation operator (which includes a certain degree of arbitrariness due to the choice of gauge function in the derivation of the Darwin interaction (Reiher and Wolf 2015)) is then given by \[\begin{equation} B^c_\text{ret}(\vb{r}_1,\vb{r}_2) % = \frac{q_1 q_2}{2} \qty( % \frac{ % \boldsymbol{\alpha}_1 % \boldsymbol{\alpha}_2 % }{r_{12}} % - \frac{ % (\vb{r}_{12} \vdot \boldsymbol{\alpha}_1) % (\vb{r}_{12} \vdot \boldsymbol{\alpha}_2) % }{r_{12}^3} % ) % \thinspace . \end{equation}\]
Towards a second-quantized four-component quasi-relativistic many-particle Hamiltonian
Let’s explore these integrals a little further by starting with the overlap between two \(4\)-spinors: \[\begin{equation} S_{\mathcal{P} \mathcal{Q}} = \braket{\phi_{\mathcal{P}}}{\phi_{\mathcal{Q}}} \thinspace , \end{equation}\] in which the expression on the right-hand side is just short-hand notation for the integral \[\begin{equation} \braket{\phi_{\mathcal{P}}}{\phi_{\mathcal{Q}}} \equiv \int \dd{\vb{r}} \phi^\dagger_{\mathcal{P}}(\vb{r}) \phi_{\mathcal{Q}}(\vb{r}) \thinspace , \end{equation}\] in which \(\phi^\dagger_{\mathcal{P}}(\vb{r}) \phi_{\mathcal{Q}}(\vb{r})\) is a Dirac distribution, like in equation \(\eqref{eq:Dirac_distribution}\). Therefore, we can write the overlap integral as \[\begin{align} S_{\mathcal{P} \mathcal{Q}} &= \braket{\phi_{\mathcal{P}}}{\phi_{\mathcal{Q}}} \\ &= \braket*{\phi^{\text{L}}_{\mathcal{P}}}{\phi^{\text{L}}_{\mathcal{Q}}} + \braket*{\phi^{\text{S}}_{\mathcal{P}}}{\phi^{\text{S}}_{\mathcal{Q}}} \\ &= \sum_u^4 \braket*{\phi^{(u)}_{\mathcal{P}}}{\phi^{(u)}_{\mathcal{Q}}} \thinspace . \end{align}\] The one-electron integrals can be written in a similar way: \[\begin{align} h^\text{D}_{\mathcal{P} \mathcal{Q}} &\equiv \matrixel{\phi_{\mathcal{P}}}{h^{c,\text{D}}}{\phi_{\mathcal{Q}}} \\ &= \int \dd{\vb{r}} \phi^\dagger_{\mathcal{P}}(\vb{r}) h^{c,\text{D}}(\vb{r}) \phi_{\mathcal{Q}}(\vb{r}) \label{eq:h_cD_integral} \thinspace . \end{align}\] We can work this out by first letting the \(1\)-particle Dirac Hamiltonian operate on its argument on the right such that we find in terms of the \(2\)-spinors: \[\begin{equation} \begin{split} h^\text{D}_{\mathcal{P} \mathcal{Q}} = &\matrixel{\phi^{\text{L}}_{\mathcal{P}}}{V^c_{\text{nuc}} + q_e \phi^c_{\text{ext}}}{\phi^{\text{L}}_{\mathcal{Q}}} + c \matrixel{\phi^{\text{L}}_{\mathcal{P}}}{\boldsymbol{\sigma} \vdot \boldsymbol{\pi}^c}{\phi^{\text{S}}_{\mathcal{Q}}} + c \matrixel{\phi^{\text{S}}_{\mathcal{P}}}{\boldsymbol{\sigma} \vdot \boldsymbol{\pi}^c}{\phi^{\text{L}}_{\mathcal{Q}}} \\ &+ \matrixel{\phi^{\text{S}}_{\mathcal{P}}}{V^c_{\text{nuc}} + q_e \phi^c_{\text{ext}} - 2 m_e c^2}{\phi^{\text{S}}_{\mathcal{Q}}} \thinspace , \end{split} \end{equation}\] and in full \(4\)-component form, we have \[\begin{align} h^\text{D}_{\mathcal{P} \mathcal{Q}} = &\sum_u^4 \matrixel*{\phi^{(u)}_{\mathcal{P}}}{V^c_{\text{nuc}}}{\phi^{(u)}_{\mathcal{Q}}} + q_e \sum_u^4 \matrixel*{\phi^{(u)}_{\mathcal{P}}}{\phi^c_{\text{ext}}}{\phi^{(u)}_{\mathcal{Q}}} + \sum_u^{3,4} \braket*{\phi^{(u)}_{\mathcal{P}}}{\phi^{(u)}_{\mathcal{Q}}} \\ &+ \matrixel*{\phi^{(1)}_{\mathcal{P}}}{\pi^c_z}{\phi^{(3)}_{\mathcal{Q}}} + \matrixel*{\phi^{(3)}_{\mathcal{P}}}{\pi^c_z}{\phi^{(1)}_{\mathcal{Q}}} + \matrixel*{\phi^{(1)}_{\mathcal{P}}}{\pi^c_x - i \pi^c_y}{\phi^{(4)}_{\mathcal{Q}}} + \matrixel*{\phi^{(3)}_{\mathcal{P}}}{\pi^c_x - i \pi^c_y}{\phi^{(2)}_{\mathcal{Q}}} \\ &+ \matrixel*{\phi^{(2)}_{\mathcal{P}}}{\pi^c_x + i \pi^c_y}{\phi^{(3)}_{\mathcal{Q}}} + \matrixel*{\phi^{(4)}_{\mathcal{P}}}{\pi^c_x + i \pi^c_y}{\phi^{(1)}_{\mathcal{Q}}} - \matrixel*{\phi^{(2)}_{\mathcal{P}}}{\pi^c_z}{\phi^{(4)}_{\mathcal{Q}}} - \matrixel*{\phi^{(4)}_{\mathcal{P}}}{\pi^c_z}{\phi^{(2)}_{\mathcal{Q}}} \thinspace . \end{align}\]
Let us now turn our attention to the 2-electron interaction integrals, for which we sometimes use physicist’s notation \[\begin{equation} g_{\mathcal{P} \mathcal{Q} \mathcal{R} \mathcal{S}} \equiv \matrixel{\phi_{\mathcal{P}} \phi_{\mathcal{R}}}{g^c}{\phi_{\mathcal{Q}} \phi_{\mathcal{S}}} \end{equation}\] but mostly use chemist’s notation: \[\begin{equation} g_{\mathcal{P} \mathcal{Q} \mathcal{R} \mathcal{S}} \equiv (\phi_{\mathcal{P}} \phi_{\mathcal{Q}} | g^c | \phi_{\mathcal{R}} \phi_{\mathcal{S}}) \thinspace . \end{equation}\] In a slightly confusing notation, we could then write \[\begin{equation} g_{\mathcal{P} \mathcal{Q} \mathcal{R} \mathcal{S}} = \int \int \dd{\vb{r}_1} \dd{\vb{r}_2} \phi^\dagger_{\mathcal{P}}(\vb{r}_1) \phi^\dagger_{\mathcal{R}}(\vb{r}_2) g^c(\vb{r}_1, \vb{r}_2) \phi_{\mathcal{Q}}(\vb{r}_1) \phi_{\mathcal{S}}(\vb{r}_2) \thinspace , \end{equation}\] but we should realize that \(\phi_{\mathcal{Q}}(\vb{r}_1) \phi_{\mathcal{S}}(\vb{r}_2)\) does not imply a matrix product (since the dimensions are not compatible). Instead, we will analyze the expression for each of the electron interactions we have previously discussed.Avoiding variational collapse
We often rely on the variational principle to determine the minimum of the total electronic energy with respect to a set of parameters. However, in relativistic theories, the unboundedness of the one-electron Dirac operator comes into play, i.e. we have a negative-energy continuum. To solve this problem, formally, the minimax principle is used: determine the minimum of the electronic energy with respect to the large component of the spinor, while guaranteeing a maximum of the energy with respect to the small component.
In practice, however, it seems that positronic negative-energy states are omitted from the density matrix, even though they are optimized in every iteration step. The proper solution algorithm then specifically targets one-electron bound states.
Current densities for many electrons
The current density for a many-electron system can be derived from Ehrenfest’s theorem, leading to the continuity equation \[\begin{equation} \pdv{t} \ev{\rho_{\vb{r}}} + \boldsymbol{\nabla} \vdot{ \vb{j}(\vb{r}) } = 0 \thinspace . \end{equation}\] The current density can then be calculated by \[\begin{equation} \vb{j}(\vb{r}) = c \sum_i^N \ev{\boldsymbol{\alpha}_i \delta(\vb{r}_i - \vb{r})}{\Psi} \end{equation}\] for the Dirac-Coulomb, the Dirac-Coulomb-Gaunt and the Dirac-Coulomb-Breit Hamiltonian since multiplicative terms vanish in the commutator.
Self-consistent field equations
Given a normalized wave function \(\ket{\Psi}\), we can write its electronic energy as \[\begin{align} E &= \ev{\hat{\mathcal{H}}_{\text{QR4C}}}{\Psi} \\ &= \sum_{\mathcal{P} \mathcal{Q}}^K h^{\text{D}}_{\mathcal{P} \mathcal{Q}} D_{\mathcal{P} \mathcal{Q}} + \frac{1}{2} \sum_{\mathcal{P} \mathcal{Q} \mathcal{R} \mathcal{S}}^K g_{\mathcal{P} \mathcal{Q} \mathcal{R} \mathcal{S}} d_{\mathcal{P} \mathcal{Q} \mathcal{R} \mathcal{S}} \thinspace , \end{align}\] in which \(D_{\mathcal{P} \mathcal{Q}}\) is the one-electron density matrix: \[\begin{equation} D_{\mathcal{P} \mathcal{Q}} = \ev{\hat{a}^\dagger_{\mathcal{P}} \hat{a}_{\mathcal{Q}}}{\Psi} \end{equation}\] and \(d_{\mathcal{P} \mathcal{Q} \mathcal{R} \mathcal{S}}\) is the two-electron density matrix: \[\begin{equation} d_{\mathcal{P} \mathcal{Q} \mathcal{R} \mathcal{S}} = \ev{\hat{a}^\dagger_{\mathcal{P}} \hat{a}^\dagger_{\mathcal{R}} \hat{a}_{\mathcal{S}} \hat{a}_{\mathcal{Q}}}{\Psi} \thinspace . \end{equation}\]
If we have a set \(\set{\phi_\mathcal{P}(\vb{r})}\) of \(K\) orthonormal 4-spinors available, we would like to minimize this electronic energy by varying these spinors under the constraint that the spinors remain orthonormal. Therefore, we set up the Lagrangian \[\begin{equation} \mathcal{L}\qty[\set{\phi_\mathcal{P}}, \set{\epsilon_{\mathcal{P} \mathcal{Q}}}] = \ev{\hat{\mathcal{H}}_{\text{QR4C}}}{\Psi} - \sum_{\mathcal{P} \mathcal{Q}}^K \epsilon_{\mathcal{P} \mathcal{Q}} ( \braket{\phi_{\mathcal{P}}}{\phi_{\mathcal{Q}}} - \delta_{\mathcal{P} \mathcal{Q}}) \thinspace , \end{equation}\] and require \(\delta{\mathcal{L}} = 0\). By doing the variational calculus, we find the self-consistent field (i.e. the Fock) equations: \[\begin{equation} \sum_{\mathcal{Q}}^K f^c_{\mathcal{P} \mathcal{Q}}(\vb{r}) \phi_\mathcal{Q}(\vb{r}) = \sum_{\mathcal{Q}}^K \epsilon_{\mathcal{P} \mathcal{Q}} \phi_{\mathcal{Q}}(\vb{r}) \thinspace , \end{equation}\] in which \(f^c_{\mathcal{P} \mathcal{Q}}(\vb{r})\) is called the Fock operator: \[\begin{equation} f^c_{\mathcal{P} \mathcal{Q}}(\vb{r}) = D_{\mathcal{P} \mathcal{Q}} h^{c, \text{D}}(\vb{r}) + \sum_{\mathcal{R} \mathcal{S}}^K d_{\mathcal{P} \mathcal{Q} \mathcal{R} \mathcal{S}} \int \dd{\vb{r}'} \phi^\dagger_\mathcal{R}(\vb{r}') g^c(\vb{r}, \vb{r}') \phi_{\mathcal{S}}(\vb{r}') \thinspace . \end{equation}\]