# Unity of Kohn-Sham Density Functional Theory and Reduced Density Matrix
Functional Theory

Neil Qiang Su<EMAIL_ADDRESS>Department of Chemistry, Key Laboratory of
Advanced Energy Materials Chemistry (Ministry of Education) and Renewable
Energy Conversion and Storage Center (RECAST), Nankai University, Tianjin
300071, China

###### Abstract

This work presents a theory to unify the two independent theoretical
frameworks of Kohn-Sham (KS) density functional theory (DFT) and reduced
density matrix functional theory (RDMFT). The generalization of the KS
orbitals to hypercomplex number systems leads to the hypercomplex KS (HCKS)
theory, which extends the search space for the density in KS-DFT to a space
that is equivalent to natural spin orbitals with fractional occupations in
RDMFT. Thereby, HCKS is able to capture the multi-reference nature of strong
correlation by dynamically varying fractional occupations. Moreover, the
potential of HCKS to overcome the fundamental limitations of KS is verified on
systems with strong correlation, including atoms of transition metals. As a
promising alternative to the realization of DFT, HCKS opens up new
possibilities for the development and application of DFT in the future.

$Introduction.$—Built upon the Hohenberg-Kohn theorem [1, 2], Kohn-Sham (KS)
density functional theory (DFT) [3, 4, 5] is a formally exact theoretical
framework toward the many-electron problem. Due to the favorable balance
between accuracy and efficiency, KS-DFT has won enormous popularity that can
manifest itself in the countless applications across physics, materials
science, chemistry, and biology [6, 7, 8]. Nonetheless, the great success of
KS-DFT, along with commonly used density functional approximations (DFAs), is
clouded by the improper treatment of strong correlation [9]. Strong
correlation represents the intractable electronic interaction stemming from
the multi-reference nature of systems, which has long posed a major challenge
to KS-DFT [10, 11, 12, 13].

It is generally recognized that intrinsic errors associated with commonly used
DFAs [14, 15, 16, 17, 18, 9, 19, 20] and the slow progress in systematically
eliminating these errors to better handle strong correlation have severely
limited the applicability of KS-DFT [18, 9, 21, 22]. The enlightening work by
Lee, Bertels, Small, and Head-Gordon [23] shows that commonly used DFAs can
better treat the strong correlation in singlet biradicals by using complex
spin-restricted orbitals in KS-DFT. This thus raises a deep question: In
addition to the errors inherent in existing DFAs, is there still any
limitation in the understanding and application of the KS-DFT framework?
Apparently, in-depth insights into this question is important for further
development and application of DFT.

Beside KS-DFT, reduced density matrix functional theory (RDMFT) [24, 25, 26,
27, 28, 29, 30, 31, 32] provides an alternative approach to the many-electron
problem. The Gilbert’s theorem [24] guarantees that the one-electron reduced
density matrix (1-RDM) instead of the density can be used as the fundamental
variable for the energy functional, which thus makes RDMFT an exact
theoretical framework independent of KS-DFT. RDMFT has proved its great
potential to overcome the fundamental limitations of KS-DFT through the
successful application in predicting dissociation energy curves [33, 34, 35],
and fundamental gaps for finite systems and extended solids [36] as well as
for Mott insulators [28]. Therefore, establishing a connection with the
theoretical framework of RDMFT will further improve our understanding of KS-
DFT and the problems of existing DFAs.

This letter seeks to come up with a theory to unify the two theoretical
frameworks of KS-DFT and RDMFT. This is achieved by generalizing the
conventional KS determinant to hypercomplex number systems. The resulting
hypercomplex KS (HCKS) theory extends the search space for the density in KS-
DFT to a space that is equivalent to natural spin orbitals with fractional
occupations in RDMFT. The potential of HCKS in capturing the physical essence
of strong correlation is demonstrated on atoms of multi-reference character,
including transition metals.

$Theory.$—The total energy $E_{\rm{tot}}[\rho_{\sigma}]$ in KS-DFT includes
the kinetic energy $T_{\rm{s}}[\rho_{\sigma}]$, the external energy
$E_{\rm{ext}}[\rho_{\sigma}]$, the Coulomb energy $E_{\rm{H}}[\rho_{\sigma}]$,
and the exchange-correlation (XC) energy $E_{\rm{XC}}[\rho_{\sigma}]$. They
are all functionals of the density that can be formed by the occupied KS
orbitals in the KS determinant ($N_{\sigma}$ being the $\sigma$-spin electron
number),

$\rho_{\sigma}({\bf r})=\sum_{k=1}^{N_{\sigma}}|\varphi_{k}^{\sigma}({\bf
r})|^{2},$ (1)

while $T_{\rm{s}}[\rho_{\sigma}]$ can be explicitly formulated as
$-\frac{1}{2}\sum_{\sigma}^{\alpha,\beta}\sum_{k=1}^{N_{\sigma}}\langle\varphi_{k}^{\sigma}|\nabla^{2}|\varphi_{k}^{\sigma}\rangle$
[3]. Therefore, the minimization of $E_{\rm{tot}}[\rho_{\sigma}]$ with respect
to $\rho_{\sigma}$ is equivalent to the minimization with respect to
$\\{\varphi_{p}^{\sigma}\\}$, subject to the orthonormalization condition

$\langle\varphi_{p}^{\sigma}|\varphi_{q}^{\sigma}\rangle=\delta_{pq}.$ (2)

In RDMFT, the four terms in the total energy are uniquely determined by 1-RDM
($\gamma_{\sigma}$), which are $T[\gamma_{\sigma}]$,
$E_{\rm{ext}}[\gamma_{\sigma}]$, $E_{\rm{H}}[\gamma_{\sigma}]$ and
$E_{\rm{XC}}[\gamma_{\sigma}]$ respectively. In terms of the natural spin
orbitals $\\{\psi_{p}^{\sigma}\\}$ and the natural occupation numbers
$\\{n_{p}^{\sigma}\\}$, $\gamma_{\sigma}$ in the spectral representation reads
$\sum_{p=1}^{K}|\psi_{p}^{\sigma}\rangle
n_{p}^{\sigma}\langle\psi_{p}^{\sigma}|$ ($K$ being the dimension of the basis
set), and the diagonal element $\langle{\bf r}|\gamma_{\sigma}|{\bf r}\rangle$
is the density,

$\rho_{\sigma}({\bf r})=\sum_{p=1}^{K}n_{p}^{\sigma}|\psi_{p}^{\sigma}({\bf
r})|^{2}.$ (3)

$T[\gamma_{\sigma}]$ has the form,
$-\frac{1}{2}\sum_{\sigma}^{\alpha,\beta}\sum_{p=1}^{K}n_{p}^{\sigma}\langle\psi_{p}^{\sigma}|\nabla^{2}|\psi_{p}^{\sigma}\rangle$.
The ground-state energy can be obtained by the minimization of
$E_{\rm{tot}}[\gamma_{\sigma}]$ with respect to $\gamma_{\sigma}$, or
equivalently, with respect to both $\\{\psi_{p}^{\sigma}\\}$ and
$\\{n_{p}^{\sigma}\\}$, subject to the orthonormalization condition,
$\langle\psi_{p}^{\sigma}|\psi_{q}^{\sigma}\rangle=\delta_{pq}$, and the Pauli
exclusion principle and the $N$-representability constraint [37],

$0\leq n_{p}^{\sigma}\leq 1,\sum_{p=1}^{K}n_{p}^{\sigma}=N_{\sigma}.$ (4)

Unlike KS-DFT, RDMFT allows dynamically varying fractional occupations to
capture the multi-reference nature of strong correlation.

The connection between KS-DFT and RDMFT is established by introducing
hypercomplex. The concept of hypercomplex [38], and the theory of Clifford
algebra that generalizes real numbers, complex numbers to quanternions,
octonions, and other hypercomplex numbers have important applications in a
variety of fields including theoretical physics [39]. Here, the HCKS orbitals
are formulated as

$\varphi_{p}^{\sigma}({\bf r})=\phi_{p}^{\sigma,0}({\bf
r})+\sum_{\mu=1}^{n}\phi_{p}^{\sigma,\mu}({\bf r})e_{\mu},$ (5)

where $\\{\phi_{p}^{\sigma,\mu}\\}$ are a set of real functions, and
$\\{e_{1},e_{2},\cdots,e_{n}\\}$ are a basis of dimension $n$ in a Clifford
algebra, such that [39]

$e_{\mu}^{2}=-1;e_{\mu}e_{\nu}=-e_{\nu}e_{\mu}.$ (6)

The conjugate hypercomplex of the HCKS orbitals are
$\bar{\varphi}_{p}^{\sigma}({\bf r})=\phi_{p}^{\sigma,0}({\bf
r})-\sum_{\mu=1}^{n}\phi_{p}^{\sigma,\mu}({\bf r})e_{\mu}$. Therefore, Eq 5
provides a general set of high dimensional orbitals for the HCKS determinant,
while complex KS orbitals are a special case for $n=1$.

Without loss of generality, $\\{\phi_{p}^{\sigma,\mu}\\}$ can be expanded on a
set of orthonormal functions $\\{\chi_{p}\\}$ and read

$\phi_{p}^{\sigma,\mu}({\bf r})=\sum_{q=1}^{K}\chi_{q}({\bf
r})V_{pq}^{\sigma,\mu}.$ (7)

Here ${\rm{V}}^{\sigma,\mu}$ is a $K\times K$ matrix associated with the
$\mu$-th component of the HCKS orbitals. The orthonormalization condition of
Eq 2 becomes [40]

$\sum_{\mu=0}^{n}{\rm{V}}^{\sigma,\mu}{\rm{V}}^{\sigma,\mu
T}=\sum_{\mu=0}^{n}{\rm{V}}^{\sigma,\mu T}{\rm{V}}^{\sigma,\mu}={\rm{I}}_{K},$
(8)

and

$\displaystyle\begin{cases}{\rm{V}}^{\sigma,\mu}{\rm{V}}^{\sigma,\nu
T}={\rm{V}}^{\sigma,\nu}{\rm{V}}^{\sigma,\mu T};\\\ {\rm{V}}^{\sigma,\mu
T}{\rm{V}}^{\sigma,\nu}={\rm{V}}^{\sigma,\nu
T}{\rm{V}}^{\sigma,\mu},\end{cases}$ (9)

where the superscript $T$ denotes the transpose, and ${\rm{I}}_{K}$ is the
$K\times K$ identity matrix. The density corresponding to the HCKS determinant
reads [40],

$\rho_{\sigma}({\bf r})=\sum_{p,q=1}^{K}\chi_{p}({\bf
r})D_{pq}^{\sigma}\chi_{q}({\bf r}),$ (10)

and ${\rm{D}}^{\sigma}$ is defined by

${\rm{D}}^{\sigma}=\sum_{\mu=0}^{n}{\rm{V}}^{\sigma,\mu
T}{\rm{I}}_{K}^{N_{\sigma}}{\rm{V}}^{\sigma,\mu},$ (11)

where ${\rm{I}}_{K}^{N_{\sigma}}$ is a $K\times K$ diagonal matrix, with the
first $N_{\sigma}$ diagonal elements being 1 and the rest being 0.
${\rm{D}}^{\sigma}$ is symmetric and can be diagonalized by an orthogonal
matrix ${\rm{U}}^{\sigma}$,

${\rm{D}}^{\sigma}={\rm{U}}^{\sigma}{\rm{\Lambda}}^{\sigma}{\rm{U}}^{\sigma
T},$ (12)

where ${\rm{\Lambda}}^{\sigma}$ is a diagonal matrix,
diag$(\lambda_{1}^{\sigma},\lambda_{2}^{\sigma},\cdots,\lambda_{K}^{\sigma})$,
and $\\{\lambda_{p}^{\sigma}\\}$ are the eigenvalues of ${\rm{D}}^{\sigma}$,
which satisfy [40]

$0\leq\lambda_{p}^{\sigma}\leq
1;\sum_{p=1}^{K}\lambda_{p}^{\sigma}=N_{\sigma}.$ (13)

By inserting Eq 12, the density of Eq 10 can be written as

$\rho_{\sigma}({\bf
r})=\sum_{p=1}^{K}{\lambda}_{p}^{\sigma}|\chi_{p}^{\sigma}({\bf r})|^{2},$
(14)

where $\chi_{p}^{\sigma}({\bf r})=\sum_{q=1}^{K}\chi_{q}({\bf
r})U_{qp}^{\sigma}$ (thereby $\\{\chi_{k}^{\sigma}\\}$ being orthonormal).
Similarly, the kinetic energy reads

$T_{\rm{s}}[\rho_{\sigma}]=-\frac{1}{2}\sum_{\sigma}^{\alpha,\beta}\sum_{p=1}^{K}\lambda_{p}^{\sigma}\langle\chi_{p}^{\sigma}|\nabla^{2}|\chi_{p}^{\sigma}\rangle.$
(15)

To further verify that $\\{\lambda_{p}^{\sigma}\\}$ in Eqs 14 and 15,
associated with any given orthonormal functions $\\{\chi_{k}^{\sigma}\\}$, can
take any values subject to Eq 13 when $K\leq n+1$, here construct a special
set of HCKS orbitals. Given a set of orthonormal $\\{\chi_{k}^{\sigma}\\}$ ,
$\\{\phi_{p}^{\sigma,\mu}\\}$ of Eq 7 are constructed with the following form

$\phi_{p}^{\sigma,\mu}=\chi_{\mu+1}^{\sigma}({\bf r})V_{p\mu+1}^{\sigma,\mu},$
(16)

i.e. the components corresponding to the same $e_{\mu}$ in all the HCKS
orbitals are formed by the same function $\chi_{\mu+1}^{\sigma}$, thereby

$\displaystyle V_{pq}^{\sigma,\mu}=\begin{cases}W_{pq}^{\sigma},&q=\mu+1\\\
0,&q\neq\mu+1,\end{cases}$ (17)

where ${\rm{W}}^{\sigma}$ is a $K\times K$ matrix. The orthonormalization
condition in terms of ${\rm{W}}^{\sigma}$ reads

${\rm{W}}^{\sigma}{\rm{W}}^{\sigma T}={\rm{I}}_{K},$ (18)

which thus requires ${\rm{W}}^{\sigma}$ to be an orthogonal matrix. Similar
derivations lead to Eqs 14 and 15 for $\rho_{\sigma}({\bf r})$ and
$T[\rho_{\sigma}]$ respectively, with
$\lambda_{p}^{\sigma}=\sum_{k=1}^{N_{\sigma}}(W_{kp}^{\sigma})^{2}$. As
${\rm{W}}^{\sigma}$ can be any orthogonal matrix, $\\{\lambda_{p}^{\sigma}\\}$
can take any values subject to Eq 13.

Therefore, the extension of KS-DFT to hypercomplex orbitals leads to a density
that has the same search space as the density of Eq 3 constructed by real
orbitals in RDMFT, with the same form of kinetic energies as well. The
resulting HCKS theory thus unifies KS-DFT and RDMFT, which extends the search
space of KS-DFT to natural spin orbitals with fractional occupations. When
$n=0$, the orbitals of Eq 5 are real and HCKS reduces the conventional KS
method; when $n=1$, the orbitals are complex, which are used in the complex,
spin-restricted KS (CRKS) method [23].

$Results.$—The numerical performance of HCKS and its potential in describing
strong correlation are evaluated with the finite basis simulation. Atoms of
multi-reference character, including transition metals, are calculated. While
two XC functionals, PBE [41] and BLYP [42, 43], are examined here, similar
results can be obtained with other commonly used functionals. All calculations
were performed using a local modified version of the NWChem package [44].

Figure 1: Density on real-space grids for the lowest singlet state of C atom.
Views along the z, y and x axes are provided respectively. (a) RKS density
(purple), (b) CRKS density (green), and (c) HCKS (or equivalently, RDMFT)
density (blue) are plotted at the isosurface value of 0.2 au. The XC
functional PBE is applied for all the calculation. The basis set used is aug-
cc-pVQZ [45, 46].

The density for the lowest singlet state of C atom is first examined. The
singlet-state C atom is of biradical nature, with two electrons in the 2p sub-
shell. Normally, the spin-restricted KS (RKS) would have one doubly occupied p
orbital and destroy the degeneracy of the p orbitals. Fig 1a shows the density
of RKS when pz is occupied, which maintains the spin symmetry but loses the
space symmetry. Different from RKS, CRKS leads to two half filled p orbitals
for the singlet C atom [23, 40]. Fig 1b shows that the CRKS density with both
px and py half filled maintains the space symmetry along the z axis, but it
cannot guarantee the apace symmetry along both x and y axes. In contrast,
HCKS, or equivalently, RDMFT, automatically converges to the spin-restricted
solution, with each p orbital 1/3 filled. Therefore, the density by HCKS can
maintain both spin and space symmetry; see Fig 1c. Similar results can be
obtained for other systems of multi-reference character such as singlet states
of O, S, and Si atoms, and even the d and f block transition metals.

Figure 2: Triplet-singlet energy gaps of C atom under varying external charged
environments. The energies are calculated with two point charges of 0.3 au
each placed equidistantly on two sides of the C atom, with the distance
between the two point charges ranging from 4 to 12 Å. KS and HCKS (or
equivalently, RDMFT) in use of the same XC functionals (PBE and BLYP) are
examined, while CCSDT [47] results are used as reference. The basis set used
is aug-cc-pVQZ. All energies are in kcal/mol.

In addition, the energies of C atom under varying external charged
environments are tested. Fig 2 shows that the triplet-singlet energy gaps by
KS seriously deviate from the results of CCSDT, especially when the point
charges are far away from the C atom. By comparison, HCKS significantly
improves the performance of KS, in use of the same XC functionals. The
occupations of the three p orbitals change gradually from $\\{1,0,0\\}$ to
$\\{1/3,1/3,1/3\\}$ as the point charges move away from the C atom. Therefore,
HCKS is able to capture the physical essence of multi-reference nature through
dynamically varying occupations.

Figure 3: Energy gaps between high and low spin-states of transition metals.
Septet-singlet, quintet-singlet, triplet-singlet energy gaps for Cr, Fe and Ni
atoms respectively are calculated. KS and HCKS (or equivalently, RDMFT) in use
of the same XC functionals (PBE and BLYP) are examined, while the experimental
data [48] are used as reference. The basis set used is def2-TZVP [49]. All
energies are in kcal/mol.

Here test also some atoms of 3d transition metals, which are generally
considered to be a big challenge to KS. Due to the partially filled 3d sub-
shell and the near degeneracy of 4s and 3d sub-shells, these atoms or systems
containing them often have a plethora of low-lying degenerate and near-
degenerate states, which make them much more difficult to correctly describe
than main-group compounds. Fig 3 shows that KS with both PBE and BLYP XC
functionals overestimates the energy gaps between high and low spin-states of
Cr, Fe and Ni, due to the lack of strong correlation in the low spin-states.
In contrast, HCKS improves the performance in use of the same XC functionals.
The fractional occupations for (nearly) degenerate orbitals render HCKS great
potential for handling strong correlation. For example, PBE@HCKS converges to
a set of spin-restricted orbitals for the singlet state of Ni, with
occupations for both $\alpha$ and $\beta$ spins being 0.916 for each 3d
orbital and 0.420 for 4s orbital. This further proves that HCKS can improve
the description of strong correlation while maintaining both spin and space
symmetry.

$Conclusions.$—This letter presented a theory to unify the two theoretical
frameworks of KS-DFT and RDMFT. This is achieved by HCKS that generalizes the
KS orbitals to hypercomplex number systems. HCKS extends the search space for
the density in KS-DFT to a space that is equivalent to natural spin orbitals
with fractional occupations in RDMFT. The test on the singlet biradical C atom
shows that HCKS can maintain both spin and space symmetry for the density with
equally fractionally occupied p orbitals, which cannot be achieved by RKS.
Besides, the calculation of energy gaps between high and low spin states, for
both C atom under varying external charged environments and atoms of
transition metals, demonstrates that HCKS is able to capture the multi-
reference nature of strong correlation by dynamically varying fractional
occupations, while KS in use of the same XC functionals cannot. Therefore,
HCKS shows great potential to overcome the fundamental limitations of KS,
which thus provides an alternative to the realization of DFT, and opens up new
channels for the development and evaluation of approximate functionals.

The Supplemental Material for this work is available on line [40]. Support
from the National Natural Science Foundation of China (Grant 22073049), the
Natural Science Foundation of Tianjin City (20JCQNJC01760), and Fundamental
Research Funds for the Central Universities (Nankai University: No. 63206008)
is appreciated. Dedicated to the 100th anniversary of Chemistry at Nankai
University.

## References

  * [1] P. Hohenberg and W. Kohn,  Phys. Rev. 136, B864 (1964).
  * [2] M. Levy,  Proc. Natl. Acad. Sci. USA 76, 6062 (1979).
  * [3] W. Kohn and L. J. Sham,  Phys. Rev. 140, A1133 (1965).
  * [4] R. G. Parr and W. Yang,  Density-Functional Theory of Atoms and Molecules,  Oxford University Press: New York, 1989.
  * [5] R. Dreizler and E. Gross,  Density Functional Theory: An Approach to the Quantum Many-Body Problem,  Springer Berlin Heidelberg, 2012.
  * [6] R. Peverati and D. G. Truhlar,  Philos. Trans. Royal. Soc. A 372, 20120476 (2014).
  * [7] R. O. Jones,  Rev. Mod. Phys. 87, 897 (2015).
  * [8] J. Sun et al.,  Nat. Chem. 8, 831 (2016).
  * [9] A. J. Cohen, P. Mori-Sánchez, and W. Yang,  Science 321, 792 (2008).
  * [10] K. Burke,  J. Chem. Phys. 136, 150901 (2012).
  * [11] A. J. Cohen, P. Mori-Sánchez, and W. Yang,  Chem. Rev. 112, 289 (2012).
  * [12] A. D. Becke,  J. Chem. Phys. 140, 18A301 (2014).
  * [13] N. Q. Su and X. Xu,  Annu. Rev. Phys. Chem. 68, 155 (2017).
  * [14] A. Ruzsinszky, J. P. Perdew, G. I. Csonka, O. A. Vydrov, and G. E. Scuseria,  J. Chem. Phys. 125, 194112 (2006).
  * [15] A. Ruzsinszky, J. P. Perdew, G. I. Csonka, O. A. Vydrov, and G. E. Scuseria,  J. Chem. Phys. 126, 104102 (2007).
  * [16] O. A. Vydrov, G. E. Scuseria, and J. P. Perdew,  J. Chem. Phys. 126, 154109 (2007).
  * [17] P. Mori-Sánchez, A. J. Cohen, and W. Yang,  Phys. Rev. Lett. 100, 146401 (2008).
  * [18] P. Mori-Sánchez, A. J. Cohen, and W. Yang,  Phys. Rev. Lett. 102, 066403 (2009).
  * [19] X. Zheng, A. J. Cohen, P. Mori-Sánchez, X. Hu, and W. Yang,  Phys. Rev. Lett. 107, 026403 (2011).
  * [20] J.-D. Chai and P.-T. Chen,  Phys. Rev. Lett. 110, 033002 (2013).
  * [21] N. Mardirossian and M. Head-Gordon,  Mol. Phys. 115, 2315 (2017).
  * [22] N. Q. Su, C. Li, and W. Yang,  Proc. Natl. Acad. Sci. USA 115, 9678 (2018).
  * [23] J. Lee, L. W. Bertels, D. W. Small, and M. Head-Gordon,  Phys. Rev. Lett. 123, 113001 (2019).
  * [24] T. L. Gilbert,  Phys. Rev. B 12, 2111 (1975).
  * [25] A. Muller,  Phys. Lett. A 105, 446 (1984).
  * [26] S. Goedecker and C. J. Umrigar,  Phys. Rev. Lett. 81, 866 (1998).
  * [27] K. Pernal,  Phys. Rev. Lett. 94, 233002 (2005).
  * [28] S. Sharma, J. K. Dewhurst, N. N. Lathiotakis, and E. K. U. Gross,  Phys. Rev. B 78, 201103(R) (2008).
  * [29] M. Piris, J. M. Matxain, X. Lopez, and J. M. Ugalde,  J. Chem. Phys. 132, 031103 (2010).
  * [30] S. Sharma, J. K. Dewhurst, S. Shallcross, and E. K. U. Gross,  Phys. Rev. Lett. 110, 116403 (2013).
  * [31] R. Schade, E. Kamil, and P. Blöchl,  Eur. Phys. J. Special Topics 226, 2677 (2017).
  * [32] C. Schilling and R. Schilling,  Phys. Rev. Lett. 122, 013001 (2019).
  * [33] O. Gritsenko, K. Pernal, and E. J. Baerends,  J. Chem. Phys. 122, 204102 (2005).
  * [34] D. R. Rohr, K. Pernal, O. V. Gritsenko, and E. J. Baerends,  J. Chem. Phys. 129, 164105 (2008).
  * [35] N. N. Lathiotakis, S. Sharma, J. K. Dewhurst, F. G. Eich, M. A. L. Marques, and E. K. U. Gross,  Phys. Rev. A 79, 040501(R) (2009).
  * [36] N. N. Lathiotakis et al.,  Z. Phys. Chem. 224, 467 (2010).
  * [37] A. J. Coleman,  Rev. Mod. Phys. 35, 668 (1963).
  * [38] I. L. Kantor and A. S. Solodovnikov,  Hypercomplex numbers,  Springer-Verlag: Berlin, 1989.
  * [39] J. Vaz Jr. and R. da Rocha Jr.,  An Introduction to Clifford Algebras and Spinors,  Oxford University Press: New York, 2019.
  * [40] See the Supplemental Material at xxx for the detailed derivation of HCKS theory.
  * [41] J. P. Perdew, K. Burke, and M. Ernzerhof,  Phys. Rev. Lett. 77, 3865 (1996).
  * [42] A. D. Becke,  Phys. Rev. A 38, 3098 (1988).
  * [43] C. Lee, W. Yang, and R. G. Parr,  Phys. Rev. B 37, 785 (1988).
  * [44] E. Aprá et al.,  J. Chem. Phys. 152, 184102 (2020).
  * [45] T. H. Dunning,  J. Chem. Phys. 90, 1007 (1989).
  * [46] R. A. Kendall, T. H. Dunning, and R. J. Harrison,  J. Chem. Phys. 96, 6796 (1992).
  * [47] I. Shavitt and R. Bartlett,  Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory,  Cambridge University Press: Cambridge, 2009.
  * [48] C. Moore,  Atomic Energy Levels as Derived from Analysis of Optical Spectra, Volume II: 24Cr-41Nb,  U.S. Government Printing Office: Washington, D.C., 1952.
  * [49] F. Weigend and R. Ahlrichs,  Phys. Chem. Chem. Phys. 7, 3297 (2005).

Supporting Information of

Unity of Kohn-Sham Density Functional Theory and Reduced Density Matrix
Functional Theory

Neil Qiang Su

_Department of Chemistry, Key Laboratory of Advanced Energy Materials
Chemistry (Ministry of Education) and Renewable Energy Conversion and Storage
Center (RECAST), Nankai University, Tianjin 300071, China
_ _Email<EMAIL_ADDRESS>
###### Contents

  1. I Reduced Density Matrix Functional Theory
  2. II Kohn-Sham Density Functional Theory
  3. III Hypercomplex numbers
  4. IV Hypercomplex Kohn-Sham Theory

##  I Reduced Density Matrix Functional Theory

In RDMFT, the electronic ground-state energy can be written as a functional of
the one-body reduced density matrix (1-RDM) $\gamma_{\sigma}$ [24, 2]

$E^{\rm{DMFT}}[\gamma_{\sigma}]=T[\gamma_{\sigma}]+E_{\rm{ext}}[\gamma_{\sigma}]+E_{\rm{H}}[\gamma_{\sigma}]+E_{\rm{XC}}[\gamma_{\sigma}],$
(S1)

where the kinetic energy is

$T[\gamma_{\sigma}]=-\frac{1}{2}\sum_{\sigma}^{\alpha,\beta}\int\nabla_{{\bf
r}}^{2}\gamma_{\sigma}({\bf r},{\bf r}^{\prime})|_{{\bf r}^{\prime}={\bf
r}}d{\bf r},$ (S2)

the external energy is

$E_{\rm{ext}}[\gamma_{\sigma}]=\sum_{\sigma}^{\alpha,\beta}\int
v_{\rm{ext}}({\bf r})\gamma_{\sigma}({\bf r},{\bf r}^{\prime})|_{{\bf
r}^{\prime}={\bf r}}d{\bf r},$ (S3)

the Coulomb energy is

$E_{\rm{H}}[\gamma_{\sigma}]=\frac{1}{2}\sum_{\sigma\sigma^{\prime}}^{\alpha,\beta}\int\frac{\gamma_{\sigma}({\bf
r}_{1},{\bf r}^{\prime}_{1})\gamma_{\sigma^{\prime}}({\bf r}_{2},{\bf
r}^{\prime}_{2})|{{}^{{\bf r}^{\prime}_{1}={\bf r}_{1}}_{{\bf
r}^{\prime}_{2}={\bf r}_{2}}}}{r_{12}}d{\bf r}_{1}d{\bf r}_{2},$ (S4)

and the exchange-correlation (XC) energy $E_{\rm{XC}}[\gamma_{\sigma}]$ is
unknown. The total 1-RDM is

$\gamma({\bf r},{\bf r}^{\prime})=\gamma_{\alpha}({\bf r},{\bf
r}^{\prime})+\gamma_{\beta}({\bf r},{\bf
r}^{\prime})=N\sum_{\sigma_{1},\cdots,\sigma_{N}}^{\alpha,\beta}\int d{\bf
r}_{2}\cdots d{\bf r}_{N}\Psi^{*}({\bf r}^{\prime}\sigma_{1},{\bf
x}_{2},\cdots,{\bf x}_{N})\Psi({\bf r}\sigma_{1},{\bf x}_{2},\cdots,{\bf
x}_{N}),$ (S5)

where ${\bf x}_{p}={\bf r}_{p}\sigma_{p}$ are combined spatial and spin
coordinates. Diagonalization of 1-RDM $\gamma_{\sigma}$ generates the natural
spin orbitals, $\\{\psi_{p}^{\sigma}\\}$, and their occupation numbers,
$\\{n_{p}^{\sigma}\\}$, thus

$\gamma_{\sigma}=\sum_{p=1}^{K}|\psi_{p}^{\sigma}\rangle
n_{p}^{\sigma}\langle\psi_{p}^{\sigma}|;\gamma_{\sigma}({\bf r},{\bf
r}^{\prime})=\langle{\bf r}|\gamma_{\sigma}|{\bf r}^{\prime}\rangle,$ (S6)

where its diagonal element is the density,

$\rho_{\sigma}({\bf r})=\langle{\bf r}|\gamma_{\sigma}|{\bf
r}\rangle=\sum_{p=1}^{K}n_{p}^{\sigma}|\psi_{p}^{\sigma}({\bf r})|^{2}.$ (S7)

The ground-state energy can be obtained by minimizing Eq. S1 with respect to
1-RDM, or equivalently, by minimizing Eq. S1 with respect to both
$\\{\psi_{p}^{\sigma}\\}$ and $\\{n_{p}^{\sigma}\\}$, subject to

$\langle\psi_{p}^{\sigma}|\psi_{q}^{\sigma}\rangle=\delta_{pq},0\leq
n_{p}^{\sigma}\leq 1,\sum_{p=1}^{K}n_{p}^{\sigma}=N_{\sigma}.$ (S8)

Here, $\sigma$ is the spin index, it can be $\alpha$ or $\beta$; $N_{\sigma}$
is the electron number of $\sigma$ spin; $K$ is the dimension of basic
functions, which should be not smaller than $N_{\sigma}$.

##  II Kohn-Sham Density Functional Theory

In KS-DFT [1, 3], the total energy is

$E^{\rm{KS}}[\rho_{\sigma}]=T_{\rm{s}}[\rho_{\sigma}]+E_{\rm{ext}}[\rho_{\sigma}]+E_{\rm{H}}[\rho_{\sigma}]+E_{\rm{XC}}[\rho_{\sigma}],$
(S9)

where the kinetic energy is

$T_{\rm{s}}[\rho_{\sigma}]=-\frac{1}{2}\sum_{\sigma}^{\alpha,\beta}\int\nabla_{{\bf
r}}^{2}\gamma_{\sigma}^{s}({\bf r},{\bf r}^{\prime})|_{{\bf r}^{\prime}={\bf
r}}d{\bf
r}=-\frac{1}{2}\sum_{\sigma}^{\alpha,\beta}\sum_{k=1}^{N_{\sigma}}\langle\varphi_{k}^{\sigma}|\nabla^{2}|\varphi_{k}^{\sigma}\rangle,$
(S10)

the external energy is

$E_{\rm{ext}}[\rho_{\sigma}]=\sum_{\sigma}^{\alpha,\beta}\int
v_{\rm{ext}}({\bf r})\rho_{\sigma}({\bf r})d{\bf r},$ (S11)

the Coulomb energy is

$E_{\rm{H}}[\rho_{\sigma}]=\frac{1}{2}\sum_{\sigma\sigma^{\prime}}^{\alpha,\beta}\int\frac{\rho_{\sigma}({\bf
r}_{1})\rho_{\sigma^{\prime}}({\bf r}_{2})}{r_{12}}d{\bf r}_{1}d{\bf r}_{2},$
(S12)

while the form of $E_{\rm{XC}}[\rho_{\sigma}]$ is unknown.

$\gamma_{\sigma}^{s}$ is 1-RDM of the noninteracting reference system, which
can be established by occupied orbitals,

$\gamma_{\sigma}^{s}=\sum_{k=1}^{N_{\sigma}}|\psi_{k}^{\sigma}\rangle\langle\psi_{k}^{\sigma}|;\gamma_{\sigma}^{s}({\bf
r},{\bf r}^{\prime})=\langle{\bf r}|\gamma_{\sigma}^{s}|{\bf
r}^{\prime}\rangle,$ (S13)

and the density is

$\rho_{\sigma}({\bf r})=\langle{\bf r}|\gamma_{\sigma}^{s}|{\bf
r}\rangle=\sum_{k=1}^{N_{\sigma}}|\varphi_{k}^{\sigma}({\bf r})|^{2}.$ (S14)

Here $\\{\varphi_{p}^{\sigma}\\}$ are the KS orbitals, which are a set of
orthonormal orbitals that span the space of dimension $K$,

$\langle\varphi_{p}^{\sigma}|\varphi_{q}^{\sigma}\rangle=\delta_{pq},$ (S15)

thereby

$1_{K}=\sum_{p=1}^{K}|\varphi_{p}^{\sigma}\rangle\langle\varphi_{p}^{\sigma}|.$
(S16)

Here $1_{K}$ represents the identity operator of the space.

##  III Hypercomplex numbers

Hypercomplex [38] and the theory of Clifford algebra [39] that generalizes
real numbers, complex numbers to quanternions, octonions, and other
hypercomplex numbers are introduced here.

Clifford algebra is a unital associative algebra generated by a vector space
with a quadratic form, which has important applications in a variety of fields
including theoretical physics. In a Clifford algebra, given a basis of
dimension $n$, $\\{e_{1},e_{2},\cdots,e_{n}\\}$, such that [39]

$e_{\mu}^{2}=-1;e_{\mu}e_{\nu}=-e_{\nu}e_{\mu},$ (S17)

further imposing closure under multiplication generates a multivector space
spanned by a basis of the $2^{n}$ distinct products of
$\\{e_{1},e_{2},\cdots,e_{n}\\}$, i.e.
$\\{1,e_{1},e_{2},\cdots,e_{1}e_{2},\cdots,e_{1}e_{2}e_{3},\cdots\\}$. These
distinct products form the basis of a hypercomplex number system. Unlike
$\\{e_{1},e_{2},\cdots,e_{n}\\}$, the remaining elements in this basis need
not anti-commute, depending on how many simple exchanges must be carried out
for the swap.

Here we consider the following hypercomplex number

$C=\sum_{\mu=0}^{n}c^{\mu}e_{\mu}=c^{0}\cdot
1+\sum_{\mu=1}^{n}c^{\mu}e_{\mu},$ (S18)

where $e_{0}=1$, and the coefficients $\\{c^{\mu}\\}$ are real. The
corresponding conjugate hypercomplex number of Eq. S18 is

$\bar{C}=c^{0}\cdot 1-\sum_{\mu=1}^{n}c^{\mu}e_{\mu}.$ (S19)

##  IV Hypercomplex Kohn-Sham Theory

Now, we generalize the KS orbitals to hypercomplex number systems,

$\varphi_{p}^{\sigma}({\bf r})=\sum_{\mu=0}^{n}\phi_{p}^{\sigma,\mu}({\bf
r})e_{\mu}.$ (S20)

The density of Eq. S14 becomes

$\rho_{\sigma}({\bf
r})=\sum_{k=1}^{N_{\sigma}}\sum_{\mu=0}^{n}[\phi_{k}^{\sigma,\mu}({\bf
r})]^{2},$ (S21)

and the kinetic energy of Eq. S10 becomes

$T_{\rm{s}}[\rho_{\sigma}]=-\frac{1}{2}\sum_{\sigma}^{\alpha,\beta}\sum_{k=1}^{N_{\sigma}}\sum_{\mu=0}^{n}\langle\phi_{k}^{\sigma,\mu}|\nabla^{2}|\phi_{k}^{\sigma,\mu}\rangle.$
(S22)

Without loss of generality, we expand $\\{\phi_{p}^{\sigma,\mu}\\}$ on a set
of orthonormal functions $\\{\chi_{p}\\}$, i.e.

$\phi_{p}^{\sigma,\mu}({\bf r})=\sum_{q=1}^{K}\chi_{q}({\bf
r})V_{pq}^{\sigma,\mu}.$ (S23)

Here $\\{{\rm{V}}^{\sigma,\mu}\\}$ are a set of $K\times K$ matrices. Eq. S15
becomes

$\sum_{\mu=0}^{n}{\rm{V}}^{\sigma,\mu}{\rm{V}}^{\sigma,\mu T}={\rm{I}}_{K},$
(S24)

and

${\rm{V}}^{\sigma,\mu}{\rm{V}}^{\sigma,\nu
T}={\rm{V}}^{\sigma,\nu}{\rm{V}}^{\sigma,\mu T}.$ (S25)

Eq. S16 leads to
$\langle\chi_{p}|\chi_{q}\rangle=\sum_{r=1}^{K}\langle\chi_{p}|\varphi_{r}^{\sigma}\rangle\langle\varphi_{r}^{\sigma}|\chi_{q}\rangle=\delta_{pq}$,
such that

$\sum_{\mu=0}^{n}{\rm{V}}^{\sigma,\mu T}{\rm{V}}^{\sigma,\mu}={\rm{I}}_{K},$
(S26)

and

${\rm{V}}^{\sigma,\mu T}{\rm{V}}^{\sigma,\nu}={\rm{V}}^{\sigma,\nu
T}{\rm{V}}^{\sigma,\mu}.$ (S27)

${\rm{I}}_{K}$ denotes the $K\times K$ identity matrix.

Here we define two $K\times K$ symmetric matrices ${\rm{D}}^{\sigma}$ and
${\rm{E}}^{\sigma}$,

${\rm{D}}^{\sigma}=\sum_{\mu=0}^{n}{\rm{V}}^{\sigma,\mu
T}{\rm{I}}_{K}^{N_{\sigma}}{\rm{V}}^{\sigma,\mu};{\rm{E}}^{\sigma}=\sum_{\mu=0}^{n}{\rm{V}}^{\sigma,\mu
T}{\rm{I}}_{K}{\rm{V}}^{\sigma,\mu T}=\sum_{\mu=0}^{n}{\rm{V}}^{\sigma,\mu
T}{\rm{V}}^{\sigma,\mu T}={\rm{I}}_{K}.$ (S28)

${\rm{I}}_{K}^{N_{\sigma}}$ is a $K\times K$ diagonal matrix, with the first
$N_{\sigma}$ diagonal elements being 1 and the rest being 0.
${\rm{D}}^{\sigma}$ can be diagonalized by an orthogonal matrix
${\rm{U}}^{\sigma}$, thus

${\rm{D}}^{\sigma}={\rm{U}}^{\sigma}{\rm{\Lambda}}^{\sigma}{\rm{U}}^{\sigma
T},$ (S29)

where ${\rm{\Lambda}}^{\sigma}$ is a diagonal matrix with diagonal elements
$\\{\lambda_{p}^{\sigma}\\}$ being the eigenvalues of ${\rm{D}}^{\sigma}$. The
positive semi-definiteness of ${\rm{D}}^{\sigma}$ guarantees that all the
eigenvalues are not negative, with the maximal eigenvalue smaller than that of
${\rm{E}}^{\sigma}$. Thereby,

$0\leq\lambda_{p}^{\sigma}\leq
1;\sum_{p=1}^{K}\lambda_{p}^{\sigma}={\rm{Tr}}\left({\rm{D}}^{\sigma}\right)={\rm{Tr}}\left(\sum_{\mu=0}^{n}{\rm{V}}^{\sigma,\mu
T}{\rm{I}}_{K}^{N_{\sigma}}{\rm{V}}^{\sigma,\mu}\right)={\rm{Tr}}\left({\rm{I}}_{K}^{N_{\sigma}}\sum_{\mu=0}^{n}{\rm{V}}^{\sigma,\mu}{\rm{V}}^{\sigma,\mu
T}\right)=N_{\sigma}.$ (S30)

The density and kinetic energy can be formulated as

$\rho_{\sigma}({\bf r})=\sum_{p,q=1}^{K}\chi_{p}({\bf
r})D_{pq}^{\sigma}\chi_{q}({\bf
r})=\sum_{p=1}^{K}\lambda_{p}^{\sigma}|\chi_{p}^{\sigma}({\bf r})|^{2},$ (S31)

and

$T[\rho_{\sigma}]=-\frac{1}{2}\sum_{\sigma}^{\alpha,\beta}\sum_{p,q=1}^{K}D_{pq}^{\sigma}\langle\chi_{p}|\nabla^{2}|\chi_{q}\rangle=-\frac{1}{2}\sum_{\sigma}^{\alpha,\beta}\sum_{p=1}^{K}\lambda_{p}^{\sigma}\langle\chi_{p}^{\sigma}|\nabla^{2}|\chi_{p}^{\sigma}\rangle,$
(S32)

where $\chi_{p}^{\sigma}({\bf r})=\sum_{q=1}^{K}\chi_{q}({\bf
r})U_{qp}^{\sigma}$.

For $K\leq n+1$, to show that $\\{\lambda_{p}^{\sigma}\\}$ can be any values
that satisfy Eq. S30, here we construct a special set of orbitals
$\\{\varphi_{p}^{\sigma}\\}$ on any given set of orthonormal functions
$\\{\chi_{p}^{\sigma}\\}$, which take the following form

$\varphi_{p}^{\sigma}({\bf
r})=\sum_{\mu=0}^{n}\phi^{\sigma,\mu}_{p}e_{\mu}=\sum_{\mu=0}^{K-1}\chi_{\mu+1}^{\sigma}V_{p\mu+1}^{\sigma,\mu}e_{\mu},$
(S33)

i.e. the components corresponding to the same $e_{\mu}$ in all the orbitals
$\\{\varphi_{p}^{\sigma}\\}$ are formed by only the function
$\chi_{\mu+1}^{\sigma}$, thereby

$\displaystyle V_{pq}^{\sigma,\mu}=\begin{cases}W_{pq}^{\sigma},&q=\mu+1\\\
0,&q\neq\mu+1,\end{cases}$ (S34)

where ${\rm{W}}^{\sigma}$ is a $K\times K$ matrix. Now the conditions of Eqs.
S24-S27 become

${\rm{W}}^{\sigma}{\rm{W}}^{\sigma T}={\rm{I}}_{K},$ (S35)
$W_{k,\mu+1}W_{l,\nu+1}\delta_{\mu\nu}=W_{k,\nu+1}W_{l,\mu+1}\delta_{\mu\nu},$
(S36) $[{\rm{W}}^{\sigma T}{\rm{W}}^{\sigma}]_{\mu+1,\mu+1}=1,$ (S37)
$[{\rm{W}}^{\sigma T}{\rm{W}}^{\sigma}]_{\mu+1,\nu+1}=[{\rm{W}}^{\sigma
T}{\rm{W}}^{\sigma}]_{\nu+1,\mu+1}.$ (S38)

Eqs. S35-S38 can be satisfied once ${\rm{W}}^{\sigma}$ is an orthogonal
matrix. Now, the density and the kinetic energy becomes

$\rho_{\sigma}({\bf
r})=\sum_{p=1}^{K}\sum_{k=1}^{N_{\sigma}}(W_{kp}^{\sigma})^{2}|\chi_{p}^{\sigma}({\bf
r})|^{2},$ (S39)

and

$T[\rho_{\sigma}]=-\frac{1}{2}\sum_{\sigma}^{\alpha,\beta}\sum_{p=1}^{K}\sum_{k=1}^{N_{\sigma}}(W_{kp}^{\sigma})^{2}\langle\chi_{p}^{\sigma}|\nabla^{2}|\chi_{p}^{\sigma}\rangle.$
(S40)

Hence $\lambda_{p}^{\sigma}=\sum_{k=1}^{N_{\sigma}}(W_{kp}^{\sigma})^{2}$.
Because ${\rm{W}}^{\sigma}$ can be any unitary matrix,
$\\{\lambda_{p}^{\sigma}\\}$ can be any values subject to

$0\leq\lambda_{p}^{\sigma}\leq
1;\sum_{p=1}^{K}\lambda_{p}^{\sigma}=N_{\sigma}.$ (S41)

When $n=0$, the orbitals of Eq. S20 are real and the hypercomplex KS (HCKS)
reduces the conventional KS method; when $n=1$, the orbitals of Eq. S20 are
complex, which are used in the complex, spin-restricted KS (CRKS) method [23].
Take the singlet state of C atom as an example to compare KS, CRKS and HCKS.
The three methods all lead to fully occupied 1s and 2s sub-shells that are
lower in energy than the 2p sub-shell. For the 2p sub-shell, the spin-
restricted KS (RKS) has one doubly occupied p orbital; CRKS leads to two half
filled p orbitals as Eq. S20 with $n=1$ has only two components for each
occupied orbitals; HCKS automatically converges to the spin-restricted
solution, with each p orbital 1/3 filled. Therefore, the density by HCKS can
maintain both spin and space symmetry. Note that in most cases especially when
the energy gap of the frontier orbitals is big, HCKS in use of commonly used
functionals gives the same results as KS.


# Theoretical derivation of diffusion-tensor coefficients for the transport of
charged particles in magnetic fields

Olivier Deligny

###### Abstract

The transport of charged particles in various astrophysical environments
permeated by magnetic fields is described in terms of a diffusion process,
which relies on diffusion-tensor parameters generally inferred from Monte-
Carlo simulations. Based on a red-noise approximation to model the two-point
correlation function of the magnetic field experienced by charged particles
between two successive times, the diffusion-tensor coefficients were
previously derived in the case of pure turbulence. In this contribution to
ECRS2022, the derivation is extended to the case of a mean field on top of the
turbulence. The results are applicable to a variety of astrophysical
environments in regimes where the Larmor radius of the particles is resonant
with the power spectrum of the turbulence wavelength (gyro-resonant regime),
or where the Larmor radius is greater than the largest turbulence wavelength
(high-rigidity regime).

In many astrophysical environments, the propagation and acceleration of high-
energy charged particles (cosmic rays) are governed by the scattering off
magnetic fields. The transport of the particles is then modelled as an
anisotropic diffusion process. Under very broad conditions, the coefficients
of the diffusion tensor can be related to the velocity correlation function of
cosmic rays, $\langle v_{0i}v_{j}(t)\rangle$, through a time integration [1],

$D_{ij}(t)=\int_{0}^{t}\mathrm{d}t^{\prime}~{}\langle
v_{0i}v_{j}(t^{\prime})\rangle,$ (1)

in the limit that $t\rightarrow\infty$. Here, $v_{0i}\equiv
v_{i}(t=0)$111Since cosmic rays are high-energy relativistic particles, the
norm of the velocity is identified to $c$ for convenience.  and
$\langle\cdot\rangle$ stands for the average quantities, taken over several
space and time correlation scales of the turbulent field. Many estimates of
these coefficients have been made from numerical simulations exploring wide
ranges of particle rigidities and turbulence levels [e.g. 2, 3, 4, 5, 6, 7, 8,
9, 10]. In this contribution, we present the main steps of a theoretical
derivation of the velocity correlation functions consistent with simulation
results in a range of rigidities gyro-resonant with the power spectrum of
turbulence. Hence, this extends results obtained in [11] in the high rigidity
regime, and those in [12] in the gyro-resonant regime limited to pure
turbulence. Without loss of generalities, the study is limited to the example
of an isotropic 3D turbulence following a Kolmogorov power spectrum without
helicity denoted as $\delta\mathbf{B}$, while the mean field is denoted as
$\mathbf{B}_{0}=B_{0}\mathbf{u}_{z}$. The turbulence level is defined as
$\eta=\delta B^{2}/(\delta B^{2}+B_{0}^{2})$.

Due to the stochastic nature of the turbulence, the velocity of the particles
is a stochastic variable as well. We are thus interested in determining the
moments of $v_{i}(t)$ from the Lorentz-Newton equation of motion for the
particles,

$\dot{v}_{i}(t)=\delta\Omega(t)~{}\epsilon_{ijk}v_{j}(t)\delta
b_{k}(t)+\Omega_{0}~{}\epsilon_{ijk}v_{j}(t)b_{0k}(t).$ (2)

Here, $\delta\Omega(t)=c^{2}Z|e|\delta B(t)/E$ is the gyrofrequency related to
the turbulence, $\Omega_{0}$ that related to the mean field, $Z|e|$ the
electric charge, $E$ the energy of the particle, and $\delta
b_{k}(t)\equiv\delta b_{k}(\mathbf{x}(t))$ the $k$-th component of the
turbulence (expressed in units of $\delta B$) at the spatial coordinate
$\mathbf{x}(t)$ of the particle at time $t$. A formal solution for $\langle
v_{i}(t)\rangle$ can be obtained by expressing the solution of Eqn. 2 as an
infinite number of Dyson series, each combining terms in powers of
$\delta\mathbf{b}$ coupled to terms in powers of $\mathbf{B}_{0}$.

Dealing with such an infinite number of Dyson series is however hardly
manageable. To circumvent this, we use the auxiliary variable introduced in
[11], $w_{i}(t)=R_{ij}(t)v_{j}(t)$, with $\mathbf{R}(t)$ the rotation matrix
of angle $\Omega_{0}t$ around $\mathbf{u}_{z}$. The equation of motion for
$\mathbf{w}$ is then

$\dot{w}_{i}(t)=\delta\Omega(t)R^{-1}_{ij}\epsilon_{jkl}\delta
b_{l}(t)R_{km}w_{m}(t),$ (3)

the formal solution of which can be expressed as a single Dyson series:

$\displaystyle\langle w_{i_{0}}(t)\rangle$ $\displaystyle=$ $\displaystyle
w_{0i_{0}}+\sum_{p=1}^{\infty}\delta\Omega^{p}~{}\epsilon_{k_{1}m_{1}n_{1}}\epsilon_{k_{2}m_{2}n_{2}}\dots\epsilon_{k_{p}m_{p}n_{p}}w_{0i_{p}}\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}\dots\int_{0}^{t_{p-1}}\mathrm{d}t_{p}$
(4) $\displaystyle\times$ $\displaystyle
R^{-1}_{i_{0}k_{1}}(t_{1})R^{-1}_{i_{1}k_{2}}(t_{2})\cdots
R^{-1}_{i_{p-1}k_{p}}(t_{p})R_{m_{1}i_{1}}(t_{1})\cdots
R_{m_{p}i_{p}}(t_{p})\langle\delta b_{n_{1}}(t_{1})\dots\delta
b_{n_{p}}(t_{p})\rangle.$

In the following, we derive the velocity correlation function parallel to the
mean field, $\langle v_{0z}v_{z}(t)\rangle$, based on this equation. The
perpendicular (and anti-symmetric) functions can be obtained in a similar way
and will be detailed in a forthcoming publication [13].

In the Gaussian approximation, the Wick theorem allows for expressing the
expectation value $\langle\delta b_{n_{1}}(t_{1})\dots\delta
b_{n_{p}}(t_{p})\rangle$ in terms of all possible permutations of products of
contractions of pairs of $\langle\delta b_{n_{i}}(t_{i})\delta
b_{n_{j}}(t_{j})\rangle$, which can be, in the case of 3D isotropic
turbulence, written as

$\langle\delta b_{n_{1}}(t_{i})\delta
b_{n_{2}}(t_{j})\rangle=\frac{\delta_{n_{1}n_{2}}}{3}\varphi(t_{i}-t_{j}).$
(5)

The correlation function $\varphi(t)$, which describes the correlation of the
turbulence experienced by a particle along its path at two different times, is
modeled in this study as a red-noise process with parameter $\tau$,

$\varphi(t)=\exp{(-t/\tau)}.$ (6)

The expression of $\tau$, similarly to that found in [3], depends on the
regime of rigidity considered. With $\rho$ the Larmor radius of the particle
expressed in units of the largest eddy scale $L_{\mathrm{max}}$ of the
turbulence, $\tau\simeq L_{\mathrm{c}}/c$ for $\rho\gtrsim\pi
L_{\mathrm{c}}/L_{\mathrm{max}}$, where $L_{\mathrm{c}}$ is the coherence
scale of the turbulence. This is because in this rigidity regime, particles
can travel over a distance $L_{\mathrm{c}}$ undergoing small deflections only.
On the other hand, for $\rho\lesssim\pi L_{\mathrm{c}}/L_{\mathrm{max}}$, the
following heuristic estimate that makes use of the kinetic energy spectrum of
the turbulence $\mathcal{E}(k)$ is observed to reproduce simulations:

$\tau\simeq\frac{2}{c}\frac{\int_{k_{\star}}^{k_{\mathrm{max}}}\mathrm{d}kk^{-1}\mathcal{E}(k)}{\int_{k_{\star}}^{k_{\mathrm{max}}}\mathrm{d}k\mathcal{E}(k)}.$
(7)

In this regime of rigidity, $\tau$ inherits a $\rho$ dependency from that of
the lower boundary in wave number
$k_{\star}(\rho)=\rho_{\star}k_{\mathrm{min}}/\rho$ with
$\rho_{\star}=2L_{\mathrm{c}}/(\pi L_{\mathrm{max}})$. The truncation in the
wavenumber integration range selects modes for which particles do not
experience spiral motions around the corresponding large-scale magnetic field
lines over several Larmor times, modes that hence prevent decorrelations from
occurring on relevant time scales.

To carry out a summation of the infinite Dyson series (Eqn. 4), we resort to a
two-step iteration procedure. First, the same partial summation scheme as in
[12] is used. The classes of “diagrams” retained in the scheme are of two
kinds: unconnected and nested ones [14]. This corresponds formally to the
Kraichnan propagator [15]. In fact, it has been shown that the summation of
the first two terms (beyond the free propagator) is sufficient to give us a
physical solution in the case of pure turbulence [12]. We follow here the same
strategy so that, making use of the properties of the Levi-Civita symbols
contracted over one index, Eqn. 4 is substituted by

$\displaystyle\langle w^{\mathrm{K}}_{z}(t)\rangle$ $\displaystyle=$
$\displaystyle
w_{0z}+\frac{\delta\Omega^{2}}{3}\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}~{}R^{-1}_{zk}(t_{1})R_{mi_{1}}(t_{1})R^{-1}_{i_{1}k}(t_{2})R_{mi_{2}}(t_{2})\varphi(t_{1}-t_{2})\langle
w^{\mathrm{K}}_{z}(t-t_{1})\rangle$ (8) $\displaystyle-$
$\displaystyle\frac{\delta\Omega^{2}}{3}\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}~{}R^{-1}_{zk_{1}}(t_{1})R_{k_{2}i_{1}}(t_{1})R^{-1}_{i_{1}k_{2}}(t_{2})R_{k_{1}i_{2}}(t_{2})\varphi(t_{1}-t_{2})\langle
w^{\mathrm{K}}_{z}(t-t_{1})\rangle$ $\displaystyle+$
$\displaystyle\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\int_{0}^{t}\mathrm{d}t_{1}\cdots\int_{0}^{t_{3}}\mathrm{d}t_{4}~{}\varphi(t_{1}-t_{4})\varphi(t_{2}-t_{3})\langle
w^{\mathrm{K}}_{z}(t-t_{1})\rangle\times$
$\displaystyle(R^{-1}_{zk_{1}}(t_{1})R_{m_{1}i_{1}}(t_{1})R^{-1}_{i_{1}k_{2}}(t_{2})R_{m_{2}i_{2}}(t_{2})R^{-1}_{i_{2}k_{2}}(t_{3})R_{m_{2}i_{3}}(t_{3})R^{-1}_{i_{3}k_{1}}(t_{4})R_{m_{1}z}(t_{4})$
$\displaystyle-$ $\displaystyle
R^{-1}_{zk_{1}}(t_{1})R_{m_{1}i_{1}}(t_{1})R^{-1}_{i_{1}k_{2}}(t_{2})R_{m_{2}i_{2}}(t_{2})R^{-1}_{i_{2}m_{2}}(t_{3})R_{k_{2}i_{3}}(t_{3})R^{-1}_{i_{3}k_{1}}(t_{4})R_{m_{1}z}(t_{4})$
$\displaystyle-$ $\displaystyle
R^{-1}_{zk_{1}}(t_{1})R_{m_{1}i_{1}}(t_{1})R^{-1}_{i_{1}k_{2}}(t_{2})R_{m_{2}i_{2}}(t_{2})R^{-1}_{i_{2}k_{2}}(t_{3})R_{m_{2}i_{3}}(t_{3})R^{-1}_{i_{3}m_{1}}(t_{4})R_{k_{1}z}(t_{4})$
$\displaystyle+$ $\displaystyle
R^{-1}_{zk_{1}}(t_{1})R_{m_{1}i_{1}}(t_{1})R^{-1}_{i_{1}k_{2}}(t_{2})R_{m_{2}i_{2}}(t_{2})R^{-1}_{i_{2}m_{2}}(t_{3})R_{k_{2}i_{3}}(t_{3})R^{-1}_{i_{3}m_{1}}(t_{4})R_{k_{1}z}(t_{4})),$

where the superscript K stands for “Kraichnan”. Next, some properties of the
rotation matrices, in particular $R^{-1}_{zk}=R_{zk}=\delta_{zk}$,
$R_{ij}(t_{1})R^{-1}_{jk}(t_{2})=R_{ik}(t_{1}-t_{2})$, and
$R_{ij}(t_{1}-t_{2})=R_{ji}(t_{2}-t_{1})$), allow us to get an explicit form
of the equation:

$\displaystyle\langle w^{\mathrm{K}}_{z}(t)\rangle$ $\displaystyle=$
$\displaystyle
w_{0z}-\frac{\delta\Omega^{2}}{3}\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t_{1}}\mathrm{d}t_{2}~{}\varphi(t_{1}-t_{2})\cos{(\Omega_{0}(t_{1}-t_{2}))}\langle
w^{\mathrm{K}}_{z}(t-t_{1})\rangle$ (9) $\displaystyle+$ $\displaystyle
4\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\int_{0}^{t}\mathrm{d}t_{1}\cdots\int_{0}^{t_{3}}\mathrm{d}t_{4}~{}\varphi(t_{1}-t_{4})\varphi(t_{2}-t_{3})\cos{(\Omega_{0}(t_{2}-t_{3}))}\cos{(\Omega_{0}(t_{1}-t_{2}+t_{3}-t_{4}))}\langle
w^{\mathrm{K}}_{z}(t-t_{1})\rangle$ $\displaystyle-$ $\displaystyle
2\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\int_{0}^{t}\mathrm{d}t_{1}\cdots\int_{0}^{t_{3}}\mathrm{d}t_{4}~{}\varphi(t_{1}-t_{4})\varphi(t_{2}-t_{3})\cos{(\Omega_{0}(t_{1}-2t_{2}+2t_{3}-t_{4}))}\langle
w^{\mathrm{K}}_{z}(t-t_{1})\rangle$ $\displaystyle+$ $\displaystyle
2\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\int_{0}^{t}\mathrm{d}t_{1}\cdots\int_{0}^{t_{3}}\mathrm{d}t_{4}~{}\varphi(t_{1}-t_{4})\varphi(t_{2}-t_{3})\cos{(\Omega_{0}(t_{1}-t_{2}+t_{3}-t_{4}))}\langle
w^{\mathrm{K}}_{z}(t-t_{1})\rangle.$

To solve this non-linear equation, we proceed with a Laplace transform. The
change of variables $t=x+x_{1}+\cdots+x_{p}$, $t_{1}=x_{1}+\cdots+x_{p}$,
$t_{2}=x_{2}+\cdots+x_{p}$,$\cdots$, $t_{p}=x_{p}$ allows for sending all
integration boundaries between 0 and $+\infty$ for the $x_{i}$ variables.
After some algebra, the equation for the Laplace transform
$\mathcal{L}\left(w(t)\right)\equiv\hat{W}(p)$ reads as

$\displaystyle\hat{W}^{\mathrm{K}}(p)$ $\displaystyle=$
$\displaystyle\frac{1}{p}-2\frac{\delta\Omega^{2}}{3}\frac{\hat{W}^{\mathrm{K}}(p)}{p}\frac{p+\tau^{-1}}{(p+\tau^{-1})^{2}+\Omega_{0}^{2}}$
(10) $\displaystyle+$ $\displaystyle
2\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\frac{\hat{W}^{\mathrm{K}}(p)}{p}\left(\frac{p+2\tau^{-1}}{(p+2\tau^{-1})^{2}+\Omega_{0}^{2}}\right)\left(\frac{(p+\tau^{-1})^{2}-\Omega_{0}^{2}}{\left((p+\tau^{-1})^{2}+\Omega_{0}^{2}\right)^{2}}\right)$
$\displaystyle+$ $\displaystyle
2\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\frac{\hat{W}^{\mathrm{K}}(p)}{p}\left(\frac{1}{p+2\tau^{-1}}\right)\left(\frac{(p+\tau^{-1})^{2}-\Omega_{0}^{2}}{\left((p+\tau^{-1})^{2}+\Omega_{0}^{2}\right)^{2}}\right)$
$\displaystyle-$ $\displaystyle
4\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\frac{\hat{W}^{\mathrm{K}}(p)}{p}\left(\frac{\Omega_{0}^{2}(p+\tau^{-1})}{\left((p+2\tau^{-1})^{2}+\Omega_{0}^{2}\right)\left((p+\tau^{-1})^{2}+\Omega_{0}^{2}\right)^{2}}\right).$

The solution for $\langle w^{\mathrm{K}}_{z}(t)\rangle$ is then obtained by
making use of the numerical Stehfest scheme of the inverse Laplace transform.

Once $\langle w^{\mathrm{K}}_{z}(t)\rangle$ is determined, the second step to
get an improved estimation of the propagator consists in including “crossed”
diagrams (that formally account for any mix of crossed and nested diagrams).
To do so, an iterative procedure is designed to carry out this summation. In
the Laplace space, the $N$th iterated $\langle w^{\mathrm{N}}_{z}(t)\rangle$
reads as [13]

$\displaystyle\hat{W}^{N}(p)$ $\displaystyle=$
$\displaystyle\frac{1}{p}-2\frac{\delta\Omega^{2}}{3}\frac{\hat{W}^{N}(p)}{p}\mathcal{L}\left(\exp{(-t/\tau)w^{N-1}(t)\cos{\Omega_{0}t}}\right)$
(11) $\displaystyle+$ $\displaystyle
2\left(\frac{\delta\Omega^{2}}{3}\right)^{2}\frac{\hat{W}^{N}(p)}{p}\mathcal{L}\left(\exp{(-2t/\tau)w^{N-1}(t)}\right)\mathcal{L}^{2}\left(\exp{(-t/\tau)w^{N-1}(t)\cos{\Omega_{0}t}}\right),$

with the initial iteration $w^{0}(t)=w^{\mathrm{K}}(t)$. In practice,
convergence is achieved after a few iterations ($N\simeq 5$).

Figure 1: Time dependence of the auto-correlation of the particle velocities
parallel to the mean field, for different values of turbulence level. A
reduced rigidity $\rho=0.03$ is chosen to illustrate the gyro-resonant regime.

The resulting picture is illustrated in Fig. 1, where the time dependence of
the auto-correlation of the particle velocities parallel to the mean field is
shown for different values of turbulence level. The reduced rigidity is chosen
to be $\rho=0.03$ so as to explore the gyro-resonant regime. The timescale of
the correlation is observed to be minimal in the case of pure turbulence
($\eta=1$) and to tend to infinity in the case of low turbulence. The various
structures beyond the exponential falloff that have been uncovered in most of
the Monte-Carlo simulations mentioned in the introduction are reproduced by
the calculation presented in this contribution.

## References

  * [1] R. Kubo, _Statistical mechanical theory of irreversible processes. 1. General theory and simple applications in magnetic and conduction problems_ , _J. Phys. Soc. Jap._ 12 (1957) 570.
  * [2] J. Giacalone and J.R. Jokipii, _The Transport of Cosmic Rays across a Turbulent Magnetic Field_ , _Astrophys. J._ 520 (1999) 204.
  * [3] F. Casse, M. Lemoine and G. Pelletier, _Transport of cosmic rays in chaotic magnetic fields_ , _Phys. Rev. D_ 65 (2001) 023002.
  * [4] J. Candia and E. Roulet, _Diffusion and drift of cosmic rays in highly turbulent magnetic fields_ , _JCAP_ 10 (2004) 007 [astro-ph/0408054].
  * [5] T. Hauff, F. Jenko, A. Shalchi and R. Schlickeiser, _Scaling Theory for Cross-Field Transport of Cosmic Rays in Turbulent Fields_ , _Astrophys. J._ 711 (2010) 997.
  * [6] F. Fraschetti and J. Giacalone, _Early-time velocity autocorrelation for charged particles diffusion and drift in static magnetic turbulence_ , _Astrophys. J._ 755 (2012) 114 [1206.6494].
  * [7] M. Fatuzzo and F. Melia, _A Numerical Assessment of Cosmic-ray Energy Diffusion through Turbulent Media_ , _Astrophys. J._ 784 (2014) 131 [1402.5469].
  * [8] A.P. Snodin, A. Shukurov, G.R. Sarson, P.J. Bushby and L.F.S. Rodrigues, _Global diffusion of cosmic rays in random magnetic fields_ , _Mon. Not. Roy. Astron. Soc._ 457 (2016) 3975 [1509.03766].
  * [9] P. Reichherzer, J. Becker Tjus, E.G. Zweibel, L. Merten and M.J. Pueschel, _Turbulence-Level Dependence of Cosmic-Ray Parallel Diffusion_ , _Mon. Not. Roy. Astron. Soc._ 498 (2020) 5051 [1910.07528].
  * [10] P. Reichherzer, L. Merten, J. Dörner, J. Becker Tjus, M.J. Pueschel and E.G. Zweibel, _Regimes of cosmic-ray diffusion in Galactic turbulence_ , _Appl. Sciences_ 4 (2022) 15 [2104.13093].
  * [11] I. Plotnikov, G. Pelletier and M. Lemoine, _Particle Transport in intense small scale magnetic turbulence with a mean field_ , _Astron. Astrophys._ 532 (2011) A68 [1105.0618].
  * [12] O. Deligny, _Transport of Charged Particles Propagating in Turbulent Magnetic Fields as a Red-noise Process_ , _Astrophys. J._ 920 (2021) 87 [2107.04391].
  * [13] O. Deligny, _in preparation_ , .
  * [14] U. Frisch, _La propagation des ondes en milieu aléatoire et les équations stochastiques_ , _Annales d’Astrophysique_ 29 (1966) 645.
  * [15] R.H. Kraichnan, _Dynamics of nonlinear stochastic systems_ , _Journal of Mathematical Physics_ 2 (1961) 124.


8 | 0.1502 | 0.4871 | 72 | 0.0001 | 0.0757 | 136 | 0.0000 | 0.0385 | 200 | 0.0001 | 0.0385
10 | 0.1172 | 0.4871 | 74 | 0.0000 | 0.0757 | 138 | 0.0000 | 0.0385 | 202 | 0.0000 | 0.0385
12 | 0.2710 | 0.4871 | 76 | 0.0002 | 0.0757 | 140 | 0.0001 | 0.0385 | 204 | 0.0001 | 0.0385
14 | 0.3076 | 0.4871 | 78 | 0.0005 | 0.0757 | 142 | 0.0001 | 0.0385 | 206 | 0.0002 | 0.0385
16 | 0.0387 | 0.2758 | 80 | 0.0000 | 0.0757 | 144 | 0.0000 | 0.0385 | 208 | 0.0000 | 0.0385
18 | 0.0293 | 0.2758 | 82 | 0.0000 | 0.0757 | 146 | 0.0000 | 0.0385 | 210 | 0.0000 | 0.0385
20 | 0.0073 | 0.2758 | 84 | 0.0000 | 0.0757 | 148 | 0.0000 | 0.0385 | 212 | 0.0000 | 0.0385
22 | 0.0146 | 0.2758 | 86 | 0.0000 | 0.0757 | 150 | 0.0000 | 0.0385 | 214 | 0.0000 | 0.0385
24 | 0.0755 | 0.2758 | 88 | 0.0001 | 0.0757 | 152 | 0.0000 | 0.0385 | 216 | 0.0000 | 0.0385
26 | 0.0586 | 0.2758 | 90 | 0.0000 | 0.0757 | 154 | 0.0000 | 0.0385 | 218 | 0.0000 | 0.0385
28 | 0.1364 | 0.2758 | 92 | 0.0001 | 0.0757 | 156 | 0.0000 | 0.0385 | 220 | 0.0001 | 0.0385
30 | 0.1556 | 0.2758 | 94 | 0.0002 | 0.0757 | 158 | 0.0001 | 0.0385 | 222 | 0.0001 | 0.0385
32 | 0.0097 | 0.1468 | 96 | 0.0049 | 0.0757 | 160 | 0.0000 | 0.0385 | 224 | 0.0024 | 0.0385
34 | 0.0073 | 0.1468 | 98 | 0.0037 | 0.0757 | 162 | 0.0000 | 0.0385 | 226 | 0.0018 | 0.0385
36 | 0.0018 | 0.1468 | 100 | 0.0009 | 0.0757 | 164 | 0.0000 | 0.0385 | 228 | 0.0005 | 0.0385
38 | 0.0037 | 0.1468 | 102 | 0.0018 | 0.0757 | 166 | 0.0000 | 0.0385 | 230 | 0.0009 | 0.0385
40 | 0.0005 | 0.1468 | 104 | 0.0002 | 0.0757 | 168 | 0.0000 | 0.0385 | 232 | 0.0001 | 0.0385
42 | 0.0000 | 0.1468 | 106 | 0.0000 | 0.0757 | 170 | 0.0000 | 0.0385 | 234 | 0.0000 | 0.0385
44 | 0.0009 | 0.1468 | 108 | 0.0005 | 0.0757 | 172 | 0.0000 | 0.0385 | 236 | 0.0002 | 0.0385
46 | 0.0018 | 0.1468 | 110 | 0.0009 | 0.0757 | 174 | 0.0000 | 0.0385 | 238 | 0.0005 | 0.0385
48 | 0.0194 | 0.1468 | 112 | 0.0097 | 0.0757 | 176 | 0.0000 | 0.0385 | 240 | 0.0048 | 0.0385
50 | 0.0146 | 0.1468 | 114 | 0.0073 | 0.0757 | 178 | 0.0000 | 0.0385 | 242 | 0.0037 | 0.0385
52 | 0.0037 | 0.1468 | 116 | 0.0018 | 0.0757 | 180 | 0.0000 | 0.0385 | 244 | 0.0009 | 0.0385
54 | 0.0073 | 0.1468 | 118 | 0.0037 | 0.0757 | 182 | 0.0000 | 0.0385 | 246 | 0.0018 | 0.0385
56 | 0.0378 | 0.1468 | 120 | 0.0189 | 0.0757 | 184 | 0.0000 | 0.0385 | 248 | 0.0095 | 0.0385
58 | 0.0293 | 0.1468 | 122 | 0.0146 | 0.0757 | 186 | 0.0000 | 0.0385 | 250 | 0.0073 | 0.0385
60 | 0.0683 | 0.1468 | 124 | 0.0342 | 0.0757 | 188 | 0.0000 | 0.0385 | 252 | 0.0171 | 0.0385
62 | 0.0781 | 0.1468 | 126 | 0.0391 | 0.0757 | 190 | 0.0000 | 0.0385 | 254 | 0.0195 | 0.0385
Table 9: Calculated sequency coefficients, $\tilde{\beta}_{\nu}$,and upper
bounds, $\tilde{B}_{\nu}$, for a $x^{4}$ potential, as shown in Fig. 7.

|  # $\delta t$ | 0.1 | 0.13 | 0.16 | 0.19 | 0.22 | 0.25 | 0.28 | 0.31 | 0.34 | 0.37 | 0.4 | 0.43 | 0.46 | 0.49 | 0.52 | 0.55
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Full | 1 | 0.971 | 0.970 | 0.968 | 0.966 | 0.964 | 0.961 | 0.958 | 0.954 | 0.950 | 0.945 | 0.940 | 0.934 | 0.928 | 0.921 | 0.914 | 0.906
2 | 0.970 | 0.968 | 0.965 | 0.962 | 0.959 | 0.955 | 0.950 | 0.946 | 0.941 | 0.936 | 0.931 | 0.926 | 0.921 | 0.916 | 0.912 | 0.907
3 | 0.969 | 0.967 | 0.964 | 0.961 | 0.958 | 0.955 | 0.952 | 0.950 | 0.947 | 0.946 | 0.944 | 0.944 | 0.945 | 0.951 | 0.960 | 0.971
4 | 0.969 | 0.967 | 0.965 | 0.963 | 0.962 | 0.961 | 0.961 | 0.961 | 0.963 | 0.968 | 0.976 | 0.986 | 0.994 | 0.997 | 0.992 | 0.978
5 | 0.970 | 0.969 | 0.968 | 0.968 | 0.968 | 0.970 | 0.973 | 0.980 | 0.989 | 0.996 | 0.998 | 0.993 | 0.982 | 0.967 | 0.950 | 0.932
6 | 0.971 | 0.971 | 0.971 | 0.973 | 0.976 | 0.982 | 0.990 | 0.997 | 0.998 | 0.992 | 0.981 | 0.966 | 0.951 | 0.934 | 0.923 | 0.922
7 | 0.972 | 0.974 | 0.976 | 0.979 | 0.986 | 0.994 | 0.999 | 0.996 | 0.987 | 0.973 | 0.960 | 0.946 | 0.937 | 0.947 | 0.966 | 0.972
8 | 0.974 | 0.977 | 0.980 | 0.987 | 0.995 | 0.999 | 0.995 | 0.984 | 0.972 | 0.959 | 0.949 | 0.955 | 0.975 | 0.986 | 0.981 | 0.964
9 | 0.976 | 0.980 | 0.986 | 0.994 | 0.999 | 0.995 | 0.985 | 0.973 | 0.962 | 0.959 | 0.971 | 0.989 | 0.992 | 0.980 | 0.949 | 0.918
10 | 0.979 | 0.984 | 0.991 | 0.998 | 0.997 | 0.988 | 0.977 | 0.967 | 0.966 | 0.981 | 0.994 | 0.991 | 0.975 | 0.941 | 0.926 | 0.936
11 | 0.981 | 0.987 | 0.996 | 0.998 | 0.991 | 0.981 | 0.971 | 0.971 | 0.985 | 0.996 | 0.990 | 0.971 | 0.942 | 0.942 | 0.960 | 0.972
12 | 0.983 | 0.991 | 0.998 | 0.996 | 0.986 | 0.976 | 0.973 | 0.986 | 0.997 | 0.991 | 0.972 | 0.948 | 0.955 | 0.975 | 0.978 | 0.942
13 | 0.986 | 0.995 | 0.998 | 0.991 | 0.982 | 0.975 | 0.984 | 0.997 | 0.993 | 0.976 | 0.954 | 0.964 | 0.982 | 0.983 | 0.943 | 0.914
14 | 0.988 | 0.997 | 0.997 | 0.988 | 0.979 | 0.981 | 0.995 | 0.996 | 0.982 | 0.960 | 0.968 | 0.987 | 0.989 | 0.947 | 0.932 | 0.948

|  # $\delta t$ | 0.1 | 0.13 | 0.16 | 0.19 | 0.22 | 0.25 | 0.28 | 0.31 | 0.34 | 0.37 | 0.4 | 0.43 | 0.46 | 0.49 | 0.52 | 0.55
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
SeqHT | 1 | 0.971 | 0.969 | 0.968 | 0.965 | 0.963 | 0.960 | 0.956 | 0.952 | 0.948 | 0.943 | 0.937 | 0.931 | 0.924 | 0.917 | 0.909 | 0.901
2 | 0.969 | 0.967 | 0.964 | 0.960 | 0.957 | 0.952 | 0.947 | 0.942 | 0.937 | 0.931 | 0.926 | 0.922 | 0.914 | 0.906 | 0.906 | 0.906
3 | 0.968 | 0.966 | 0.963 | 0.959 | 0.956 | 0.951 | 0.947 | 0.945 | 0.943 | 0.940 | 0.936 | 0.938 | 0.941 | 0.941 | 0.957 | 0.965
4 | 0.968 | 0.966 | 0.963 | 0.959 | 0.959 | 0.958 | 0.956 | 0.955 | 0.956 | 0.963 | 0.971 | 0.982 | 0.989 | 0.982 | 0.983 | 0.968
5 | 0.968 | 0.967 | 0.965 | 0.963 | 0.964 | 0.965 | 0.967 | 0.977 | 0.987 | 0.990 | 0.988 | 0.985 | 0.976 | 0.945 | 0.935 | 0.925
6 | 0.968 | 0.969 | 0.968 | 0.967 | 0.972 | 0.978 | 0.985 | 0.990 | 0.993 | 0.986 | 0.968 | 0.960 | 0.940 | 0.908 | 0.916 | 0.916
7 | 0.969 | 0.971 | 0.972 | 0.974 | 0.983 | 0.991 | 0.992 | 0.988 | 0.983 | 0.965 | 0.938 | 0.944 | 0.931 | 0.912 | 0.955 | 0.962
8 | 0.971 | 0.974 | 0.977 | 0.981 | 0.992 | 0.998 | 0.989 | 0.972 | 0.962 | 0.947 | 0.928 | 0.949 | 0.968 | 0.945 | 0.963 | 0.963
9 | 0.972 | 0.977 | 0.983 | 0.988 | 0.996 | 0.992 | 0.977 | 0.963 | 0.954 | 0.945 | 0.951 | 0.975 | 0.984 | 0.926 | 0.938 | 0.923
10 | 0.974 | 0.981 | 0.989 | 0.994 | 0.993 | 0.984 | 0.967 | 0.958 | 0.961 | 0.971 | 0.973 | 0.978 | 0.962 | 0.877 | 0.939 | 0.932
11 | 0.976 | 0.985 | 0.994 | 0.993 | 0.988 | 0.977 | 0.964 | 0.961 | 0.980 | 0.985 | 0.966 | 0.954 | 0.932 | 0.884 | 0.970 | 0.956
12 | 0.978 | 0.989 | 0.997 | 0.990 | 0.981 | 0.973 | 0.962 | 0.977 | 0.988 | 0.977 | 0.937 | 0.935 | 0.943 | 0.902 | 0.950 | 0.944
13 | 0.980 | 0.993 | 0.997 | 0.984 | 0.976 | 0.972 | 0.975 | 0.980 | 0.982 | 0.954 | 0.929 | 0.955 | 0.969 | 0.889 | 0.928 | 0.905
14 | 0.983 | 0.996 | 0.996 | 0.979 | 0.973 | 0.978 | 0.984 | 0.984 | 0.972 | 0.946 | 0.950 | 0.972 | 0.964 | 0.865 | 0.929 | 0.938

Table 10: The fidelity of the adiabatically prepared $\lambda\phi^{4}$ ground
state as a function of $\delta t$ and the number of adiabatic steps, for
$n_{q}=5$ with $\phi_{\rm max}=4$ and $\lambda=10$. Each adiabatic step
consists of one first-order Trotter step. Results are shown in the top panels
of Fig. 9.

# $\delta t$ | 0.1 | 0.13 | 0.16 | 0.19 | 0.22 | 0.25 | 0.28 | 0.31 | 0.34 | 0.37 | 0.4 | 0.43 | 0.46 | 0.49 | 0.52 | 0.55 | 0.58 | 0.61 | 0.64
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
F | 1 | 0.973 | 0.974 | 0.975 | 0.975 | 0.976 | 0.977 | 0.978 | 0.979 | 0.980 | 0.982 | 0.983 | 0.984 | 0.985 | 0.987 | 0.987 | 0.988 | 0.988 | 0.988 | 0.988
2 | 0.975 | 0.976 | 0.977 | 0.979 | 0.980 | 0.982 | 0.984 | 0.986 | 0.988 | 0.989 | 0.991 | 0.992 | 0.993 | 0.994 | 0.994 | 0.995 | 0.994 | 0.992 | 0.989
3 | 0.976 | 0.978 | 0.980 | 0.982 | 0.985 | 0.987 | 0.989 | 0.991 | 0.993 | 0.995 | 0.996 | 0.997 | 0.998 | 0.997 | 0.995 | 0.992 | 0.988 | 0.983 | 0.979
4 | 0.978 | 0.980 | 0.983 | 0.986 | 0.989 | 0.991 | 0.993 | 0.995 | 0.997 | 0.998 | 0.998 | 0.996 | 0.995 | 0.992 | 0.991 | 0.989 | 0.985 | 0.982 | 0.976
5 | 0.979 | 0.983 | 0.986 | 0.989 | 0.992 | 0.995 | 0.997 | 0.998 | 0.998 | 0.997 | 0.996 | 0.995 | 0.994 | 0.992 | 0.990 | 0.988 | 0.987 | 0.984 | 0.976
6 | 0.981 | 0.985 | 0.989 | 0.992 | 0.995 | 0.997 | 0.998 | 0.997 | 0.996 | 0.996 | 0.996 | 0.995 | 0.994 | 0.993 | 0.992 | 0.989 | 0.987 | 0.983 | 0.968
7 | 0.983 | 0.987 | 0.991 | 0.995 | 0.997 | 0.997 | 0.997 | 0.996 | 0.996 | 0.996 | 0.996 | 0.996 | 0.996 | 0.994 | 0.991 | 0.988 | 0.986 | 0.979 | 0.961
8 | 0.985 | 0.989 | 0.993 | 0.996 | 0.997 | 0.997 | 0.996 | 0.996 | 0.996 | 0.997 | 0.998 | 0.998 | 0.995 | 0.992 | 0.988 | 0.986 | 0.986 | 0.977 | 0.961
9 | 0.987 | 0.991 | 0.995 | 0.997 | 0.997 | 0.997 | 0.996 | 0.997 | 0.997 | 0.999 | 0.999 | 0.997 | 0.994 | 0.990 | 0.988 | 0.987 | 0.984 | 0.975 | 0.952
10 | 0.988 | 0.993 | 0.997 | 0.997 | 0.997 | 0.997 | 0.997 | 0.998 | 0.999 | 0.999 | 0.998 | 0.996 | 0.992 | 0.990 | 0.988 | 0.989 | 0.987 | 0.973 | 0.948

# $\delta t$ | 0.1 | 0.13 | 0.16 | 0.19 | 0.22 | 0.25 | 0.28 | 0.31 | 0.34 | 0.37 | 0.4 | 0.43 | 0.46 | 0.49 | 0.52 | 0.55 | 0.58 | 0.61 | 0.64
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
T | 1 | 0.973 | 0.974 | 0.974 | 0.975 | 0.976 | 0.977 | 0.978 | 0.979 | 0.980 | 0.981 | 0.982 | 0.984 | 0.985 | 0.984 | 0.985 | 0.987 | 0.986 | 0.985 | 0.986
2 | 0.974 | 0.975 | 0.977 | 0.978 | 0.980 | 0.981 | 0.982 | 0.984 | 0.986 | 0.987 | 0.986 | 0.989 | 0.991 | 0.987 | 0.988 | 0.993 | 0.990 | 0.980 | 0.985
3 | 0.975 | 0.977 | 0.979 | 0.980 | 0.983 | 0.986 | 0.986 | 0.987 | 0.989 | 0.990 | 0.990 | 0.995 | 0.994 | 0.984 | 0.988 | 0.992 | 0.984 | 0.959 | 0.971
4 | 0.976 | 0.979 | 0.982 | 0.983 | 0.986 | 0.987 | 0.988 | 0.992 | 0.996 | 0.993 | 0.989 | 0.994 | 0.991 | 0.975 | 0.986 | 0.988 | 0.980 | 0.939 | 0.961
5 | 0.977 | 0.981 | 0.984 | 0.985 | 0.989 | 0.992 | 0.991 | 0.992 | 0.992 | 0.991 | 0.987 | 0.991 | 0.989 | 0.967 | 0.986 | 0.988 | 0.981 | 0.918 | 0.952
6 | 0.979 | 0.983 | 0.986 | 0.988 | 0.992 | 0.994 | 0.992 | 0.992 | 0.994 | 0.990 | 0.981 | 0.990 | 0.987 | 0.957 | 0.984 | 0.990 | 0.977 | 0.894 | 0.942
7 | 0.980 | 0.985 | 0.988 | 0.989 | 0.994 | 0.996 | 0.993 | 0.990 | 0.991 | 0.987 | 0.982 | 0.991 | 0.990 | 0.958 | 0.990 | 0.986 | 0.976 | 0.863 | 0.933
8 | 0.981 | 0.987 | 0.990 | 0.991 | 0.995 | 0.994 | 0.991 | 0.991 | 0.993 | 0.985 | 0.981 | 0.990 | 0.986 | 0.944 | 0.989 | 0.987 | 0.972 | 0.836 | 0.916
9 | 0.982 | 0.989 | 0.993 | 0.993 | 0.994 | 0.995 | 0.989 | 0.988 | 0.992 | 0.988 | 0.979 | 0.988 | 0.985 | 0.937 | 0.982 | 0.985 | 0.970 | 0.827 | 0.898
10 | 0.984 | 0.991 | 0.994 | 0.991 | 0.995 | 0.995 | 0.990 | 0.985 | 0.991 | 0.986 | 0.974 | 0.990 | 0.982 | 0.924 | 0.986 | 0.983 | 0.968 | 0.781 | 0.879

Table 11: The fidelity of the adiabatically prepared $\lambda\phi^{4}$ ground
state as a function of $\delta t$ and the number of adiabatic steps, for
$n_{q}=5$ with $\phi_{\rm max}=4$ and $\lambda=10$. Each adiabatic step
consists of one second-order Trotter step. Results are shown in the bottom
panels of Fig. 9.

|  # $\delta t$ | 0.1 | 0.15 | 0.2 | 0.25 | 0.3 | 0.35 | 0.4 | 0.45 | 0.5 | 0.55
---|---|---|---|---|---|---|---|---|---|---|---
Full | 1 | 0.973 | 0.974 | 0.975 | 0.977 | 0.978 | 0.980 | 0.982 | 0.984 | 0.985 | 0.986
2 | 0.974 | 0.976 | 0.978 | 0.980 | 0.983 | 0.985 | 0.986 | 0.987 | 0.988 | 0.990
3 | 0.975 | 0.978 | 0.981 | 0.984 | 0.986 | 0.989 | 0.990 | 0.990 | 0.988 | 0.986
4 | 0.977 | 0.980 | 0.983 | 0.986 | 0.989 | 0.992 | 0.992 | 0.985 | 0.982 | 0.982
5 | 0.978 | 0.982 | 0.985 | 0.989 | 0.992 | 0.991 | 0.986 | 0.981 | 0.979 | 0.980
6 | 0.979 | 0.984 | 0.988 | 0.992 | 0.992 | 0.989 | 0.987 | 0.977 | 0.976 | 0.979
7 | 0.981 | 0.986 | 0.990 | 0.992 | 0.991 | 0.987 | 0.983 | 0.973 | 0.977 | 0.977
8 | 0.982 | 0.988 | 0.992 | 0.992 | 0.990 | 0.985 | 0.985 | 0.971 | 0.975 | 0.968
9 | 0.983 | 0.990 | 0.993 | 0.992 | 0.988 | 0.984 | 0.982 | 0.968 | 0.968 | 0.972
10 | 0.985 | 0.992 | 0.993 | 0.990 | 0.987 | 0.984 | 0.980 | 0.962 | 0.971 | 0.968
11 | 0.986 | 0.994 | 0.993 | 0.989 | 0.987 | 0.982 | 0.976 | 0.956 | 0.970 | 0.962
12 | 0.987 | 0.995 | 0.992 | 0.988 | 0.987 | 0.980 | 0.976 | 0.960 | 0.963 | 0.960
13 | 0.989 | 0.995 | 0.991 | 0.987 | 0.986 | 0.979 | 0.974 | 0.950 | 0.965 | 0.963
14 | 0.990 | 0.995 | 0.990 | 0.987 | 0.985 | 0.976 | 0.972 | 0.949 | 0.962 | 0.955

|  # $\delta t$ | 0.1 | 0.15 | 0.2 | 0.25 | 0.3 | 0.35 | 0.4 | 0.45 | 0.5 | 0.55
---|---|---|---|---|---|---|---|---|---|---|---
SeqHT | 1 | 0.973 | 0.974 | 0.975 | 0.976 | 0.978 | 0.980 | 0.981 | 0.983 | 0.984 | 0.985
2 | 0.974 | 0.976 | 0.978 | 0.980 | 0.983 | 0.984 | 0.985 | 0.985 | 0.986 | 0.987
3 | 0.975 | 0.978 | 0.980 | 0.983 | 0.986 | 0.986 | 0.988 | 0.987 | 0.985 | 0.982
4 | 0.977 | 0.980 | 0.982 | 0.985 | 0.989 | 0.989 | 0.990 | 0.981 | 0.978 | 0.977
5 | 0.978 | 0.982 | 0.984 | 0.988 | 0.991 | 0.986 | 0.982 | 0.976 | 0.973 | 0.973
6 | 0.979 | 0.983 | 0.986 | 0.991 | 0.991 | 0.983 | 0.982 | 0.972 | 0.968 | 0.971
7 | 0.980 | 0.985 | 0.989 | 0.991 | 0.990 | 0.979 | 0.977 | 0.967 | 0.973 | 0.968
8 | 0.982 | 0.987 | 0.990 | 0.991 | 0.989 | 0.976 | 0.979 | 0.963 | 0.965 | 0.957
9 | 0.983 | 0.989 | 0.991 | 0.990 | 0.987 | 0.973 | 0.974 | 0.960 | 0.957 | 0.960
10 | 0.984 | 0.991 | 0.991 | 0.989 | 0.986 | 0.972 | 0.971 | 0.953 | 0.962 | 0.953
11 | 0.986 | 0.992 | 0.990 | 0.987 | 0.986 | 0.969 | 0.969 | 0.944 | 0.961 | 0.948
12 | 0.987 | 0.993 | 0.989 | 0.987 | 0.985 | 0.967 | 0.965 | 0.951 | 0.949 | 0.942
13 | 0.988 | 0.993 | 0.988 | 0.986 | 0.985 | 0.964 | 0.964 | 0.935 | 0.953 | 0.944
14 | 0.990 | 0.993 | 0.987 | 0.985 | 0.984 | 0.961 | 0.961 | 0.938 | 0.949 | 0.935

Table 12: The fidelity of the adiabatically prepared $\lambda\phi^{4}$ ground
state as a function of $\delta t$ and the number of adiabatic steps, for
$n_{q}=12$ with $\phi_{\rm max}=4$ and $\lambda=10$. Each adiabatic step
consists of one second-order Trotter step. Results are shown in Fig. 17.