Title: ETH_Sem5 / AdvQC / 11_Decoupling_Negative-Energy States_FW_DKH

Q: Which ``disadvantages'' do the small components of the spinor have? What is the main idea to circumvent these problems?
A: Small components of the spinor are the reason for:
unboundedness of Dirac Hamiltonian \(\rightarrow\) variational instability
increased computational cost due to the increased basis-set size.
Idea: Block-diagonalize Dirac Hamiltonian by unitary transformation \(U\),
\[f_{bd} = U\,f\,U^\dagger = \begin{pmatrix}f_{+} & 0 \\ 0 & f_{-}\end{pmatrix}\]
i.e., 2 new operators take care of electronic and positronic states separately. One is only interested in \(f_{+}\), which contains the energies of the bound states.
Then, the positive- and negative-energy states are completely decoupled, variational collapse can be avoided and a basis set is only needed for the upper two large components of each atomic/molecular spinor.
Similarly to the Hamiltonian, the spinor is also reduced by the unitary transformation and has only one non-vanishing 2-spinor component
\[\begin{aligned} & \text{class-I solution: }\ \tilde{\psi}_{+} = U \psi = \begin{pmatrix}\tilde{\psi}^L \\ 0 \end{pmatrix} \\ & \text { and } \quad\quad \text{class-II solution: }\ \tilde{\psi}_{-} = \begin{pmatrix} 0 \\ \tilde{\psi}^S \end{pmatrix} \end{aligned}\]

Q: Which main techniques exist to decouple the bound states from the positronic ones?
A: Transformation Techniques: \[f_{bd} = U\,f\,U^\dagger = \begin{pmatrix}f_{+} & 0 \\ 0 & f_{-}\end{pmatrix}\]
- FW: exact and simple for a free particle (fpFW), however in a potential, no exact decoupling is possible \(\rightarrow\) ``Solution'': sequence of unitary transformations to damp off-diagonal elements step by step in order of \(1/c\) (or \(1/m\) in different unit system). BUT: does not converge \(\rightarrow\) only for perturbation theory and not variational approach
- DKH: same as FW but order according to powers in \(V\) \(\rightarrow\) works! Standard approach is DKH2
Elimination Techniques: \[\Psi^S = f(\Psi^L)\] e.g., kinetic balance. Insert this expression for \(\Psi^S\) in equation in split notation to solve for \(\Psi^L\) \(\rightarrow\) only one basis is needed!
- Pauli or Breit-Pauli: also in order of \(1/c\) \(\rightarrow\) only for perturbation theory and not variational approach. Historically important but not used anymore.
- Regular approximations (ZORA): standard approach to be used

Q: What is the formal solution to find the unitary transformation which decouples the positronic states?
A: Formalize relation between large and small component by \(X\)-operator,
\[\psi^S = X\psi^L\]
It can be shown (see exercise) that \(X\) fulfills a nonlinear operator identity,
\[X = \frac{1}{2m_ec^2}\left\{c\bm{\sigma}\cdot\bm{p} - \comm{X}{V} - Xc\bm{\sigma}\cdot\bm{p}X\right\}\]
whose solution, however, is as complicated as solution of Dirac equation itself.
Next, it can be shown that the most general unitary transformation can be constructed in terms of the \(X\)-operator (see exercise)
\[U = U(X) = \begin{pmatrix} (1+X^\dagger X)^{-1/2} & (1+X^\dagger X)^{-1/2}X^\dagger \\ -\mathrm{e}^{\imag\varphi}(1+XX^\dagger)^{-1/2}X & \mathrm{e}^{\imag\varphi}(1+XX^\dagger)^{-1/2} \end{pmatrix}\]
and therefore the Hamiltonian for the electronic states reads
\[f_{+} = \frac{1}{\sqrt{1+X^\dagger X}}\Bigg\{ V + c\bm{\sigma}\cdot\bm{p}X + X^\dagger c\bm{\sigma}\cdot\bm{p} + X^\dagger(V-2m_ec^2)X\Bigg\}\frac{1}{\sqrt{1+X^\dagger X}}\]

Q: What is the Foldy-Wouthuysen transformation?
A: Historically, it was the first attempt to achieve the block-diagonalization of the Dirac Hamiltonian \(h^D\) by Foldy and Wouthuysen (1950).
They derived a closed-form expressions for both the unitary transformation and the decoupled Hamiltonian for the case of a free particle.
\[U_0 = A_p \begin{pmatrix} 1 & \bm{\sigma}\cdot\bm{P}_p \\ -\bm{\sigma}\cdot\bm{P}_p & 1\end{pmatrix}\]
\[\begin{aligned} A_p & \equiv \sqrt{\frac{E_p+m_ec^2}{2E_p}}, \\ \bm{P}_p & \equiv\frac{c\bm{p}}{E_p+m_ec^2}, \\ R_p & \equiv \bm{\alpha}\cdot\bm{P}_p = \frac{c\bm{\alpha}\cdot\bm{p}}{E_p+m_ec^2}\end{aligned}\]
are the kinematic factors. \(E_p\) is the positive square-root of the energy-momentum dispersion relation,
\[E_p = \sqrt{c^2p^2+m_e^2c^4}\]
which contains a squared momentum eigenvalue (not operator).
Clearly, this can be directly calculated in momentum space. For this reason, the theory is often (but unnecessarily) derived in momentum space.
All what is required for the evaluation of \(E_p\) and therefore of \(A_p\) and \(R_p\) is a one-electron basis, in which the matrix representation of \(\bm{p}^2\) is diagonal (then, the \(\bm{p}^2\) eigenvalues realized in this basis are found on the diagonal).

Q: Which trick can be used for the FW transformation in position space?
A: Trick for position-space formulations (Hess, 1986):
Diagonalize the non-relativistic kinetic energy matrix
\[ \bm{T} = \{T_{\mu\nu}\} = \left\{\biggl<{\lambda_\mu}\bigg|{\frac{\mom^2}{2m_e}}\bigg|{\lambda_\nu}\biggr>\right\}\]
which directly yields \(\mom^2\) after multiplication of its eigenvalues by \(2m_e\).
Here, \(\lambda_\mu\) are position-space one-electron functions (such as Guassian functions). The benefit is that the non-relativistic kinetic energy operator \(-\hbar^2 \nabla^2/(2m_e)\) contains the squared momentum operator \(\mom^2 = -\hbar^2 \nabla^2\). An operator in a position-space basis can then be transformed into one in the \(\mom^2\)-basis by transforming it with the eigenvectors of the kinetic energy matrix \(\bm{T}\) scaled by \(2m_e\). As a fortunate consequence, an explicit momentum-space representation of all operators is not necessary!
(Reiher. M. Sequential Decooupling of Negative-Energy States in Douglas-Kroll-Hess Theory. In Handbook of Relativistic Quantum Chemistry Liu, W.; Ed.; Springer: Berlin, Heidelberg, 2015)

Q: Is it possible to perform the free-particle Foldy-Wouthuysen transformation also in the presence of a scalar potential \(V\) of any form, i.e., \(f = h^D + V\)?
A: Yes, it can still be performed in closed form
\[f_1 = U_0\left(h^D+V\right)U_0^\dagger = \mathcal{E}_0 + \mathcal{E}_1 + \mathcal{O}_1\]
\[\begin{aligned}\text{with } \qquad &\mathcal{E}_0 = \beta E_p-m_ec^2, \\ & \mathcal{E}_1 = A_p\,V\,A_p + A_p\,R_p\,V\,R_p\,A_p, \\ & \mathcal{O}_1 = \beta A_p\,\comm{R_p}{V}\,A_p\end{aligned}\]
The subscripts of \(\mathcal{E}_0\), \(\mathcal{E}_1\) and \(\mathcal{O}_1\) denote the order in the scalar potential \(V\). The notation reflects the block-diagonal super-structure of the resulting operators: \(\mathcal{E}\) denotes a block-diagonal operator with \(LL\) and \(SS\) blocks, which is called even, while \(\mathcal{O}\) denotes an operator with entries on the two off-diagonal \(LS\) and \(SL\) blocks, which is called odd.
NB1: \(R_p\) is linear in \(\bm{\alpha}\) and each of its elements is odd. Hence, \(\mathcal{O}_1\) is odd.
NB2: \(A_p\,R_p\,V\,R_p\,A_p\) is bilinear in \(R_p\) and therefore even because the
product of two odd matrices is even!
NB3: the presence of a potential energy operator \(\bm{V}\) leads to an odd operator and exact decoupling is not achieved by the free-particle FW (fpFW)!

Q: What was the main idea of Foldy and Wouthuysen to resolve the problem that exact decoupling is not possible in the presence of a potential energy operator \(V\)?
A: Idea by F&W: eliminate these odd operators term by term through a sequence of unitary transformations:
\[ U = \cdots U_4U_3U_2U_1U_0 = \prod_{i=0}^\infty U_i \]
The idea is thus to damp the off-diagonal elements step by step. For exact unitary transformations \(U\), i.e., infinite power series, the energy eigenvalues of the transformed operator will be exactly equal to the original ones.
In actual applications the power series expansions of the Fock operator and therefore also of the matrices \(U_i\) may always be truncated after a finite number of terms to produce approximately block-diagonal, but efficient, one-electron operators.
The problem is that it does not converge! Therefore, it can only be used for a perturbation theory approach but not in a variational approach.

Q: But how to choose any of these unitary transformations in the sequence from Foldy Wouthuysen?
A: Most general unitary transformation \(U\) can be parametrized by antihermitean operator \(W\)
\[W^\dagger = W^{\star,T} = -W\]
Well-known example in quantum chemistry: exponential parametrization
\[\begin{aligned}U&=\exp(W) \\ \thus U^\dagger&=\exp(-W) \\\thus U\cdot U^\dagger&=\exp(W)\cdot\exp(-W)=1 \end{aligned}\]
Other examples for analytic forms for parametrizations of \(U\) in terms of \(W\) are
\[\begin{aligned}\text{Douglas&Kroll:}\, &U=\sqrt{1+W^2}+W \ ; \\ \text{Cayley-type:}\, &U=\frac{a+W}{a-W} \ ; \\ \text{McWeeny:}\, & U = \frac{1+W}{\sqrt{1-W^2}}\end{aligned}\]
Expansion of \(U\) in a Taylor series is natural in a term by term elimination of odd operators and is key in the derivation of explicit analytic expressions
for \(f_{+}\). We can cover all different choices of parametrized unitary transformations by making the general ansatz,
\[\begin{aligned}& U=a_0 \mathbf{1}+a_1 W+a_2 W^2+a_3 W^3+\cdots=a_0 \mathbf{1}+\sum_{k=1}^{\infty} a_k W^k \\ & \thus U^\dagger = a_0 \mathbf{1}+\sum_{k=1}^{\infty}(-1)^k a_k W^k \end{aligned}\]
In order to guarantee \(UU^\dagger = \bm{1}\), the coefficients \(a_k\) have to satisfy a set of constraints. Odd coefficients can be chosen arbitrarily; even coefficients are functions of the lower odd ones.

Q: What is the difference between Foldy-Wouthuysen (FW) and Douglas-Kroll-Hess (DKH) expansions? What is the problem with FW?
A: Every transformation method aiming at decoupling the Dirac Hamiltonian
\[\begin{aligned} h^D+V &=(\beta-1) m_e c^2+V+c \bm{\alpha} \cdot \bm{p} \\ &=\mathcal{E}_0^{(0)}+\mathcal{E}_1^{(0)}+\mathcal{O}_0^{(0)} \\ &=\mathcal{E}_{[-2]}^{(0)}+\mathcal{E}_{[0]}^{(0)}+\mathcal{O}_{[-1]}^{(0)} \end{aligned}\]
has necessarily to eliminate the odd term \(\mathcal{O}_{[-1]}^{(0)} = \mathcal{O}_{0}^{(0)}\) in the first step.
The two notations refer to the only two order parameters available, namely \(V\) and \(1/c\)!
Historically, \(1/c\) expansions, which allow one to approach the well-known nonrelativistic limit via \(c\rightarrow\infty\), were the only ones seriously considered for decades. Therefore, we strictly refer to \(\)1/c expansions as Foldy-Wouthuysen expansions. This is in accordance with the original intention of Foldy and Wouthuysen (1950).
Expansions in terms of the scalar potential \(V\) are called Douglas-Kroll-Hess (DKH) expansions.
All \(1/c\)-expansions are ill-defined and therefore only DKH expansions are possible for Dirac electrons in external fields!

Q: Explain the procedure of the FW expansion considering the first step!
A: Consider the first step of the FW transformation, which eliminates the first odd operator of order \(c^{-1}\). The most general ansatz for the initial free-particle transformation \(U_0\) then reads
\[U_{V=0}=U_0=U_0\left(W_{[0]}\right)=a_{0,0} \mathbf{1}+\sum_{k=1}^{\infty} a_{0, k} W_{[0]}^k\]
Transformation of the one-electron operator \(f_0=\left(h^D+V\right)\) now yields
\[\begin{aligned} f_1 &= U_0f_0U_0^{\dagger} = \left[a_{0,0} \mathbf{1}+\sum_{k=1}^{\infty} a_{0, k} W_{[0]}^k\right] \left(\mathcal{E}_{[-2]}+\mathcal{E}_{[0]}+{O}_{[-1]}^{(0)}\right) \left[a_{0,0} \mathbf{1}+\sum_{k=1}^{\infty} (-1)^k a_{0, k} W_{[0]}^k\right] \\ &= \mathcal{E}_{[-2]}+\mathcal{E}_{[0]}+\mathcal{E}_{[2]}+\underbrace{\mathcal{O}_{[-1]}^{(0)}+a_{0,0} a_{0,1}\left[W_{[0]}, \mathcal{E}_{[-2]}\right]}_{\text {odd }} + \sum_{k=2}^{\infty} \mathcal{E}_{[2 k]}^{(1)}+\sum_{k=1}^{\infty} \mathcal{O}_{[2 k-1]}^{(1)} \end{aligned}\]
Condition on \(W_{[0]}\) to eliminate lowest-order odd term in \(f_1\) then reads
\[\begin{aligned} \mathcal{O}_{[-1]}^{(1)} \equiv \mathcal{O}_{[-1]}^{(0)}+a_{0,0} a_{0,1}\left[W_{[0]}, \mathcal{E}_{[-2]}\right] \stackrel{!}{=} 0 \\ \thus W_{[0]}=\frac{a_{0,0}}{a_{0,1}} \beta \frac{c \bm{\alpha} \cdot \bm{p}}{2 m_e c^2}=\frac{a_{0,0}}{a_{0,1}} \beta \frac{\bm{\alpha} \cdot \bm{p}}{2 m_e c} \end{aligned}\]
Note: \(W_{[0]}\) is odd, antihermitean and of first order in \(1/c\)!
The leading terms of \(f_1\) are independent of the chosen parametrization of \(U_0\) \(\rightarrow\) independent of expansion coeffs. \(a_{0, k}\).

Q: Which operators are recovered by the first step of the FW transformation?
A: The explicit form of the lowest-order even terms is
\[\begin{aligned} \mathcal{E}_{[-2]} &=(\beta-1) m_e c^2 \\ \mathcal{E}_{[0]} &=V+\beta \frac{\bm{p}^2}{2 m_e} \\ \mathcal{E}_{[2]} &=-\beta \frac{\bm{p}^4}{8 m_e^3 c^2}-\frac{[\bm{\alpha} \cdot \bm{p},[\bm{\alpha} \cdot \bm{p}, V]]}{8 m_e^2 c^2}\end{aligned}\]
We therefore recover the (corrected) rest energy \(\mathcal{E}_{[-2]}\), the non-relativistic Schrödinger operator \(\mathcal{E}_{[0]}\), and the lowest-order relativistic corrections (see Pauli elimination and classical mechanics).
Note that these operators are still four-dimensional. The relevant part needed for the construction of the \(f_{+}\)operator is obtained when \(\beta\) is replaced by \(\mathbf{1}_2\) (and for \(f_{-}\) if \(\beta\) is replaced by \(-\mathbf{1}_2\)).

Q: Why is the Foldy-Wouthuysen expansion ill-defined?
A: Since the square root occurring in the relativistic energy-momentum relation \(E_p\) possesses branching points at \(x \equiv p / m_e c=\pm \mathrm{i}\), any series expansion of \(E_p\) around the static nonrelativistic limit \(x=0\) is only related to the exact expression for \(E_p\) for non-ultrarelativistic values of the momentum, i.e., \(|x|<1\). This is most easily seen by rewriting \(E_p\) as
\[E_p=m_e c^2 \sqrt{1+x^2} \equiv m_e c^2 \exp \left[\frac{1}{2} \ln \left(1+x^2\right)\right]\]
From this equation it is evident that \(E_p\) has singularities at \(\pm\mathrm{i}\). Thus only within the shaded disk of radius 1 series expansions of \(E_p\) around \(x=0\) are legitimate, and only this simply connected domain contains the whole real axis, i.e., all physically relevant values for the momentum \(x=p/m_ec\).
For larger momenta, which inevitably arise in a complete basis-set
expansion, the series expansion up to any arbitrary or even infinite order in
\(p/m_ec\) does not represent the original function \(E_p\). For these momenta, this series does not even converge at all!
The singular behavior of the series expansions of \(E_p\), and hence of \(A_p\) and \(R_p\), becomes therefore the worse the more terms of the expansions are taken into account.

Q: Explain the main principle and procedure of the Douglas-Kroll-Hess transformation! Which problem arises?
A: While an expansion in powers of \(1/c\) (FW) is ill-defined, only one expansion parameter remains, which is the external potential that leads to the Douglas-Kroll-Hess (DKH) transformation. In order to eliminate the odd operator \(\mathcal{O}_0\) of the Dirac Hamiltonian order by order in the scalar potential \(V\), an odd operator that depends on \(\bm V\) must be generated. Only the special closed-form free-particle Foldy-Wouthuysen transformation produces an operator \(\mathcal{O}_1\) linear in \(V\) as indicated by the subscript. This is the mandatory starting point for subsequent transformation steps.
After the initial transformation \(U_0\) (see fpFW), the new odd part \(\mathcal{O}_1\) in \(f_1\) needs to be deleted by the next transformation \(U_1\), and all that we require is that new odd terms are of higher order in the scalar potential \(V\).
Then the procedure is basically the same as in the FW expansion! It is important to note that upon each transformation two new even terms with increasing power in \(V\) are obtained, i.e.,
\[\begin{aligned} \text{after }U_0\text{: } h_\mathrm{DKH1} = \sum_{k=0}^1\mathcal{E}_k \\ \text{after }U_1\text{: } h_\mathrm{DKH3} = \sum_{k=0}^3\mathcal{E}_k \\ \text{after }U_2\text{: } h_\mathrm{DKH5} = \sum_{k=0}^5\mathcal{E}_k \end{aligned}\]
All terms contributing to \(h_\mathrm{DKH5}\) are independent of the coefficients \(a_{1,k}\), i.e., they are invariant under an arbitrary change of the parametrization of \(U_1\)! Terms contributing to the fifth- and higher-order DKH corrections, however, depend on \(a_{1,k}\)! The dependence of numerical results on the expansion coefficients is negligible and diminishes the better the higher orders approximate infinite order. But in physical sense not nice!

Q: Sketch the energy of one electron atoms \(E_{\mathrm{1s}}/Z^2\) as a function of the nuclear charge \(Z^2\) for the non-relativistic Schrödinger case, DKH1 to DKH5 and the analytical Dirac eigenvalue! Furthermore, plot \(E_{\mathrm{1s}}\) for different DKH orders! What do you observe?
A: Observations from these results:
- DKH1 is no reliable approximation, but already DKH2 is very accurate.
- Odd DKH-orders are below the reference, even orders are above \(\rightarrow\) Oscillatory convergence.
- Parameter dependence is much smaller than the change of energy order by order.

Q: Outline the general procedure to incorporate relativistic DKH calculations into non-relativistic standard procedures!
A: Only the additional \(\braket{{\chi_i}|{\bm{p}\cdot V\cdot\bm{p}}|{\chi_j}}\) integral are needed which are easy to evaluate from the already calculated \(\braket{{\chi_i}|{V}|{\chi_j}}\) integrals as \(\bm{p}\ket{\chi^l} = c_{(l-1)}\ket{\chi^{(l-1)}} + c_{(l+1)}\ket{\chi^{(l+1)}}\).

Q: Which algebraic exact decoupling methods exist?
A: infinite-order DKH feasible with polynomial cost algorithm (D. Peng, K. Hirao, J. Chem. Phys. 130 2009 044102)
Barysz-Sadlej-Snijders (BSS): two-step protocol: fpFW transformation + solution of an \(R\)-operator equation identical to the \(X\)-operator equation
Exact Two-Component Decoupling (X2C): Direct solution of \(U(X)\); requires, however, solution of Dirac equation first; still efficient within a one-electron approximation such as the one-electron DKH scheme
\[\begin{aligned} & \begin{pmatrix}F_{LL} & F_{LS} \\ F_{SL} & F_{SS}\end{pmatrix} \begin{pmatrix}C_L^{(+)}\\C_S^{(+)}\end{pmatrix} = \epsilon^{(+)}\begin{pmatrix}C_L^{(+)}\\ C_S^{(+)}\end{pmatrix} \ \ \\ & X = C_S^{(+)}\left(C_L^{(+)}\right)^{-1} \ \ \\ & U_X = \begin{pmatrix} \Omega & \Omega X^\dagger \\-\Omega^\prime X & \Omega^\prime \end{pmatrix} \end{aligned}\]
Comparison of DKH, BSS, and X2C:\(N =\) number of basis functions.
\(N\) | DKH2 | DKH4 | DKH6 | X2C | DKH8 | BSS | DKH10 | \(t_N/t_{100}\)
\(100\) | 1.8 | 2.6 | 4.4 | 7.0 | 7.4 | 7.5 | 12.6 |
\(200\) | 12.6 | 19.3 | 33.2 | 54.2 | 57.1 | 58.3 | 97.8 | \(\sim7\)
\(300\) | 41.5 | 63.9 | 110.8 | 181.1 | 191.0 | 194.5 | 326.0 | \(\sim25\)
In practical calculations: Atomic construction of U matrix: Accurate local variants do not require matrix transformations in full primitive basis (D. Peng, M. Reiher, Local Relativistic Exact Decoupling, J. Chem. Phys. 136 2012 244108)

Q: Derive the modified Dirac equation! Why is it useful?
A: 1. Replace small component \(\psi^S\) by pseudo-large component \(\varphi^L\) \[c\bm{\sigma\cdot\mom}\varphi^L = 2m_ec^2\psi^S \thus \psi^S = \frac{\bm{\sigma\cdot\mom}}{2m_ec}\varphi^L\]
2. Insert into split form of the four-component Fock equation yields: \[\begin{aligned} V\psi^L + \underbrace{c\bm{\sigma\cdot\mom}\frac{\bm{\sigma\cdot\mom}}{2m_ec}}_{=(\bm{\sigma\cdot\mom})(\bm{\sigma\cdot\mom})/2m_e = \mom^2\bm{1}_2/2m_e \equiv T}\varphi^L &= \epsilon\psi^L \\ \bm{\sigma\cdot\mom}\psi^L + (V-2m_ec^2)\frac{\bm{\sigma\cdot\mom}}{2m_ec}\varphi^L &= \epsilon\frac{\bm{\sigma\cdot\mom}}{2m_ec}\varphi^L \end{aligned}\]
3. Multiplying the second equation by \((\bm{\sigma\cdot\mom})/2m_e\) we arrive at the modified Dirac equation\[\begin{aligned} (V-\epsilon)\psi^L + \operatorname{T}\varphi^L &= 0 \\ \operatorname{T}\psi^L + \left[\frac{1}{4m_e^2c^2}(\bm{\sigma\cdot\mom})(V-\epsilon)(\bm{\sigma\cdot\mom})-\operatorname{T}\right]\varphi^L &= 0 \end{aligned}\]
Essential ingredient of X2C: diagonalize the modified Dirac equation.
most important point: \(\varphi^L\) and \(\psi^L\) can now be expanded into the same basis set!