# Exact solutions to the quantum many-body problem using the geminal density
matrix

Nicholas Cox  ICFO–Institut de Ciencies Fotoniques, The Barcelona Institute of
Science and Technology, 08860 Castelldefels, Barcelona, Spain  The College of
Optics and Photonics (CREOL), University of Central Florida, Orlando, Florida
32816, USA<EMAIL_ADDRESS>
###### Abstract

It is virtually impossible to directly solve the Schrödinger equation for a
many-electron wave function due to the exponential growth in degrees of
freedom with increasing particle number. The two-body reduced density matrix
(2-RDM) formalism reduces this coordinate dependence to that of four particles
irrespective of the wave function’s dimensionality, providing a promising path
to solve the many-body problem. Unfortunately, errors arise in this approach
because the 2-RDM cannot practically be constrained to guarantee that it
corresponds to a valid wave function. Here we approach this so-called
$N$-representability problem by expanding the 2-RDM in a complete basis of
two-electron wave functions and studying the matrix formed by the expansion
coefficients. This quantity, which we call the geminal density matrix (GDM),
is found to evolve in time by a unitary transformation that preserves
$N$-representability. This evolution law enables us to calculate eigenstates
of strongly correlated systems by a fictitious adiabatic evolution in which
the electron-electron interaction is slowly switched on. We show how this
technique is used to diagonalize atomic Hamiltonians, finding that the problem
reduces to the solution of $\sim N(N-1)/2$ two-electron eigenstates of the
Helium atom on a grid of electron-electron interaction scaling factors.

##  I Introduction

In 1955, Löwdin [1, 2, 3] and Mayer [4] presented similar methods to express
the ground state energy of a many-electron quantum system as a functional of
the two-body reduced density matrix ($2$-RDM). Their work inspired belief in
the feasibility of solving complex many-body problems with an effective two-
particle analysis [5]. However, early calculations significantly
underestimated experimental ground state energies because the $2$-RDM was not
adequately constrained to ensure that it represents a valid many-body wave
function. Full determination of these constraints, known as the
$N$-representability conditions [6], would indeed yield a method to reduce the
many-body problem to an effective two-particle system. The development of new
constraints continues to improve the accuracy of 2-RDM calculations, but the
general problem remains unsolved.

As detailed in Ref. [7], modern $2$-RDM analysis primarily employs one of two
methods. The first begins by deriving a contracted Schrödinger equation (CSE)
that computes the energy as a function of the two and four-body reduced
density matrices ($4$-RDM) [8, 9, 10, 11]. Approximating the $4$-RDM as a
function of the $2$-RDM allows one to solve the CSE for a set of candidate
eigenstates from which the physically valid states are selected by imposing
$N$-representability conditions. The second technique, density matrix
variational theory (DMVT), aims to directly minimize the energy as a
functional of the $2$-RDM [12, 13, 14, 15, 16]. The most successful
applications of DMVT use a convex optimization scheme called semi-definite
programming [17, 15] in which the $N$-representability conditions are included
by a set of positivity conditions that restrict the search space for
solutions.

Here we take a conceptually simpler approach that begins by expanding the
$2$-RDM in a basis of two-electron geminal [18] eigenstates. The resulting
expansion coefficients are collected into a quantity we call the geminal
density matrix (GDM) that can be used to compute many-body observables from
effective two-body operators. The technique was first introduced by Bopp in an
attempt to calculate the ground state energy of selected ions [19]. Although
his method was exact, his results differed quite significantly from
experimental ground state energies. These errors were later attributed to the
non-$N$-representability of the assumed ground state matrix [6]. Until now,
very little work has been done to advance this matrix-based approach.

We will detail how the GDM formalism enables us to calculate the stationary
states of a general many-body Hamiltonian with two-body (electron-electron)
interactions. To do so, we need to place Bopp’s work on a solid theoretical
foundation to make sense of $N$-representability in the context of GDMs. Most
importantly, we find that the GDM must evolve unitarily in time by the
Liouville-Von Neumann equation in order to produce the same expectation values
as the time-dependent $N$-electron wave function. Since the wave function is
only useful insofar as it generates observables, any matrix that reproduces
these quantities is clearly a faithful representation of the quantum state.

Because the equation of motion preserves $N$-representability, we find it
useful to examine a Hamiltonian that slowly switches on the electron-electron
interaction by some time-dependent scaling. We show by the adiabatic theorem
that eigenstates of the non-interacting (initial) system evolve into those of
the interacting (final) system. In this way, we find it possible to construct
$N$-electron stationary states using $\sim N(N-1)/2$ eigenstates of an
effective two-electron Hamiltonian computed on a grid of interaction scaling
strengths.

As an example, we show that the effective Hamiltonian for an arbitrary atom or
ion reduces to a coordinate-scaled Helium Hamiltonian with some specific
electron-electron interaction strength. The result is that all atomic electron
energy eigenstates can be found strictly from the solution to this Helium atom
problem. While atoms provide the simplest use case, the formalism presented
here applies equally well to molecular and solid state systems.

The paper is organized as follows. Section II introduces the $2$-RDM then
expands it in a two-electron basis to define the GDM. Section III gives
examples of valid density matrices that serve as the starting point for
eigenstate calculations. Section IV derives the unitary time evolution law of
the GDM, which Section V applies to the cases in which the electron-electron
interaction is switched on suddenly or slowly. The latter leads to the
adiabatic theorem, which is exploited in section VI to solve the many-body
Schrödinger equation.

Appendix A derives the four necessary $N$-representability constraints listed
in Section II. Appendix B defines the matrix transformation imposed by a
change of geminal basis, and Appendix C provides an alternate derivation of
the GDM equation of motion found in Section IV.

Hartree atomic units with $\hbar=e=m_{0}=1/(4\pi\epsilon_{0})=1$ will be used
throughout.

##  II The geminal density matrix

This section presents the GDM as the basic quantity needed to calculate the
observables of a many-body quantum system. It begins with an introductory
derivation of the $2$-RDM and proceeds to define the GDM by expanding the
$2$-RDM in a two-electron basis. The section ends by presenting a short list
of necessary $N$-representability conditions that are derived in Appendix A.

###  II.1 The 2-RDM

The state of a given $N$-electron system is fully described by its wave
function $\Psi(\mathbf{x}_{1},\dots,\mathbf{x}_{N})$, where we combined
spatial and spin degrees of freedom into the symbol

$\mathbf{x}_{i}=\mathbf{r}_{i}\sigma_{i}.$ (1)

The wave function is the position representation
$\Psi(\mathbf{x}_{1},\dots,\mathbf{x}_{N})=\braket{\mathbf{x}_{1},\dots,\mathbf{x}_{N}}{\alpha}$
of the abstract vector $\ket{\alpha}$, and it must be anti-symmetric under
coordinate exchange $\mathbf{x}_{i}\leftrightarrow\mathbf{x}_{j}$ for any pair
$i$ and $j$. We impose the normalization condition

$\int|\Psi(\mathbf{x}_{1},\dots,\mathbf{x}_{N})|^{2}\prod_{i=1}^{N}d\mathbf{x}_{i}=1,$
(2)

with the integral over $d\mathbf{x}$ including a sum over spin coordinates by

$\int f(\mathbf{x})d\mathbf{x}=\sum_{\sigma}\int
f(\mathbf{r}\mathbf{\sigma})d\mathbf{r}.$ (3)

Observables are represented by $N$-body linear operators with position
representation
$A(\mathbf{x}_{1},\dots,\mathbf{x}_{N})=\braket{\mathbf{x}_{1},\dots,\mathbf{x}_{N}}{\hat{A}}{\mathbf{x}_{1},\dots,\mathbf{x}_{N}}$.
An operator must be symmetric under particle exchange so that it maps anti-
symmetrized wave functions to anti-symmetrized wave functions. In general, $A$
contains one-body components acting on the coordinates of single particles and
two-body components acting on pairs. Denoting one-body contributions by lower
case letters and two-body by upper case, $A$ is expressed in the position
representation as

$A(\mathbf{x}_{1},\dots,\mathbf{x}_{N})=\sum_{k}a_{1}(\mathbf{x}_{k})+\sum_{i,j>i}A_{2}(\mathbf{x}_{i},\mathbf{x}_{j}).$
(4)

We will find it convenient to combine Eq. 4 into the single sum over pairs

$A(\mathbf{x}_{1},\dots,\mathbf{x}_{N})=\sum_{i,j>i}A(\mathbf{x}_{i},\mathbf{x}_{j}),$
(5)

with the summand $A(\mathbf{x}_{i},\mathbf{x}_{j})$ defined by

$A(\mathbf{x}_{i},\mathbf{x}_{j})=A_{1}(\mathbf{x}_{i},\mathbf{x}_{j})+A_{2}(\mathbf{x}_{i},\mathbf{x}_{j}).$
(6)

The term $A_{1}(\mathbf{x}_{i},\mathbf{x}_{j})$ is the promotion of a one-body
operator to act on pairs

$A_{1}(\mathbf{x}_{i},\mathbf{x}_{j})=\frac{a_{1}(\mathbf{x}_{i})+a_{1}(\mathbf{x}_{j})}{N-1},$
(7)

where the denominator divides out the overcounting that occurs in the pair
sum.

Using Eq. 5, we calculate the expectation value of $A$ by the inner product

$\displaystyle\braket{A}$
$\displaystyle=\sum_{i,j>i}\int\Psi^{*}(\mathbf{x}_{1},\dots,\mathbf{x}_{N})$
$\displaystyle\times
A(\mathbf{x}_{i},\mathbf{x}_{j})\Psi(\mathbf{x}_{1},\dots,\mathbf{x}_{N})\prod_{i=1}^{N}d\mathbf{x}_{i}.$
(8)

Each integral in the pair sum is found to be equivalent by first swapping the
names $\mathbf{x}_{i}\leftrightarrow\mathbf{x}_{1}$ and
$\mathbf{x}_{j}\leftrightarrow\mathbf{x}_{2}$ then permuting the same pairs
within the argument lists of $\Psi$ and $\Psi^{*}$. Since each wave function
changes sign under coordinate permutation, the argument swapping leaves the
integrand unchanged.

Combining like variables by the shorthand
$\mathbf{X}=\mathbf{x}_{1},\mathbf{x}_{2}$ and
$\mathbf{Y}=\mathbf{x}_{3},\dots,\mathbf{x}_{N}$, we find that

$\braket{A}=\begin{pmatrix}N\\\
2\end{pmatrix}\int\Psi^{*}(\mathbf{X},\mathbf{Y})A(\mathbf{X})\Psi(\mathbf{X},\mathbf{Y})d\mathbf{X}d\mathbf{Y},$
(9)

where $d\mathbf{X}=d\mathbf{x}_{1}d\mathbf{x}_{2}$ and
$d\mathbf{Y}=d\mathbf{x}_{3}\cdots d\mathbf{x}_{N}$. The prefactor, equal to
$N(N-1)/2$, is the number of equivalent integrals in the pair sum. The
operator $A(\mathbf{X})$ is exactly the contribution of the single pair in Eq.
6,

$\displaystyle
A(\mathbf{X})=\frac{a_{1}(\mathbf{x}_{1})+a_{2}(\mathbf{x}_{2})}{N-1}+A_{2}(\mathbf{X}).$
(10)

Noting that $A(\mathbf{X})$ does not depend on $\mathbf{Y}$, we would like to
remove it from the integral over $d\mathbf{Y}$. We cannot factor out
$A(\mathbf{X})$ directly because it may contain derivatives that act on the
integral by the chain rule. We circumvent this problem by introducing a set of
primed coordinates $\mathbf{X}^{\prime}$ upon which the operator
$A(\mathbf{X}^{\prime})$ is taken to act. This definition permits the re-
formulation of Eq. 9 as

$\displaystyle\braket{A}=$ $\displaystyle\begin{pmatrix}N\\\
2\end{pmatrix}\int
d\mathbf{X}d\mathbf{X}^{\prime}\delta(\mathbf{X}^{\prime}-\mathbf{X})$
$\displaystyle\times
A(\mathbf{X}^{\prime})\int\Psi^{*}(\mathbf{X},\mathbf{Y})\Psi(\mathbf{X}^{\prime},\mathbf{Y})d\mathbf{Y}.$
(11)

In Eq. II.1 we defined the delta function

$\delta(\mathbf{X}^{\prime}-\mathbf{X})=\prod_{i=1}^{N}\delta(\mathbf{r}_{i}^{\prime}-\mathbf{r}_{i})\delta_{\sigma_{i}^{\prime},\sigma_{i}}$
(12)

that includes Dirac delta functions for position and Kronecker delta functions
for spin polarization.

Finally, we define the $2$-RDM

$\rho(\mathbf{X},\mathbf{X}^{\prime})=\begin{pmatrix}N\\\
2\end{pmatrix}\int\Psi^{*}(\mathbf{X},\mathbf{Y})\Psi(\mathbf{X}^{\prime},\mathbf{Y})d\mathbf{Y},$
(13)

so that Eq. II.1 simplifies to

$\braket{A}=\int
d\mathbf{X}d\mathbf{X}^{\prime}\delta(\mathbf{X}^{\prime}-\mathbf{X})A(\mathbf{X}^{\prime})\rho(\mathbf{X},\mathbf{X}^{\prime}).$
(14)

Eqs. 13 and 14 provide an enormous complexity reduction compared to Eq. II.1,
as they depend only on $\mathbf{X}$ and $\mathbf{X}^{\prime}$ regardless of
the number of particles under investigation.

###  II.2 From the 2-RDM to the GDM

We will now represent the $2$-RDM in matrix form, starting by expressing it as
the position representation of a non-local two-body linear operator $\hat{D}$:

$\rho(\mathbf{X},\mathbf{X}^{\prime})=\braket{\mathbf{X}^{\prime}}{\hat{D}}{\mathbf{X}}.$
(15)

Introducing a complete set of geminal basis states $\ket{i}$ with wave
functions $\psi_{i}(\mathbf{X})=\braket{\mathbf{X}}{i}$, we expand $\hat{D}$
by inserting two resolutions of the identity $1=\sum_{i}\ket{i}\bra{i}$ to
find

$\hat{D}=\sum_{mn}\ket{m}\braket{m}{\hat{D}}{n}\bra{n}.$ (16)

Defining $D_{mn}=\braket{m}{\hat{D}}{n}$ and taking the position
representation by pre-multiplying $\bra{\mathbf{X}^{\prime}}$ and post-
multiplying $\ket{\mathbf{X}}$, the $2$-RDM as expressed in Eq. 15 takes the
form

$\rho(\mathbf{X},\mathbf{X}^{\prime})=\sum_{mn}D_{mn}\psi_{n}^{*}(\mathbf{X})\psi_{m}(\mathbf{X}^{\prime}).$
(17)

We can prove the validity of this expansion by construction; applying
$\int\psi_{n^{\prime}}(\mathbf{X})\psi_{m^{\prime}}^{*}(\mathbf{X}^{\prime})d\mathbf{X}d\mathbf{X}^{\prime}$
to either side of Eq. 17 isolates element $D_{m^{\prime}n^{\prime}}$ as a
functional of the $2$-RDM. See Eq. 129 in Appendix A for more detail.

Inserting Eq. 17 into Eq. 14 yields

$\braket{A}=\sum_{mn}D_{mn}A_{nm},$ (18)

where

$A_{nm}=\int\psi_{n}^{*}(\mathbf{X})A(\mathbf{X})\psi_{m}(\mathbf{X})d\mathbf{X}.$
(19)

We now define the GDM as the matrix $\mathbf{D}$ of coefficients $D_{mn}$.
Similarly defining $\mathbf{A}$ to have coefficients $A_{mn}$, Eq. 18 reduces
to the matrix form

$\braket{A}=\operatorname{Tr}[\mathbf{D}\mathbf{A}].$ (20)

Eq. 20 was derived without approximation so it exactly reproduces any
observable quantity given by the many-body wave function. Thus, knowledge of
the matrix $\mathbf{D}$ is equivalent to knowledge of
$\Psi(\mathbf{x}_{1},\dots,\mathbf{x}_{N})$.

Although we can always generate a unique GDM from the many-body wave function,
the process can not be reversed to construct the wave function from the GDM.
As a result, it is necessary to constrain $\mathbf{D}$ so that it is
guaranteed to satisfy $N$-representability. In Appendix A we derive the
following four necessary conditions:

$\displaystyle\mathbf{D}$ $\displaystyle=\mathbf{D}^{\dagger}$ (21a)
$\displaystyle 0\leq$ $\displaystyle D_{nn}\leq 1$ (21b)
$\displaystyle\operatorname{Tr}[\mathbf{D}]$
$\displaystyle=\begin{pmatrix}N\\\ 2\end{pmatrix}$ (21c) $\displaystyle
0\leq\operatorname{Tr}[\mathbf{D}^{2}]$ $\displaystyle\leq\begin{pmatrix}N\\\
2\end{pmatrix}.$ (21d)

Condition (21a) says the matrix must be Hermitian. We can interpret rule (21b)
as the maximum occupation number for a given geminal and (21c) as the total
number of electron pairs present in the wave function. The final expression,
(21d), is derivable from the first three and provides a way to distinguish
states in a manner that is invariant under unitary matrix transformation. In
fact, we will find that all matrices of interest in this work satisfy the
strict equality $\operatorname{Tr}[\mathbf{D}^{2}]=N(N-1)/2$.

In Section III we find that Eq. 21 is an insufficient set of constraints
because it is possible to define a non-$N$-representable GDM that obeys these
rules.

##  III Matrix Examples

This section begins with examples of simple GDMs followed by a demonstration
of the $N$-representability problem. Subsection III.2 explores the application
of these matrices to solve for the stationary states of interacting many-body
Hamiltonians. It is found that solving such systems may be possible through a
time-dependent analysis.

###  III.1 Matrices with a non-interacting geminal basis

Using an orthonormal basis of one-electron wave functions
$\phi_{i}(\mathbf{x})$, we can build $N$-electron Slater determinants by

$\Psi_{\\{\alpha\\}}(\mathbf{x}_{1},\dots,\mathbf{x}_{N})=\hat{S}_{-}\prod_{i=1}^{N}\phi_{\alpha_{i}}(\mathbf{x}_{i}).$
(22)

In Eq. 22, we defined a configuration $\\{\alpha\\}$ to be an ordered
collection of integers $\alpha_{i}$ that specify the single-particle
eigenstates included in a given product. The operator $\hat{S}_{-}$ transforms
the product into an anti-symmetrized wave function by the determinant operator
in Eqs. 134 and 135.

Per Eq. 17, the GDM is defined with respect to a complete basis of two-
electron eigenstates. For now we take this basis to be the two-particle Slater
determinants

$\psi_{\mathbf{n}}(\mathbf{X})=\frac{1}{\sqrt{2}}\left(\phi_{n_{1}}(\mathbf{x}_{1})\phi_{n_{2}}(\mathbf{x}_{2})-\phi_{n_{1}}(\mathbf{x}_{2})\phi_{n_{2}}(\mathbf{x}_{1})\right)$
(23)

labeled by the pair of integers $\mathbf{n}=\\{n_{1},n_{2}\\}$. For a many-
electron state defined by a single configuration $\\{\alpha\\}$, we find in
Appendix A that the $2$-RDM can be represented as a rank four tensor indexed
by the pairs

$D_{\mathbf{mn}}=\begin{cases}1,&\text{if }\mathbf{m}=\mathbf{n}\text{ and
}n_{1},n_{2}\in\\{\alpha\\}\\\ 0,&\text{otherwise}.\end{cases}$ (24)

To transform Eq. 24 into a matrix with indices $D_{mn}$, we must assign each
pair $\mathbf{n}$ to a single integer. For this purpose we choose the mapping
shown in Table 1.

Table 1: Geminal index map $\mathbf{n}$ | $(1,2)$ | $(1,3)$ | $(2,3)$ | $(1,4)$ | $(2,4)$ | $(3,4)$ | …
---|---|---|---|---|---|---|---
$n$ | $1$ | $2$ | $3$ | $4$ | $5$ | $6$ | $\dots$

We now have all the necessary tools to give an example; suppose our wave
function is the three-electron Slater determinant
$\Psi(\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{x}_{3})=\hat{S}_{-}\phi_{1}(\mathbf{x}_{1})\phi_{2}(\mathbf{x}_{2})\phi_{3}(\mathbf{x}_{3})$.
From Eq. 24 it is clear that $D_{\mathbf{nn}}=1$ for
$\mathbf{n}=\\{1,2\\},\\{1,3\\}$ and $\\{2,3\\}$. Converting to a matrix by
Table 1, we find that $D_{nn}=1$ for $n=1,2,3$. This example generalizes
trivially to $N$-electron Slater determinants, which are represented by
matrices with $N(N-1)/2$ ones placed along the diagonal. Such matrices obey
the strict equality of Eq. 21d

$\operatorname{Tr}[\mathbf{D}^{2}]=\begin{pmatrix}N\\\ 2\end{pmatrix}$ (25)

and the equivalent idempotence condition

$\mathbf{D}^{2}=\mathbf{D}.$ (26)

It is tempting to conclude that we can arbitrarily place ones on the diagonal,
but this assumption fails immediately for the matrix $\mathbf{D}^{\prime}$
with $D^{\prime}_{nn}=1$ for $n=1,2,4$. According to Table 1, the Slater
determinant which generates this matrix must contains pairs
$\\{1,2\\},\\{1,3\\},$ and $\\{1,4\\}$. It is impossible for a single three-
electron product to contain all four basis functions
$\phi_{1}(\mathbf{x}),\phi_{2}(\mathbf{x}),\phi_{3}(\mathbf{x})$ and
$\phi_{4}(\mathbf{x})$, so the matrix is not $N$-representable. We conclude
that the restrictions in Eq. 21 are insufficient to ensure
$N$-representability because $\mathbf{D}^{\prime}$ satisfies all the rules and
fails to correspond to a valid wave function.

We may also create valid states that do not satisfy Eqs. 25 and 26 by taking a
linear superposition of configurations

$\Psi(\mathbf{x}_{1},\dots,\mathbf{x}_{N})=\sum_{\\{\alpha\\}}C_{\\{\alpha\\}}\Psi_{\\{\alpha\\}}(\mathbf{x}_{1},\dots,\mathbf{x}_{N}),$
(27)

whose diagonal elements are (Eq. 142)

$\displaystyle
D_{\mathbf{n}\mathbf{n}}=\sum_{\\{\alpha\\}\ni\mathbf{n}}\left|C_{\\{\alpha\\}}\right|^{2}.$
(28)

The remaining elements $D_{\mathbf{mn}}$ with $\mathbf{m}\neq\mathbf{n}$ are
only non-zero when two configurations share all but two basis functions (Eq.
141). Take for example the constant superposition of $M$ disjoint anti-
symmetrized configurations

$\displaystyle\Psi$
$\displaystyle(\mathbf{x}_{1},\dots,\mathbf{x}_{N})=\frac{1}{\sqrt{M}}\hat{S}_{-}\Big{[}\Big{(}\phi_{1}(\mathbf{x}_{1})\cdots\phi_{N}(\mathbf{x}_{N})\Big{)}$
$\displaystyle+\Big{(}\phi_{N+1}(\mathbf{x}_{1})\cdots\phi_{2N}(\mathbf{x}_{N})\Big{)}$
$\displaystyle+\dots+\Big{(}\phi_{(M-1)N+1}(\mathbf{x}_{1})\cdots\phi_{MN}(\mathbf{x}_{N})\Big{)}\Big{]},$
(29)

where each $\phi_{i}(\mathbf{x})\neq\phi_{j}(\mathbf{x})$. Clearly, the off-
diagonal elements of Eq. III.1 are all zero. We then find from Eq. 28 that
$D_{\mathbf{nn}}=1/M$ because each pair in the expansion appears in a single
configuration with coefficient $1/\sqrt{M}$. It is straightforward to show
that $\operatorname{Tr}[\mathbf{D}^{2}]=N(N-1)/(2M^{2})$, which decreases as
the number $M$ of disjoint configurations increases.

###  III.2 Matrices representing a real system

The previous examples were defined without regard to their relationship to a
physical system. We now aim to make use of the GDM to solve the eigenvalue
relation

$H(\mathbf{x}_{1},\dots,\mathbf{x}_{N})\Psi_{i}(\mathbf{x}_{1},\dots,\mathbf{x}_{N})=\mathcal{E}_{i}\Psi_{i}(\mathbf{x}_{1},\dots,\mathbf{x}_{N})$
(30)

for an $N$-electron Hamiltonian $H(\mathbf{x}_{1},\dots,\mathbf{x}_{N})$
containing one and two-body terms as in Eq. 4. We discovered in Section II
that any wave function $\Psi_{i}(\mathbf{x}_{1},\dots,\mathbf{x}_{N})$
corresponds to a matrix $\mathbf{D}_{i}$ that produces the eigenvalues
$\mathcal{E}_{i}$ of Eq. 30 by

$\mathcal{E}_{i}=\operatorname{Tr}[\mathbf{D}_{i}\mathbf{H}].$ (31)

The matrix $\mathbf{H}$ has elements

$H_{mn}=\int\psi_{m}^{*}(\mathbf{X})H(\mathbf{X})\psi_{n}(\mathbf{X})d\mathbf{X},$
(32)

where $H(\mathbf{X})$ is the effective two-particle Hamiltonian (Eq. 10)
acting on the geminal basis functions $\psi_{i}(\mathbf{X})$. Instead of using
the two-electron Slater determinants of Eq. 23, it will be advantageous to
choose the basis that diagonalizes $H(\mathbf{X})$ by

$H(\mathbf{X})\psi_{j}(\mathbf{X})=E_{j}\psi_{j}(\mathbf{X}).$ (33)

Suppose that we have solved for at least the $N(N-1)/2$ lowest energy
eigenstates of Eq. 31 and computed the diagonal matrix $\mathbf{H}$ by Eq. 32.
We may then discover the ground state by finding the $N$-representable matrix
$\mathbf{D}_{i}$ that yields the minimum possible energy in Eq. 31. Dropping
subscript $i$ from the GDM, the minimum energy state that satisfies the rules
of Eq. 21 is

$\displaystyle D_{mn}$ $\displaystyle=\begin{cases}1&\text{if }m=n\text{ and
}n\leq\begin{pmatrix}N\\\ 2\end{pmatrix}\\\ 0&\text{otherwise}.\end{cases}$
(34)

Unfortunately, the matrix in Eq. 34 is not guaranteed to be $N$-representable
because it was discovered through a minimization subject to an insufficient
set of constraints. It is nonetheless worthwhile to introduce Eq. 34 for two
reasons, the first being that it is precisely the form of the ground state
matrices postulated by Bopp [19] for his atomic calculations. The second use
is to gain insight into the $N$-representability problem and the path toward
its resolution.

We may wonder if the correct ground state GDM may be non-diagonal in contrast
to Eq. 34. However, borrowing intuition from single-particle mixed-case
density matrices, we expect off-diagonal elements to introduce temporal
density oscillations that render the state non-stationary. Another alternative
is that the eigenstate matrices are diagonal but the elements may be any real
number between 0 and 1. In this case, the GDM formalism seems not to reduce
difficulty of the energy minimization because the optimization occurs over a
possibly infinite set of matrix coefficients.

There is a hint following from a physical argument that stationary states do
not exhibit such non-integral occupation levels. Suppose $N$ electrons begin
in a Slater determinant in a system governed by a Hamiltonian with electron-
electron interactions. The initial GDM is trivially $N$-representable as it is
constructed directly from the wave function, but it is clearly not stationary.
Given enough time, however, we expect the system to relax to the ground state
by processes like the emission of radiation. Using single particle density
matrices as a guide again, we may conjecture that the electron-radiation
interaction proceeds as a unitary transformation of the GDM.

Since $\operatorname{Tr}[\mathbf{D}^{2}]$ is conserved under unitary
transformation, the ground state must also satisfy
$\operatorname{Tr}[\mathbf{D}^{2}]=N(N-1)/2$ and $\mathbf{D}^{2}=\mathbf{D}$.
These conditions are inconsistent with the presence of non-integer occupation
numbers that decrease the trace of the squared GDM. Accordingly, we may
predict that all stationary states are represented by $N(N-1)/2$ ones along
the diagonal but the remaining difficulty is to determine which collection of
occupied states defines an $N$-representable GDM.

The preceding discussion suggests that Bopp’s model may be closer to correct
than previously thought. After all, his method calculated the ground state of
O5+ with a remarkably low error of $0.017\%$ (See tables I and II of Ref
[19]). Although the method fared worse for Be+ with an error of $0.86\%$, the
matrix designated as the first excited state curiously had energy within
$0.040\%$ of the experimental ground state. While it is always possible that
these anomalies are the result of coincidence, it deserves to be investigated
whether Bopp’s errors originated from a simple improper accounting of two-
electron eigenstates.

Clearly, it is not immediately obvious that the GDM should actually evolve by
unitary matrix transformation. For this reason, the next section carefully
derives the evolution law.

##  IV The time evolution equation

Subsection III.2 emphasized the importance of modeling the time evolution of
the GDM for the discovery of many-body eigenstates. In this section we derive
the governing equation, finding that the GDM indeed evolves in unitary fashion
by the Liouville-Von Neumann equation.

We begin by generalizing the equations of Section II to apply at arbitrary
times. Most importantly, our wave function will be the position representation
of the time-dependent abstract operator $\alpha(t)$

$\Psi(\mathbf{x}_{1},\dots,\mathbf{x}_{N}|t)=\braket{\mathbf{x}_{1},\dots,\mathbf{x}_{N}}{\alpha(t)}.$
(35)

The expectation value of a generally time-dependent linear operator
$\hat{A}(t)$ for the state in Eq. 35 is

$\braket{A}(t)=\braket{\alpha(t)}{\hat{A}(t)}{\alpha(t)}.$ (36)

We can also calculate observables by the equivalent $2$-RDM formulation,
defining the time-dependent $2$-RDM

$\rho(\mathbf{X},\mathbf{X}^{\prime}|t)=\int\Psi^{*}(\mathbf{X},\mathbf{Y}|t)\Psi(\mathbf{X}^{\prime},\mathbf{Y}|t)d\mathbf{Y}$
(37)

following Eq. 13. The expectation value is then given by

$\braket{A}(t)=\int\delta(\mathbf{X}-\mathbf{X}^{\prime})A(\mathbf{X}^{\prime}|t)\rho(\mathbf{X},\mathbf{X}^{\prime}|t)d\mathbf{X}d\mathbf{X}^{\prime}.$
(38)

The question that we need to answer is: How must
$\rho(\mathbf{X},\mathbf{X}^{\prime}|t)$ evolve in time so that the
expectation values computed by Eq. 38 match those given in Eq. 36 by the many-
body wave function? The ability to compute the same observable quantities for
a given time means that $\rho(\mathbf{X},\mathbf{X}^{\prime}|t)$ and its
corresponding $\mathbf{D}(t)$ furnish a complete representation of the quantum
state.

The first step in answering the posed question is to differentiate Eq. 38 to
find

$\displaystyle\frac{d}{dt}\braket{A}(t)$
$\displaystyle=\Braket{\frac{dA}{dt}}$ $\displaystyle+\int
d\mathbf{X}d\mathbf{X}^{\prime}\delta(\mathbf{X}-\mathbf{X}^{\prime})A(\mathbf{X}^{\prime})\dot{\rho}(\mathbf{X},\mathbf{X}^{\prime}|t).$
(39)

We must find a matching expression for Eq. 37 in the hope of finding an
equation that connects $\dot{\rho}(\mathbf{X},\mathbf{X}^{\prime}|t)$ to known
quantities. For brevity, we will group all $N$ coordinate sets together into
the symbol

$\overline{\mathbf{X}}=\mathbf{x}_{1},\dots,\mathbf{x}_{N}$ (40)

so that we can express the $N$-electron wave function at some initial time
$t_{0}$ as $\Psi(\overline{\mathbf{X}}|t_{0})$. This wave function evolves
according to the unitary time evolution operator $U(t,t_{0})$ by

$\Psi(\overline{\mathbf{X}}|t)=U(t,t_{0})\Psi(\overline{\mathbf{X}}|t_{0}).$
(41)

Expressing the time-dependent operator $A(\overline{\mathbf{X}}|t)$ as a pair
sum following Eq. 5, we find the position representation of the inner product
in Eq. 36

$\displaystyle\braket{A}(t)$
$\displaystyle=\sum_{i,j>i}\int\Psi^{*}(\overline{\mathbf{X}}|t_{0})U^{\dagger}(t,t_{0})$
$\displaystyle\times
A(\mathbf{x}_{i},\mathbf{x}_{j}|t)U(t,t_{0})\Psi(\overline{\mathbf{X}}|t_{0})d\overline{\mathbf{X}}.$
(42)

Recalling the notation of Section II, we have defined the differential
$d\overline{\mathbf{X}}=d\mathbf{x}_{1}\cdots d\mathbf{x}_{N}$ which includes
spin sums by Eq. 3.

We proceed to differentiate Eq. IV, distributing derivatives by the chain rule
to $U^{\dagger}(t,t_{0})$, $A(\mathbf{x}_{i},\mathbf{x}_{j})$ and
$U(t,t_{0})$. Applying the time-dependent Schrödinger equations
$i\partial_{t}U(t)=H(t)U(t)$ and
$-i\partial_{t}U^{\dagger}(t)=U^{\dagger}(t)H(t)$, we find that

$\displaystyle\frac{d}{dt}\braket{A}(t)=\Braket{\frac{dA}{dt}}-iK(t),$ (43)

where $K(t)$ is the expectation value of the many-body commutator

$\displaystyle K(t)$ $\displaystyle=\int
d\overline{\mathbf{X}}\Psi^{*}(\overline{\mathbf{X}}|t)$
$\displaystyle\times\sum_{\begin{subarray}{c}i,j>i\\\
k,l>k\end{subarray}}\left[H(\mathbf{x}_{i},\mathbf{x}_{j}|t),A(\mathbf{x}_{k},\mathbf{x}_{l}|t)\right]\Psi(\overline{\mathbf{X}}|t)d\overline{\mathbf{X}}.$
(44)

In Eq. IV we expressed the many-body Hamiltonian $H(\overline{\mathbf{X}}|t)$
as a pair sum per Eq. 5.

Simplify Eq. IV by splitting $H(\mathbf{x}_{i},\mathbf{x}_{j}|t)$ and
$A(\mathbf{x}_{k},\mathbf{x}_{l}|t)$ into their one and two body components
following Eq. 6. Defining $K_{\alpha\beta}(t)$ to be the integral involving
the commutator between $\alpha$-body terms of the Hamiltonian and $\beta$-body
terms of $A$, $K(t)$ decomposes into the sum

$K(t)=K_{11}(t)+K_{12}(t)+K_{21}(t)+K_{22}(t).$ (45)

We will now calculate $K_{1\beta}(t)$, which deals with the commutator between
one-body Hamiltonian terms and both one and two-body portions of $A$. By
swapping each $\mathbf{x}_{k}$ to $\mathbf{x}_{1}$ and $\mathbf{x}_{l}$ to
$\mathbf{x}_{2}$, the sum over $k$ and $l>k$ in IV reduces to $N(N-1)/2$
identical integrals. Continuing to expand
$H_{1}(\mathbf{x}_{i},\mathbf{x}_{j})$ by Eq. 7, we can separate
$\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ contributions to find

$\displaystyle K_{1\beta}(t)$ $\displaystyle=\begin{pmatrix}N\\\
2\end{pmatrix}\int d\overline{\mathbf{X}}\Psi^{*}(\overline{\mathbf{X}}|t)$
$\displaystyle\times\Bigg{(}$
$\displaystyle\sum_{i,j>i}\left[\frac{h_{1}(\mathbf{x}_{i}|t)}{N-1},A_{\beta}(\mathbf{x}_{1},\mathbf{x}_{2}|t)\right]$
$\displaystyle+$
$\displaystyle\sum_{i,j>i}\left[\frac{h_{1}(\mathbf{x}_{j}|t)}{N-1},A_{\beta}(\mathbf{x}_{1},\mathbf{x}_{2}|t)\right]\Bigg{)}\Psi(\overline{\mathbf{X}}|t).$
(46)

Most terms in Eq. IV cancel because operators acting on different coordinates
always commute. The only surviving contributions are those in which
$\mathbf{x}_{i}$ or $\mathbf{x}_{j}$ are equal to $\mathbf{x}_{1}$ or
$\mathbf{x}_{2}$. From the second line we pick up $N-1$ copies of the
commutator $[h_{1}(\mathbf{x}_{1}),A_{\beta}(\mathbf{x}_{1},\mathbf{x}_{2})]$
by fixing $i=1$ and running over all $N-1$ elements in the $j$ sum. We also
find from this expression a single instance of the commutator
$[h_{1}(\mathbf{x}_{2}),A_{\beta}(\mathbf{x}_{1},\mathbf{x}_{2})]$. In the
third line we find $N-2$ more copies of the commutator
$[h_{1}(\mathbf{x}_{2}),A_{\beta}(\mathbf{x}_{1},\mathbf{x}_{2})]$ by fixing
$i=2$ and summing over the $N-2$ values of $j$, so that each portion appears a
total of $N-1$ times.

Noting that the above argument applies identically to the computation of
$K_{\alpha 1}$, we find that

$\displaystyle K_{\alpha\beta}(t)$ $\displaystyle=\begin{pmatrix}N\\\
2\end{pmatrix}(N-1)\int\Psi^{*}(\overline{\mathbf{X}}|t)$
$\displaystyle\times\left[H_{\alpha}(\mathbf{X}|t),A_{\beta}(\mathbf{X}|t)\right]\Psi(\overline{\mathbf{X}}|t)d\overline{\mathbf{X}}$
(47)

for $\alpha\beta=11,12,21$. Eq. IV does not immediately apply to $K_{22}(t)$
because, unlike in Eq. IV, two-body operators cannot be separated into single-
coordinate expressions. As a result, 3-coordinate terms such as
$[H(\mathbf{x}_{1},\mathbf{x}_{2}),A(\mathbf{x}_{1},\mathbf{x}_{3})]$ remain
in the equation for $K_{22}(t)$. By restricting two-body terms of the
Hamiltonian and operator $A$ to depend only on position in the form
$\sum_{i,j>i}f(|\mathbf{r}_{i}-\mathbf{r}_{j}||t)$, two body terms will always
commute so that $[H_{2}(\mathbf{X}|t),A_{2}(\mathbf{X}|t)]=0$ and Eq. IV
applies trivially to $\alpha\beta=22$. The consequences of this requirement
are discussed in the Conclusion.

Absorb the prefactor $(N-1)$ of Eq. IV into the Hamiltonian by defining

$H^{\prime}(\mathbf{X}|t)=(N-1)H(\mathbf{X}|t).$ (48)

Summing the $K_{\alpha\beta}$ by Eq. 45 allows us to express $K(t)$ as

$\displaystyle K(t)$ $\displaystyle=\begin{pmatrix}N\\\ 2\end{pmatrix}\int
d\mathbf{X}d\mathbf{Y}\Big{\\{}$
$\displaystyle\Big{[}H^{\prime}(\mathbf{X}|t)\Psi^{*}(\mathbf{X},\mathbf{Y}|t)\Big{]}A(\mathbf{X}|t)\Psi(\mathbf{X},\mathbf{Y}|t)$
$\displaystyle-\Psi^{*}(\mathbf{X},\mathbf{Y}|t)H^{\prime}(\mathbf{X}|t)A(\mathbf{X}|t)\Psi(\mathbf{X},\mathbf{Y}|t)\Big{\\}}.$
(49)

In Eq. IV, we separated the terms of the commutator and chose one Hermitian
operator $H^{\prime}(\mathbf{X}|t)$ to act on the left copy of the wave
function. We can finally substitute Eq. 37 into Eq. IV to find the derivative
of the expectation value by Eq. 43:

$\displaystyle\frac{d}{dt}\braket{A}(t)$
$\displaystyle=\Braket{\frac{dA}{dt}}-i\int
d\mathbf{X}d\mathbf{X}^{\prime}\delta(\mathbf{X}-\mathbf{X}^{\prime})$ (50)
$\displaystyle\times
A(\mathbf{X}^{\prime})[H^{\prime}(\mathbf{X})-H^{\prime}(\mathbf{X}^{\prime})]\rho(\mathbf{X},\mathbf{X}^{\prime}|t).$
(51)

Note that one copy of the Hamiltonian depends on coordinates $\mathbf{X}$ as a
consequence of its acting to the left in the inner product. Comparing Eq. 50
to Eq. IV, we find the desired expression for
$\dot{\rho}(\mathbf{X},\mathbf{X}^{\prime}|t)$:

$\dot{\rho}(\mathbf{X},\mathbf{X}^{\prime}|t)=-i[H(\mathbf{X})-H(\mathbf{X}^{\prime})]\rho(\mathbf{X},\mathbf{X}^{\prime}|t).$
(52)

We proceed to derive a more convenient matrix representation of Eq. 52
starting with a time-dependent geminal expansion of the $2$-RDM

$\rho(\mathbf{X},\mathbf{X}^{\prime}|t)=\sum_{mn}D_{mn}(t)\psi_{n}^{*}(\mathbf{X}|t)\psi_{m}(\mathbf{X}^{\prime}|t).$
(53)

Note that Eq. 53 represents the most general case in which both the matrix
elements $D_{mn}(t)$ and the geminal basis functions may vary in time.
Plugging into Eq. 52 gives

$\displaystyle i\frac{d}{dt}$
$\displaystyle\braket{\mathbf{X}^{\prime}}{\hat{D}}{\mathbf{X}}=\sum_{mn}D_{mn}(t)\psi_{m}(\mathbf{X}^{\prime}|t)\left[H^{\prime}(\mathbf{X}|t)\psi_{n}^{*}(\mathbf{X}|t)\right]$
$\displaystyle-\sum_{mn}D_{mn}(t)\left[H^{\prime}(\mathbf{X}^{\prime}|t)\psi_{m}(\mathbf{X}^{\prime}|t)\right]\psi_{n}^{*}(\mathbf{X}|t),$
(54)

where we have rearranged the $\psi_{i}(\mathbf{X})$ into the most convenient
order. Substitute into Eq. IV the following identities

$\displaystyle\psi_{m}(\mathbf{X}^{\prime}|t)$
$\displaystyle=\braket{\mathbf{X}^{\prime}}{m(t)}$ (55)
$\displaystyle\psi_{n}^{*}(\mathbf{X}|t)$
$\displaystyle=\braket{n(t)}{\mathbf{X}}$ (56) $\displaystyle
H^{\prime}(\mathbf{X}^{\prime}|t)\psi_{m}(\mathbf{X}^{\prime}|t)$
$\displaystyle=\braket{\mathbf{X}^{\prime}}{\hat{H}^{\prime}(t)}{m(t)}$ (57)
$\displaystyle H^{\prime}(\mathbf{X}|t)\psi_{n}^{*}(\mathbf{X}|t)$
$\displaystyle=\braket{n(t)}{\hat{H}^{\prime}(t)}{\mathbf{X}},$ (58)

so that everything on the right-hand side sits between
$\bra{\mathbf{X}^{\prime}}$ and $\ket{\mathbf{X}}$. Recalling the definition
of $\hat{D}(t)$ (Eq. 16), we extract from Eq. IV the abstract operator
equation

$\frac{d}{dt}\hat{D}(t)=-i[\hat{H}^{\prime}(t),\hat{D}(t)].$ (59)

Eq. 59 is the familiar Liouville-Von Neumann equation.

We can specialize Eq. 59 to a matrix equation by choosing a time-independent
basis $\ket{i}$ with which to expand
$\hat{D}(t)=\sum_{mn}D_{mn}\ket{m}\bra{n}$ and
$\hat{H}(t)=H_{mn}\ket{m}\bra{n}$. The result is that

$\dot{\mathbf{D}}(t)=-i[\mathbf{H}^{\prime}(t),\mathbf{D}(t)],$ (60)

which we find by an alternate derivation in Appendix C. Eq. 60 can be shown to
evolve through matrix transformation by the operator $\mathbf{U}(t,t_{0})$ as

$\mathbf{D}(t)=\mathbf{U}(t,t_{0})\mathbf{D}(t_{0})\mathbf{U}^{\dagger}(t,t_{0}).$
(61)

Plugging Eq. 61 into Eq. 60 and matching expressions on either side we find
that $\mathbf{U}(t)$ must satisfy

$\dot{\mathbf{U}}(t,t_{0})=-i\mathbf{H}(t)\mathbf{U}(t,t_{0}),$ (62)

and it is a routine procedure to successively integrate Eq. 62 to infinite
order for

$\mathbf{U}(t,t_{0})=\mathcal{T}_{t}\exp\left(-i\int_{t_{0}}^{t}d\tau\mathbf{H}^{\prime}(\tau)\right).$
(63)

$\mathbf{U}(t)$, being the exponential of a skew-Hermitian matrix, is clearly
unitary. As mentioned at the end of Section II, we now see that the unitary
transformation-invariant property

$\mathbf{D}^{2}=\mathbf{D}$ (64)

holds for any state connected to a Slater determinant by a time-dependent
Hamiltonian. Since states can only change by such a unitary transformation,
Eq. 64 must be true for any GDM with fixed particle number.

##  V Switching on the electron-electron interaction

Having confirmed unitary evolution of the GDM, we will now exploit this
property to solve the many-body Schrödinger equation. We do so by studying an
$N$-electron Hamiltonian in which the Coulomb potential is switched on by a
temporal function $\lambda(t)$. Using Hartree atomic units with
$\hbar=m_{0}=|e|=1/(4\pi\epsilon_{0})=1$, the Hamiltonian is written

$H(\overline{\mathbf{X}}|t)=\sum_{i,j>i}H_{1}(\mathbf{x}_{i},\mathbf{x}_{j})+\frac{\lambda(t)}{|\mathbf{r}_{i}-\mathbf{r}_{j}|},$
(65)

with $\lambda(0)=0$ and $\lambda(t)=1$ for $t$ larger than some switching time
$T$.

As we saw in Eq. 10, Eq. 65 reduces in the $2$-RDM formalism to the effective
two-electron form

$H(\mathbf{X}|t)=H_{1}(\mathbf{X})+\frac{\lambda(t)}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}.$
(66)

We introduce two natural basis sets for the study of Eq. 66, the first being
the non-interacting basis that diagonalizes the Hamiltonian at $t=0$ by

$H(\mathbf{X}|0)\psi_{i}(\mathbf{X})=E_{i}\psi_{i}(\mathbf{X}).$ (67)

Similarly, we define the interacting basis by fixing $t=T$ and solving

$\displaystyle H(\mathbf{X}|T)\psi_{i}(\mathbf{X})$
$\displaystyle=\left(H_{1}(\mathbf{X})+\frac{1}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}\right)\psi_{i}(\mathbf{X})$
$\displaystyle=E_{i}\psi_{i}(\mathbf{X}).$ (68)

Using the subscript $I$ and $N$ for the interacting and non-interacting bases,
respectively, we can change between the two representations by the unitary
transformation (see Appendix B):

$\mathbf{D}_{I}(t)=\mathbf{U}_{I}^{N}\mathbf{D}_{N}(t)\big{(}\mathbf{U}_{I}^{N}\big{)}^{\dagger}.$
(69)

In the following, we observe the time evolution of a Slater determinant under
the influence of the Hamiltonian in Eqs. 65 and 66. Section V.1 treats the
case in which the Coulomb interaction is quickly switched on and Section V.2
details a slow adiabatic change which allows us to construct fully-interacting
many-body solutions.

###  V.1 The sudden approximation

Beginning with the electrons in a Slater determinant state, we instantaneously
turn on the Coulomb interaction by the step function $\lambda(t)=u(t-T)$. The
sudden approximation posits that the electron gas remains unchanged during
switching but starts to evolve according to the Hamiltonian
$H(\overline{\mathbf{X}}|t)$ for $t>T$. Since the Hamiltonian is constant for
$t>T$, the time-dependent GDM is found from Eqs. 61 and 63 to be

$\mathbf{D}(t)=e^{i\mathbf{H}^{\prime}(T)(t-T)}\mathbf{D}(T)e^{-i\mathbf{H}^{\prime}(T)(t-T)}.$
(70)

Selecting the interacting basis, we have that the Hamiltonian $\mathbf{H}$ is
diagonal while $\mathbf{D}(T)$ is non-diagonal by the transformation in Eq.
69. This choice of basis simplifies the matrix equation in Eq. 70 to

$D_{mn}(t)=e^{iE_{nm}^{\prime}(t-T)}D_{mn}(T),$ (71)

with $E^{\prime}_{mn}=(N-1)(E_{m}-E_{n})$ accounting for the multiplicative
constant attached to $H^{\prime}(\mathbf{X})$ in Eq. 48.

With the time-dependent GDM given by Eq. 71, we can calculate the observable
quantity described by the operator with position representation

$\rho(\mathbf{x};\mathbf{X})=\frac{\delta(\mathbf{x}_{1}-\mathbf{x})+\delta(\mathbf{x}_{2}-\mathbf{x})}{N-1}.$
(72)

Recalling that $\mathbf{x}=\mathbf{r}\sigma$, the expectation value
$\braket{\rho(\mathbf{x})}(t)$ gives the electron density at position
$\mathbf{r}$ with spin $\sigma$. For simplicity, we will define the symbol

$\rho(\mathbf{x},t)=\braket{\rho(\mathbf{x})}(t),$ (73)

which we calculate by the trace relation

$\rho(\mathbf{x},t)=\operatorname{Tr}[\mathbf{D}(t)\boldsymbol{\rho}(\mathbf{x})].$
(74)

The matrix elements of $\boldsymbol{\rho}(\mathbf{x})$ take the simplified
form

$\rho_{nm}(\mathbf{x})=\frac{2}{N-1}\int\psi_{n}^{*}(\mathbf{x},\mathbf{x}_{2})\psi_{m}(\mathbf{x},\mathbf{x}_{2})d\mathbf{x}_{2}.$
(75)

Computing the trace in Eq. 74 finally yields the electron density

$\rho(\mathbf{x},t)=\sum_{mn}D_{mn}(T)e^{-iE^{\prime}_{mn}(t-T)}\rho_{nm}(\mathbf{x}).$
(76)

We perform a simple check that particle number is conserved; integrating Eq.
76 over $\mathbf{x}$, we find from orthonormality and Eq. 21c that
$\int\rho(\mathbf{x},t)d\mathbf{x}=N$ as required. In the next subsection we
delve deeper into the implications of the electron density equation.

####  V.1.1 The origin of incoherent quantum fluctuations

We can extract a surprising amount of insight from the simple relation in Eq.
76. Separating the stationary and oscillating terms allows us to express the
density as

$\displaystyle\rho(\mathbf{x},t)$
$\displaystyle=\sum_{l}D_{ll}\rho_{ll}(\mathbf{x})+\sum_{\begin{subarray}{c}m,n\\\
E_{m}^{\prime}=E_{n}^{\prime}\end{subarray}}\text{Re}\\{P_{mn}(\mathbf{x})\\}$
$\displaystyle+\sum_{\begin{subarray}{c}p,q\\\ E_{p}^{\prime}\neq
E_{q}^{\prime}\end{subarray}}|P_{pq}(\mathbf{x})|\cos\left(E_{pq}^{\prime}t+\theta_{pq}(\mathbf{x})\right),$
(77)

with

$\displaystyle P_{ij}(\mathbf{x})$
$\displaystyle=D_{ij}(T)\rho_{ji}(\mathbf{x})$ (78)
$\displaystyle\theta_{ij}(\mathbf{x})$
$\displaystyle=\arg\left(P_{ij}(\mathbf{x})\right).$ (79)

The first sum in Eq. V.1.1 gives the contribution to the density from diagonal
matrix elements, and the second and third are from the degenerate and non-
degenerate off-diagonals, respectively.

Defining the time average functional for $t>T$ by

$\overline{f(\mathbf{x},t)}=\lim_{\tau\rightarrow\infty}\frac{1}{\tau}\int_{T}^{T+\tau}f(\mathbf{x},t^{\prime})dt^{\prime},$
(80)

we find the time-averaged density at $\mathbf{x}$ to be determined entirely by
the diagonal and degenerate terms

$\overline{\rho(\mathbf{x},t)}=\sum_{l}D_{ll}\rho_{ll}(\mathbf{x})+\sum_{\begin{subarray}{c}m,n\\\
E_{m}^{\prime}=E_{n}^{\prime}\end{subarray}}\text{Re}\\{P_{mn}(\mathbf{x})\\}.$
(81)

On the other hand, the non-degenerate off-diagonals introduce temporal density
fluctuations with variance

$\sigma^{2}(\mathbf{x})=\overline{\rho^{2}(\mathbf{x},t)}-\overline{\rho(\mathbf{x},t)}^{2}=\frac{1}{2}\sum_{\begin{subarray}{c}p,q\\\
E_{p}^{\prime}\neq E_{q}^{\prime}\end{subarray}}|P_{pq}(\mathbf{x})|^{2}.$
(82)

The fluctuations increase in magnitude as the number of off-diagonal terms
increases, an idea which can also be understood directly from Eq. V.1.1
whereby summing an increasing number of out-of-phase cosines leads to peaks in
the density that decrease in duration but increase in intensity. A similar
argument can be made for the spatial extent of these quasi-random density
spikes.

Eqs. V.1.1–82 are best understood in the context of a physical example.
Suppose we have a system of electrons experiencing the potential of a set of
nuclei at positions $\mathbf{R}_{i}(t)$ that are free to move in time. The
resulting effective electronic Hamiltonian is

$\displaystyle H^{\prime}(\mathbf{X}|t)$
$\displaystyle=\sum_{i=1}^{2}\left(-\frac{\nabla_{i}^{2}}{2}-\sum_{j}\frac{1}{|\mathbf{r}_{i}-\mathbf{R}_{j}(t)|}\right)$
$\displaystyle+(N-1)\frac{\lambda(t)}{|\mathbf{r}_{1}-\mathbf{r}_{2}|},$ (83)

noting that we used $H^{\prime}(\mathbf{X}|t)=(N-1)H(\mathbf{X}|t)$. Once
again the presence of $\lambda(t)=u(t-T)$ indicates that we abruptly switch on
the Coulomb potential at time $t=T$.

The electron-nuclei system evolves according to the coupled equations

$\displaystyle\dot{\mathbf{D}}(t)$
$\displaystyle=-i[\mathbf{H}^{\prime}(t),\mathbf{D}(t)]$ (84) $\displaystyle
M_{i}\ddot{\mathbf{R}}_{i}(t)$ $\displaystyle=\sum_{j\neq
i}\frac{1}{|\mathbf{R}_{i}(t)-\mathbf{R}_{j}(t)|^{2}}-\sum_{\sigma}\int\frac{\rho(\mathbf{r}\sigma,t)}{|\mathbf{r}-\mathbf{R}_{i}|^{2}}d\mathbf{r},$
(85)

where Eq. 85 is the classical non-relativistic equation of motion for nucleus
$i$ with mass $M_{i}$. Each nucleus feels a repulsive Coulomb force induced by
the other nuclei and an attractive Coulomb force from the electron gas. The
sum over electrons spins is explicitly written following Eq. 3. We take the
system to be at rest for $t<T$, meaning that all $\mathbf{R}_{i}(t)$ are fixed
in place and the electrons are in an eigenstate represented by a diagonal
matrix.

When $t>T$, the electron density fluctuates by a series of spikes localized in
space and time. Each peak causes an abrupt change in the electron-nucleus
Coulomb force in Eq. 85 which causes a near-instantaneous scattering of the
nuclei. The resulting motion of $\mathbf{R}_{i}(t)$ affects the Hamiltonian in
Eq. V.1.1 and subsequently the GDM by Eq. 84. The net effect is the transfer
of energy from the electrons to the nuclei. This process continues until an
equilibrium is reached wherein net energy ceases to flow between the
subsystems. There is an alternate picture that expands Eq. 85 as a sum of
normal modes of motion excited by the oscillatory terms of Eq. V.1.1. Summing
the out-of-phase normal mode oscillation amplitudes yields the same
thermalized motion as the scattering picture.

Instead of evolving by the non-physical Hamiltonian in Eq. V.1.1, we could
have fixed $\lambda=1$ and included the effect of a time dependent
electromagnetic potential $A(\mathbf{r},t)$. In this case we would see that
the field excites electron density fluctuations by introducing non-zero off-
diagonal elements to $\mathbf{D}$, which induces motion in the nuclei leading
again to thermalization. This example and the previous one illustrate a major
strength of the GDM formalism. By properly treating all electrons together, we
begin to see the emergence of classical behavior in a quantum system.

###  V.2 The degenerate adiabatic theorem

Section IV derived the equation of motion (Eq. 59)

$\frac{d}{dt}\hat{D}(t)=-i[\hat{H}^{\prime}(t),\hat{D}(t)]$ (86)

for the abstract density operator $\hat{D}(t)$. We found the corresponding
matrix equation (Eq. 60) after expanding $\hat{D}(t)$ in a time-independent
geminal basis.

Here we derive the adiabatic theorem, beginning by determining the matrix
equation of motion for a GDM expanded in a time-dependent basis in which the
Hamiltonian is always diagonal. That is, the basis functions satisfy

$\hat{H}(t)\ket{i(t)}=E_{i}(t)\ket{i(t)},$ (87)

where $E_{i}(t)$ is the instantaneous eigenvalue of state $\ket{i(t)}$ at time
$t$. Labeling matrix elements by $\mathcal{D}_{pq}(t)$, we have that

$\hat{D}(t)=\sum_{pq}\mathcal{D}_{pq}(t)\ket{p(t)}\bra{q(t)}.$ (88)

Inserting Eq. 88 into Eq. 86 and picking out the $mn$ component by pre-
multiplying $\bra{m(t)}$ and post-multiplying $\ket{n(t)}$, we find that the
right-hand side evaluates to $-iE_{mn}^{\prime}(t)D_{mn}(t)$ after applying
Eq. 87. Repeating for the left-hand side, by the chain rule we have that

$\displaystyle\Braket{m(t)}{\frac{d}{dt}\hat{D}(t)}{n(t)}=\dot{\mathcal{D}}_{mn}(t)$
$\displaystyle+\sum_{p}\mathcal{D}_{pn}(t)\braket{m(t)}{\dot{p}(t)}+\sum_{q}\mathcal{D}_{mq}(t)\braket{\dot{q}(t)}{n(t)}.$
(89)

Because $(d/dt)\braket{q(t)}{n(t)}=0$ by orthonormality at time $t$, we can
simplify Eq. V.2 by substituting
$\braket{\dot{q}(t)}{n(t)}\rightarrow-\braket{q(t)}{\dot{n}(t)}$ and taking
these terms to be the elements of the skew-Hermitian matrix $\mathbf{M}$ with
coefficients

$M_{ij}(t)=\Braket{i(t)}{\frac{d}{dt}}{j(t)}.$ (90)

After doing so, the final equation of motion for the $mn$ element is

$\displaystyle\dot{\mathcal{D}}_{mn}(t)=-iE_{mn}^{\prime}(t)\mathcal{D}_{mn}(t)$
$\displaystyle-\sum_{p}\mathcal{D}_{pn}(t)M_{mp}(t)+\sum_{q}\mathcal{D}_{mq}(t)M_{qn}(t).$
(91)

We will now simplify Eq. V.2 for a Hamiltonian that varies slowly over a long
time interval $T$. In the context of this work, $T$ will be the switching time
that appears in the Coulomb potential scaling $\lambda(t)$ of Eq. 66. Closely
following Ref. [20], we will discover how $\boldsymbol{\mathcal{D}}(t)$
evolves when $T\rightarrow\infty$. Defining a natural time

$s(t)=\frac{t}{T}$ (92)

scaled by the switching duration, we use the fact that $(d/dt)=T(d/ds)$ to see
that

$\displaystyle\dot{\mathcal{D}}_{mn}(s)$
$\displaystyle=-iTE_{mn}^{\prime}(s)\mathcal{D}_{pn}(s)$ $\displaystyle+$
$\displaystyle\sum_{q}\mathcal{D}_{mq}(s)M_{qn}(s)-\sum_{p}M_{mp}(s)\mathcal{D}_{pn}(s).$
(93)

In Eq. V.2 we used that $M_{ij}(s)=(1/T)M_{ij}(t)$ and chose a more convenient
ordering for the sum terms. The first term on the right hand side represents
the dynamical phase, which we factor out by defining $\tilde{\mathcal{D}}(t)$
as

$\mathcal{D}_{mn}(s)=\tilde{\mathcal{D}}_{mn}(s)e^{-iT\int_{0}^{s}E_{mn}^{\prime}(s^{\prime})ds^{\prime}},$
(94)

so that

$\displaystyle\frac{d}{ds}\tilde{\mathcal{D}}_{mn}(s)=\sum_{\begin{subarray}{c}q\\\
E_{q}=E_{n}\end{subarray}}\tilde{\mathcal{D}}_{mq}(s)M_{qn}(s)-\sum_{\begin{subarray}{c}p\\\
E_{p}=E_{m}\end{subarray}}M_{mp}(s)\tilde{\mathcal{D}}_{pn}(s)$
$\displaystyle+\sum_{\begin{subarray}{c}q\\\ E_{q}\neq
E_{n}\end{subarray}}\frac{d}{ds}\int_{0}^{s}ds^{\prime}\tilde{\mathcal{D}}_{mq}(s^{\prime})M_{qn}(s^{\prime})e^{-iT\int_{0}^{s^{\prime}}E_{nq}^{\prime}(s^{\prime\prime})ds^{\prime\prime}}$
$\displaystyle-\sum_{\begin{subarray}{c}p\\\ E_{p}\neq
E_{m}\end{subarray}}\frac{d}{ds}\int_{0}^{s}ds^{\prime}M_{mp}(s^{\prime})\tilde{\mathcal{D}}_{pn}(s^{\prime})e^{-iT\int_{0}^{s^{\prime}}E_{pm}^{\prime}(s^{\prime\prime})ds^{\prime\prime}}.$
(95)

In Eq. V.2 we separated the $p$ and $q$ sums into terms that accumulate
dynamical phase and those that do not. To the latter we have applied the
identity operator $\hat{1}=(d/ds)\int_{0}^{s}ds^{\prime}$.

We will now evaluate the last line of Eq. V.2 to understand the implications
of the dynamical phase. To begin, we simplify by defining

$F_{p}(s^{\prime})=M_{mp}(s^{\prime})\mathcal{D}_{pn}(s^{\prime})$ (96)

and

$g_{pm}(s^{\prime})=\int_{0}^{s^{\prime}}E_{pm}^{\prime}(s^{\prime\prime})ds^{\prime\prime}$
(97)

so that the integral over $ds^{\prime}$ (which we call $I_{p}(s)$) takes the
form

$I_{p}(s)=\int_{0}^{s}F_{p}(s^{\prime})e^{-iTg_{pm}(s^{\prime})}ds^{\prime}.$
(98)

Using that $\dot{g}_{pm}(s^{\prime})=E^{\prime}_{pm}(s^{\prime})$ by Eq. 97,
we multiply the integrand by the identity
$1=iT\dot{g}_{pm}(s^{\prime})/(iTE^{\prime}_{pm}(s^{\prime}))$ so that

$I_{p}(s)=\int_{0}^{s}\frac{F_{p}(s^{\prime})}{iTE_{pm}^{\prime}(s^{\prime})}\left(iT\dot{g}_{pm}(s^{\prime})e^{iTg_{pm}(s^{\prime})}\right)ds^{\prime}.$
(99)

Because
$iT\dot{g}_{pm}(s^{\prime})e^{iTg_{pm}(s^{\prime})}=(d/ds^{\prime})e^{iTg_{pm}(s^{\prime})}$,
we are able to integrate Eq. 99 parts to finally find

$\displaystyle I_{p}(s)$
$\displaystyle=\frac{1}{iT}\left[\frac{F_{p}(s^{\prime})}{E_{pm}^{\prime}(s^{\prime})}e^{iTg_{pm}(s^{\prime})}\right]_{0}^{s}$
$\displaystyle-\frac{1}{iT}\int_{0}^{s}\frac{d}{ds^{\prime}}\left(\frac{F_{p}(s^{\prime})}{E^{\prime}_{pm}(s^{\prime})}\right)e^{iTg_{pm}(s^{\prime})}ds^{\prime}.$
(100)

As long as $F_{p}(s^{\prime})$ is differentiable, we can take the adiabatic
approximation by letting $T\rightarrow\infty$ so that $I_{p}(s)\rightarrow 0$.
Since the integral in the second line proceeds identically, we conclude that
any contribution that accumulates dynamical phase will evaluate to zero.

Canceling these integral terms and reverting to the unscaled time $t$
simplifies Eq. V.2 to

$\frac{d}{dt}\tilde{\mathcal{D}}_{mn}(t)=\sum_{\begin{subarray}{c}q\\\
E_{q}=E_{n}\end{subarray}}\tilde{\mathcal{D}}_{mq}(t)M_{qn}(t)-\sum_{\begin{subarray}{c}p\\\
E_{p}=E_{m}\end{subarray}}M_{mp}(t)\tilde{\mathcal{D}}_{pn}(t),$ (101)

which we study in two separate cases. First suppose that states $m$ and $n$
are non-degenerate. Under this condition, the sum restriction $E_{q}=E_{n}$
implies that $E_{q}\neq E_{m}$ because the contrary would contradict the
assumption of non-degeneracy ($E_{m}\neq E_{n}$). The same applies for the $p$
sum so that only terms depending on $\mathcal{D}_{ij}(t)$ with $E_{i}(s)\neq
E_{j}$ appear on the right hand side. If we choose matrices with $D_{ij}(0)=0$
for all states $i,j$ with $E_{i}\neq E_{j}$, these terms remain $D_{ij}(t)=0$
at all times because these elements only intermix with each other.

What remains is to evaluate Eq. 101 when $m$ and $n$ are degenerate. In this
case, $p$ and $q$ both run over all states degenerate with $m$ and $n$ so that
the right-hand side reduces to a single sum. We see from Eq. 94 that the
degeneracy condition implies that
$\tilde{\boldsymbol{\mathcal{D}}}(t)=\boldsymbol{\mathcal{D}}(t)$. Taking the
basis to be ordered by increasing $E_{i}(t)$, we can define submatrices
$\boldsymbol{\mathcal{D}}_{\mu}(t)$ for each degenerate subspace $\mu$. Noting
that the terms in Eq. 101 are in the form of matrix multiplications, we find
that the equation of motion separates into the commutation relations

$\dot{\boldsymbol{\mathcal{D}}}_{\mu}(t)=-[\mathbf{M}_{\mu}(t),\boldsymbol{\mathcal{D}}_{\mu}(t)]$
(102)

for each subspace $\mu$. This expression details the accumulation of geometric
phase within a given degenerate subspace as time progresses. The final GDM is
simply the block diagonal direct sum

$\boldsymbol{\mathcal{D}}(t)=\text{diag}[\boldsymbol{\mathcal{D}}_{1}(t),\boldsymbol{\mathcal{D}}_{2}(t),\dots,\boldsymbol{\mathcal{D}}_{\mu}(t),\dots],$
(103)

valid for all times including $t>T$. As a final remark, note that Eq. 102
describes a unitary transformation because the matrix $\mathbf{M}(t)$ is skew-
Hermitian. We may also choose to define a Hermitian matrix
$\tilde{\mathbf{M}}(t)=-i\mathbf{M}(t)$ so that the expression more closely
resembles Eq. 60.

####  V.2.1 Practical considerations of the degenerate adiabatic theorem

The application of Eqs. 102 and 103 requires knowledge of the matrix
$\mathbf{M}(t)$. Using the relation $(d/dt)=\dot{\lambda}(d/d\lambda)$, we can
calculate the matrix elements from Eq. 90 in the position representation

$M_{ij}(t)=-i\dot{\lambda}(t)\int\psi_{i}^{*}(\mathbf{X}|\lambda)\frac{\partial}{\partial\lambda}\psi_{j}(\mathbf{X}|\lambda)d\mathbf{X}.$
(104)

However, our only prescription for the instantaneous eigenstates was that they
diagonalized the Hamiltonian at time $t$. We have not fixed the connection
between eigenstates $\psi_{i}(\mathbf{X}|\lambda)$ and
$\psi_{i}(\mathbf{X}|\lambda+\Delta\lambda)$. In fact, an eigenstate solver
may output a different arbitrary rotation for each value of $\lambda$ so that
$\psi_{j}(\mathbf{X}|\lambda)$ is discontinuous in $\lambda$ and the
$T\rightarrow\infty$ limit of Eqs. V.2 is no longer valid. Thus, we must
require that $\psi_{j}(\mathbf{X}|\lambda)$ is smooth as a function of
$\lambda$. In Section VI we will perturb the Hamiltonian with an asymmetric
potential to lift the degeneracy so that the connection between time-adjacent
eigenstates is obvious.

The final issue pertains to the $\lambda$-dependent spectrum of geminal
eigenstates. As eigenvalues approach each other with increasing $\lambda$,
they may cross or anti-cross depending on symmetry. Anti-crossings are handled
perfectly fine by Eqs. 102 and 103, but the derivation of the adiabatic
theorem breaks down at a crossing point. In a sense, the point acts as a pole
where the degenerate submatrices intermix by instantaneously picking up
geometric phase within the expanded degenerate subspace at the intersection.

We determine the behavior around these accidental degeneracies by a physical
argument. First note that the total energy as computed by the many-body wave
function must be differentiable, which we see upon differentiating the energy
expectation value by the Hellman-Feynman theorem:

$\frac{d\mathcal{E}}{dt}=\dot{\lambda}(t)\sum_{i,j>i}\int\frac{|\Psi(\overline{\mathbf{X}}|t)|^{2}}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}d\overline{\mathbf{X}}.$
(105)

It follows that $\mathcal{E}(t)$ is differentiable because the right hand side
of Eq. 105 is finite at all times. The GDM that represents the wave function
computes the energy as

$\mathcal{E}(t)=\sum_{\mu}\operatorname{Tr}[\boldsymbol{\mathcal{D}}_{\mu}(t)]E_{\mu}(\lambda(t)),$
(106)

where $\operatorname{Tr}[\boldsymbol{\mathcal{D}}_{\mu}(t)]$ is the population
in a given degenerate subspace with energy $E_{\mu}(\lambda(t))$. An abrupt
re-distribution of population at a crossing point $t_{0}$ would yield an
$\mathcal{E}(t)$ that is not differentiable at $t_{0}$. Since we know by Eq.
105 that $\mathcal{E}(t)$ must be differentiable, we conclude that no
instantaneous population interchange can occur at a level crossing.

##  VI Computing many-body eigenstates using the adiabatic theorem

We are now equipped to solve the many-electron Schrödinger equation

$\left(\sum_{ij}H(\mathbf{x}_{i},\mathbf{x}_{j})\right)\Psi_{n}(\overline{\mathbf{X}})=\mathcal{E}_{n}\Psi_{n}(\overline{\mathbf{X}})$
(107)

for the matrices $\mathbf{D}_{n}$ that represent wave functions
$\Psi_{n}(\overline{\mathbf{X}})$. The Hamiltonian of interest has terms

$\displaystyle
H(\mathbf{x}_{i},\mathbf{x}_{j})=H_{1}(\mathbf{x}_{i},\mathbf{x}_{j})+\frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|},$
(108)

containing the one-body contribution

$H_{1}(\mathbf{x}_{i},\mathbf{x}_{j})=\frac{1}{N-1}\left[-\frac{\nabla_{i}^{2}}{2}-\frac{\nabla_{j}^{2}}{2}+v(\mathbf{x}_{i})+v(\mathbf{x}_{j})\right]$
(109)

along with the Coulomb potential $1/|\mathbf{r}_{i}-\mathbf{r}_{j}|$. For
reasons that will be explained shortly, we treat Eqs. 107 and 108 as a special
case of the two-parameter Hamiltonian

$\displaystyle H(\epsilon,\lambda)$
$\displaystyle=\sum_{i,j>i}\Bigg{[}H_{1}(\mathbf{x}_{i},\mathbf{x}_{j})$
$\displaystyle+\epsilon\left(\frac{v_{p}(\mathbf{x}_{i})+v_{p}(\mathbf{x}_{j})}{N-1}\right)+\lambda\left(\frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}\right)\Bigg{]}$
(110)

that reduces to the original operator when $\epsilon=0$ and $\lambda=1$. The
one-body perturbing potential $v_{p}(\mathbf{x})$ is chosen to have symmetry
group containing only the identity so that the spectrum of the auxiliary many-
body system

$H(\overline{\mathbf{X}}|\epsilon,\lambda)\Psi_{n}(\overline{\mathbf{X}}|\epsilon,\lambda)=\mathcal{E}_{n}(\epsilon,\lambda)\Psi_{n}(\overline{\mathbf{X}}|\epsilon,\lambda)$
(111)

is nondegenerate when $\epsilon\neq 0$.

We represent a given eigenstate
$\Psi_{n}(\overline{\mathbf{X}}|\epsilon,\lambda)$ of Eq. 111 by its GDM
$\boldsymbol{\mathcal{D}}(\epsilon,\lambda)$, dropping the subscript $n$. The
electronic energy is then calculated by

$\mathcal{E}_{n}(\epsilon,\lambda)=\sum_{i}\mathcal{D}_{ii}(\epsilon,\lambda)E_{i}(\epsilon,\lambda),$
(112)

where each $E_{i}(\epsilon,\lambda)$ is an eigenvalue of the effective two-
particle Hamiltonian $H(\mathbf{X}|\epsilon,\lambda)$ satisfying the
Schrödinger equation

$H(\mathbf{X}|\epsilon,\lambda)\psi_{i}(\mathbf{X}|\epsilon,\lambda)=E_{i}(\epsilon,\lambda)\psi_{i}(\mathbf{X}|\epsilon,\lambda)$
(113)

for the geminal eigenstate $\psi_{i}(\mathbf{X}|\epsilon,\lambda)$.

As mentioned, the analysis takes parameters $\epsilon$ and $\lambda$ to be
functions of time. The dependence is chosen in such a way that the Hamiltonian
slowly changes from non-interacting to fully interacting over some time
interval. Beginning at $t=0$ with $\epsilon(0)=\lambda(0)=0$, we slowly ramp
the symmetry-breaking parameter $\epsilon(t)$ over time $T_{1}$ so that
$\epsilon(T_{1})=1$ and $\lambda(T_{1})=0$. The nondegenerate eigenstates at
$T_{1}$ remain Slater determinants due to the absence of electron-electron
interactions.

The crucial next step is to smoothly switch on the Coulomb interaction by
increasing $\lambda(t)$ from $0$ at $T_{1}$ to $1$ at $T_{2}$ with
$\epsilon=1$. Assuming adiabatic evolution in which
$T_{2}-T_{1}\rightarrow\infty$, we have by Eq. 102 that the GDM follows the
trivial relation $\dot{\boldsymbol{\mathcal{D}}}(t)=0$ because the two-
electron spectrum is nondegenerate. In simple terms, all elements of the GDM
stay fixed while the geminal eigenstates evolve to satisfy Eq. 113 for each
$\lambda(t)$ with $\epsilon=1$. Fig. 1 is an illustration of how these
parameter-dependent spectra may look. The curves do not represent the
diagonalization of a real Hamiltonian, but they provide insight into how the
technique works.

Figure 1: Example illustration of the eigenvalues of $H(\mathbf{X}|1,\lambda)$
versus $\lambda$. $\lambda=0$ and $\lambda=1$ eigenstate indices are marked by
numbers on the vertical axes. The curves do not represent the diagonalization
of a real Hamiltonian.

For the simplest non-trivial example, consider a 3-electron state constructed
from the geminal eigenstates of Fig. 1. Suppose that the GDM with geminals 1–3
occupied represents a valid Slater determinant at $\lambda(T_{1})=0$. As time
progresses, the matrix elements remain unchanged as each eigenstate evolves to
its $\lambda(T_{2})=1$ counterpart with $E_{i}(1,\lambda)$ varying smoothly
with $\lambda$. The smoothness condition is imposed to define the behavior at
level crossings; the total energy calculated by Eq. 112 should be
differentiable at all times by Eq. 105. As a consequence of the level crossing
between eigenvalues 3 and 4, the $\lambda=1$ (interacting) GDM with the 3
lowest energy eigenstates occupied is not $N$-representable if the
configuration $\\{1,2,4\\}$ is not a valid $\lambda=0$ Slater determinant.

Generalizing to $N$ electrons, the principal task is to compute at least
$N(N-1)/2$ geminal eigenstates on a grid of $\lambda$ points to generate
curves like those in Fig. 1. This number of states is a lower bound due to
constraints on the initial Slater determinant as well as the existence of
level crossings. With sufficiently fine $\lambda$ sampling, it is possible to
resolve all such crossings to connect the $\lambda=0$ and $\lambda=1$
geminals. This two-electron diagonalization can be performed with standard
techniques like the configuration interaction, noting that the adiabatic
theorem only serves to select which states are $N$-representable by virtue of
their being the result of evolution under a time-dependent Hamiltonian.

After determining valid symmetry-broken states by the above procedure, we must
finally ramp $\epsilon$ to zero at time $T_{3}$ to discover solutions to the
target Hamiltonian with $\epsilon=0$ and $\lambda=1$. During this process, the
degenerate subspaces recombine as $\epsilon$ decreases. We can calculate the
exact total energy at $T_{3}$ by

$\mathcal{E}=\sum_{\mu}\operatorname{Tr}[\boldsymbol{\mathcal{D}}_{\mu}(0)]E_{\mu}(0,1),$
(114)

where $\operatorname{Tr}[\boldsymbol{\mathcal{D}}_{\mu}(0)]$ is the time-
invariant population of subspace $\mu$. Even though the Slater determinants
were chosen at $T_{1}$, we may simply back-propagate to determine the
population at the initial time. Repeating for various $t=T_{1}$ Slater
determinants, the ground state of the interacting system is that which yields
the minimum energy by Eq. 114.

For observables other than the energy, the best we can do is calculate the
$\epsilon$-dependent expectation value by

$\braket{A}(\epsilon)=\operatorname{Tr}[\boldsymbol{\mathcal{D}}(T_{1})\mathbf{A}(\epsilon)].$
(115)

The matrix $\mathbf{A}$ can be determined to arbitrary precision by computing
it from the geminal eigenstates with sufficiently small $\epsilon$.
Conceptually, this limiting procedure of switching on and off the perturbation
$\epsilon$ can be considered to account for the geometric phase that would be
picked up in a degenerate adiabatic evolution.

###  VI.1 Calculating atomic energy eigenstates

This section details a simple but powerful example of the calculation of many-
body eigenstates using the GDM. We will find that the many-body Schrödinger
equation for an arbitrary atom or ion is analyzable strictly through the
solution of an appropriately-scaled Helium atom problem.

Begin with the Hamiltonian $H_{Z,N}$ for a central potential with nuclear
charge $Z$ and $N$ electrons. The fully-interacting system obeys the effective
two-particle Hamiltonian of Eq. 66 with $\lambda=1$:

$\displaystyle H_{Z,N}(\mathbf{r}_{1},\mathbf{r}_{2}|1)$
$\displaystyle=\frac{1}{N-1}\Bigg{[}-\frac{\nabla_{1}^{2}}{2}-\frac{\nabla_{2}^{2}}{2}$
$\displaystyle-Z\left(\frac{1}{|\mathbf{r}_{1}|}+\frac{1}{|\mathbf{r}_{2}|}\right)+\frac{N-1}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}\Bigg{]}.$
(116)

It will be convenient to reduce the bracketed expression into a form that
resembles the Helium atom Hamiltonian

$H_{2,2}(\mathbf{r}_{1},\mathbf{r}_{2}|\lambda)=-\frac{\nabla_{1}^{2}}{2}-\frac{\nabla_{1}^{2}}{2}-\frac{2}{|\mathbf{r}_{1}|}-\frac{2}{|\mathbf{r}_{2}|}+\frac{\lambda}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}.$
(117)

for some choice of $\lambda$. We find this simplification by transforming the
coordinates into a yet-undetermined natural scale

$\overline{\mathbf{r}}=a\mathbf{r}.$ (118)

This transformation changes Eq. VI.1 to

$\displaystyle H_{Z,N}(\overline{\mathbf{r}}_{1},\overline{\mathbf{r}}_{2}|1)$
$\displaystyle=\frac{a^{2}}{N-1}\Bigg{[}-\frac{\overline{\nabla}_{1}^{2}}{2}-\frac{\overline{\nabla}_{2}^{2}}{2}$
$\displaystyle-\frac{Z}{a}\left(\frac{1}{|\overline{\mathbf{r}}_{1}|}+\frac{1}{|\overline{\mathbf{r}}_{2}|}\right)+\frac{N-1}{a}\frac{1}{|\overline{\mathbf{r}}_{1}-\overline{\mathbf{r}}_{2}|}\Bigg{]}.$
(119)

Eq. VI.1 reduces to a scaled copy of the Helium Hamiltonian in Eq. 117 if we
choose

$a=\frac{Z}{2},$ (120)

from which it follows that $\lambda=2(N-1)/Z$. In the end we relate the $Z,N$
and $2,2$ (Helium) Hamiltonians by

$H_{Z,N}(\overline{\mathbf{r}}_{1},\overline{\mathbf{r}}_{2}|1)=\left(\frac{Z}{2}\right)^{2}\frac{1}{N-1}H_{2,2}\Bigg{(}\overline{\mathbf{r}}_{1},\overline{\mathbf{r}}_{2}\Bigg{|}\frac{2(N-1)}{Z}\Bigg{)}.$
(121)

Eq. 121 allows us to diagonalize any atomic Hamiltonian by solving for
eigenstates of the Helium atom with the appropriate values of $\lambda$.
Accounting for coordinate scaling, the geminal basis functions are

$\psi_{Z,N}(\mathbf{r}_{1},\mathbf{r}_{2}|1)=\psi_{2,2}\left(\frac{Z\mathbf{r}_{1}}{2},\frac{Z\mathbf{r}_{2}}{2}\right.\left|\frac{2(Z-1)}{Z}\right).$
(122)

Eqs. 121 and 122 can be used in conjunction with the two-parameter solution
method outlined in Section VI. Finally, we observe that for any neutral atom
with $Z=N$ it suffices to calculate $\sim N(N-1)/2$ eigenstates of
$H_{2,2}(\overline{\mathbf{r}}_{1},\overline{\mathbf{r}}_{1}|\lambda)$ for
$\lambda$ between $0$ and $2$.

##  VII Discussion

We introduced the geminal density matrix as the basic variable to describe a
many-body quantum state. The GDM was found to faithfully represent the
$N$-electron wave function because it reproduces the expectation value of an
arbitrary observable. The GDM is remarkably simple because it derives directly
from the exchange symmetries of the $N$-electron wave function and a given
operator $A$ without postulating the existence of statistical ensembles.
Furthermore, it is analyzed strictly by matrix operations without needing
higher order tensors.

The key behavior of the GDM is that it evolves by unitary matrix
transformation, which we derived by assuming that the two-body contributions
of the Hamiltonian and operator $A$ depended only on position. This likely
imposes no limitations on the applicability of the GDM because it is unclear
whether two-body interactions with derivatives exist. We can certainly solve
for the eigenstates of a system with Coulomb repulsion along with more exotic
interactions like the attractive potential mediated by lattice vibrations in
the theory of superconductivity. In the case that two-body derivative terms
exist, we must restrict the valid operators for which the expectation value is
computable to include only one-body operators and the Hamiltonian itself.

The simple time evolution law allowed us to recover the classical intuition of
electron-nuclei thermalization from an exact quantum-mechanical treatment of
the electron gas. We continued to derive a degenerate adiabatic theorem which
we exploited to calculate the stationary states of an arbitrary many-body
Hamiltonian. We found that the problem of directly minimizing the energy
functional for an $N$-electron wave function reduces to the calculation of
around $N(N-1)/2$ eigenstates of an effective two-electron Hamiltonian on a
grid of electron-electron interaction scaling strengths.

We finally displayed the power of this diagonalization method by applying it
to atomic Hamiltonians, which reduced to an analysis of an effective Helium
atom. So long as we know $\sim Z(Z-1)/2$ eigenstates of the Helium atom with
coulomb interaction scaled by $\lambda$ from $0$ to $2$, we are able to
compute the exact eigenstates of any atom. While the solution to the two-body
problem is not trivial, it provides an incredible speedup for the solution of
the many-body problem.

###### Acknowledgements.

I thank Jens Biegert for the fruitful conversations that inspired this work
along with his guidance and support throughout its completion. I also thank
Eric Van Stryland and David Hagan for critical readings of the manuscript and
acknowledge funding from the US Fulbright Student Program and the Air Force
Office of Scientific Research grant FA9550-20-1-0322.

##  Appendix A Matrix Constraints

This appendix derives constraints on the form of $\mathbf{D}$ that follow
directly from the definitions (Eq. 13)

$\rho(\mathbf{X},\mathbf{X}^{\prime})=\int\Psi^{*}(\mathbf{X},\mathbf{Y})\Psi(\mathbf{X}^{\prime},\mathbf{Y})d\mathbf{Y}$
(123)

and (Eq. 17)

$\rho(\mathbf{X},\mathbf{X}^{\prime})=\sum_{mn}D_{mn}\psi^{*}_{n}(\mathbf{X})\psi_{m}(\mathbf{X}^{\prime}).$
(124)

These restrictions will serve as necessary $N$-representability conditions for
the GDM to represent a valid $N$-electron wave function.

The first rule arises from Eq. 123 which has the property
$\rho(\mathbf{X}^{\prime},\mathbf{X})=\rho^{*}(\mathbf{X},\mathbf{X}^{\prime})$.
Applying this transformation directly to the expansion in Eq. 124 yields

$\displaystyle\sum_{mn}D_{mn}\psi_{n}^{*}(\mathbf{X}^{\prime})\psi_{m}(\mathbf{X})=\sum_{mn}D_{mn}^{*}\psi_{n}(\mathbf{X})\psi_{m}^{*}(\mathbf{X}^{\prime}).$

After swapping the sum index labels on the right hand side, we find the
symmetry to be satisfied when $D_{nm}^{*}=D_{mn}$. This relationship implies
the matrix identity

$\mathbf{D}^{\dagger}=\mathbf{D}.$ (125)

We uncover another $N$-representability requirement by fixing
$\mathbf{X}^{\prime}=\mathbf{X}$ in Eq. 123 and integrating both sides over
the remaining free coordinates $\mathbf{X}$. The integral on the right-hand
side reduces to unity by the normalization condition of the wave function.
Choosing a convenient representation for the left hand side gives

$\int\delta(\mathbf{X}-\mathbf{X}^{\prime})\rho(\mathbf{X},\mathbf{X}^{\prime})d\mathbf{X}d\mathbf{X}^{\prime}=\begin{pmatrix}N\\\
2\end{pmatrix}.$ (126)

Once again expanding $\rho(\mathbf{X},\mathbf{X}^{\prime})$ by Eq. 124, we
find using the orthonormality of the geminal basis that

$\displaystyle\operatorname{Tr}[\mathbf{D}]=\begin{pmatrix}N\\\
2\end{pmatrix}.$ (127)

The requirement for antisymmetry under the exchange
$\mathbf{x}_{1}\leftrightarrow\mathbf{x}_{2}$ (or
$\mathbf{x}_{1}^{\prime}\leftrightarrow\mathbf{x}_{2}^{\prime}$) follows from
that of the many-body wave function, so that
$\rho(\mathbf{x}_{2},\mathbf{x}_{1},\mathbf{X}^{\prime})=-\rho(\mathbf{X},\mathbf{X}^{\prime})$.
Swapping these coordinates in the two-body expansion of Eq. 124 gives

$\displaystyle\rho(\mathbf{x}_{2},\mathbf{x}_{1},\mathbf{X}^{\prime})$
$\displaystyle=\sum_{mn}D_{mn}\psi_{n}^{*}(\mathbf{x}_{2},\mathbf{x}_{1})\psi_{m}(\mathbf{X}^{\prime})$
$\displaystyle=-\rho(\mathbf{X},\mathbf{X}^{\prime}),$ (128)

indicating that the property is inherited from the anti-symmetry of the basis
functions and does not further restrict $\mathbf{D}$.

Unfortunately, we have now found all the $N$-representability conditions that
follow directly from Eqs. 123 and 124. To further understand the
$N$-representability problem we must choose a basis and derive expressions for
the matrix elements $D_{mn}$. We compute these matrix elements by pre-
multiplying Eq. 124 by $\psi_{n}(\mathbf{X})\psi_{m}^{*}(\mathbf{X}^{\prime})$
and integrating over $d\mathbf{X}$ and $d\mathbf{X}^{\prime}$ to find

$D_{mn}=\int\psi_{n}(\mathbf{X})\psi_{m}^{*}(\mathbf{X}^{\prime})\rho(\mathbf{X},\mathbf{X}^{\prime})d\mathbf{X}d\mathbf{X}^{\prime}.$
(129)

Continuing to substitute Eq. 123 into Eq. 129 yields the simple equation

$D_{mn}=\int\Theta_{n}^{*}(\mathbf{Y})\Theta_{m}(\mathbf{Y})d\mathbf{Y},$
(130)

with overlap functions $\Theta_{m}(\mathbf{Y})$ defined to be

$\Theta_{m}(\mathbf{Y})=\int\psi_{m}^{*}(\mathbf{X})\Psi(\mathbf{X},\mathbf{Y})d\mathbf{X}.$
(131)

We choose the geminal basis functions $\psi_{i}(\mathbf{X})$ to be those
formed by the anti-symmetrized product of two single particle spinors
$\phi_{i}(\mathbf{x})$. Grouping the two index labels into the symbol
$\mathbf{n}=\\{n_{1},n_{2}\\}$, the basis functions take the form

$\psi_{\mathbf{n}}(\mathbf{X})=\frac{1}{\sqrt{2}}\left[\phi_{n_{1}}(\mathbf{x}_{1})\phi_{n_{2}}(\mathbf{x}_{2})-\phi_{n_{2}}(\mathbf{x}_{1})\phi_{n_{1}}(\mathbf{x}_{2})\right].$
(132)

Because our basis functions are labeled by two integers, our $2$-RDM expansion
will temporarily take the form of a rank four tensor with components
$D_{\mathbf{mn}}$. It will eventually be necessary to map each $\mathbf{n}$ to
a single integer index to flatten this tensor into a matrix (see the Table 1
and the surrounding discussion).

The most general wave function $\Psi(\mathbf{X},\mathbf{Y})$ that will appear
in Eq. 131 is a possibly infinite linear superposition of $N$-body Slater
determinants

$\Psi(\mathbf{X},\mathbf{Y})=\sum_{\\{\alpha\\}}C_{\\{\alpha\\}}\Psi_{\\{\alpha\\}}(\mathbf{X},\mathbf{Y}).$
(133)

As in the main body of this text, we defined a configuration $\\{\alpha\\}$ to
be an ordered list of single-particle spinors present in a given product of
states. The Slater determinants are formed by the antisymmetrization operator

$\Psi(\mathbf{x}_{1},\dots,\mathbf{x}_{N})=\hat{S}_{-}\prod_{i=1}^{N}\phi_{\alpha_{i}}(\mathbf{x}_{i})$
(134)

equivalent to the determinant expression

$\Psi_{\\{\alpha\\}}(\mathbf{X},\mathbf{Y})=\frac{1}{\sqrt{N!}}\left|\begin{matrix}\phi_{\alpha_{1}}(\mathbf{x}_{1})&\phi_{\alpha_{2}}(\mathbf{x}_{1})&\dots&\phi_{\alpha_{N}}(\mathbf{x}_{1})\\\
\phi_{\alpha_{1}}(\mathbf{x}_{2})&\phi_{\alpha_{2}}(\mathbf{x}_{2})&\dots&\phi_{\alpha_{N}}(\mathbf{x}_{2})\\\
\vdots&\vdots&\ddots&\vdots\\\
\phi_{\alpha_{1}}(\mathbf{x}_{N})&\phi_{\alpha_{2}}(\mathbf{x}_{N})&\dots&\phi_{\alpha_{N}}(\mathbf{x}_{N})\end{matrix}\right|.$
(135)

We now expand each $\Psi_{\\{\alpha\\}}(\mathbf{X},\mathbf{Y})$ in Eq. 133
along minors of the top two rows of Eq. 135 to isolate the $\mathbf{x}_{1}$
and $\mathbf{x}_{2}$ dependence. The result is

$\displaystyle\Psi(\mathbf{X},\mathbf{Y})=\frac{1}{\sqrt{N(N-1)}}\sum_{\\{\alpha\\}}C_{\\{\alpha\\}}\sum_{i,j>i}(-1)^{i+j-1}\left[\phi_{\alpha_{i}}(\mathbf{x}_{1})\phi_{\alpha_{j}}(\mathbf{x}_{2})-\phi_{\alpha_{j}}(\mathbf{x}_{1})\phi_{\alpha_{i}}(\mathbf{x}_{2})\right]\Psi_{\\{\alpha\\}_{ij}}(\mathbf{Y}),$
(136)

where the reduced configuration
$\\{\alpha\\}_{ij}=\\{\alpha\\}\setminus\\{\alpha_{i},\alpha_{j}\\}$ is the
set subtraction of $\alpha_{i}$ and $\alpha_{j}$ from the original list of
states. The state $\Psi_{\\{\alpha\\}_{ij}}(\mathbf{Y})$ is the determinant of
the matrix formed by removing rows $1$ and $2$ and columns $i$ and $j$ from
Eq. 135. As the normalization of this $N-2$ electron state requires the
prefactor $1/\sqrt{(N-2)!}$, we multiplied by its inverse which partially
canceled with the $1/\sqrt{N!}$ prefactor.

Continuing to normalize the $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$ dependence
into a two-electron Slater determinant
$\psi_{\alpha_{i}\alpha_{j}}(\mathbf{X})$, we finally have

$\displaystyle\Psi(\mathbf{X},\mathbf{Y})=\sqrt{\frac{2}{N(N-1)}}$
$\displaystyle\sum_{\\{\alpha\\}}C_{\\{\alpha\\}}\sum_{i,j>i}(-1)^{i+j-1}$
$\displaystyle\times\psi_{\alpha_{i}\alpha_{j}}(\mathbf{X})\Psi_{{\\{\alpha\\}}_{ij}}(\mathbf{Y}).$
(137)

Plugging Eq. A into Eq. 131 for the overlap $\Theta_{\mathbf{m}}(\mathbf{Y})$,
the integration over $d\mathbf{X}$ reduces the two-electron wave functions to
$\delta_{\alpha_{i},m_{1}}\delta_{\alpha_{j},m_{2}}$ by orthonormality. Thus,
a given configuration that does not contain $\mathbf{m}=\\{m_{1},m_{2}\\}$
will not contribute to $D_{\mathbf{mn}}$. Consequently, we may reduce the
$\\{\alpha\\}$ (configuration) sum into one over
$\\{\alpha\\}\ni\\{m_{1},m_{2}\\}$. The remaining sum is reduced to the single
term with $(i,j)=\\{m_{1},m_{2}\\}$ so that

$\displaystyle\Theta_{\mathbf{m}}(\mathbf{Y})=\sqrt{\frac{2}{N(N-1)}}\sum_{\\{\alpha\\}\ni\mathbf{m}}$
$\displaystyle
C_{\\{\alpha\\}}\mathcal{S}_{\alpha}[\mathbf{m}]\Psi_{{\\{\alpha\\}}_{\mathbf{m}}}(\mathbf{Y}),$
(138)

where
$\\{\alpha\\}_{\mathbf{m}}=\\{\alpha\\}_{m_{1}m_{2}}=\\{\alpha\\}\setminus\\{m_{1},m_{2}\\}$.
The symbol $S_{\alpha}[\mathbf{m}]$ is the sign function

$S_{\alpha}[\mathbf{m}]=(-1)^{I_{\alpha}[m_{1}]+I_{\alpha}[m_{2}]-1}$ (139)

with $I_{\alpha}[p]$ the index of basis function $p$ in configuration
$\\{\alpha\\}$. We absorb this sign into the expansion coefficient by defining
$\mathcal{C}_{\\{\alpha\\}}=C_{\\{\alpha\\}}S_{\alpha}[\mathbf{m}]$.

Finally, using that $\Theta_{n}^{*}(\mathbf{Y})$ is the complex conjugate of
Eq. 138, we compute $D_{\mathbf{mn}}$ by Eq. 130:

$\displaystyle D_{\mathbf{mn}}$
$\displaystyle=\sum_{\begin{subarray}{c}\\{\alpha\\}\ni\mathbf{m}\\\
\\{\beta\\}\ni\mathbf{n}\end{subarray}}\mathcal{C}^{*}_{\\{\beta\\}}\mathcal{C}_{\\{\alpha\\}}\int\Psi^{*}_{{\\{\beta\\}}_{\mathbf{n}}}(\mathbf{Y})\Psi_{{\\{\alpha\\}}_{\mathbf{m}}}(\mathbf{Y})d\mathbf{Y}.$
(140)

The integral, being the inner product between orthonormal $N-2$ electron
Slater determinants, equals one when
$\\{\alpha\\}_{\mathbf{n}}=\\{\beta\\}_{\mathbf{m}}$ and zero otherwise.
Therefore,

$\displaystyle D_{\mathbf{m}\mathbf{n}}$
$\displaystyle=\sum_{\begin{subarray}{c}\\{\alpha\\}\ni\mathbf{n}\\\
\\{\beta\\}\ni\mathbf{m}\end{subarray}}\mathcal{C}^{*}_{\\{\beta\\}}\mathcal{C}_{\\{\alpha\\}}\delta_{\\{\alpha\\}_{\mathbf{n}},\\{\beta\\}_{\mathbf{m}}}.$
(141)

The diagonal matrix elements are then found from Eq. 141 to take the simple
form

$\displaystyle
D_{\mathbf{n}\mathbf{n}}=\sum_{\\{\alpha\\}\ni\mathbf{n}}\left|C_{\\{\alpha\\}}\right|^{2}.$
(142)

Because Eq. 142 is a sum over the magnitude squared of all expansion
coefficients of configurations containing $\mathbf{n}$, the overall
normalization condition $\sum_{\\{\alpha\\}}|C_{\\{\alpha\\}}|^{2}=1$ implies
that

$0\leq D_{\mathbf{nn}}\leq 1.$ (143)

The maximum diagonal value, $D_{\mathbf{nn}}=1$, occurs when every
configuration $\\{\alpha\\}$ contains $\mathbf{n}$. In this case we encounter
the additional rule that

$D_{\mathbf{nn}}=1\implies\forall\mathbf{m}\neq\mathbf{n},D_{\mathbf{mn}}=D_{\mathbf{nm}}=0,$
(144)

meaning that $1$ on the diagonal in position $\mathbf{n}$ forces all other
elements in the column and row $\mathbf{n}$ to zero.

The proof of Eq. 144 proceeds as follows. Per Eq. 142, any non-zero term must
simultaneously satisfy the conditions $\\{\alpha\\}\ni\\{m_{1},m_{2}\\}$,
$\\{\beta\\}\ni\\{n_{1},n_{2}\\}$ and
$\\{\alpha\\}\setminus\\{m_{1},m_{2}\\}=\\{\beta\\}\setminus\\{n_{1},n_{2}\\}$.
By this equality, $\\{\alpha\\}\setminus\\{m_{1},m_{2}\\}$ does not contain
$\\{n_{1},n_{2}\\}$. Since $\\{\alpha\\}$ is formed by the set addition of
some $\\{m_{1},m_{2}\\}\neq\\{n_{1},n_{2}\\}$, we have that
$\\{\alpha\\}\not\ni\\{n_{1},n_{2}\\}$. Supposing now that
$D_{\mathbf{mn}}\neq 0$ implies existence of some $\\{\alpha\\}$ in the state
expansion that does not contain $\\{n_{1},n_{2}\\}$. This contradicts the
requirement that must be met for $D_{\mathbf{nn}}=1$ so we conclude that the
existence of $1$ on a diagonal implies all other elements in that row and
column are $0$.

By the Hermiticity (Eq. 125) of $\mathbf{D}$, it can always be transformed
into diagonal form by a unitary basis transformation (see Appendix B). The
resulting diagonal matrix obeying Eq. 143 has the property

$0\leq\operatorname{Tr}[\mathbf{D}^{2}]\leq\begin{pmatrix}N\\\
2\end{pmatrix},$ (145)

which follows trivially from the fact that $a^{2}\leq a$ for a number $a\leq
1$. $\operatorname{Tr}[\mathbf{D}^{2}]$ is a basis-indepedent quantity as it
is invariant under unitary transformation by the cyclic property of the trace.

We finally summarize the necessary $N$-representability constraints on the
matrix $\mathbf{D}$:

$\displaystyle\mathbf{D}$ $\displaystyle=\mathbf{D}^{\dagger}$ (146a)
$\displaystyle 0\leq$ $\displaystyle D_{nn}\leq 1$ (146b)
$\displaystyle\operatorname{Tr}[\mathbf{D}]$
$\displaystyle=\begin{pmatrix}N\\\ 2\end{pmatrix}$ (146c) $\displaystyle
0\leq\operatorname{Tr}[\mathbf{D}^{2}]$ $\displaystyle\leq\begin{pmatrix}N\\\
2\end{pmatrix}.$ (146d)

##  Appendix B Change of basis

The formula for a change of basis is identical to the transformation for any
density matrix, but we re-derive it here for completeness. We begin with the
two-body density matrix

$\displaystyle\mathbf{D}=\sum_{mn}D_{mn}\psi^{*}_{n}(\mathbf{X})\psi_{m}(\mathbf{X})$
(147)

and introduce a new orthonormal basis with wave functions
$\phi_{i}(\mathbf{X})$. Expanding the initial states $\psi_{n}(\mathbf{X})$ in
terms of the new, we find

$\displaystyle\mathbf{D}$
$\displaystyle=\sum_{mn}D_{mn}\left(\sum_{j}U^{*}_{jn}\phi^{*}_{j}(\mathbf{X})\right)\left(\sum_{i}U_{im}\phi_{i}(\mathbf{X}^{\prime})\right)$
$\displaystyle=\sum_{ij}\left(\sum_{mn}U_{im}D_{mn}U^{*}_{jn}\right)\phi^{*}_{j}(\mathbf{X})\phi_{i}(\mathbf{X}^{\prime}).$
(148)

The parenthetical term gives the expression for $D_{ij}$, which we take to be
the coefficients of matrix $\mathbf{D}^{\prime}$. We finally find the matrix
form for the change of basis

$\mathbf{D}^{\prime}=\mathbf{U}\mathbf{D}\mathbf{U}^{\dagger},$ (149)

where $\mathbf{U}$ is the unitary matrix of coefficients $U_{ij}$.

##  Appendix C Alternate derivation of the GDM evolution equation

We found the matrix Liouville-von Neumann equation (Eq. 60) by deriving its
most general operator equivalent then specializing to the case of a time-
independent basis. We can perform an alternate derivation by starting with the
time-independent expansion of the $2$-RDM

$\rho(\mathbf{X},\mathbf{X}^{\prime}|t)=\sum_{mn}D_{mn}(t)\psi_{n}^{*}(\mathbf{X})\psi_{m}(\mathbf{X}^{\prime}).$
(150)

We further assume that our dummy operator $A(\mathbf{X})$ is time independent.
Noting that $(d/dt)\braket{A}(t)=-iK(t)$ for a time-independent
$A(\mathbf{X})$, we restart from Eq. IV,

$\frac{d}{dt}\braket{A}(t)=-i\int\psi_{n}^{*}(\overline{\mathbf{X}}|t)[H(\mathbf{X}|t),A(\mathbf{X})]\psi_{m}(\overline{\mathbf{X}}|t)d\overline{\mathbf{X}}.$
(151)

Expressing Eq. 151 in terms of the $2$-RDM as

$\displaystyle\frac{d}{dt}\braket{A}(t)$ $\displaystyle=\int
d\mathbf{X}d\mathbf{X}^{\prime}\delta(\mathbf{X}-\mathbf{X}^{\prime})$
$\displaystyle\times\left[H^{\prime}(\mathbf{X}|t),A(\mathbf{X})\right]\rho(\mathbf{X},\mathbf{X}^{\prime}|t),$
(152)

we can apply the expansion in Eq. 150 to Eq. C. In terms of the abstract
effective two-electron operators, the result is

$\frac{d}{dt}\braket{A}(t)=\sum_{mn}D_{mn}[\braket{m}{\hat{H}^{\prime}(t)\hat{A}}{n}-\braket{m}{\hat{A}\hat{H}^{\prime}(t)}{n}],$
(153)

into which we can insert the identity $1=\sum_{i}\ket{i}\bra{i}$ to form the
matrix equation

$\frac{d}{dt}\braket{A}(t)=-i\operatorname{Tr}\left[\mathbf{D}(t)[\mathbf{H}^{\prime}(t),\mathbf{A}]\right].$
(154)

For comparison, directly differentiating the trace relation
$\braket{A}(t)=\operatorname{Tr}[\mathbf{D}(t)\mathbf{A}]$ gives

$\frac{d}{dt}\braket{A}(t)=\operatorname{Tr}[\dot{\mathbf{D}}(t)\mathbf{A}],$
(155)

whose result must agree with Eq. 154.

Using the cyclic property and linearity of the trace, we can rearrange Eq. 154
as
$\operatorname{Tr}[\mathbf{D}\mathbf{H}^{\prime}\mathbf{A}-\mathbf{D}\mathbf{A}\mathbf{H}^{\prime}]=\operatorname{Tr}[\mathbf{D}\mathbf{H}^{\prime}\mathbf{A}-\mathbf{H}^{\prime}\mathbf{D}\mathbf{A}]$
so that Eq. 154 takes the form

$\frac{d}{dt}\braket{A}(t)=\operatorname{Tr}[-i[\mathbf{D}(t),\mathbf{H}^{\prime}(t)]\mathbf{A}(t)].$
(156)

Comparing Eq. 156 to Eq. 155 yields the condition

$\dot{\mathbf{D}}(t)=-i[\mathbf{D}(t),\mathbf{H}^{\prime}(t)]$ (157)

to ensure the GDM gives the same expectation values as the wave function.

## References

  * Löwdin [1955a] P. O. Löwdin, Quantum theory of many-particle systems. I. Physical interpretations by means of density matrices, natural spin-orbitals, and convergence problems in the method of configurational interaction, Physical Review 97, 1474 (1955a).
  * Löwdin [1955b] P.-O. Löwdin, Quantum Theory of Many-Particle Systems. II. Study of the Ordinary Hartree-Fock Approximation, Physical Review 97, 1490 (1955b).
  * Löwdin [1955c] P. O. Löwdin, Quantum theory of many-particle systems. III. Extension of the Hartree-Fock scheme to include degenerate systems and correlation effects, Physical Review 97, 1509 (1955c).
  * Mayer [1955] J. E. Mayer, Electron Correlation, Physical Review 100, 1579 (1955).
  * Coulson [1960] C. A. Coulson, Present state of molecular structure calculations, Reviews of Modern Physics 32, 170 (1960).
  * Coleman [1963] A. J. Coleman, Structure of fermion density matrices, Reviews of Modern Physics 35, 668 (1963).
  * Mazziotti [2012] D. A. Mazziotti, Two-Electron Reduced Density Matrix as the Basic Variable in Many-Electron Quantum Chemistry and Physics, Chemical Reviews 112, 244 (2012).
  * Nakatsuji [1976] H. Nakatsuji, Equation for the direct determination of the density matrix, Physical Review A 14, 41 (1976).
  * Nakatsuji and Yasuda [1996] H. Nakatsuji and K. Yasuda, Direct Determination of the Quantum-Mechanical Density Matrix Using the Density Equation, Phys. Rev. Lett. 76, 1039 (1996).
  * Mazziotti [1998] D. A. Mazziotti, Contracted Schrödinger equation: Determining quantum energies and two-particle density matrices without wave functions, Phys. Rev. A 57, 4219 (1998).
  * Mazziotti [1999] D. A. Mazziotti, Comparison of contracted Schrödinger and coupled-cluster theories, Phys. Rev. A 60, 4396 (1999).
  * Garrod and Percus [1964] C. Garrod and J. K. Percus, Reduction of the N ‐Particle Variational Problem, Journal of Mathematical Physics 5, 1756 (1964).
  * Mihailović and Rosina [1975] M. Mihailović and M. Rosina, The variational approach to the density matrix for light nuclei, Nuclear Physics A 237, 221 (1975).
  * Mazziotti and Erdahl [2001] D. A. Mazziotti and R. M. Erdahl, Uncertainty relations and reduced density matrices: Mapping many-body quantum mechanics onto four particles, Physical Review A - Atomic, Molecular, and Optical Physics 63, 1 (2001).
  * Nakata  _et al._ [2001] M. Nakata, H. Nakatsuji, M. Ehara, M. Fukuda, K. Nakata, and K. Fujisawa, Variational calculations of fermion second-order reduced density matrices by semidefinite programming algorithm, The Journal of Chemical Physics 114, 8282 (2001).
  * Mazziotti [2004] D. A. Mazziotti, Realization of Quantum Chemistry without Wave Functions through First-Order Semidefinite Programming, Phys. Rev. Lett. 93, 213001 (2004).
  * Vandenberghe and Boyd [1996] L. Vandenberghe and S. Boyd, Semidefinite Programming, SIAM Review 38, 49 (1996), publisher: Society for Industrial and Applied Mathematics.
  * Surján [1999] P. R. Surján, An Introduction to the Theory of Geminals, in  _Correlation and Localization_, Vol. 203, edited by A. de Meijere, K. N. Houk, H. Kessler, J.-M. Lehn, S. V. Ley, S. L. Schreiber, J. Thiem, B. M. Trost, F. Vögtle, H. Yamamoto, P. R. Surján, R. J. Bartlett, F. Bogár, D. L. Cooper, B. Kirtman, W. Klopper, W. Kutzelnigg, N. H. March, P. G. Mezey, H. Müller, J. Noga, J. Paldus, J. Pipek, M. Raimondi, I. Røeggen, J. Q. Sun, P. R. Surján, C. Valdemoro, and S. Vogtner (Springer Berlin Heidelberg, Berlin, Heidelberg, 1999) pp. 63–88, series Title: Topics in Current Chemistry.
  * Bopp [1959] F. Bopp, Ableitung der Bindungsenergie von N-Teilchen-Systemen aus 2-Teilchen-Dichtematrizen, Zeitschrift für Physik 156, 348 (1959).
  * Aguiar Pinto  _et al._ [2002] A. Aguiar Pinto, K. Fonseca Romero, and M. Thomaz, Adiabatic approximation in the density matrix approach: non-degenerate systems, Physica A: Statistical Mechanics and its Applications 311, 169 (2002).