text:
\section{Introduction}\label{sec:Intro}
In 1955, L\"owdin  \cite{lowdinQuantumTheoryManyparticle1955a,lowdinQuantumTheoryManyParticle1955,lowdinQuantumTheoryManyparticle1955b} and Mayer \cite{mayerElectronCorrelation1955} 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 \cite{coulsonPresentStateMolecular1960}. 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 \cite{colemanStructureFermionDensity1963}, 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. \cite{mazziottiTwoElectronReducedDensity2012}, modern $2$-RDM analysis primarily employs one of two methods. The first begins by deriving a contracted Schr\"odinger equation (CSE) that computes the energy as a function of the two and four-body reduced density matrices ($4$-RDM) \cite{nakatsujiEquationDirectDetermination1976,nakatsujiDirectDeterminationQuantumMechanical1996,mazziottiContractedSchrodingerEquation1998,mazziottiComparisonContractedSchrodinger1999}. 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 \cite{garrodReductionParticleVariational1964,mihailovicVariationalApproachDensity1975,mazziottiUncertaintyRelationsReduced2001,nakataVariationalCalculationsFermion2001,mazziottiRealizationQuantumChemistry2004}. The most successful applications of DMVT use a convex optimization scheme called semi-definite programming \cite{vandenbergheSemidefiniteProgramming1996,nakataVariationalCalculationsFermion2001} 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 \textit{geminal} \cite{surjanIntroductionTheoryGeminals1999} 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 \cite{boppAbleitungBindungsenergieNTeilchenSystemen1959}. 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 \cite{colemanStructureFermionDensity1963}. 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 \ref{sec:Density} introduces the $2$-RDM then expands it in a two-electron basis to define the GDM. Section \ref{sec:States} gives examples of valid density matrices that serve as the starting point for eigenstate calculations. Section \ref{sec:Time} derives the unitary time evolution law of the GDM, which Section \ref{sec:Evolution} 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 \ref{sec:Solve} to solve the many-body Schr\"odinger equation.

Appendix \ref{sec:Properties} derives the four necessary $N$-representability constraints listed in Section \ref{sec:Density}. Appendix \ref{sec:Basis} defines the matrix transformation imposed by a change of geminal basis, and Appendix \ref{sec:Alternate} provides an alternate derivation of the GDM equation of motion found in Section \ref{sec:Time}.

Hartree atomic units with $\hbar = e = m_0 = 1/(4\pi \epsilon_0) = 1$ will be used throughout.

\section{The geminal density matrix}\label{sec:Density}
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 \ref{sec:Properties}.

\subsection{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
\begin{equation}
\mathbf{x}_i = \mathbf{r}_i\sigma_i.
\end{equation}
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
\begin{equation}\label{eq:norm}
	\int |\Psi(\mathbf{x}_1,\dots,\mathbf{x}_N)|^2 \prod_{i=1}^N d\mathbf{x}_i = 1,
\end{equation}
with the integral over $d\mathbf{x}$ including a sum over spin coordinates by
\begin{equation}\label{eq:spinsum}
	\int f(\mathbf{x}) d\mathbf{x} = \sum_{\sigma} \int f(\mathbf{r}\mathbf{\sigma}) d\mathbf{r}.
\end{equation}

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
\begin{equation}\label{eq:sop}
	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).
\end{equation}
We will find it convenient to combine Eq. \ref{eq:sop} into the single sum over pairs
\begin{equation}\label{eq:Atot}
	A(\mathbf{x}_1,\dots,\mathbf{x}_N) = \sum_{i,j>i} A(\mathbf{x}_i,\mathbf{x}_j),
\end{equation}
with the summand $A(\mathbf{x}_i,\mathbf{x}_j)$ defined by
\begin{equation}\label{eq:Asep1}
	A(\mathbf{x}_i,\mathbf{x}_j) = A_1(\mathbf{x}_i,\mathbf{x}_j) + A_2(\mathbf{x}_i,\mathbf{x}_j).
\end{equation}
The term $A_1(\mathbf{x}_i,\mathbf{x}_j)$ is the promotion of a one-body operator to act on pairs
\begin{equation}\label{eq:Aeff1}
	A_1(\mathbf{x}_i,\mathbf{x}_j) =  \frac{a_1(\mathbf{x}_i) + a_1(\mathbf{x}_j)}{N-1},
\end{equation}
where the denominator divides out the overcounting that occurs in the pair sum.

Using Eq. \ref{eq:Atot}, we calculate the expectation value of $A$ by the inner product
\begin{align}\label{eq:AA1}
	\braket{A} &= \sum_{i,j>i}\int \Psi^*(\mathbf{x}_1,\dots,\mathbf{x}_N) \nonumber \\
	&\times A(\mathbf{x}_i,\mathbf{x}_j)\Psi(\mathbf{x}_1,\dots,\mathbf{x}_N)\prod_{i=1}^N d\mathbf{x}_i.
\end{align}
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
\begin{equation}\label{eq:A11}
	\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},
\end{equation}
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. \ref{eq:Asep1},
\begin{align}\label{eq:A2}
	A(\mathbf{X}) = \frac{a_1(\mathbf{x}_1) + a_2(\mathbf{x}_2)}{N-1} + A_2(\mathbf{X}).
\end{align}

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}'$ upon which the operator $A(\mathbf{X}')$ is taken to act. This definition permits the re-formulation of Eq. \ref{eq:A11} as
\begin{align}\label{eq:A1}
	\braket{A} = &\begin{pmatrix}N\\2\end{pmatrix}\int d\mathbf{X} d\mathbf{X}' \delta(\mathbf{X}' - \mathbf{X}) \nonumber \\
	&\times A(\mathbf{X}')\int \Psi^*(\mathbf{X},\mathbf{Y})\Psi(\mathbf{X}',\mathbf{Y})d\mathbf{Y}.
\end{align}
In Eq. \ref{eq:A1} we defined the delta function
\begin{equation}
	\delta(\mathbf{X}' - \mathbf{X}) = \prod_{i=1}^N \delta(\mathbf{r}_i' - \mathbf{r}_i)\delta_{\sigma_i',\sigma_{i}}
\end{equation}
that includes Dirac delta functions for position and Kronecker delta functions for spin polarization.

Finally, we define the $2$-RDM
\begin{equation}\label{eq:rho1}
	\rho(\mathbf{X},\mathbf{X}') = \begin{pmatrix}N\\2\end{pmatrix}\int \Psi^*(\mathbf{X},\mathbf{Y})\Psi(\mathbf{X}',\mathbf{Y})d\mathbf{Y},
\end{equation}
so that Eq. \ref{eq:A1} simplifies to
\begin{equation}\label{eq:At1}
	\braket{A} = \int d\mathbf{X} d\mathbf{X}' \delta(\mathbf{X}' - \mathbf{X}) A(\mathbf{X}')\rho(\mathbf{X},\mathbf{X}').
\end{equation}
Eqs. \ref{eq:rho1} and \ref{eq:At1} provide an enormous complexity reduction compared to Eq. \ref{eq:AA1}, as they depend only on $\mathbf{X}$ and $\mathbf{X}'$ regardless of the number of particles under investigation.

\subsection{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}$:
\begin{equation}\label{eq:rhopos}
	\rho(\mathbf{X},\mathbf{X}') = \braket{\mathbf{X}'|\hat{D}|\mathbf{X}}.
\end{equation}
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
\begin{equation}\label{eq:Dab}
	\hat{D} = \sum_{mn}\ket{m}\braket{m|\hat{D}|n}\bra{n}.
\end{equation}
Defining $D_{mn} = \braket{m|\hat{D}|n}$ and taking the position representation by pre-multiplying $\bra{\mathbf{X}'}$ and post-multiplying $\ket{\mathbf{X}}$, the $2$-RDM as expressed in Eq. \ref{eq:rhopos} takes the form
\begin{equation}\label{eq:rhoexp}
	\rho(\mathbf{X},\mathbf{X}') = \sum_{mn}D_{mn}\psi_n^*(\mathbf{X})\psi_m(\mathbf{X}').
\end{equation}
We can prove the validity of this expansion by construction; applying $\int \psi_{n'}(\mathbf{X}) \psi_{m'}^*(\mathbf{X}') d\mathbf{X} d\mathbf{X}'$ to either side of Eq. \ref{eq:rhoexp} isolates element $D_{m'n'}$ as a functional of the $2$-RDM. See Eq. \ref{eq:Dmn1} in Appendix \ref{sec:Properties} for more detail.

Inserting Eq. \ref{eq:rhoexp} into Eq. \ref{eq:At1} yields
\begin{equation} \label{eq:Ata}
	\braket{A} =  \sum_{mn}D_{mn} A_{nm},
\end{equation}
where
\begin{equation}
	A_{nm} = \int \psi_n^*(\mathbf{X})A(\mathbf{X})\psi_m(\mathbf{X}) d\mathbf{X}.
\end{equation}
We now define the GDM as the matrix $\mathbf{D}$ of coefficients $D_{mn}$. Similarly defining $\mathbf{A}$ to have coefficients $A_{mn}$, Eq. \ref{eq:Ata} reduces to the matrix form
\begin{equation}\label{eq:tr1}
	\braket{A} = \Tr[\mathbf{D}\mathbf{A}].
\end{equation}

Eq. \ref{eq:tr1} 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 \ref{sec:Properties} we derive the following four necessary conditions:
\begin{subequations}
\label{eq:ruless}
\begin{align}
    \mathbf{D} &= \mathbf{D}^\dagger \label{eq:hermitian0}\\
    0\leq &D_{nn} \leq 1 \label{eq:occupation0}\\
    \Tr[\mathbf{D}] &= \begin{pmatrix}
    N \\ 2
    \end{pmatrix}  \label{eq:trace0}\\
     0\leq \Tr[\mathbf{D}^2] &\leq \begin{pmatrix}
    N \\ 2
    \end{pmatrix}. \label{eq:tracesquared0}
\end{align}
\end{subequations}
Condition (\ref{eq:hermitian0}) says the matrix must be Hermitian. We can interpret rule (\ref{eq:occupation0}) as the maximum occupation number for a given geminal and (\ref{eq:trace0}) as the total number of electron pairs present in the wave function. The final expression, (\ref{eq:tracesquared0}), 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 $\Tr[\mathbf{D}^2] = N(N-1)/2$.

In Section \ref{sec:States} we find that Eq. \ref{eq:ruless} is an insufficient set of constraints because it is possible to define a non-$N$-representable GDM that obeys these rules.

\section{Matrix Examples}\label{sec:States}
This section begins with examples of simple GDMs followed by a demonstration of the $N$-representability problem. Subsection \ref{sec:real} 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.

\subsection{Matrices with a non-interacting geminal basis}\label{sec:slater}

Using an orthonormal basis of one-electron wave functions $\phi_i(\mathbf{x})$, we can build $N$-electron Slater determinants by
\begin{equation}\label{eq:basis1}
	\Psi_{\conf{\alpha}}(\mathbf{x}_1,\dots,\mathbf{x}_N) = \hat{S}_- \prod_{i=1}^N \phi_{\alpha_i}(\mathbf{x}_{i}).
\end{equation}
In Eq. \ref{eq:basis1}, we defined a configuration $\conf{\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. \ref{eq:symop1} and \ref{eq:slater}.

Per Eq. \ref{eq:rhoexp}, 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
\begin{equation}\label{eq:nibasis}
	\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)
\end{equation}
labeled by the pair of integers $\mathbf{n} = \conf{n_1,n_2}$. For a many-electron state defined by a single configuration $\conf{\alpha}$, we find in Appendix \ref{sec:Properties} that the $2$-RDM can be represented as a rank four tensor indexed by the pairs
\begin{equation}\label{eq:Drulenn}
	D_{\mathbf{mn}} = \begin{cases}
		1, & \text{if } \mathbf{m}=\mathbf{n}\text{ and } n_1,n_2 \in \conf{\alpha} \\
		0, & \text{otherwise}.
	\end{cases}
\end{equation}
To transform Eq. \ref{eq:Drulenn} 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 \ref{tab:map}.
\begin{table}[h]
	\centering
	\caption{Geminal index map}
	\label{tab:map}
	\begin{tabular}[t]{lccccccc}
		\toprule
		$\mathbf{n}$ & $(1,2)$ & $(1,3)$ &$(2,3)$ & $(1,4)$ & $(2,4)$& $(3,4)$&\dots\\
		$n$ & $1$ & $2$ & $3$ & $4$ & $5$ & $6$ &$\dots$ \\
		\bottomrule
	\end{tabular}
\end{table}

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. \ref{eq:Drulenn} it is clear that $D_{\mathbf{nn}} = 1$ for $\mathbf{n} = \conf{1,2},\conf{1,3}$ and $\conf{2,3}$. Converting to a matrix by Table \ref{tab:map}, 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. \ref{eq:tracesquared0}
\begin{equation}\label{eq:trtr1}
	\Tr[\mathbf{D}^2] = \begin{pmatrix}N\\2\end{pmatrix}
\end{equation}
and the equivalent idempotence condition
\begin{equation}\label{eq:DD1}
	\mathbf{D}^2 = \mathbf{D}.
\end{equation}

It is tempting to conclude that we can arbitrarily place ones on the diagonal, but this assumption fails immediately for the matrix $\mathbf{D}'$ with $D'_{nn} = 1$ for $n = 1,2,4$. According to Table \ref{tab:map}, the Slater determinant which generates this matrix must contains pairs $\conf{1,2},\conf{1,3},$ and $\conf{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. \ref{eq:ruless} are insufficient to ensure $N$-representability because $\mathbf{D}'$ satisfies all the rules and fails to correspond to a valid wave function.

We may also create valid states that do not satisfy Eqs. \ref{eq:trtr1} and \ref{eq:DD1} by taking a linear superposition of configurations
\begin{equation}
	\Psi(\mathbf{x}_1,\dots,\mathbf{x}_N)= \sum_{\conf{\alpha}}C_{\conf{\alpha}}\Psi_{\conf{\alpha}}(\mathbf{x}_1,\dots,\mathbf{x}_N),
\end{equation}
whose diagonal elements are (Eq. \ref{eq:Dmnfinal})
\begin{align}\label{eq:Ddisj}
	D_{\mathbf{n}\mathbf{n}} = \sum_{\conf{\alpha}\ni \mathbf{n}}\left|C_{\conf{\alpha}}\right|^2.
\end{align}
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. \ref{eq:Dndiag}). Take for example the constant superposition of $M$ disjoint anti-symmetrized configurations
\begin{align}\label{eq:Psidis}
	\Psi&(\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) \nonumber \\
	&+\Big(\phi_{N+1}(\mathbf{x}_1)\cdots\phi_{2N}(\mathbf{x}_N)\Big) \nonumber \\
	&+ \dots + \Big(\phi_{(M-1)N+1}(\mathbf{x}_1)\cdots\phi_{MN}(\mathbf{x}_N)\Big)\Big],
\end{align}
where each $\phi_i(\mathbf{x}) \neq \phi_j(\mathbf
x)$. Clearly, the off-diagonal elements of Eq. \ref{eq:Psidis} are all zero. We then find from Eq. \ref{eq:Ddisj} 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 $\Tr[\mathbf{D}^2] = N(N-1)/(2M^2)$, which decreases as the number $M$ of disjoint configurations increases.

\subsection{Matrices representing a real system}\label{sec:real}
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
\begin{equation}\label{eq:Hn}
	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)
\end{equation}
for an $N$-electron Hamiltonian $H(\mathbf{x}_1,\dots,\mathbf{x}_N)$ containing one and two-body terms as in Eq. \ref{eq:sop}. We discovered in Section \ref{sec:Density} 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. \ref{eq:Hn} by
\begin{equation}\label{eq:Eibasis}
	\mathcal{E}_i = \Tr[\mathbf{D}_i\mathbf{H}].
\end{equation}
The matrix $\mathbf{H}$ has elements
\begin{equation}\label{eq:Hmn}
	H_{mn} = \int \psi_m^*(\mathbf{X})H(\mathbf{X})\psi_n(\mathbf{X})d\mathbf{X},
\end{equation}
where $H(\mathbf{X})$ is the effective two-particle Hamiltonian (Eq. \ref{eq:A2}) acting on the geminal basis functions $\psi_i(\mathbf{X})$. Instead of using the two-electron Slater determinants of Eq. \ref{eq:nibasis}, it will be advantageous to choose the basis that diagonalizes $H(\mathbf{X})$ by
\begin{equation}\label{eq:H2basis}
H(\mathbf{X})\psi_j(\mathbf{X}) = E_j\psi_j(\mathbf{X}).
\end{equation}

Suppose that we have solved for at least the $N(N-1)/2$ lowest energy eigenstates of Eq. \ref{eq:Eibasis} and computed the diagonal matrix $\mathbf{H}$ by Eq. \ref{eq:Hmn}. We may then discover the ground state by finding the $N$-representable matrix $\mathbf{D}_i$ that yields the minimum possible energy in Eq. \ref{eq:Eibasis}. Dropping subscript $i$ from the GDM, the minimum energy state that satisfies the rules of Eq. \ref{eq:ruless} is
\begin{align}\label{eq:Drule}
	D_{mn} &= \begin{cases}
		1  & \text{if } m=n \text{ and } n \leq \begin{pmatrix}N\\2\end{pmatrix} \\
		0 & \text{otherwise}.
	\end{cases}
\end{align}
Unfortunately, the matrix in Eq. \ref{eq:Drule} 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. \ref{eq:Drule} for two reasons, the first being that it is precisely the form of the ground state matrices postulated by Bopp \cite{boppAbleitungBindungsenergieNTeilchenSystemen1959} 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. \ref{eq:Drule}. 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 $\Tr[\mathbf{D}^2]$ is conserved under unitary transformation, the ground state must also satisfy $\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 O$^{5+}$ with a remarkably low error of $0.017\%$ (See tables I and II of Ref \cite{boppAbleitungBindungsenergieNTeilchenSystemen1959}). 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.

\section{The time evolution equation}\label{sec:Time}
Subsection \ref{sec:real} 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 \ref{sec:Density} to apply at arbitrary times. Most importantly, our wave function will be the position representation of the time-dependent abstract operator $\alpha(t)$
\begin{equation}\label{eq:psiabs}
	\Psi(\mathbf{x}_1,\dots,\mathbf{x}_N|t) = \braket{\mathbf{x}_1,\dots,\mathbf{x}_N|\alpha(t)}.
\end{equation}
The expectation value of a generally time-dependent linear operator $\hat{A}(t)$ for the state in Eq. \ref{eq:psiabs} is
\begin{equation}\label{eq:Awave}
	\braket{A}(t) = \braket{\alpha(t)|\hat{A}(t)|\alpha(t)}.
\end{equation}
We can also calculate observables by the equivalent $2$-RDM formulation, defining the time-dependent $2$-RDM
\begin{equation}\label{eq:rt}
	\rho(\mathbf{X},\mathbf{X}'|t) = \int \Psi^*(\mathbf{X},\mathbf{Y}|t)\Psi(\mathbf{X}',\mathbf{Y}|t)d\mathbf{Y}
\end{equation}
following Eq. \ref{eq:rho1}. The expectation value is then given by
\begin{equation}\label{eq:rhoev}
	\braket{A}(t) = \int \delta(\mathbf{X}-\mathbf{X}') A(\mathbf{X}'|t)\rho(\mathbf{X},\mathbf{X}'|t)d\mathbf{X} d\mathbf{X}'.
\end{equation}
The question that we need to answer is: How must $\rho(\mathbf{X},\mathbf{X}'|t)$ evolve in time so that the expectation values computed by Eq. \ref{eq:rhoev} match those given in Eq. \ref{eq:Awave} 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}'|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. \ref{eq:rhoev} to find
\begin{align}\label{eq:rhocomp}
	\frac{d}{dt}\braket{A}(t) &=  \Braket{\frac{dA}{dt}} \nonumber \\
	&+ \int  d\mathbf{X} d\mathbf{X}'\delta(\mathbf{X} - \mathbf{X}')A(\mathbf{X}')\dot{\rho}(\mathbf{X},\mathbf{X}'|t).
\end{align}
We must find a matching expression for Eq. \ref{eq:rt} in the hope of finding an equation that connects $\dot{\rho}(\mathbf{X},\mathbf{X}'|t)$ to known quantities. For brevity, we will group all $N$ coordinate sets together into the symbol
\begin{equation}
	\overline{\mathbf{X}} = \mathbf{x}_1,\dots,\mathbf{x}_N
\end{equation}
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
\begin{equation}
	\Psi(\overline{\mathbf{X}}|t) = U(t,t_0)\Psi(\overline{\mathbf{X}}|t_0).
\end{equation}
Expressing the time-dependent operator $A(\overline{\mathbf{X}}|t)$ as a pair sum following Eq. \ref{eq:Atot}, we find the position representation of the inner product in Eq. \ref{eq:Awave}
\begin{align}\label{eq:At0}
    \braket{A}(t) &= \sum_{i,j>i}\int \Psi^*(\overline{\mathbf{X}}|t_0)U^\dagger(t,t_0) \nonumber \\
     &\times A(\mathbf{x}_i,\mathbf{x}_j|t) U(t,t_0) \Psi(\overline{\mathbf{X}}|t_0)d\overline{\mathbf{X}}.
\end{align}
Recalling the notation of Section \ref{sec:Density}, we have defined the differential $d\overline{\mathbf{X}} = d\mathbf{x}_1\cdots d\mathbf{x}_N$ which includes spin sums by Eq. \ref{eq:spinsum}.

We proceed to differentiate Eq. \ref{eq:At0}, 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{\"o}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
\begin{align}\label{eq:Along1}
	\frac{d}{dt}\braket{A}(t) = \Braket{\frac{dA}{dt}} -iK(t),
\end{align}
where $K(t)$ is the expectation value of the many-body commutator
\begin{align}\label{eq:Along2}
	K(t) &=  \int d\overline{\mathbf{X}}\Psi^*(\overline{\mathbf{X}}|t) \nonumber \\
	&\times\sum_{\substack{i,j>i \\ k,l>k}}
	\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}}.
\end{align}
In Eq. \ref{eq:Along2} we expressed the many-body Hamiltonian $H(\overline{\mathbf{X}}|t)$ as a pair sum per Eq. \ref{eq:Atot}.

Simplify Eq. \ref{eq:Along2} 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. \ref{eq:Asep1}. 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
\begin{equation}\label{eq:Ktot}
	K(t) = K_{11}(t) + K_{12}(t) + K_{21}(t) + K_{22}(t).
\end{equation}

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 \ref{eq:Along2} reduces to $N(N-1)/2$ identical integrals. Continuing to expand $H_1(\mathbf{x}_i,\mathbf{x}_j)$ by Eq. \ref{eq:Aeff1}, we can separate $\mathbf{x}_i$ and $\mathbf{x}_j$ contributions to find
\begin{align}\label{eq:K1b}
	K_{1\beta}(t) &= \begin{pmatrix}N\\2\end{pmatrix} \int d\overline{\mathbf{X}} \Psi^*(\overline{\mathbf{X}}|t) \nonumber \\
	\times\Bigg(&\sum_{i,j>i}\left[\frac{h_1(\mathbf{x}_i|t)}{N-1},A_{\beta}(\mathbf{x}_1,\mathbf{x}_2|t)\right] \nonumber \\
	+&\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).
\end{align}

Most terms in Eq. \ref{eq:K1b} 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
\begin{align}\label{eq:Kab}
	K_{\alpha \beta}(t) &= \begin{pmatrix}N\\2\end{pmatrix} (N-1)\int  \Psi^*(\overline{\mathbf{X}}|t) \nonumber \\
	&\times\left[H_{\alpha}(\mathbf{X}|t),A_{\beta}(\mathbf{X}|t)\right]\Psi(\overline{\mathbf{X}}|t) d\overline{\mathbf{X}}
\end{align}
for $\alpha \beta = 11,12,21$. Eq. \ref{eq:Kab} does not immediately apply to $K_{22}(t)$ because, unlike in Eq. \ref{eq:K1b}, 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. \ref{eq:Kab} applies trivially to $\alpha \beta = 22$. The consequences of this requirement are discussed in the Conclusion.

Absorb the prefactor $(N-1)$ of Eq. \ref{eq:Kab} into the Hamiltonian by defining
\begin{equation}\label{eq:Hprime}
	H'(\mathbf{X}|t) = (N-1)H(\mathbf{X}|t).
\end{equation}
Summing the $K_{\alpha\beta}$ by Eq. \ref{eq:Ktot} allows us to express $K(t)$ as
\begin{align}\label{eq:Kwf}
	K(t) &=  \begin{pmatrix}N\\2\end{pmatrix} \int d\mathbf{X} d\mathbf{Y} \Big\{ \nonumber \\
	&  \Big[H'(\mathbf{X}|t)\Psi^*(\mathbf{X},\mathbf{Y}|t)\Big]
	A(\mathbf{X}|t)\Psi(\mathbf{X},\mathbf{Y}|t) \nonumber \\
	&- \Psi^*(\mathbf{X},\mathbf{Y}|t)H'(\mathbf{X}|t)
	A(\mathbf{X}|t)\Psi(\mathbf{X},\mathbf{Y}|t)\Big\}.
\end{align}
In Eq. \ref{eq:Kwf}, we separated the terms of the commutator and chose one Hermitian operator $H'(\mathbf{X}|t)$ to act on the left copy of the wave function. We can finally substitute Eq. \ref{eq:rt} into Eq. \ref{eq:Kwf} to find the derivative of the expectation value by Eq. \ref{eq:Along1}:
\begin{align}\label{eq:Afromr}
	\frac{d}{dt}\braket{A}(t) &= \Braket{\frac{dA}{dt}}  -i\int d\mathbf{X} d\mathbf{X}' \delta(\mathbf{X} - \mathbf{X}') \\
	&\times A(\mathbf{X}')[H'(\mathbf{X}) - H'(\mathbf{X}')]\rho(\mathbf{X},\mathbf{X}'|t).
\end{align}
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. \ref{eq:Afromr} to Eq. \ref{eq:rhocomp}, we find the desired expression for $\dot{\rho}(\mathbf{X},\mathbf{X}'|t)$:
\begin{equation}\label{eq:rhodot}
	\dot{\rho}(\mathbf{X},\mathbf{X}'|t) = -i[H(\mathbf{X}) - H(\mathbf{X}')]\rho(\mathbf{X},\mathbf{X}'|t).
\end{equation}

We proceed to derive a more convenient matrix representation of Eq. \ref{eq:rhodot} starting with a time-dependent geminal expansion of the $2$-RDM
\begin{equation}\label{eq:rhoegen}
	\rho(\mathbf{X},\mathbf{X}'|t) = \sum_{mn}D_{mn}(t) \psi_n^*(\mathbf{X}|t)\psi_m(\mathbf{X}'|t).
\end{equation}
Note that Eq. \ref{eq:rhoegen} 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. \ref{eq:rhodot} gives
\begin{align}\label{eq:derivative}
	i\frac{d}{dt}&\braket{\mathbf{X}'|\hat{D}|\mathbf{X}} =  \sum_{mn}D_{mn}(t) \psi_m(\mathbf{X}'|t)\left[H'(\mathbf{X}|t)\psi_n^*(\mathbf{X}|t)\right] \nonumber \\
	&- \sum_{mn}D_{mn}(t) \left[H'(\mathbf{X}'|t)\psi_m(\mathbf{X}'|t)\right]\psi_n^*(\mathbf{X}|t),
\end{align}
where we have rearranged the $\psi_i(\mathbf{X})$ into the most convenient order. Substitute into Eq. \ref{eq:derivative} the following identities
\begin{align}\label{eq:Dcom}
	\psi_m(\mathbf{X}'|t) &= \braket{\mathbf{X}'|m(t)} \\
	\psi_n^*(\mathbf{X}|t) &= \braket{n(t)|\mathbf{X}} \\
	H'(\mathbf{X}'|t)\psi_m(\mathbf{X}'|t) &= \braket{\mathbf{X}'|\hat{H}'(t)|m(t)} \\
	H'(\mathbf{X}|t)\psi_{n}^*(\mathbf{X}|t) &= \braket{n(t)|\hat{H}'(t)|\mathbf{X}},
\end{align}
so that everything on the right-hand side sits between $\bra{\mathbf{X}'}$ and $\ket{\mathbf{X}}$. Recalling the definition of $\hat{D}(t)$ (Eq. \ref{eq:Dab}), we extract from Eq. \ref{eq:derivative} the abstract operator equation
\begin{equation}\label{eq:comabs}
	\frac{d}{dt}\hat{D}(t) = -i[\hat{H}'(t),\hat{D}(t)].
\end{equation}
Eq. \ref{eq:comabs} is the familiar Liouville-Von Neumann equation.

We can specialize Eq. \ref{eq:comabs} 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
\begin{equation}\label{eq:Lvn}
    \dot{\mathbf{D}}(t) = -i[\mathbf{H}'(t),\mathbf{D}(t)],
\end{equation}
which we find by an alternate derivation in Appendix \ref{sec:Alternate}. Eq. \ref{eq:Lvn} can be shown to evolve through matrix transformation by the operator $\mathbf{U}(t,t_0)$ as
\begin{equation}\label{eq:Dopu}
    \mathbf{D}(t) = \mathbf{U}(t,t_0)\mathbf{D}(t_0)\mathbf{U}^\dagger(t,t_0).
\end{equation}
Plugging Eq. \ref{eq:Dopu} into Eq. \ref{eq:Lvn} and matching expressions on either side we find that $\mathbf{U}(t)$ must satisfy
\begin{equation}\label{eq:UUU}
	\dot{\mathbf{U}}(t,t_0) = -i \mathbf{H}(t)\mathbf{U}(t,t_0),
\end{equation}
and it is a routine procedure to successively integrate Eq. \ref{eq:UUU} to infinite order for
\begin{equation}\label{eq:U}
    \mathbf{U}(t,t_0) = \mathcal{T}_t \exp\left(-i\int_{t_0}^t d\tau \mathbf{H}'(\tau)\right).
\end{equation}
$\mathbf{U}(t)$, being the exponential of a skew-Hermitian matrix, is clearly unitary. As mentioned at the end of Section \ref{sec:Density}, we now see that the unitary transformation-invariant property
\begin{equation}\label{eq:idemp}
	\mathbf{D}^2 = \mathbf{D}
\end{equation}
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. \ref{eq:idemp} must be true for any GDM with fixed particle number.

\section{Switching on the electron-electron interaction}\label{sec:Evolution}
Having confirmed unitary evolution of the GDM, we will now exploit this property to solve the many-body Schr\"odinger 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
\begin{equation}\label{eq:Hamfull}
	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|},
\end{equation}
with $\lambda(0) = 0 $ and $\lambda(t) = 1$ for $t$ larger than some switching time $T$.

As we saw in Eq. \ref{eq:A2}, Eq. \ref{eq:Hamfull} reduces in the $2$-RDM formalism to the effective two-electron form
\begin{equation}\label{eq:Ht}
	H(\mathbf{X}|t) = H_1(\mathbf{X}) + \frac{\lambda(t)}{|\mathbf{r}_1 - \mathbf{r}_2|}.
\end{equation}

We introduce two natural basis sets for the study of Eq. \ref{eq:Ht}, the first being the non-interacting basis that  diagonalizes the Hamiltonian at $t=0$ by
\begin{equation}
	H(\mathbf{X}|0)\psi_i(\mathbf{X}) = E_i \psi_i(\mathbf{X}).
\end{equation}
Similarly, we define the interacting basis by fixing $t=T$ and solving
\begin{align}
	H(\mathbf{X}|T)\psi_i(\mathbf{X}) &= \left(H_1(\mathbf{X}) + \frac{1}{|\mathbf{r}_1 - \mathbf{r}_2|}\right)\psi_{i}(\mathbf{X}) \nonumber \\
	&= E_i \psi_i(\mathbf{X}).
\end{align}
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 \ref{sec:Basis}):
\begin{equation}\label{eq:baschange}
	\mathbf{D}_I(t) = \mathbf{U}_{I}^N\mathbf{D}_N(t)\big(\mathbf{U}_{I}^N\big)^\dagger.
\end{equation}

In the following, we observe the time evolution of a Slater determinant under the influence of the Hamiltonian in Eqs. \ref{eq:Hamfull} and \ref{eq:Ht}. Section \ref{sec:sudden} treats the case in which the Coulomb interaction is quickly switched on and Section \ref{sec:adiabatic} details a slow adiabatic change which allows us to construct fully-interacting many-body solutions.

\subsection{The sudden approximation}\label{sec:sudden}

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. \ref{eq:Dopu} and \ref{eq:U} to be
\begin{equation}\label{eq:bt}
	\mathbf{D}(t) = e^{i\mathbf{H}'(T)(t-T)}\mathbf{D}(T)e^{-i\mathbf{H}'(T)(t-T)}.
\end{equation}
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. \ref{eq:baschange}. This choice of basis simplifies the matrix equation in Eq. \ref{eq:bt} to
\begin{equation}\label{eq:Dtt}
	D_{mn}(t) = e^{iE_{nm}'(t-T)}D_{mn}(T),
\end{equation}
with $E'_{mn} = (N-1)(E_m - E_n)$ accounting for the multiplicative constant attached to $H'(\mathbf{X})$ in Eq. \ref{eq:Hprime}.

With the time-dependent GDM given by Eq. \ref{eq:Dtt}, we can calculate the observable quantity described by the operator with position representation
\begin{equation}\label{eq:density}
	\rho(\mathbf{x};\mathbf{X}) = \frac{\delta(\mathbf{x}_1 - \mathbf{x}) + \delta(\mathbf{x}_2 - \mathbf{x})}{N-1}.
\end{equation}
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
\begin{equation}
	\rho(\mathbf{x},t) = \braket{\rho(\mathbf{x})}(t),
\end{equation}
which we calculate by the trace relation
\begin{equation}\label{eq:rhotrace}
	\rho(\mathbf{x},t) = \Tr[\mathbf{D}(t)\boldsymbol{\rho}(\mathbf{x})].
\end{equation}
The matrix elements of $\boldsymbol{\rho}(\mathbf{x})$ take the simplified form
\begin{equation}\label{eq:rhot2}
	\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.
\end{equation}

Computing the trace in Eq. \ref{eq:rhotrace} finally yields the electron density
\begin{equation}\label{eq:rhot}
    \rho(\mathbf{x},t) = \sum_{mn}D_{mn}(T)e^{-iE'_{mn}(t-T)}\rho_{nm}(\mathbf{x}).
\end{equation}
We perform a simple check that particle number is conserved; integrating Eq. \ref{eq:rhot} over $\mathbf{x}$, we find from orthonormality and Eq. \ref{eq:trace0} 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.
\subsubsection{The origin of incoherent quantum fluctuations}
We can extract a surprising amount of insight from the simple relation in Eq. \ref{eq:rhot}. Separating the stationary and oscillating terms allows us to express the density as
\begin{align}\label{eq:rhoo}
    \rho(\mathbf{x},t) &= \sum_l D_{ll}\rho_{ll}(\mathbf{x}) + \sum_{\substack{m,n\\E_{m}'=E_n'}}\text{Re}\{P_{mn}(\mathbf{x})\} \nonumber \\
    &+ \sum_{\substack{p,q\\E_{p}'\neq E_q'}}|P_{pq}(\mathbf{x})|\cos\left(E_{pq}'t + \theta_{pq}(\mathbf{x})\right),
\end{align}
with
\begin{align}
	P_{ij}(\mathbf{x}) &= D_{ij}(T)\rho_{ji}(\mathbf{x}) \\
		\theta_{ij}(\mathbf{x}) &= \arg \left( P_{ij}(\mathbf{x})\right).
\end{align}
The first sum in Eq. \ref{eq:rhoo} 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
\begin{equation}
	\overline{f(\mathbf{x},t)} = \lim_{\tau \rightarrow \infty}\frac{1}{\tau}\int_{T}^{T+\tau}f(\mathbf{x},t')dt',
\end{equation}
we find the time-averaged density at $\mathbf{x}$ to be determined entirely by the diagonal and degenerate terms
\begin{equation}
	\overline{\rho(\mathbf{x},t)} = \sum_l D_{ll}\rho_{ll}(\mathbf{x}) + \sum_{\substack{m,n\\E_{m}'=E_n'}}\text{Re}\{P_{mn}(\mathbf{x})\}.
\end{equation}
On the other hand, the non-degenerate off-diagonals introduce temporal density fluctuations with variance
\begin{equation}\label{eq:sigma}
    \sigma^2(\mathbf{x}) = \overline{\rho^2(\mathbf{x},t)} - \overline{\rho(\mathbf{x},t)}^2  = \frac{1}{2}\sum_{\substack{p,q\\E_{p}'\neq E_q'}}|P_{pq}(\mathbf{x})|^2.
\end{equation}
The fluctuations increase in magnitude as the number of off-diagonal terms increases, an idea which can also be understood directly from Eq. \ref{eq:rhoo} 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. \ref{eq:rhoo}--\ref{eq:sigma} 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
\begin{align}\label{eq:Hpr}
	H'(\mathbf{X}|t) &= \sum_{i=1}^2\left(-\frac{\nabla_i^2}{2} - \sum_j\frac{1}{|\mathbf{r}_i - \mathbf{R}_j(t)|}\right) \nonumber \\
	&+ (N-1)\frac{\lambda(t)}{|\mathbf{r}_1 - \mathbf{r}_2|},
\end{align}
noting that we used $H'(\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
\begin{align}
	\dot{\mathbf{D}}(t) &= -i[\mathbf{H}'(t),\mathbf{D}(t)]\label{eq:Dcou} \\
	M_i\ddot{\mathbf{R}}_i(t) &= \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} \label{eq:Mc},
\end{align}
where Eq. \ref{eq:Mc} 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. \ref{eq:spinsum}. 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. \ref{eq:Mc} which causes a near-instantaneous scattering of the nuclei. The resulting motion of $\mathbf{R}_i(t)$ affects the Hamiltonian in Eq. \ref{eq:Hpr} and subsequently the GDM by Eq. \ref{eq:Dcou}. 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. \ref{eq:Mc} as a sum of normal modes of motion excited by the oscillatory terms of Eq. \ref{eq:rhoo}. 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. \ref{eq:Hpr}, 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.

\subsection{The degenerate adiabatic theorem}\label{sec:adiabatic}
Section \ref{sec:Time} derived the equation of motion (Eq. \ref{eq:comabs})
\begin{equation}\label{eq:em2}
	\frac{d}{dt}\hat{D}(t) = -i[\hat{H}'(t),\hat{D}(t)]
\end{equation}
for the abstract density operator $\hat{D}(t)$. We found the corresponding matrix equation (Eq. \ref{eq:Lvn}) 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
\begin{equation}\label{eq:evec}
	\hat{H}(t)\ket{i(t)} = E_i(t)\ket{i(t)},
\end{equation}
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
\begin{equation}\label{eq:adex}
	\hat{D}(t) = \sum_{pq}\mathcal{D}_{pq}(t)\ket{p(t)}\bra{q(t)}.
\end{equation}
Inserting Eq. \ref{eq:adex} into Eq. \ref{eq:em2} 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}'(t)D_{mn}(t)$ after applying Eq. \ref{eq:evec}. Repeating for the left-hand side, by the chain rule we have that
\begin{align}\label{eq:emb}
	&\Braket{m(t)|\frac{d}{dt}\hat{D}(t)|n(t)} = \dot{\mathcal{D}}_{mn}(t) \nonumber \\
	&+\sum_{p}\mathcal{D}_{pn}(t)\braket{m(t)|\dot{p}(t)} + \sum_{q}\mathcal{D}_{mq}(t)\braket{\dot{q}(t)|n(t)}.
\end{align}
Because $(d/dt)\braket{q(t)|n(t)} = 0$ by orthonormality at time $t$, we can simplify Eq. \ref{eq:emb} 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
\begin{equation}\label{eq:Mel}
	M_{ij}(t) = \Braket{i(t)|\frac{d}{dt}|j(t)}.
\end{equation}
After doing so, the final equation of motion for the $mn$ element is
\begin{align}\label{eq:Dmnexpansion}
	&\dot{\mathcal{D}}_{mn}(t) = -iE_{mn}'(t)\mathcal{D}_{mn}(t) \nonumber \\
	 &-\sum_{p}\mathcal{D}_{pn}(t)M_{mp}(t) + \sum_{q}\mathcal{D}_{mq}(t)M_{qn}(t).
\end{align}

We will now simplify Eq. \ref{eq:Dmnexpansion} 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. \ref{eq:Ht}. Closely following Ref. \cite{aguiarpintoAdiabaticApproximationDensity2002}, we will discover how $\boldsymbol{\mathcal{D}}(t)$ evolves when $T \rightarrow \infty$. Defining a natural time
\begin{equation}
	s(t) = \frac{t}{T}
\end{equation}
scaled by the switching duration, we use the fact that $(d/dt) = T(d/ds)$ to see that
\begin{align}\label{eq:qqq}
	\dot{\mathcal{D}}_{mn}(s) &= - iTE_{mn}'(s) \mathcal{D}_{pn}(s) \nonumber \\
	+&\sum_{q}\mathcal{D}_{mq}(s)M_{qn}(s) - \sum_{p}M_{mp}(s)\mathcal{D}_{pn}(s).
\end{align}
In Eq. \ref{eq:qqq} 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
\begin{equation}\label{eq:dynphase}
	\mathcal{D}_{mn}(s) =  \tilde{\mathcal{D}}_{mn}(s)e^{-iT\int_{0}^s E_{mn}'(s') ds'},
\end{equation}
so that
\begin{align}\label{eq:Dtilde}
	&\frac{d}{ds}\tilde{\mathcal{D}}_{mn}(s) = \hspace{-4pt}\sum_{\substack{q\\E_q=E_n}}\hspace{-4pt}\tilde{\mathcal{D}}_{mq}(s)M_{qn}(s) -\hspace{-4pt}\sum_{\substack{p\\E_p=E_m}}\hspace{-4pt}M_{mp}(s)\tilde{\mathcal{D}}_{pn}(s) \nonumber \\
	&+ \sum_{\substack{q\\E_q\neq E_n}}\frac{d}{ds}\int_0^s ds' \tilde{\mathcal{D}}_{mq}(s')M_{qn}(s')e^{-iT\int_0^{s'} E_{nq}'(s'')ds''} \nonumber \\
	&-\sum_{\substack{p\\E_p\neq E_m}} \frac{d}{ds}\int_0^s ds' M_{mp}(s')\tilde{\mathcal{D}}_{pn}(s')e^{-iT\int_0^{s'} E_{pm}'(s'')ds''}.
\end{align}
In Eq. \ref{eq:Dtilde} 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'$.

We will now evaluate the last line of Eq. \ref{eq:Dtilde} to understand the implications of the dynamical phase. To begin, we simplify by defining
\begin{equation}\label{eq:Fp}
	F_p(s') = M_{mp}(s')\mathcal{D}_{pn}(s')
\end{equation}
and
\begin{equation}\label{eq:gmp}
	g_{pm}(s') = \int_0^{s'} E_{pm}'(s'') ds''
\end{equation}
so that the integral over $ds'$ (which we call $I_p(s)$) takes the form
\begin{equation}
	I_{p}(s) = \int_{0}^s  F_p(s')e^{-iT g_{pm}(s')}ds'.
\end{equation}
Using that $\dot{g}_{pm}(s') = E'_{pm}(s')$ by Eq. \ref{eq:gmp}, we multiply the integrand by the identity $1 = iT\dot{g}_{pm}(s')/(iTE'_{pm}(s'))$ so that
\begin{equation}\label{eq:Ipa}
	I_{p}(s) = \int_{0}^s  \frac{F_p(s')}{iTE_{pm}'(s')}\left(iT\dot{g}_{pm}(s')e^{iT g_{pm}(s')}\right)ds'.
\end{equation}
Because $iT \dot{g}_{pm}(s')e^{iT g_{pm}(s')} = (d/ds')e^{iT g_{pm}(s')}$, we are able to integrate Eq. \ref{eq:Ipa} parts to finally find
\begin{align}\label{eq:Ip}
	I_p(s) &= \frac{1}{iT}\left[\frac{F_{p}(s')}{E_{pm}'(s')}e^{iTg_{pm}(s')}\right]_{0}^s \nonumber \\
	&- \frac{1}{iT}\int_0^s \frac{d}{ds'}\left(\frac{F_{p}(s')}{E'_{pm}(s')}\right)e^{iTg_{pm}(s')}ds'.
\end{align}
As long as $F_p(s')$ 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. \ref{eq:Dtilde} to
\begin{equation}\label{eq:adif}
	\frac{d}{dt}\tilde{\mathcal{D}}_{mn}(t) = \hspace{-5pt}\sum_{\substack{q\\E_q = E_n}}\tilde{\mathcal{D}}_{mq}(t)M_{qn}(t)
	- \hspace{-5pt}\sum_{\substack{p\\E_p = E_m}}M_{mp}(t)\tilde{\mathcal{D}}_{pn}(t),
\end{equation}
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. \ref{eq:adif} 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. \ref{eq:dynphase} 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 $\bs{\mathcal{D}}_\mu(t)$ for each degenerate subspace $\mu$. Noting that the terms in Eq. \ref{eq:adif} are in the form of matrix multiplications, we find that the equation of motion separates into the commutation relations
\begin{equation}\label{eq:Ddot}
	\dot{\boldsymbol{\mathcal{D}}}_\mu(t) = -[\mathbf{M}_\mu(t),\boldsymbol{\mathcal{D}}_\mu(t)]
\end{equation}
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
\begin{equation}\label{eq:Dsum}
	\boldsymbol{\mathcal{D}}(t) = \text{diag}[\boldsymbol{\mathcal{D}}_1(t),\boldsymbol{\mathcal{D}}_2(t),\dots,\boldsymbol{\mathcal{D}}_\mu (t),\dots],
\end{equation}
valid for all times including $t>T$. As a final remark, note that Eq. \ref{eq:Ddot} 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. \ref{eq:Lvn}.

\subsubsection{Practical considerations of the degenerate adiabatic theorem}

The application of Eqs. \ref{eq:Ddot} and \ref{eq:Dsum} 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. \ref{eq:Mel} in the position representation
\begin{equation}\label{eq:connection}
	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}.
\end{equation}
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. \ref{eq:Ip} is no longer valid. Thus, we must require that $\psi_j(\mathbf{X}|\lambda)$ is smooth as a function of $\lambda$. In Section \ref{sec:Solve} 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. \ref{eq:Ddot} and \ref{eq:Dsum}, 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:
\begin{equation}\label{eq:finite}
	\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}}.
\end{equation}
It follows that $\mathcal{E}(t)$ is differentiable because the right hand side of Eq. \ref{eq:finite} is finite at all times. The GDM that represents the wave function computes the energy as
\begin{equation}\label{eq:Etime}
	\mathcal{E}(t) = \sum_{\mu}\Tr[\boldsymbol{\mathcal{D}}_{\mu}(t)]E_\mu(\lambda(t)),
\end{equation}
where $\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. \ref{eq:finite} that $\mathcal{E}(t)$ must be differentiable, we conclude that no instantaneous population interchange can occur at a level crossing.

\section{Computing many-body eigenstates using the adiabatic theorem}\label{sec:Solve}
We are now equipped to solve the many-electron Schr\"odinger equation
\begin{equation}\label{eq:Hfin}
	\left(\sum_{ij}H(\mathbf{x}_i,\mathbf{x}_j)\right)\Psi_n(\overline{\mathbf{X}}) = \mathcal{E}_n \Psi_n(\overline{\mathbf{X}})
\end{equation}
for the matrices $\mathbf{D}_n$ that represent wave functions $\Psi_n(\overline{\mathbf{X}})$. The Hamiltonian of interest has terms
\begin{align}\label{eq:Hijf}
	H(\mathbf{x}_i,\mathbf{x}_j) =  H_1(\mathbf{x}_i,\mathbf{x}_j) + \frac{1}{|\mathbf{r}_i - \mathbf{r}_j|},
\end{align}
containing the one-body contribution
\begin{equation}
	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]
\end{equation}
along with the Coulomb potential $1/|\mathbf{r}_i - \mathbf{r}_j|$. For reasons that will be explained shortly, we treat Eqs. \ref{eq:Hfin} and \ref{eq:Hijf} as a special case of the two-parameter Hamiltonian
\begin{align}\label{eq:Vpert}
	H(\epsilon,\lambda) &= \sum_{i,j>i}\Bigg[ H_1(\mathbf{x}_i,\mathbf{x}_j) \nonumber \\
	&+ \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]
\end{align}
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
\begin{equation}\label{eq:H2p}
	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)
\end{equation}
is nondegenerate when $\epsilon \neq 0$.

We represent a given eigenstate $\Psi_n(\overline{\mathbf{X}}|\epsilon,\lambda)$ of Eq. \ref{eq:H2p} by its GDM $\boldsymbol{\mathcal{D}}(\epsilon,\lambda)$, dropping the subscript $n$. The electronic energy is then calculated by
\begin{equation}\label{eq:Esumup}
	\mathcal{E}_n(\epsilon,\lambda) = \sum_{i}\mathcal{D}_{ii}(\epsilon,\lambda)E_i(\epsilon,\lambda),
\end{equation}
where each $E_i(\epsilon,\lambda)$ is an eigenvalue of the effective two-particle Hamiltonian $H(\mathbf{X}|\epsilon,\lambda)$ satisfying the Schr\"odinger equation
\begin{equation}\label{eq:g2e}
	H(\mathbf{X}|\epsilon,\lambda)\psi_i(\mathbf{X}|\epsilon,\lambda) = E_i(\epsilon,\lambda)\psi_i(\mathbf{X}|\epsilon,\lambda)
\end{equation}
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. \ref{eq:Ddot} 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. \ref{eq:g2e} for each $\lambda(t)$ with $\epsilon = 1$. Fig. \ref{fig:adi} 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.

\begin{figure}[h]
	\includegraphics{adiabatic_paper.pdf}
	\caption{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.}
	\label{fig:adi}
\end{figure}

For the simplest non-trivial example, consider a 3-electron state constructed from the geminal eigenstates of Fig. \ref{fig:adi}. 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. \ref{eq:Esumup} should be differentiable at all times by Eq. \ref{eq:finite}. 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 $\conf{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. \ref{fig:adi}. 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
\begin{equation}\label{eq:Efind}
	\mathcal{E} = \sum_\mu \Tr[\boldsymbol{\mathcal{D}}_\mu(0)]E_\mu(0,1),
\end{equation}
where $\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. \ref{eq:Efind}.

For observables other than the energy, the best we can do is calculate the $\epsilon$-dependent expectation value by
\begin{equation}
	\braket{A}(\epsilon) = \Tr[\boldsymbol{\mathcal{D}}(T_1) \mathbf{A}(\epsilon)].
\end{equation}
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.

\subsection{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\"odinger 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. \ref{eq:Ht} with $\lambda = 1$:
\begin{align}\label{eq:Hzn}
	H_{Z,N}(\mathbf{r}_1,\mathbf{r}_2|1) &= \frac{1}{N-1}\Bigg[-\frac{\nabla_1^2}{2} - \frac{\nabla_2^2}{2} \nonumber \\
	&- Z\left(\frac{1}{|\mathbf{r}_1|} + \frac{1}{|\mathbf{r}_2|}\right) + \frac{N-1}{|\mathbf{r}_1 - \mathbf{r}_2|}\Bigg].
\end{align}
It will be convenient to reduce the bracketed expression into a form that resembles the Helium atom Hamiltonian
\begin{equation}\label{eq:H22}
	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|}.
\end{equation}
for some choice of $\lambda$. We find this simplification by transforming the coordinates into a yet-undetermined natural scale
\begin{equation}
	\overline{\mathbf{r}} = a \mathbf{r}.
\end{equation}
This transformation changes Eq. \ref{eq:Hzn} to
\begin{align}\label{eq:Hznt}
	H_{Z,N}(\overline{\mathbf{r}}_1,\overline{\mathbf{r}}_2|1) &= \frac{a^2}{N-1}\Bigg[-\frac{\overline{\nabla}_1^2}{2} - \frac{\overline{\nabla}_2^2}{2} \nonumber \\
	&-\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].
\end{align}
Eq. \ref{eq:Hznt} reduces to a scaled copy of the Helium Hamiltonian in Eq. \ref{eq:H22} if we choose
\begin{equation}
	a = \frac{Z}{2},
\end{equation}
from which it follows that $\lambda = 2(N-1)/Z$. In the end we relate the $Z,N$ and $2,2$ (Helium) Hamiltonians by
\begin{equation}\label{eq:Hzn22}
	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).
\end{equation}
Eq. \ref{eq:Hzn22} 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
\begin{equation}\label{eq:psizn}
	\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).
\end{equation}
Eqs. \ref{eq:Hzn22} and \ref{eq:psizn} can be used in conjunction with the two-parameter solution method outlined in Section \ref{sec:Solve}. 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$.

\section{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.

\acknowledgments
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.

redpajama_set_name: RedPajamaArXiv


text:
<?xml version="1.0" encoding="utf-8"?>
<metadata>
	<layout title="COM_REGOFFICES_COUNTRY">
		<message>COM_REGOFFICES_COUNTRY_DESC</message>
	</layout>
	<!-- Add fields to the request variables for the layout. -->
	<fields name="request">
		<fieldset name="request"
			addfieldpath="/administrator/components/com_content/models/fields">

			<field name="id" type="modal_article"
				label="COM_REGOFFICES_FIELD_SELECT_COUTNRY_LABEL"
				required="true"
				edit="true"
				clear="false"
				description="COM_REGOFFICES_FIELD_SELECT_COUTNRY_DESC"
			/>
		</fieldset>
	</fields>

</metadata>

redpajama_set_name: RedPajamaGithub


text:
package com.datazuul.webapps.cmslight.datasource;

import org.apache.tapestry.IRequestCycle;

import com.datazuul.webapps.cmslight.beans.IUser;
import com.datazuul.webapps.cmslight.beans.IWebsite;

/**
 * Interface for a service which provides {@link IUser} instances, and each
 * users {@link IWebsite} instances.
 *
 * @author Ralf Eichinger
 */
public interface IUserSource {
	/**
	 * Name ("pixocmslight.user-source") of an application
	 * extension that implements {@link IUserSource}.
	 */
	public static final String USER_SOURCE_EXTENSION_NAME = "pixocmslight.user-source";

	/**
	 * Invoked to ensure the source is initialized. This must always be invoked
	 * before any other methods.
	 */
	public void initialize(IRequestCycle cycle);

	/**
	 * Create and return a new {@link IUser}defined in this user database.
	 *
	 * @param username
	 *            Username of the new user
	 *
	 * @exception IllegalArgumentExceptionif
	 *                the specified username is not unique
	 */
	public IUser createUser(String username);

	/**
	 * Finalize access to the underlying persistence layer.
	 *
	 * @exception Exception
	 *                if a database access error occurs
	 */
	public void close() throws Exception;

	/**
	 * Return the existing {@link IUser}with the specified username, if any;
	 * otherwise return null.
	 *
	 * @param username
	 *            Username of the user to retrieve
	 */
	public IUser findUser(String username);

	/**
	 * Return the set of {@link IUser}s defined in this user database.
	 */
	public IUser[] findUsers();

	/**
	 * Initiate access to the underlying persistence layer.
	 *
	 * @exception Exception
	 *                if a database access error occurs
	 */
	public void open() throws Exception;

	/**
	 * Remove the specified {@link IUser}from this database.
	 *
	 * @param user
	 *            IUser to be removed
	 *
	 * @exception IllegalArgumentException
	 *                if the specified user is not associated with this database
	 */
	public void removeUser(IUser user);

	/**
	 * Save any pending changes to the underlying persistence layer.
	 *
	 * @exception Exception
	 *                if a database access error occurs
	 */
	public void save() throws Exception;
}


text:
Vitsa Kassotou (née Schtakleff) (b Jerusalem 14 Mar 1905 - d Nicosia 5 Jul 1992)
29 Dec 2021 / Marina Parisinou / 8 Comments
Before succumbing to the juggernaut of consumerism, the Christmas season was filled with the smells of baked goods wafting from home ovens. None evokes Christmas more vividly in my mind than the spicy, rich aroma of Yiayia Vitsa's German Biscuits. (Biscuits in Brit speak, cookies in American lingo.)
Yiayiá in Greek means grandmother. Yiayia Vitsa, my maternal grandmother, was an amazing baker and cook who had made a name for herself in her adopted hometown of Nicosia, Cyprus. Her business activities were a means of supplementing the family income and eventually hers alone, as well as an outlet for her creativity. I'm quite certain, however, that Yiayia would never have imagined that this is what life had in store for her. In Jerusalem both her father's family and the one she built afterwards with her husband were of comfortable means. How to conceive that they'd end up destitute in a foreign land?
Paraskevi (Vitsa) Schtakleff
Paraskevi Schtakleff was born in Jerusalem on 14 March 1905. Her father, John (Hanna) Schtakleff, was an ethnic Bulgarian/Macedonian whose father had emigrated to Palestine from the town of Tetovo in what today is the Republic of North Macedonia (not to be confused with the Tetovo of Bulgaria). Paraskevi's mother, Eugenie Agathopoulou (Ευγενία Αγαθοπούλου), was a Greek from Istanbul (or Smyrna?). Eugenie's father operated an itinerant theatre troupe in which he and his two daughters were the principal actors. Their travels took them to Jerusalem where great-grandfather John, upon seeing great-grandmother Eugenie perform, fell for her and hastened to ask her father for her hand which he was granted forthwith.
Eugenie bore him six children, the first being Paraskevi. The nickname Vitsa derives from the diminutive of her name: Paraskevitsa. She wasn't fond of that nickname which gave her lots of trouble in childhood because in Greek it's also a word meaning stick or whip. She impressed upon her three daughters that if they ever had daughters of their own, they could give them any name they wished except for hers. (As her only female grandchild, I'm particularly grateful for that act of humility and mercy!)
Eugenie with her five children: Vitsa on the far left. Next to her, Cocas (standing), Nandos (seated).
On the right, Marika (standing), Colias (seated) – Jerusalem, ca 1918
At the time of her birth, Jerusalem was part of the Ottoman Empire. The family was quite well off: Great-grandfather John ran a flour mill and bakery which supplied the Ottoman army with bread. Vitsa was raised in a comfortable home and went to school at the Notre Dame de Sion convent, receiving her education mostly in French. They spoke Greek at home and Arabic in the street.
Before reaching her teens, the Ottoman Empire collapsed and the fortunes of the family gradually followed suit as the British, the new overlords of Palestine, did not bring their business to John, possibly out of mistrust resulting from his close connection to the Ottomans.
In the meantime Vitsa had been taking piano lessons and had become quite accomplished. They had a piano at home where she was frequently called upon to entertain visitors. British visitors called her Freda. (Paraskevi is the Greek name for Friday.) Her piano teacher, a man by the name of Vamvoudakis, had wanted to marry her but Vitsa gave her heart instead to a Greek man who had emigrated from the island of Samos, probably around the time of her birth.
Paraskevi (Vitsa) Schtakleff, aged 16 – Jerusalem, 1921
Emmanuel (Manolis) Kassotis was a nephew of the Greek Orthodox Patriarch of Jerusalem Damianos I who took him under his wing and brought him to Jerusalem when Manolis was about 12 years old. (More about my grandfather in a later post.) He subsequently took over the running of a vast estate that belonged to the Patriarchate, located some 30 km east of Jerusalem at the village of Al-Bureij (Breij). In Greek they called it Brets (Μπρετς).
Vitsa and Manolis's marriage in July 1922 put an end to my grandmother's modern city life. My grandfather was a stern and pious man and although he loved my grandmother, he was also quite jealous. It's not difficult to see how a man of the church would have a hard time reconciling with the lifestyle of a cosmopolitan young woman a dozen years his junior. So Vitsa's sleeves and hems were lowered. Her new life in Breij was a far cry from the hustle and bustle of Jerusalem and she spent several lonely years there. Eventually my grandfather bought her a piano which was a great comfort to her, as was the arrival of the first of three daughters in 1924.
Five years later: Vitsa as a wife and mother
with eldest daughter, Vasso (aged 2) and Manolis – Palestine, 18 Aug 1926
I have been told that when my grandmother married, she didn't know how to roast a chicken and my grandfather had to show her. Hardly surprising given her dedication to the piano. Was it in Breij that she learnt to make her delicious hummus and babaganoush, her tabbouleh and Palestinian lahmajoun (which looked very different from the Armenian version known in Cyprus)? Perhaps she was taught by the fellahahs who worked on the estate. I suspect keeping busy in the kitchen, in addition to being a requirement of her new role as married woman and mother, also provided some relief from loneliness. And perhaps having been initiated to the culinary arts, she expanded her repertoire to reflect her cultural background. Her scrumptious Soutzoukakia Smyrneika (Smyrna meatballs) (yes, I was a carnivore back in the day) must have come from her mother's side of the family whereas the Tlatchen Piper (a Bulgarian pepper dip – which a cousin has recognised as a Macedonian dip by a different name) she was surely taught by her father's Slav relatives.
When the first daughter reached schooling age, the family moved to Jerusalem where two more girls followed, the middle one being my mother. Breij remained my grandfather's workplace and the family's vacation getaway.
Life in Jerusalem was full. The girls went to school; Vitsa was a busy homemaker and also on the board of the Greek community's charitable organisation. Like the rest of the large Greek community of Jerusalem, the Kassotis participated in Greek religious and national celebrations, and frequented the Greek Club in the Greek Colony where my grandmother's hands were often found hopping around on the piano keyboard. On Sundays their lunch table frequently hosted guests and featured a roast chicken which Aunt Feely – Mum's cousin and close friend – still drools thinking about.
"… Mrs Cassotou at the piano…"
The Palestine Post – 31 Mar 1939
With the partition of Palestine in May 1948 the family became refugees although at first they didn't know it. They fled to Cyprus hoping to wait out the war so they could go back to their home and life in Jerusalem. But they were never allowed to return and found themselves destitute.
Baking for Family and Business
The main evidence I have for the beginning of Yiayia's baking activities is based on her order books, the first of which is dated 1961. It's a decrepit note book filled with a mishmash of recipes, cooking notes, cost calculations and orders, in a mishmash of Greek and English. In stark contrast with later years' books in which each month's orders filled a whole page and often spanned two or more, in 1961 all of her January-November orders took just one page, including two cakes for my October christening! A second page listed the Christmas cakes.
Jan-Nov 1961 orders
And then there were the wedding cakes for which she developed quite a reputation in Nicosia. She baked her first for her youngest daughter's wedding: When Aunt Mary married journalist Alex Efthyvoulou in May 1961, the five-storey fruitcake covered in white hard icing was made by her mother. Two more wedding cakes followed that same year, both five-storeyed, one weighing 100 lb (45+ kg), the second 140 lb (60.5 kg). These were the first "commercial" ones, priced at CY£50 and CY£62 respectively – not a trivial amount at the time.
By 1977 she had made 51 wedding cakes, ie around three a year. It was a major production which she used to stress over. Not only did it require intensive preparation and work, but she would also have to ensure that the whole construction could be safely transported. In some instances she would go to the venue herself to oversee the assembly of the floors. Some of the weddings (and cakes) made the newspapers. By late 1970s she started slowing down and only made a few more, the last one being in 1984 for Evi, eldest daughter of Ellie Savvides (née Louisidis), another Jerusalemite friend of the family.
Wedding Cakes by Vitsa Kassotou
On the right, the wedding of Minister of Interior & Defence Polykarpos Giorkadjis
and Photini Michaelidou – Jun 1967
Then there were the birthday cakes. So many birthday cakes! Dozens of Nicosia kids got a thrill on their birthdays from what my brother and I, as well as our cousins, took for granted every year: a custom-made cake. The process would begin with a visit to Yiayia's flat where we'd thumb through the pages, sticky from all the dough Yiayia's fingers had deposited on them, of dozens of magazines, many of them sent to her from the US by her brother Nando. Disney characters, fairy tales, trains, swimming pools, speed boats, rockets, and so many other designs would be transformed by Yiayia's hands into cakes that looked amazing and tasted equally good.
The Efthyvoulos family and Yiayia on cousin Dimitris's 13th birthday
Skateboard cake!
Cake decorating techniques have evolved a great deal since my grandmother's time and some of her cakes may seem outdated. But Yiayia was a pioneer for her time and a master of sugar, frosting and the piping bag.
Her pièce de resistance was her swan. With the graceful bend of its neck and its angel-like wings, it was placed on a mirror which Yiayia decorated with sugar lilies. I had that for my sixteenth birthday but the swan's beauty proved its undoing. We didn't have the heart to cut it straight away, we thought we'd look at it for a few more days. And then for a few days more. And then one day we realised that the cake had gone bad.
Swan on lake with lilies
The largest birthday cake Yiayia baked must have been the one ordered by the US Embassy for America's bicentennial, which, if I remember correctly, had the map of America drawn on it.
Birthday Cakes Slideshow:
By far the most memorable cake was the one ordered by the wife of longtime mayor of Nicosia Lellos Demetriades. Lellos (as we all referred to him affectionately) had worked in partnership with Mustafa Akıncı, mayor of occupied north Nicosia, to build a unified sewerage system for the city. It was the beginning of what came to be known as the Nicosia Master Plan, a bicommunal plan for developing the two sides of the divided capital in a manner that would allow it to operate as one, should it ever be reunified. What Mrs Demetriades and Yiayia designed, based on a cartoon that appeared in the press, can safely be characterised as unique: the two mayors sitting across each other on… toilets, while at the edge of the cake were placed some… ahem… number-twos! (Not real ones, made of sugar of course!)
"Long live the sewerage system!"
Lellos on the right, Akıncı on the left.
One of the many images that come to mind when I think of Yiayia is of her sitting in our living room during one of her regular visits, in her favourite armchair next to the fireplace, with her Slavic colours – fair skin, rosy cheeks, blue eyes – her long white hair braided like a halo around her head, and her long arthritic fingers tapping nervously as she pondered a structural issue for an upcoming creation. It could have been in one of those afternoons that she figured out how to turn an ice cream cone into a swan's beak.
Yiayia by the fireplace – Nicosia, 1976
Then there were the everyday cakes – "Geography" (aka marble cake) and the like – and food for family, friends and many others. The names of the three Kassotou sisters appear in practically every month's order list. In other words: we were the regular beneficiaries of Yiayia's cooking and baking.
There were of course the expected seasonal ups and downs but, come Christmas, her kitchen and her hands went into overdrive. What came out of Yiayia's oven was a reflection of who she was: a Greek-Palestinian whose adult years in her native land were lived under a British Mandate. Her Christmas repertoire included Greek melomakarona, Palestinian ma'amoul and English Christmas cakes. And, to break with the theme of origin and tradition, she also baked German Biscuits.
Christmas cakes – early 1960s
(Photos by Jules Parisinos)
German Biscuits
We don't know how the German Biscuits sneaked into Yiayia's cookery book. Most likely she came across the recipe in one of her many baking magazines. I also don't know if she started making them in Jerusalem or Nicosia, although I suspect the latter. Regardless, Christmas in our family was not Christmas without them.
I missed them sorely after her death in 1992. A few years later I came to the conclusion that Christmas without Yiayia's German Biscuits was no longer sustainable. So I baked a batch. I had made them once before, when I still lived in Cyprus and Yiayia was still very much alive to give me some pointers. I kept some and sent the rest home. The response I got from the family was nothing short of elation: it was as if I had given them Christmas back.
And so the tradition was passed on to me, Yiayia's only granddaughter. The next year I baked more, the year after that even more. And just like my grandmother, I distributed them to the family. I would also give them to friends both here in the US and in Cyprus in lieu of purchased gifts and they loved them and came to expect them every year. And they were not shy to voice their displeasure if I missed a year. The packages I dispatched – express, of course – from San Francisco to Cyprus, London and the East Coast ran to at least a few hundred dollars at the post office.
My baking crop – Oakland, CA, Dec 2006
For several years, Pat, a good friend from work, was my most able assistant and egg-washer par excellence, and as such became part of the tradition. We baked together for several years until she moved to Atlanta – and even one year after that, when she came back to San Francisco especially. Other friends also became gracious volunteers at times. The last time we baked with Pat, the dough was over 20 kg (about 46 lb) producing 830 biscuits.
Oakland, CA – 2001
San Francisco, CA – 2009
Pat and I had amazing synergy. We worked – and chatted – non-stop for hours!
As I would recover from this exhausting – for me – feat, I would reflect in amazement on Yiayia's own output. Not only did she produce more biscuits than I, but they were only one amongst a whole array of other Christmas goodies, in equally impressive quantities. One little old lady who walked with a cane, in her home kitchen! And there was me, fit in my 30s and 40s, utterly knocked out by a mere 20 kg of dough!
Inevitably the recipe has evolved in my hands over the years. Shopping for the ingredients that first time presented some minor challenges. I remember asking at my local supermarket in Berkeley, CA, for 'mixed peel' and drawing blank looks. Luckily I spotted a package myself and learnt that in the US it's referred to as 'candied fruit'. One glance at the ingredients list, full of preservatives and colourings, is all it took for me to come to terms with the fact that the recipe would have to be adjusted. Our knowledge of nutrition has come a long way since Yiayia's days and the biscuits would have to keep up with the times.
So I substituted dried fruit (organic, of course) for the mixed peel. And out went the Spry shortening (with its trans fat), in came the butter (having been restored to the world's good graces). The flour I used was mixed 50%-50% white and wholemeal from the beginning. Then I discovered Whole White Wheat flour – flour from a white berry with its wholemeal texture being very close to the refined flour from the more customary red wheatberries – and thus the biscuits became 100% wholemeal. In later years I even made a couple of batches with canola oil rather than butter when I got word that a cousin in London was on a vegan diet for medical reasons. What remains unchanged, however, is the sliver of blanched almond placed in the centre of each biscuit, like a heart. Yiayia wouldn't have it any other way!
Baking – San Francisco, Dec 2013
German friends have identified the biscuits as (at least akin to) the traditional Lebkuchen found in German Christmas markets. Most everyone thinks they are ginger biscuits although there is no ginger in them. There's cloves and cinnamon, citrus zest, and almonds and honey… (No, sorry, the recipe is not available to the public! The three Kassotou sisters have stipulated that it remain strictly in the family, at least while they're still on this earth. Two of them are sadly gone but thankfully the youngest – and strictest about this rule! – is still with us so I'm keeping this recipe card close to my chest for now and hopefully for a while still.)
I have not abandoned the tradition – at least not in theory. In practice I haven't baked for several years now, what with one thing and another. Even though I contemplated the possibility, it wasn't meant to be this year either. Next year, inshallah!
But the holidays are here and I've missed the German Biscuits and I'm sorry I haven't been able to share them with family and friends so we can remember my beloved grandmother. What I can do though is honour Yiayia Vitsa and her baking through this post — with much love.
PS – When the pandemic started with its lockdowns, I began going through the material I have on Yiayia Vitsa – her order books, diaries, etc – thinking I'd write a book about her and her kitchen. But then the Semiramis book took precedence. In drafting this post, I realised how much more there is to share about Yiayia so a book about her has to happen at some point.
Greek ClubJerusalemNicosiaSchtakleffsYiayia Vitsa
29 Dec 2021 at 12.51pm
wonderful! I remember baking w/you one year! we had a great result!
Thanks for the reminder and memories! Stay well! hope to see you soon!
Marina Parisinou (Post author)
29 Dec 2021 at 1.49pm
Thank you, Chère. Yes, you were one my gracious volunteers!
Corinna & Feely Antipas
Dearest cousin Marina- I hung onto your every word . This is truly priceless as I am mesmerized by Thea Vitsa's talent and her story. Yes, please do write the book- your post has me wanting more!
Ever so grateful for you sharing these treasures! Thank you . And next time you bake the biscuits, please send a care package to Montréal : ))))
Many thanks for reading and your sweet words, Corinna! I'm glad Yiayia's story touched you. With a little bit of luck, there'll be a package coming your way next holiday season! 😉
Vicky Aloneftis
30 Dec 2021 at 9.26am
What a great memory of your grandmother. Baking is a unique talent, nothing to do with cooking.
30 Dec 2021 at 10.26am
Thank you, Vicky!
Munir e Zacharia
10 Jan 2022 at 1.33am
Thanks for sharing your memories, My mother also used to bake German cakes and cookies, I imagine she learned that at Schnellar.
10 Jan 2022 at 1.19pm
Thank you for reading! Another person mentioned the German cookies from Schneller so I'm beginning to suspect that that's where my grandmother discovered them.

redpajama_set_name: RedPajamaCommonCrawl