# Dynamical quasi-condensation in the weakly interacting Fermi-Hubbard model

Iva Březinová<EMAIL_ADDRESS>Institute for Theoretical Physics,
Vienna University of Technology, Wiedner Hauptstraße 8-10/136, 1040 Vienna,
Austria, EU  Markus Stimpfle  Institute for Theoretical Physics, Vienna
University of Technology, Wiedner Hauptstraße 8-10/136, 1040 Vienna, Austria,
EU  Stefan Donsa  Institute for Theoretical Physics, Vienna University of
Technology, Wiedner Hauptstraße 8-10/136, 1040 Vienna, Austria, EU  Angel
Rubio  Max Planck Institute for the Structure and Dynamics of Matter, Hamburg,
Germany, EU  Center for Computational Quantum Physics (CCQ), Flatiron
Institute, New York, NY, USA

(August 31, 2024)

###### Abstract

We study dynamical (quasi)-condensation in the Fermi-Hubbard model starting
from a completely uncorrelated initial state of adjacent doubly occupied
sites. We show that upon expansion of the system in one dimension, dynamical
(quasi)-condensation occurs not only for large interactions via the
condensation of doublons, but also for small interactions. The behavior of the
system is distinctly different in the two parameter regimes, underlining a
different mechanism at work. We address the question whether the dynamical
(quasi-)condensation effect persists in the thermodynamic limit. For this
purpose, we use the two-particle reduced density matrix method, which allows
the extension to large system sizes, long propagation times, and two-
dimensional (2D) systems. Our results indicate that the effect vanishes in the
thermodynamic limit. However, especially in 2D, further investigation beyond
numerically tractable system sizes calls for the use of quantum simulators,
for which we show that the described effect can be investigated by probing
density fluctuations.

##  I Introduction

Superconductivity is one of the most intriguing collective electronic
phenomena in solid state physics. The search for an explanation of high-
temperature superconductivity in correlated materials still drives a
considerable amount of experimental and theoretical work. Traditionally,
superconductivity has been viewed as an effect of a system at (or near)
equilibrium [1], but pioneering experiments (see e.g. [2, 3, 4, 5, 6, 7, 8, 9,
10]) have shown that states with signatures of superconductivity can be
induced by external driving. This has stimulated new interest in dynamical
condensation effects in fermionic systems (see e.g. [11, 12, 13, 14, 15, 16,
17]).
The workhorse for theoretical investigations of correlated phenomena in solid
state physics is the Fermi-Hubbard model [18, 19, 20, 21, 22]. It has been
shown analytically by Yang [23] that the Fermi-Hubbard model in arbitrary
dimensions has a special symmetry, the so-called $\eta$-symmetry, which gives
rise to the appearance of a fermionic condensate in the excitation spectrum of
the Fermi-Hubbard model. This so-called $\eta$-condensate saturates the
maximally allowed occupation number for pair states, whose value was also
derived by Yang [24], and shows off-diagonal long range order [23]. Several
recent theoretical works have been devoted to showing that ground states of
the Fermi-Hubbard model can be driven to metastable non-equilibrium states
with large overlap with the $\eta$-condensate inheriting its properties (see
e.g. [12, 13]).
Concomitantly, it was found in numerical studies [25, 26] and subsequently
measured experimentally [27] that initially uncorrelated states of hard-core
bosons form a quasi-condensate in one-dimension (1D) upon free expansion. This
effect has been explained through the physics of emergent Hamiltonians where
the dynamical system inherits the properties of a ground state condensate of
the emergent non-equilibrium Hamiltonian [28]. This study is highly relevant
also for fermions since in the case of large interactions the fermions form
pairs (doublons), which to a high degree of accuracy can be described as hard-
core bosons. This mapping has been used to study dynamical quasi-condensation
of fermions in 1D at large interactions in e.g. [29, 14].
In this work we investigate dynamical (quasi)-condensation upon a quench with
similar initial conditions, but in the opposite limit of weak electron
interactions $U\ll 1J$, where $J$ is the hopping matrix element. We show that
this dynamical (quasi-)condensation effect exhibits a markedly different
behavior than its counterpart at large $U$. We exploit our newly developed
time-dependent two-particle reduced density matrix (TD2RDM) method [30, 31,
32] to reach the required long expansion times for systems sizes of several
tens of sites, while laying the ground based on small systems and exact
calculations. We compare the results within the TD2RDM method with exact
results for small systems to extrapolate its accuracy for larger systems where
exact results are not available. Furthermore, the TD2RDM method allows to
extend the investigations to two-dimensional (2D) systems enabling us to
address the question whether the Mermin-Wagner-Hohenberg theorem [33, 34],
which prohibits any condensation in 1D in equilibrium in the thermodynamic
limit, is also valid in this dynamical situation.
The paper is structured as follows: We introduce the system under
investigation in Sec. II and review the essential building blocks of the
TD2RDM method in Sec. III emphasizing necessary extensions to incorporate the
$\eta$-symmetry. We investigate dynamical (quasi-)condensation in 1D in Sec.
IV. We first discuss small systems that are treatable exactly. These small
systems show already signatures of the effect and serve as benchmarks for the
TD2RDM method. Extensions to substantially larger systems are performed within
the TD2RDM method. In Sec. IV.3 we extend our research to 2D and conclude in
Sec. V. As units we use $\hbar=m=e=1$ unless otherwise stated.

##  II Out-of-equilibrium Fermi-Hubbard model

The system under investigation is the Fermi-Hubbard model in 1D and 2D (see
Fig. 1) given by

$\hat{H}=-J\sum_{\langle
i,j\rangle}\sum_{\sigma}\hat{a}_{i\sigma}^{\dagger}\hat{a}_{j\sigma}+U\sum_{i}\hat{n}_{i}^{\uparrow}\hat{n}_{i}^{\downarrow},$
(1)

where $\langle i,j\rangle$ denotes nearest-neighbor hopping on a 1D or 2D
lattice, and
$\hat{n}_{i}^{\uparrow(\downarrow)}=\hat{a}_{i\uparrow(\downarrow)}^{\dagger}\hat{a}_{i\uparrow(\downarrow)}$.

Figure 1: Fermi-Hubbard model in (a) 1D, and narrow strips in 2D with (b) two
or (c) three transversal sites. The total number of sites is $M_{s}$. The
orange bar represents the investigated initial condition of a block of
adjacent doubly occupied sites either placed in the center (as shown in the
images) or left aligned (not shown). The upper panel shows exemplary a few
snapshots of the ensuing density fluctuations.

The number of sites is $M_{s}$. As boundary conditions we use either periodic
boundary conditions or hard-wall boundary conditions. In the case of hard-wall
boundary conditions, the particle density is strictly zero at the entire
boundary of the system, which extends beyond the sites plotted in Fig. 1. For
example, in 1D, where the sites are enumerated from $1$ to $M_{s}$, the hard-
wall boundary is positioned at the sites $0$ and $M_{s}+1$.
In 1D we use half-filling, i.e. the number of particles $N$ is given by
$N=M_{s}$ and the number of spin-up and spin-down particles is equal. As
initial states $|\Psi(0)\rangle$ we use $M_{s}/2$ adjacent doubly occupied
sites, either centered in the middle of the system as shown in Fig. 1 (a) or
aligned to the left edge. Similar uncorrelated initial conditions have been
used extensively as benchmarks for approximate methods to solve the time-
dependent multi-particle Schrödinger equation, see e.g. [35, 36].
To explore the effect of dimensionality, we extend the 1D system by
transversal degrees of freedom as depicted in Fig. 1 (b) and (c). In this case
the filling corresponds to $N=M_{s}/2$ in Fig. 1 (b) and $N=M_{s}/3$ in Fig. 1
(c), and the initial conditions used in these 2D cases are those depicted in
the figure.
All these states are initially entirely uncorrelated, i.e. correspond to
Hartree-Fock states. The correlation, however, strongly increases as a
function of time as will be discussed in Sec. IV, which is key to the
emergence of the (quasi-)condensate. The initial states are highly excited,
i.e. many states in the entire excitation spectrum have a non-negligible
overlap. One might think that an ensuing dynamical (quasi-)condensation effect
requires a particularly large overlap with the $\eta$-condensate
$|\eta\rangle$, which is an exact excited eigenstate of the Hamiltonian Eq. 1.
The overlap between our initial condition and the $\eta$-condensate is,
however, small as can be calculated analytically using the construction
$|\eta\rangle=\left(\hat{\eta}^{+}\right)^{N/2}|0\rangle$, where $|0\rangle$
is the vacuum state and

$\hat{\eta}^{+}=\sum_{j}(-1)^{j}\hat{a}_{j\downarrow}^{\dagger}\hat{a}_{j\uparrow}^{\dagger}.$
(2)

Using simple combinatorial arguments, we obtain in 1D in the case of half-
filling where $N=M_{s}$ that

$|\langle\eta|\Psi(0)\rangle|^{2}=\binom{M_{s}}{M_{s}/2}^{-1},$ (3)

and remains so for $t>0$. For large $M_{s}$ Eq. 3 goes like $\sim 2^{-M_{s}}$
and is thus exponentially suppressed. Therefore, although the dynamical
(quasi-)condensate described in this paper has the same energy expectation
value and similar properties as the $\eta$ condensate upon its dynamical
formation, it cannot be explained by an actual proximity to it.
To investigate the dynamics of the system in the following, we solve the
Schrödinger equation with the Hamiltonian Eq. 1 either by using exact
diagonalization or by employing the TD2RDM method as discussed below.

##  III TD2RDM method and $\eta$-symmetry

Within the TD2RDM method, we avoid the propagation of the wavefunction, and
thus the problem of exponential scaling altogether, and instead resort to the
propagation of the two-particle reduce density matrix (2RDM). The price we pay
is a partial neglect of three-particle correlations, as briefly described in
the following (for more details see [30, 31, 32]).
The 2RDM is obtained from a wavefunction $|\Psi(t)\rangle$ by tracing out
$N-2$ particles as

$D_{12}(t)=N(N-1)\text{Tr}_{3\ldots N}|\Psi(t)\rangle\langle\Psi(t)|.$ (4)

Similarly, the three-particle reduced density matrix (3RDM), required for the
propagation of the equations of motion of the 2RDM, is given by

$D_{123}(t)=N(N-1)(N-2)\text{Tr}_{4\ldots N}|\Psi(t)\rangle\langle\Psi(t)|,$
(5)

and relates to the 2RDM via $D_{12}=1/(N-2)\text{Tr}_{3}D_{123}$. For later
reference we also introduce the one-particle reduced density matrix (1RDM)
given by

$D_{1}(t)=1/(N-1)\text{Tr}_{2}D_{12}=N\text{Tr}_{2\ldots
N}|\Psi(t)\rangle\langle\Psi(t)|.$ (6)

In a given single-particle basis the elements of $D_{12}(t)$ can be written as

$D_{j_{1}\sigma_{1}j_{2}\sigma_{2}}^{i_{1}\sigma^{\prime}_{1}i_{2}\sigma^{\prime}_{2}}=\langle\Psi(t)|\hat{a}_{i_{1}\sigma^{\prime}_{1}}^{\dagger}\hat{a}_{i_{2}\sigma^{\prime}_{2}}^{\dagger}\hat{a}_{j_{2}\sigma_{2}}\hat{a}_{j_{1}\sigma_{1}}|\Psi(t)\rangle.$
(7)

We omit the numerical index (e.g. $12$ for the 2RDM) whenever it is obvious
from the number of indices in a specific single-particle basis which RDM we
consider (as in Eq. 7). In the present case, a suitable basis corresponds to
the individual sites of the Fermi-Hubbard model. In this basis, the
uncorrelated initial conditions considered in this paper can be easily
constructed as

$D_{i\uparrow j\downarrow}^{i\uparrow j\downarrow}=1$ (8)

for all doubly occupied sites $i$ and $j$ and zero otherwise. In the present
case of a total spin singlet case it turns out that the construction of the
entire $D_{12}$ is not necessary because the spin block
$D_{12}^{\uparrow\downarrow}$ contains all the information [30, 32]. The
equations of motion for this spin-block within the Fermi-Hubbard model are
then given by

$\displaystyle i\partial_{t}D_{j_{1}\uparrow j_{2}\downarrow}^{i_{1}\uparrow
i_{2}\downarrow}$ $\displaystyle=\sum_{n}h_{n}^{i_{1}}D_{j_{1}\uparrow
j_{2}\downarrow}^{n\uparrow
i_{2}\downarrow}+\sum_{n}h_{n}^{i_{2}}D_{j_{1}\uparrow
j_{2}\downarrow}^{i_{1}\uparrow n\downarrow}$
$\displaystyle+U\delta^{i_{1},i_{2}}D_{j_{1}\uparrow
j_{2}\downarrow}^{i_{1}\uparrow i_{2}\downarrow}$
$\displaystyle-\sum_{n}h_{j_{1}}^{n}D_{n\uparrow
j_{2}\downarrow}^{i_{1}\uparrow
i_{2}\downarrow}-\sum_{n}h_{j_{2}}^{n}D_{j_{1}\uparrow
n\downarrow}^{i_{1}\uparrow i_{2}\downarrow}$
$\displaystyle-U\delta_{j_{1},j_{2}}D_{j_{1}\uparrow
j_{2}\downarrow}^{i_{1}\uparrow i_{2}\downarrow}$
$\displaystyle+UD_{i_{1}\uparrow j_{2}\uparrow j_{1}\downarrow}^{i_{1}\uparrow
i_{2}\uparrow i_{1}\downarrow}+UD_{j_{1}\uparrow i_{2}\uparrow
j_{2}\downarrow}^{i_{1}\uparrow i_{2}\uparrow i_{2}\downarrow},$
$\displaystyle-UD_{j_{1}\uparrow j_{2}\uparrow j_{1}\downarrow}^{j_{1}\uparrow
i_{2}\uparrow i_{1}\downarrow}-UD_{j_{1}\uparrow j_{2}\uparrow
j_{2}\downarrow}^{i_{1}\uparrow j_{2}\uparrow i_{2}\downarrow}$ (9)

with

$h_{j}^{i}=-J\delta_{j}^{i+1}-J\delta_{j}^{i-1}.$ (10)

These equations scale like $M_{s}^{4}$ with the number of sites, which is key
to the efficiency of the TD2RDM method. The closure of the equations of motion
requires a reconstruction of the 3RDM via the 2RDM using an approximate stable
reconstruction functional [30, 31, 32]. Briefly, our reconstruction
functionals are based on the cumulant expansion of the 3RDM [37]

$D_{123}=\hat{A}D_{1}D_{2}D_{3}+\hat{A}\Delta_{12}D_{3}+\Delta_{123},$ (11)

which represents a separation into elements with different levels of particle
correlations. $\hat{A}$ is an anti-symmetrization operator that creates only
permutations that give non-equivalent terms,

$\Delta_{12}=D_{12}-\hat{A}D_{1}D_{2}$ (12)

is the two-particle cumulant representing two-particle correlations, and
$\Delta_{123}$ is the three-particle cumulant. Using the cumulant expansion,
an approximate reconstruction functional of $D_{123}$ boils down to finding
physically motivated reconstruction functionals for the three-particle
cumulant $\Delta_{123}$. We have shown previously that a stable and accurate
propagation of Eq. 9 requires that the reconstructed 3RDM correctly contracts
into the two-particle space [30]. Only then are conservation of energy and
spin symmetries guaranteed at all times during time propagation. This
contraction consistency [30] can be achieved employing the unitary
decomposition (see e.g. [38, 39, 40]), which allows to decompose a tensor into
its trace-free kernel and the orthogonal component (with respect to the
Hilbert-Schmidt inner product), which carries traces. The orthogonal component
is an exact functional of the traces of the tensor and can thus be easily
determined. We employ here the reconstruction functional by Valdemoro and
coworkers [41] and enforce construction consistency [30, 32]. Enforcing
contraction consistency leads to a scaling with $M_{s}^{5}$. Due to the
linearity of the unitary decomposition, this reconstruction of the 3RDM
amounts to neglecting the kernel (i.e. the trace free component) of
$\Delta_{123}$ in Eq. 11. We do not use the Nakatsuji-Yasuda reconstruction
[42, 32] of the three-particle cumulant here because the build up of the
(quasi)-condensate leads to large two-particle cumulants as a function of time
without concomitant increases in the three-particle cumulants such that the
Nakatsuji-Yasuda reconstruction, while being initially more accurate, leads to
overall larger reconstruction errors as time progresses.
The $\eta$-symmetry of the Hamiltonian leads to a further constant of motion
that has to be considered, i.e.
$[\hat{H},\hat{\eta}^{+}\hat{\eta}^{-}]=[\hat{H},\hat{\eta}_{z}]=0$ with
$\eta^{+}$ given in Eq. 2, $\eta^{-}$ being its hermitian conjugate, and
$\hat{\eta}_{z}=\frac{1}{2}\sum_{j}\left(\hat{n}_{j\uparrow}+\hat{n}_{j\downarrow}-1\right)$
[23]. While the conservation of the expectation value of
$\langle\eta_{z}\rangle$ is guaranteed by the conservation of the particle
number, the conservation of $\langle\eta^{+}\eta^{-}\rangle$ requires further
scrutiny. Using the 2RDM the expectation value
$\langle\eta^{+}\eta^{-}\rangle$ can easily be calculated as

$\langle\hat{\eta}^{+}\hat{\eta}^{-}\rangle=\sum_{i,j}(-1)^{i-j}D^{i\uparrow
i\downarrow}_{j\uparrow j\downarrow},$ (13)

and its behavior within the TD2RDM can be investigated by calculating its
time-derivative. The fact that the time-derivative of Eq. 13 has to be zero,
however, does not lead to further constraints on the 3RDM reconstruction
functional as it requires only exchange symmetry. (This is in contrast to e.g.
the conservation of energy [35] which requires contraction consistency [30]).
However, the TD2RDM suffers from the N-representability problem, which
necessitates the application of purification [40, 31, 43, 32]. We, therefore,
have to adapt our purification procedure, which so far took only into account
contraction consistency between the 2RDM and the 1RDM as well as energy
conservation [43, 32]. The modified purification procedure, which guarantees
the conservation of all constants of motion including
$\langle\eta^{+}\eta^{-}\rangle$, is described in App. A.

##  IV Dynamical (quasi-)condensation

Our figure of merit to asses the presence of a ferimionic (quasi-)condensate
is the largest eigenvalue of the 2RDM $g_{1}$ (i.e. the highest geminal
occupation number) obtained from diagonalization of the 2RDM

$D_{12}(t)=\sum_{j=1}^{r(r-1)/2}g_{j}(t)|g_{j}(t)\rangle\langle g_{j}(t)|,$
(14)

where $|g_{j}(t)\rangle$ are the geminal states (pair-states) and $r$ is the
number of single-particle basis states (orbitals). In our case $r=2M_{s}$. In
our convention, the geminal occupation numbers are ordered in descending order
and $g_{1}$ is the largest occupation number. The geminals can be grouped
according to their spin symmetry into spin-singlet and spin-tripled states. In
all cases reported below the largest geminal occupation number $g_{1}$ belongs
to a spin-singlet state $|g_{1}\rangle$.
Yang showed [24] that $g_{j}$ is bounded from above by

$g_{\text{max}}=\frac{N(r-N+2)}{r}$ (15)

for a system of $N$ fermions in $r$ modes and with the normalization of the
2RDM as given in Eq. 4. Clearly, $g_{\text{max}}$ is macroscopic, i.e. $O(N)$,
showing that the smallest reduced density matrix allowing macroscopic
occupations in fermionic systems is the 2RDM, in contrast to bosonic
condensation, which manifests itself by a macroscopic occupation of an
eigenstate of the 1RDM according to Penrose and Onsager [44]. The maximal
geminal occupation number for bosonic systems thus scales as $O(N^{2})$.
In case of half-filling, $N=r/2=M_{s}$, we obtain

$g_{\text{max}}=\frac{N+2}{2},$ (16)

and for the other cases considered here, i.e. $N=r/4=M_{s}/2$ and
$N=r/6=M_{s}/3$, we obtain $g_{\text{max}}=(3N+2)/4$ and
$g_{\text{max}}=(5N+2)/6$, respectively.
When $g_{\text{max}}$ is reached, all other $r(r-1)/2-1$ states have equal
weight given by [45]

$\bar{g}=\frac{2N(N-2)}{r(r-2)}.$ (17)

It has been shown that states reaching $g_{\text{max}}$ correspond to extreme
antisymmetric geminal power (AGP) states [45], which can be constructed by
antisymmetrization of a product ansatz of one pair-state. The
$\eta$-condensate is an example for the realization of these states.
In the search for the dynamical (quasi-)condensation we will monitor the time
evolution of $g_{1}(t)$, whose proximity to $g_{\text{max}}$ will signify the
presence of a fermionic (quasi-)condensate.

Figure 2: Distribution of geminal occupation numbers for several distinct
states. Blue: uncorrelated (Hartree-Fock) state with entropy $S_{2}^{0}$, red:
extreme AGP state that reaches the predicted upper bound of $g_{\text{max}}$
(Eq. 15) and has an entropy of $S_{2}^{*}$, and green: a state with the
highest entropy $S_{2}^{\text{max}}$ (Eq. 18).

When $g_{1}(t)$ approaches $g_{\text{max}}$, all other eigenvalues should
become close to equal. This behavior can be analyzed by means of a single
quantity, i.e. the entropy of the distribution of geminal occupation numbers
(see Fig. 2)

$S_{2}(t)=-\sum_{i=1}^{r(r-1)/2}g_{i}(t)\ln{g_{i}(t)}.$ (18)

We evaluate $S_{2}$ by renormalizing $\sum_{i=1}^{r(r-1)/2}g_{i}=1$. $S_{2}$
has been used as a measure of entanglement and correlations in the context of
fermionic systems in e.g. [46, 47]. We denote the 2RDM entropy for an extreme
AGP as $S_{2}^{*}$. Note that $S_{2}^{*}$ is close to but not equal to the
maximum $S_{2}^{\text{max}}$ of $S_{2}$ given by an equal distribution of
$g_{i}$, see Fig. 2.

###  IV.1 Small systems in 1D

We start our investigation with $M_{s}=8$ sites and exact results.

Figure 3: Exact results for the $M_{s}=8$ site Fermi-Hubbard model in 1D for
different $U$, boundary conditions, and alignments of the initial state. (a)
and (b) $U=0.1J$, (c) and (d) $U=4J$. Plotted is the maximal geminal
occupation number as a function of time, $g_{1}(t)$, in (a) and (c), and the
particle density at site one as a function of time, $\rho_{1}(t)$, in (b) and
(d). In all sub-figures solid lines correspond to hard wall boundary
conditions and a left aligned initial state, dashed lines correspond to hard
wall boundary conditions and a centered initial state [as in Fig. 1 (a) for
$M_{s}=8$], and dashed-dotted lines corresponds to periodic boundary
conditions and left aligned initial state.

Tuning $U$ from $U\ll 1J$ to $U\gg 1J$ we observe significant differences in
the system’s behavior, see Fig. 3. For $U=4J$ and periodic boundary
conditions, $g_{1}(t)$ fluctuates around $g_{1}(t)\approx 2$ as a function of
time regardless of the position of the initial state [see Fig. 3 (c) for the
example of a left-aligned initial state]. The system thermalizes within
$t\approx 5J^{-1}$, soon after the particle density completes a full circle
around the system [Fig. 3 (d)]. This behavior is also observed for hard wall
boundary and a centered initial condition, with thermalization occurring
before any significant increase in $g_{1}(t)$.
Significantly different behavior is observed when employing a left-aligned
initial state and hard wall boundary conditions. In this scenario, more space
(i.e. half of the system) is available for the system to expand freely until
reaching the boundary. During this period, a noticeable enhancement of the
pair-occupation number occurs, until $g_{1}$ reaches at the time
$T_{\text{max}}$ a maximum of about $g_{1}(T_{\text{max}})\approx 4$. Here and
in the following, $T_{\text{max}}$ represents the time of the first
significant maximum in $g_{1}(t)$ [see Fig. 3 (c)]. Overall, the process of
thermalization is slower, allowing for the observation of a second local
maximum during the ”back-and-forth” movement of the density, compare Fig. 3
(c) and (d). These observations are in line with the physics of emergent
Hamiltonians [28], elucidating the phenomenon of dynamical quasicondensation
in the context of hard-core bosons upon free expansion [25]: The large
interaction parameter $U$ induces a pairing between the fermions, facilitating
a mapping to hard-core bosons, also known as doublons (see e.g. [14]). If the
system can expand freely for long enough time, quasicondensation emerges.
However, the interaction with the boundary ultimately destroys the validity of
the emergent Hamiltonian and with it the quasicondensate.
In our analysis of the dynamical (quasi-)condensation effect, we focus on
evaluating $T_{\text{max}}$, and the associated amplitude of the geminal
occupation number $g_{1}(T_{\text{max}}$). For $U\geq 1J$ we observe that
$T_{\text{max}}$ scales linearly with $U$, i.e. $T_{\text{max}}\propto U$
[Fig. 4 (b)]. This behavior can be attributed to the rescaled hopping matrix
element of the doublons, which is proportional to $J^{2}/U$ (see e.g. [14]).
The dynamical (quasi-)condensation effect is least pronounced in the range of
$U\approx 1J$, where we observe a substantially smaller enhancement of
$g_{1}(T_{\text{max}})$ to values around $g_{1}(T_{\text{max}})\gtrsim 3.5$
(Fig. 5). As $U$ increases, $g_{1}(T_{\text{max}})$ shows a steady rise,
reaching saturation at $g_{1}(T_{\text{max}})\approx 4$.

For small $U<1J$, the behavior of the system is strikingly different. We
observe the emergence of a broad maximum, with $g_{1}(t)$ approaching
$g_{\text{max}}$ independently of the initial state’s position and the chosen
boundary conditions [Fig. 3 (a)].

Figure 4: The time of the first maximum in $g_{1}(t)$, $T_{\text{max}}$, (see
Fig. 3) as a function of $U$ for the 1D Fermi-Hubbard model with $M_{s}=8$,
hard-wall boundary conditions, and a left aligned initial state. Solid lines
are TD2RDM results, dashed lines correspond to exact results. The gray dotted
lines represent linear fits. (a) Interval $U\in[0.02,0.5]J$ with inverted
axis, (b) $U\in[1,25]J$. The inset in (a) shows the full $U$ interval in
linear scale. The gray dotted lines in each figure correspond to linear fits.

The specific quantitative value of $T_{\text{max}}$ is unaffected by the
initial state’s position and exhibits only a mild dependence on the selected
boundary conditions. Notably, the particle density undergoes multiple
fluctuations between the systems boundaries during the development of the
broad maximum in $g_{1}(t)$, without thermalizing the system. This observation
leads us to the conclusion that this phenomenon represents a novel dynamical
(quasi-)condensation effect, distinct from those captured by an emergent
Hamiltonian. We present a detailed analysis that further supports this
statement in App. B, where we show that the Hamiltonian Eq. 1 for $U<1$ does
not lead to an emergent time-depedent Hamiltonian that is (approximately)
conserved during time evolution as required by [28].
Accordingly, in contrast to the linear scaling with $U$ for $U>1J$, we observe
an inverse scaling behavior of $T_{\text{max}}$ for $U<1J$, i.e.
$T_{\text{max}}\propto 1/U$, see Fig. 4 (a).

Figure 5: The maximal value of the largest occupied geminal at
$T_{\text{max}}$, $g_{1}(T_{\text{max}})$, for different $U$ for the 1D Fermi-
Hubbard model with $M_{s}=8$, hard-wall boundary conditions, and a left
aligned initial state. Solid lines are TD2RDM results, dashed lines correspond
to exact results. The black dashed horizontal line marks the maximal possible
two-state occupation number $g_{\text{max}}=5$ for $M_{s}=8$. The inset shows
the full analyzed $U$ interval (the exact results extend towards larger
values).

$g_{1}(T_{\text{max}})$ in turn, increases with decreasing $U$ reaching again
values close to $g_{1}(T_{\text{max}})\approx 4$.
We now proceed to analyze whether these observations are captured by our
approximate TD2RDM method, as an underpinning for our studies with larger
systems, where the effect is more pronounced. For $U>1J$ we observe that the
TD2RDM method is in excellent agreement with the exact results, accurately
predicting the linear dependence of $T_{\text{max}}$ on $U$ and the slope. For
$U<1J$ the TD2RDM method correctly predicts the $1/U$ behavior but slightly
overestimates the proportionality constant. Simultaneously, the amplitude at
$T_{\text{max}}$ is overestimated but the overall behavior is again captured
very well, see Fig. 5.
To gain a deeper understanding of these discrepancies, we further scrutinize
the dynamics of all geminal occupation numbers as predicted by the approximate
TD2RDM method and compare them to the exact results.

Figure 6: Comparison between exact and TD2RDM results for the dynamics of all
geminal occupations numbers $g_{i}(t)$ for the 1D Fermi-Hubbard model with
$M_{s}=8$, hard-wall boundary conditions and a left aligned initial state. The
largest geminal occupation number $g_{1}(t)$ is drawn with a thicker line and
darker color. (a) Exact results for $U=0.1J$, (b) TD2RDM results for $U=0.1J$,
(c) exact results for $U=4J$, and (d) TD2RDM results for $U=4J$.

We observe that as $g_{1}(t)$ evolves towards its maximum, all other geminal
occupation numbers approach each other, consistent with the expectation that
when $g_{\text{max}}$ is reached, all other geminals are occupied by
$\bar{g}$, Eq. 17 [see Fig. 6 (a)]. However, the exact results show a small
time shift between the maximum of $g_{1}(t)$ and the time at which all other
$g_{i>1}(t)$ are closest. This leads to $g_{1}(T_{\text{max}})$ being
noticeably smaller than $g_{\text{max}}$. Furthermore, the exact results show
a small revival of the maximum at $t\approx 160J^{-1}$. In contrast, the
TD2RDM results lack this shift, resulting in an overestimation of the
amplitude $g_{1}(T_{\text{max}})$. Additionally, the TD2RDM results exhibit
prominent periodic revivals, indicating that relaxation effects are
underestimated. Notably, the TD2RDM prediction for $U=4J$ is remarkably
accurate, even capturing the third revival of $g_{1}(t)$.
This analysis can be complemented by means of the entropy $S_{2}(t)$, which
concentrates the information on all geminal occupation numbers into one
quantity. We observe that $S_{2}(t)$ increases as a function of time, reaching
$S_{2}^{*}$ for an extreme AGP state at a time close to $T_{\text{max}}$ [Fig.
7 (a)]. The TD2RDM results closely follow the exact results initially, but
then $S_{2}(t)$ remains large within the exact results while spuriously
fluctuating within the TD2RDM results. Another quantity that gives insights
into pair-correlations is the two-particle cumulant Eq. 12. Evaluating the
Frobenius norms of the two-particle cumulants for spin-polarized pairs,
$|\Delta_{12}^{\uparrow\uparrow}|_{2}$, and spin-unpolarized pairs
$|\Delta_{12}^{\uparrow\downarrow}|_{2}$ as a measure for two-particle
correlations reveals that the TD2RDM overestimates the pairing of spin-
unpolarized pairs while underestimating the production of spin-polarized
pairs, see Fig. 7 (b) and (c).
We, therefore, focus in our further analysis on the scaling behavior of
$T_{\text{max}}$ and $g_{1}(T_{\text{max}})$, which is overall very well
captured by the TD2RDM method, see Figs. 4 and 5.

Figure 7: Comparison between exact (dashed lines) and TD2RDM results (solid
lines) for (a) the dynamics of the entropy $S_{2}$, and the Frobenius norm of
two-particle cumulants (b) $|\Delta_{12}^{\uparrow\downarrow}|_{2}$ and (c)
$|\Delta_{12}^{\uparrow\uparrow}|_{2}$. In (a) the solid horizontal line marks
the maximal value of $S_{2}$, $S_{2}^{\text{max}}$, and the dotted horizontal
line denotes the value of $S_{2}$ obtained for a perfect condensate,
$S_{2}^{*}$. Red and orange lines are for $U=0.1J$, green and blues lines are
for $U=4J$. The system corresponds to the 1D Fermi-Hubbard model with
$M_{s}=8$ sites, hard wall boundary conditions and a left aligned initial
state.

We would like to point out that our previous results [32] indicate that the
TD2RDM method is more accurate for larger systems with $M_{s}\gtrsim 20$ than
for systems as small as $M_{s}=8$. While for $M_{s}=8$ the TD2RDM method
overestimates the density fluctuations for $t>60J^{-1}$ (not shown) consistent
with the deviations obtained for the geminal occupation numbers, we have
observed almost perfect agreement for the density fluctuations in case of
$M_{s}=18$ in [32]. The predictions for $T_{\text{max}}$ and the amplitude
$g_{1}(T_{\text{max}})$ within TD2RDM might, therefore, be even more accurate
than a straight forward extrapolation of the results for $M_{s}=8$ allows to
judge.

###  IV.2 Larger systems in 1D

We now turn to the fundamental question of whether the observed effect
exhibits characteristics of a quasicondensate that vanishes in the
thermodynamic limit, or whether this dynamical setting enables a circumvention
of the Mermin-Wagner-Hohenberg theorem which prohibits condensation in 1D in
equilibrium [34]. To address this question, we extend our investigation using
the TD2RDM method to explore significantly larger system sizes for which there
are no exact benchmarks available, bearing in mind that the TD2RDM might
slightly overestimate both $T_{\text{max}}$ as well as $g_{\text{max}}$ while
maintaining a reliable prediction of the overall behavior and scaling.
We start our investigation with a system of $M_{s}=20$ for which we make
detailed comparison to the system with $M_{s}=8$.

Figure 8: TD2RDM results for the Fermi-Hubbard model in 1D with $M_{s}=20$,
$U=0.1J$, hard-wall boundary conditions, and a centered initial condition. (a)
Geminal occupation numbers, (b) entropies $S_{2}$ (left y-axis) and $S_{1}$
(dashed line, right y-axis), and (c) particle-density fluctuations as a
function of time. The vertical lines mark the time $T_{\text{max}}$. In (b)
the upper solid horizontal line corresponds to $S_{2}^{\text{max}}$, the lower
corresponds to $S_{1}^{\text{max}}$, and the dotted horizontal line
corresponds to $S_{2}^{*}$.

Fig. 8 shows one particular example of the behavior for $M_{s}=20$. Similar to
the case of $M_{s}=8$, the largest geminal occupation number $g_{1}(t)$ rises
to a pronounced maximum of $g_{1}(T_{\text{max}})\approx 10$ close to
$g_{\text{max}}=11$ and several revivals appear [Fig. 8 (a)]. These strong
revivals are again most likely spurious.
It is instructive to explore in parallel the density fluctuations of the
system [Fig. 8 (c)]. We observe strong fluctuations in between the boundaries
of the system and several reflections from the hard walls except for a short
interval in time around $T_{\text{max}}$, where the particle density becomes
homogeneous. This effect is accompanied by a local maximum in the single-
particle entropy $S_{1}$ obtained from the diagonalization of the 1RDM (Eq. 6)

$D_{1}(t)=\sum_{j=1}^{r}n_{j}(t)|n_{j}(t)\rangle\langle n_{j}(t)|,$ (19)

given by

$S_{1}(t)=-\sum_{j=1}^{r}n_{j}(t)\ln{n_{j}(t)}.$ (20)

We have again renormalized $\sum_{i=1}^{r}n_{i}=1$ to calculate the entropy
$S_{1}$. The maximum at $T_{\text{max}}$ is given by
$S_{1}(T_{\text{max}})=3.68$, which is very close to the maximal value of
$S_{1}^{\text{max}}=3.69$ for equally distributed natural occupation numbers
for $r=2M_{s}=40$ [see Fig. 8 (b)]. Overall, $S_{1}(t)$ follows the curve of
the two-particle entropy $S_{2}(t)$.
The homogeneous distribution of the particle density over the entire system
follows from the fact that an extreme AGP state belongs to one of the quantum
many-body states that maximize $S_{1}$, and is a direct consequence of the
fact that the observed state is close to an extreme AGP state, with natural
orbitals distributed over the entire system and their occupation numbers being
almost equal. This striking effect of the particle density becoming
homogeneous during a small but finite time interval $\Delta T$ around
$T_{\text{max}}$ opens the door for an experimental study of these effects
within the platforms of ultracold quantum simulators. While the increase of
the pair-state occupation number is not easily accessible experimentally, the
monitoring of the density fluctuations of the system has become an
experimental routine on these platforms (see e.g. [48, 49, 50, 51, 52], for a
current experimental realization of a 1D periodic system see [53]) and would
allow to measure $T_{\text{max}}$ for even larger systems and other
geometries. $\Delta T$ depends on the width of the maximum in $g_{1}(t)$,
which increases with $T_{\text{max}}$, i.e. with $1/U$ [see Fig. 9 (a)]. In
Fig. 9 (a) we plot $\Delta T$ exemplary for hard wall boundary conditions by
determining the time interval for which the weight of the Fourier components
with $k>0$ for the density fluctuations in Fig. 8 (c) fall below a certain
limit.
When considering different boundary conditions and positions of the initial
state, we observe again only a weak dependence of the dynamical
(quasi-)condensation effect on them. $T_{\text{max}}$ as a function of $U$ is
practically equal for both boundary conditions, see Fig. 9 (a).

Figure 9: (a) $T_{\text{max}}$ and (b) the amplitude of the maximum of the
largest geminal occupation number $g_{1}(T_{\text{max}})$ as a function of $U$
for $M_{s}=20$ sites and different boundary conditions. In (a) on the second
y-axis we plot with crosses exemplary for hard wall boundary conditions
$\Delta T$ (i.e. the width of the time window, where the density fluctuations
become homogeneous). Figure 10: Absolute square of the momentum distribution
of the state $|g_{1}(T_{\text{max}})\rangle$ for $M_{s}=20$, $U=0.1J$ and (a)
periodic boundary conditions, and (b) hard-wall boundary. (c) Absolute square
of the momentum distribution of the $\eta$-condensate for $M_{s}=20$ and
periodic boundary conditions.

$g_{1}(T_{\text{max}})$ behaves similarly for the two different boundary
conditions but overall the decay with $U$ is slower in case of periodic
boundary conditions and less monotonic. A further analysis of the state
$|g_{1}(T_{\text{max}})\rangle$ reveals another interesting consequence of the
boundary conditions. When comparing the state $|g_{1}(T_{\text{max}})\rangle$
to the $\eta$-condensate in momentum space we observe for periodic boundary
conditions that $|g_{1}(T_{\text{max}})\rangle$ shows the typical
$\eta$-pairing of $k$ and $k-\pi$ for particles of different spin, see Fig.
10. $|g_{1}(T_{\text{max}})\rangle$ shows only small deviations from the
$\eta$-condensate [see Fig. 10 (c)] most notably a larger spread and a local
minimum across the $k_{2}=\pi-k_{1}$ lines in momentum space [Fig. 10 (a)]. In
contrast, for the system with hard-wall boundary conditions, the reflexions at
the boundaries lead to the emergence of additional pairings, most notably
Cooper-pair like pairs with $k_{1}=-k_{2}$, see Fig. 10 (b). Note that the
absolute square of the momentum distribution for the $\eta$-state is
practically equal for periodic and hard-wall boundary conditions, the only
difference being a different phase of the amplitudes and two more data points
in Fourier space due to explicitly taking into account the vanishing density
at site number $0$ and $M_{s}+1$ for hard-wall boundary conditions. Despite
these differences in the properties of the (quasi-)condensate state, the
geminal occupation numbers show similar behavior for both boundary conditions
[Fig. 9].
We now turn to analyzing the behavior of the (quasi)-condensation effect with
system size.

Figure 11: Dynamical quasi-condensation as functions of system size $M_{s}$.
(a) $T_{\text{max}}$, (b) $g_{1}(T_{\text{max}})$, and (c)
$g_{1}(T_{\text{max}})$ relative to the maximum $g_{\text{max}}$ (Eq. 15) for
hard-wall boundary conditions (solid), periodic boundary conditions (dashed),
$U=0.1J$ (dots), and $U=0.2J$ (squares). The gray dotted line in (b)
corresponds to a fit to a square root function.

Since parameter scans are numerically increasingly demanding with system size
we focus on two parameters for the on-site interaction, i.e. $U=0.1J$ and
$U=0.2J$. We observe for both values of $U$ that $T_{\text{max}}$ increases
with $M_{s}$ in a non-linear way [Fig. 11 (a)]. The influence of the boundary
conditions is negligible. In agreement with the previous investigations for
$M_{s}=8$ we observe that $T_{\text{max}}$ is larger for $U=0.1J$ than for
$U=0.2J$ for all investigated $M_{s}$. The amplitude of the first maximum
$g_{1}(T_{\text{max}})$ increases monotonically with $M_{s}$ [Fig. 11 (b)].
Most importantly, in the given available interval, the increase is
proportional to $\sqrt{M_{s}}$. Consequently, we observe a monotonic decrease
of $g_{1}(T_{\text{max}})$ with respect to $g_{\text{max}}$. Since
$g_{\text{max}}$ grows linearly with $M_{s}$ (Eq. 15) our results indicate
that $g_{1}(T_{\text{max}})/g_{\text{max}}\propto 1/\sqrt{M_{s}}$ in the
thermodynamic limit. In other words, if the $\sqrt{M_{s}}$-dependence of
$g_{1}(T_{\text{max}})$ is preserved in the thermodynamic limit of
$M_{s}\rightarrow\infty$, then in the thermodynamic limit,
$g_{1}(T_{\text{max}})/g_{\text{max}}\rightarrow 0$. Our results thus predict
that the observed effect is a dynamical quasi-condensate.
In the following section, we will investigate whether we might obtain a true
condensation in 2D, i.e. in the presence of transverse degrees of motion, by
extending the system to narrow 2D stripes.

###  IV.3 Extension to 2D systems

We restrict our study here to hard-wall boundary conditions as these are more
easily realizable experimentally in 2D systems. The geometries we investigate
here are depicted in Fig. 1 (b) and (c). For these systems with large $M_{s}$
systematic scans with $U$ and $M_{s}$ are increasingly expensive even within
TD2RDM such that we focus here on the question whether we get an indication of
a convergence to finite values of $g_{1}(T_{\text{max}})/g_{\text{max}}$ in
the thermodynamic limit. For this purpose we stick to $U=0.1J$ and $U=0.2J$
for different geometries and numbers of sites $M_{s}$.
When extending the 1D system by one transversal site along the entire length
as in Fig. 1 (b) we observe a build-up of a maximum in $g_{1}(t)$ over time,
see Fig. 12 (a) for different $M_{s}$ and a filling of $1/4$ (i.e.
$N=M_{s}/2=r/4$), and Fig. 12 (b) for $g_{1}(t)$ relative to $g_{\text{max}}$.

Figure 12: (a) The largest geminal occupation number $g_{1}(t)$, and (b)
$g_{1}(t)$ relative to $g_{\text{max}}$ as a function of time for $U=0.1J$ for
different 2D stripe geometries with $M_{s}=8\times 2$ and $N=8$,
$M_{s}=20\times 2$ and $N=20$, and $M_{s}=30\times 2$ and $N=30$. (c) The
largest value $g_{1}(T_{\text{max}})$ as a function of the number of sites
$M_{s}$ for the systems in (a). The gray dashed line represents a fit to
square root function of $M_{s}$.

In this case, the prediction for $g_{\text{max}}$ (Eq. 15) is equal to
$g_{\text{max}}=(3N+2)/4$, i.e. larger than in the case of half-filling, Eq.
16. It is striking that the maxima are much broader than in 1D extending over
longer periods of time. For $g_{1}(T_{\text{max}})$ we observe a sub-linear
(square root) increase with $M_{s}$ [Fig. 12 (c)], while convergence of
$g_{1}(T_{\text{max}})/g_{\text{max}}$ towards finite values would require a
linear increase. The conclusion from these results is thus that the
condensation effect will vanish in the thermodynamic limit.
Increasing the system further by another transverse degree of freedom as in
Fig. 1 changes the picture quite strongly, see Fig. 13. For the system with
$M_{s}=20\times 3$ sites we observe a fast increase in $g_{1}(t)$ initially
but then a flattening and a further increase at a much smaller pace. Until a
time of around $t=100J^{-1}$, which is already numerically quite costly, the
largest geminal occupation number $g_{1}$ does not come close to
$g_{\text{max}}=(5N+2)/6=17$ for this case of sixth filling. It is difficult
to estimate based on these observations whether the system develops a
(quasi-)condensate at some later point in time but the curves obtained so far
indicate that the effect does not persist in 2D.

Figure 13: The largest geminal occupation number $g_{1}(t)$ as a function of
time for $U=0.2J$ and number of particles $N=20$ and different stripe
geometries with $M_{s}=20\times 1$ (1D), $M_{s}=20\times 2$, and
$M_{s}=20\times 3$ sites amounting to different fillings.

##  V Conclusions

We have numerically investigated dynamical quasi-condensation in the Fermi-
Hubbard model starting from a completely uncorrelated initial state of
adjacent doubly occupied sites. In 1D we have shown that upon expansion of the
system the largest eigenvalue of the two-particle reduced density matrix
(2RDM), $g_{1}(t)$, develops a local maximum that comes close to the
theoretical upper limit $g_{\text{max}}$ predicted by Yang [24], signaling the
appearance of fermionic pair condensation. This dynamical quasi-condensation
is accompanied by strong two-particle correlations as measured by the two-
particle cumulants and the two-particle entropy obtained from the eigenstates
of the 2RDM. This condensation effect appears for all investigated values of
the interaction $U$, but shows a distinctly different behavior for small
$U<1J$ as compared to large $U>1J$, where $J$ is the hopping matrix element.
In the case of $U>1J$, the quasi-condensation effect has been explained by the
physics of an emergent Hamiltonian [28], and requires a completely free
expansion. In contrast, for small $U<1J$ we observe that the system undergoes
many interactions with the boundaries of the system on the characteristic time
scale $T_{\text{max}}$ during which the dynamical quasi-condensate emerges.
The effect is only weakly dependent on whether we chose hard-wall or periodic
boundary conditions. Interestingly, the two-particle quasi-condensate state
does dependent on the boundary conditions and features the typical
$\eta$-condensate pairing in case of periodic boundary conditions and both an
$\eta$-condensate pairing and Cooper-pair-like pairing in case of hard-wall
boundary conditions. Moreover, the scaling of both $T_{\text{max}}$ and the
amplitude $g_{1}(T_{\text{max}})$ with $U$ is markedly different for $U<1J$
and $U>1J$. For $U<1J$ we observe that $T_{\text{max}}\propto 1/U$, while for
$U>1J$ one obtains $T_{\text{max}}\propto U$.
It is so far an open question whether it is possible to circumvent the Mermin-
Wagner-Hohenberg theorem [33, 34] in a dynamical setting, which prohibits the
appearance of condensation in 1D in the thermodynamic limit in equilibrium. We
address this question for our weakly interacting quasi-condensation effect by
expanding to system sizes of up to $60$ sites. To propagate these systems over
the required long periods of time we employ our newly developed time-dependent
2RDM (TD2RDM) method. Based on comparisons with exact results for small system
sizes we have demonstrated that the TD2RDM method accurately predicts the
essential features of the effect such that extrapolation to large systems,
where exact benchmarks do not exist, can be made. The analysis of the effect
with increasing system size reveals that the maximum of $g_{1}$ grows
proportionally to the square root of the system size while $g_{\text{max}}$
grows linearly. Our results thus indicate that
$g_{1}(T_{\text{max}})/g_{\text{max}}$ vanishes in the thermodynamic limit, a
defining feature of a quasi-condensate.
To further expand on this question, we have extended our system to narrow 2D
stripes allowing for transversal motion during expansion. For systems with two
transverse sites, we observe broad maxima in $g_{1}$ over longer periods of
time compared to the 1D case. The scaling of $g_{1}(T_{\text{max}})$ with the
system size, however, again indicates a vanishing condensation effect in the
thermodynamic limit. A further increase, of the transverse degrees of freedom
seems to lead to a strong deterioration of the effect. In this case, $g_{1}$
increases strongly initially but then flattens, featuring a slow but steady
increase over the entire investigated time interval during which, however,
$g_{1}$ does not come as close to $g_{\text{max}}$ as in the previous cases.
Further increase in the number of transversal degrees of freedom to approach a
sufficiently large square lattice is planned in future.
Our results open the door to further scrutinize this effect in the platform of
experimental quantum simulators with single-site resolution [48, 49, 50, 51,
52] where substantially larger system sizes both in 1D and 2D could be probed.
We have shown that the appearance of the quasi-condensation effect is
accompanied with spatial particle density fluctuations becoming homogeneous
over the entire system for a short but finite period of time around
$T_{\text{max}}$. Strong density fluctuations reappear shortly after
$T_{\text{max}}$. This effect can be traced back to the single-particle
entropy obtained from the eigenvalues of the one-particle reduced density
matrix showing a local maximum. Since probing the particle density
distribution is an experimental routine nowadays, observing a homogeneous
density over finite periods of time could serve as a strong indication of the
dynamical (quasi-)condensation effect.
Finally, our results highlight the potential of the TD2RDM method to deliver
predictions for system sizes and time scales not reachable by any other
method. We would like to point out that our present code, while exploiting
several symmetries of the system, is not fully optimized and heavily
parallelized yet. Further improvement of its optimization and employing multi-
node parallelization should allow to compute the dynamics in systems with a
number of sites of about a factor of two larger than the largest systems in
the present study over similar time scales. With this extension,
investigations of condensation effects on a square lattice in 2D should become
possible. The TD2RDM method could thus provide benchmarks and guidance for
studies on non-equilibrium systems on quantum simulators for large systems.

## Acknowledgements

We thank Slavatore R. Manmana and Joseph Tindall for helpful discussions and
hints on literature. IB thanks the Simons Foundation for the great hospitality
and support during her research visit at the CCQ of the Flatiron Institute,
where parts of this research were conducted. The Flatiron Institute is a
division of the Simons Foundation. We acknowledge support from the Max Planck-
New York City Center for Non-Equilibrium Quantum Phenomena, Cluster of
Excellence ‘CUI: Advanced Imaging of Matter’- EXC 2056 - project ID. This
research was funded by the Austrian Science Fund (FWF) grant P 35539-N.
Calculations were performed on the Vienna Scientific Cluster (VSC4).

##  Appendix A Conservation of $\eta$-symmetry within purification and details
on the numerical implementation

In order to preserve N-representability of the 2RDM at least partially during
time propagation, we apply an a-posteriori purification procedure after
propagation time steps [30, 31, 43, 32]. Without purification the TD2RDM
method tends to produce 2RDMs with negative eigenvalues that ultimately may
lead to instabilities [35]. Purification entails removing iteratively the
defective part from the 2RDM while preserving its contraction to the 1RDM.
This can be facilitated by applying the unitary decomposition of the 2RDM
[40],

$D_{12}=D_{12;\perp}\\{D_{1}\\}+D_{12;K},$ (21)

where the kernel $D_{12;K}$ has vanishing traces, and all the trace
information is contained in the orthogonal component, $D_{12,\perp}$.
$D_{12,\perp}\\{D_{1}\\}$ indicates that the orthogonal component is a
functional of the 1RDM. We have empirically shown [30, 31, 32] that restoring
the positive-semidefinitness of the 2RDM and the two-hole RDM is sufficient to
stabilize the equations of motion in most cases [32]. Accordingly, we
determine the defective part of the 2RDM

$D_{12}^{<}=\sum_{g_{i}<0}g_{i}|g_{i}\rangle\langle g_{i}|$ (22)

and similarly the defective part of the two-hole RDM, $Q_{12}^{<}$.
Subtracting the kernel of the defective components of both the 2RDM and the
two-hole RDM from the 2RDM iteratively leads to a positive semi-definite 2RDM,
while preserving contraction consistency to the 1RDM. In general, however,
other symmetries like e.g. energy conservation are broken in this process
unless further amendments are applied. It has been shown [43] that energy
conservation can be restored by removing from $D_{12}^{<}$ and $Q_{12}^{<}$
the contributions that would lead to violations of energy conservation, before
calculating the kernel. Similarly, to enforce conservation of
$\langle\hat{\eta}^{+}\hat{\eta}-\rangle$, Eq. 13, we set ${D^{<}}_{j\uparrow
j\downarrow}^{i\uparrow i\downarrow}$ and ${Q^{<}}_{j\uparrow
j\downarrow}^{i\uparrow i\downarrow}$to zero for all $i$ and $j$. We then
determine the kernel of the correspondingly corrected defective components
$D_{12}^{<}$ and $Q_{12}^{<}$, and subtract them from the 2RDM. This procedure
is iterative repeated until the required threshold for the smallest geminal
occupation number is reached, or until the selected maximal number of
iterations is reached.
For the studies in the main text we have applied a threshold of $0$ for the
smallest geminal occupation number and at most $10$ iteration steps within the
purification for all systems except for $M_{s}=60$, where only one step is
applied due to its scaling with $M_{s}^{6}$ within the current implementation
using full diagonalization of the 2RDM. Purification is applied after each
global time step of $dt=0.02J^{-1}$ within which the equations of motion, Eq.
9, are solved with adaptive time steps using a Runge-Kutta-Fehlberg propagator
of $4^{th}$ and $5^{th}$ order. The different orders are used to adapt the
time steps to achieve a prescribed accuracy. We have checked that changing the
global time step to $dt=0.01J^{-1}$, i.e. applying up to twice as many
purification steps, leads to the same results to high degree of accuracy.

##  Appendix B Formalism of emergent Hamiltonians

We explore the applicability of the formalism of emergent Hamiltonians to the
dynamical quasi-condensation effect for $U<1$. Following [28], we define a
Hamiltonian $\hat{H}_{0}$, such that the initial state is an eigenstate of
$\hat{H}_{0}$. One particular choice following Ref. [28] for the present
initial condition is

$\hat{H}_{0}=\frac{1}{M_{s}}\sum_{j=1}^{M_{s}}j(\hat{n}_{j}^{\uparrow}+\hat{n}_{j}^{\downarrow}).$
(23)

There are, however, other possible choices such as
$\hat{H}_{0}=U\sum_{j=1}^{M_{s}}\hat{n}_{j}^{\uparrow}\hat{n}_{j}^{\downarrow}$,
but they lead to essentially the same conclusions, see below. The system
evolves under the influence of $\hat{H}$, Eq. 1, such that
$\hat{P}=\hat{H}_{0}-\hat{H}$ can be regarded as a quench. Starting from the
initial condition $(\hat{H}_{0}-\lambda)|\Psi(0)\rangle=0$ one applies the
time evolution operator to obtain

$\left(e^{-i\hat{H}t}\hat{H}_{0}e^{i\hat{H}t}-\lambda\right)|\Psi(t)\rangle:=\hat{M}(t)|\Psi(t)\rangle.$
(24)

The operator $\hat{M}(t)$ is given by

$\begin{split}\hat{M}(t)&=\hat{H}-\lambda+\hat{P}-it[\hat{H},\hat{P}]+\frac{(-it)^{2}}{2}\left[\hat{H},[\hat{H},\hat{P}]\right]+...\\\
&=\hat{H}-\lambda+\hat{P}-it\hat{Q}+\sum_{n=1}^{\infty}\frac{(-it)^{n+1}}{(n+1)!}\mathcal{\hat{H}}_{n},\end{split}$
(25)

with the definition $\hat{Q}=[\hat{H},P]$, and $\mathcal{\hat{H}}_{n}$
describing nested commutators starting with
$\mathcal{\hat{H}}_{1}=[\hat{H},\hat{Q}]$. According to [28] the emergent
Hamiltonian is applicable if $\hat{M}(t)$ is a local operator. One particular
family studied corresponds to $[\hat{H},\hat{Q}]=0$ up to boundary terms,
where the time propagated state is exponentially close to the eigenstate of
the emergent Hamiltonian for times proportional to the system size [28].
For the following discussion it is convenient to separate our Hamiltonian Eq.
1 as $\hat{H}=\hat{T}+\hat{W}$, where $\hat{T}$ stands for the hopping
operator proportional to $J$ and $\hat{W}$ is the on-site interaction
proportional to $U$. For our system with hard wall boundary conditions
(similar results hold for periodic boundary conditions) we obtain

$\begin{split}\hat{Q}=[\hat{H},\hat{P}]=[\hat{T},\hat{P}]&=\frac{J}{M_{s}}\sum_{j=1}^{M_{s}-1}\sum_{\sigma}\hat{a}_{j+1\sigma}^{\dagger}\hat{a}_{j\sigma}\\\
&-\frac{J}{M_{s}}\sum_{j=1}^{M_{s}-1}\sum_{\sigma}\hat{a}_{j\sigma}^{\dagger}\hat{a}_{j+1\sigma}\end{split}$
(26)

(note that with the prefactor $(-i)$ from Eq. 25 this contribution to
$\hat{M}(t)$ is Hermitian as it should be), and

$\begin{split}\mathcal{\hat{H}}_{1}=[\hat{H},\hat{Q}]=\frac{2J^{2}}{M_{s}}&\left(n_{M_{s}}^{\uparrow}+n_{M_{s}}^{\downarrow}-n_{1}^{\uparrow}-n_{1}^{\downarrow}\right)\\\
+\frac{UJ}{M_{s}}&\Bigg{[}\sum_{j=1}^{M_{s}-1}\hat{a}_{j+1\uparrow}^{\dagger}a_{j\uparrow}\left(\hat{n}_{j+1}^{\downarrow}-\hat{n}_{j}^{\downarrow}\right)+\text{H.c.}\\\
&+\sum_{j=1}^{M_{s}-1}\hat{a}_{j+1\downarrow}^{\dagger}a_{j\downarrow}\left(\hat{n}_{j+1}^{\uparrow}-\hat{n}_{j}^{\uparrow}\right)+\text{H.c.}\Bigg{]}.\end{split}$
(27)

The first contribution comes from the commutator between $\hat{T}$ and
$\hat{Q}$ (and corresponds to the one found in [28] for the non-interacting
system), while the second term comes from the commutator between $\hat{W}$ and
$\hat{Q}$. Clearly, the commutator $[\hat{H},\hat{Q}]$ does not vanish and
leads to a term growing quadratically in time in $\hat{M}(t)$, Eq. 25. Higher
orders $\mathcal{\hat{H}}_{n}$ contain elements of the form

$\begin{split}\frac{U^{n}J}{M_{s}}&\Bigg{[}\sum_{j=1}^{M_{s}-1}\hat{a}_{j+1\uparrow}^{\dagger}a_{j\uparrow}\left(\hat{n}_{j+1}^{\downarrow}-\hat{n}_{j}^{\downarrow}\right)^{n}+(-1)^{n+1}\text{H.c.}\\\
&+\sum_{j=1}^{M_{s}-1}\hat{a}_{j+1\downarrow}^{\dagger}a_{j\downarrow}\left(\hat{n}_{j+1}^{\uparrow}-\hat{n}_{j}^{\uparrow}\right)^{n}+(-1)^{n+1}\text{H.c.}\Bigg{]},\end{split}$
(28)

which originate from the evaluation of the $n$ nested commutators
$[\hat{W},[\hat{W},\ldots,[\hat{W},\hat{Q}]]]$. In addition,
$\mathcal{\hat{H}}_{n}$ contains all lower orders in $U$, i.e. $U^{j}$ with
$j\in[0,n-1]$ originating from different sequences of $\hat{T}$ and $\hat{W}$
in the nested commutators and leading to products of $2(j+1)$ creation and
annihilation operators. These products contain elements of the form
$a_{j}^{\dagger}a_{j+m}$ with $m\leq n$ (up to $m=M_{s}-1$, i.e. up to the
span over the entire system). The contribution independent of $U$ originates
from repeated application of the commutator with $\hat{T}$ and contains only
boundary terms. $\mathcal{M}(t)$ therefore contains arbitrary high orders of
products of creation and annihilation operators with an increasing span over
the entire system, albeit being suppressed by $U^{m}$ with $U<1$ in our case
and $m$ the number of commutators containing $\hat{W}$. While the
contributions from Eq. 28 to Eq. 25 can be evaluated to all orders leading to
an extensive sum of local exponential operators, the other terms amount to an
infinite series of non-local operators and thus cannot be easily treated.
These results underline our conclusion that the present effect is not captured
by an emergent Hamiltonian.

## References

## References

  * Abrikosov  _et al._ [1965] A. A. Abrikosov, L. P. Gor’kov, and I. Y. Dzyaloshinskii,  _Quantum Field Theoretical Methods in Statistical Physics_ (Pergamon Press Ltd., Oxford, 1965).
  * Fausti  _et al._ [2011] D. Fausti, R. I. Tobey, N. Dean, S. Kaiser, A. Dienst, M. C. Hoffmann, S. Pyon, T. Takayama, H. Takagi, and A. Cavalleri, Light-Induced Superconductivity in a Stripe-Ordered Cuprate, Science 331, 189 (2011).
  * Nicoletti  _et al._ [2014] D. Nicoletti, E. Casandruc, Y. Laplace, V. Khanna, C. R. Hunt, S. Kaiser, S. S. Dhesi, G. D. Gu, J. P. Hill, and A. Cavalleri, Physical Review B 90, 100503(R) (2014).
  * Cremin  _et al._ [2019] K. A. Cremin, J. Zhang, C. C. Homes, G. D. Gu, Z. Sun, M. M. Fogler, A. J. Millis, D. N. Basov, and R. D. Averitt, Photoenhanced metastable c-axis electrodynamics in stripe-ordered cuprate La ${}_{\textrm{1.885}}$ Ba ${}_{\textrm{0.115}}$ CuO ${}_{\textrm{4}}$, Proceedings of the National Academy of Sciences 116, 19875 (2019).
  * Rajasekaran  _et al._ [2018] S. Rajasekaran, J. Okamoto, L. Mathey, M. Fechner, V. Thampy, G. D. Gu, and A. Cavalleri, Probing optically silent superfluid stripes in cuprates, Science 359, 575 (2018).
  * Suzuki  _et al._ [2019] T. Suzuki, T. Someya, T. Hashimoto, S. Michimae, M. Watanabe, M. Fujisawa, T. Kanai, N. Ishii, J. Itatani, S. Kasahara, Y. Matsuda, T. Shibauchi, K. Okazaki, and S. Shin, Photoinduced possible superconducting state with long-lived disproportionate band filling in FeSe, Communications Physics 2, 115 (2019).
  * Mitrano  _et al._ [2016] M. Mitrano, A. Cantaluppi, D. Nicoletti, S. Kaiser, A. Perucchi, S. Lupi, P. Di Pietro, D. Pontiroli, M. Riccò, S. R. Clark, D. Jaksch, and A. Cavalleri, Possible light-induced superconductivity in K3C60 at high temperature, Nature 530, 461 (2016).
  * Budden  _et al._ [2021] M. Budden, T. Gebert, M. Buzzi, G. Jotzu, E. Wang, T. Matsuyama, G. Meier, Y. Laplace, D. Pontiroli, M. Riccò, F. Schlawin, D. Jaksch, and A. Cavalleri, Evidence for metastable photo-induced superconductivity in K3C60, Nature Physics 17, 611 (2021).
  * Buzzi  _et al._ [2020] M. Buzzi, D. Nicoletti, M. Fechner, N. Tancogne-Dejean, M. Sentef, A. Georges, T. Biesner, E. Uykur, M. Dressel, A. Henderson, T. Siegrist, J. Schlueter, K. Miyagawa, K. Kanoda, M.-S. Nam, A. Ardavan, J. Coulthard, J. Tindall, F. Schlawin, D. Jaksch, and A. Cavalleri, Photomolecular High-Temperature Superconductivity, Physical Review X 10, 031028 (2020).
  * Bloch  _et al._ [2022] J. Bloch, A. Cavalleri, V. Galitski, M. Hafezi, and A. Rubio, Strongly correlated electron–photon systems, Nature 606, 41 (2022).
  * Sentef  _et al._ [2016] M. A. Sentef, A. F. Kemper, A. Georges, and C. Kollath, Theory of light-enhanced phonon-mediated superconductivity, Phys. Rev. B 93, 144506 (2016).
  * Tindall  _et al._ [2019] J. Tindall, B. Buča, J. Coulthard, and D. Jaksch, Heating-Induced Long-Range $\eta$ Pairing in the Hubbard Model, Physical Review Letters 123, 030603 (2019).
  * Kaneko  _et al._ [2019] T. Kaneko, T. Shirakawa, S. Sorella, and S. Yunoki, Photoinduced $\eta$ Pairing in the Hubbard Model, Physical Review Letters 122, 077002 (2019).
  * Cook and Clark [2020] M. W. Cook and S. R. Clark, Controllable finite-momenta dynamical quasicondensation in the periodically driven one-dimensional Fermi-Hubbard model, Physical Review A 101, 033604 (2020).
  * Paeckel  _et al._ [2020] S. Paeckel, B. Fauseweh, A. Osterkorn, T. Köhler, D. Manske, and S. R. Manmana, Detecting superconductivity out of equilibrium, Physical Review B 101, 180507(R) (2020).
  * Buzzi  _et al._ [2021] M. Buzzi, G. Jotzu, A. Cavalleri, J. I. Cirac, E. A. Demler, B. I. Halperin, M. D. Lukin, T. Shi, Y. Wang, and D. Podolsky, Higgs-mediated optical amplification in a nonequilibrium superconductor, Phys. Rev. X 11, 011055 (2021).
  * Dolgirev  _et al._ [2022] P. E. Dolgirev, A. Zong, M. H. Michael, J. B. Curtis, D. Podolsky, A. Cavalleri, and E. Demler, Periodic dynamics in superconductors induced by an impulsive optical quench, Communications Physics 5, 234 (2022).
  * Hubbard [1963] J. Hubbard, Electron correlations in narrow energy bands, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 276, 238 (1963).
  * Hubbard [1964] J. Hubbard, Electron correlations in narrow energy bands III. An improved solution, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 281, 401 (1964).
  * Gutzwiller [1963] M. C. Gutzwiller, Effect of Correlation on the Ferromagnetism of Transition Metals, Physical Review Letters 10, 159 (1963).
  * Kanamori [1963] J. Kanamori, Electron Correlation and Ferromagnetism of Transition Metals, Progress of Theoretical Physics 30, 275 (1963).
  * Qin  _et al._ [2022] M. Qin, T. Schäfer, S. Andergassen, P. Corboz, and E. Gull, The Hubbard Model: A Computational Perspective, Annual Review of Condensed Matter Physics 13, 275 (2022).
  * Yang [1989] C. N. Yang, $\eta$ pairing and off-diagonal long-range order in a Hubbard model, Physical Review Letters 63, 2144 (1989).
  * Yang [1962] C. N. Yang, Concept of Off-Diagonal Long-Range Order and the Quantum Phases of Liquid He and of Superconductors, Reviews of Modern Physics 34, 694 (1962).
  * Rigol and Muramatsu [2004] M. Rigol and A. Muramatsu, Emergence of Quasicondensates of Hard-Core Bosons at Finite Momentum, Physical Review Letters 93, 230404 (2004).
  * Rodriguez  _et al._ [2006] K. Rodriguez, S. R. Manmana, M. Rigol, R. M. Noack, and A. Muramatsu, Coherent matter waves emerging from Mott-insulators, New Journal of Physics 8, 169 (2006).
  * Vidmar  _et al._ [2015] L. Vidmar, J. P. Ronzheimer, M. Schreiber, S. Braun, S. S. Hodgman, S. Langer, F. Heidrich-Meisner, I. Bloch, and U. Schneider, Dynamical quasicondensation of hard-core bosons at finite momenta, Phys. Rev. Lett. 115, 175301 (2015).
  * Vidmar  _et al._ [2017] L. Vidmar, D. Iyer, and M. Rigol, Emergent Eigenstate Solution to Quantum Dynamics Far from Equilibrium, Physical Review X 7, 021012 (2017).
  * Heidrich-Meisner  _et al._ [2008] F. Heidrich-Meisner, M. Rigol, A. Muramatsu, A. E. Feiguin, and E. Dagotto, Ground-state reference systems for expanding correlated fermions in one dimension, Physical Review A 78, 013620 (2008).
  * Lackner  _et al._ [2015] F. Lackner, I. Březinová, T. Sato, K. L. Ishikawa, and J. Burgdörfer, Propagating two-particle reduced density matrices without wave functions, Physical Review A 91, 023412 (2015).
  * Lackner  _et al._ [2017] F. Lackner, I. Březinová, T. Sato, K. L. Ishikawa, and J. Burgdörfer, High-harmonic spectra from time-dependent two-particle reduced-density-matrix theory, Physical Review A 95, 033414 (2017).
  * Donsa  _et al._ [2023] S. Donsa, F. Lackner, J. Burgdörfer, M. Bonitz, B. Kloss, A. Rubio, and I. Březinová, Nonequilibrium correlation dynamics in the one-dimensional Fermi-Hubbard model: A testbed for the two-particle reduced density matrix theory, Physical Review Research 5, 033022 (2023).
  * Mermin and Wagner [1966] N. D. Mermin and H. Wagner, Absence of Ferromagnetism or Antiferromagnetism in One- or Two-Dimensional Isotropic Heisenberg Models, Physical Review Letters 17, 1133 (1966).
  * Hohenberg [1967] P. C. Hohenberg, Existence of Long-Range Order in One and Two Dimensions, Physical Review 158, 383 (1967).
  * Akbari  _et al._ [2012] A. Akbari, M. J. Hashemi, A. Rubio, R. M. Nieminen, and R. van Leeuwen, Challenges in truncating the hierarchy of time-dependent reduced density matrices equations, Physical Review B 85, 235121 (2012).
  * Schlünzen  _et al._ [2017] N. Schlünzen, J.-P. Joost, F. Heidrich-Meisner, and M. Bonitz, Nonequilibrium dynamics in the one-dimensional Fermi-Hubbard model: Comparison of the nonequilibrium Green-functions approach and the density matrix renormalization group method, Physical Review B 95, 165139 (2017).
  * Mazziotti [1998] D. A. Mazziotti, Approximate solution for electron correlation through the use of Schwinger probes, Chemical Physics Letters 289, 419 (1998).
  * Absar and Coleman [1976] I. Absar and A. Coleman, One electron orbitals intrinsic to the reduced hamiltonian, Chemical Physics Letters 39, 609 (1976).
  * Harriman [1978] J. E. Harriman, Geometry of density matrices. ii. reduced density matrices and n-representability, Phys. Rev. A 17, 1257 (1978).
  * Mazziotti [2002] D. A. Mazziotti, Purification of correlated reduced density matrices, Phys. Rev. E 65, 026704 (2002).
  * Colmenero  _et al._ [1993] F. Colmenero, C. Pérez del Valle, and C. Valdemoro, Approximating q-order reduced density matrices in terms of the lower-order ones. I. General relations, Physical Review A 47, 971 (1993).
  * Yasuda and Nakatsuji [1997] K. Yasuda and H. Nakatsuji, Direct determination of the quantum-mechanical density matrix using the density equation. II., Physical Review A 56, 2648 (1997).
  * Joost  _et al._ [2022] J.-P. Joost, N. Schlünzen, H. Ohldag, M. Bonitz, F. Lackner, and I. Březinová, Dynamically screened ladder approximation: Simultaneous treatment of strong electronic correlations and dynamical screening out of equilibrium, Physical Review B 105, 165155 (2022).
  * Penrose and Onsager [1956] O. Penrose and L. Onsager, Bose-einstein condensation and liquid helium, Phys. Rev. 104, 576 (1956).
  * Coleman and Yukalov [2000] A. J. Coleman and V. I. Yukalov,  _Reduced Density Matrices, Coulson’s Challenge_ (Springer-Verlag Berlin Heidelberg, 2000).
  * Gigena  _et al._ [2021] N. Gigena, M. Di Tullio, and R. Rossignoli, Many-body entanglement in fermion systems, Physical Review A 103, 052424 (2021).
  * Ferreira  _et al._ [2022] D. L. B. Ferreira, T. O. Maciel, R. O. Vianna, and F. Iemini, Quantum correlations, entanglement spectrum, and coherence of the two-particle reduced density matrix in the extended Hubbard model, Physical Review B 105, 115145 (2022).
  * Cheuk  _et al._ [2015] L. W. Cheuk, M. A. Nichols, M. Okan, T. Gersdorf, V. V. Ramasesh, W. S. Bakr, T. Lompe, and M. W. Zwierlein, Quantum-Gas Microscope for Fermionic Atoms, Physical Review Letters 114, 193001 (2015).
  * Parsons  _et al._ [2015] M. F. Parsons, F. Huber, A. Mazurenko, C. S. Chiu, W. Setiawan, K. Wooley-Brown, S. Blatt, and M. Greiner, Site-Resolved Imaging of Fermionic Li 6 in an Optical Lattice, Physical Review Letters 114, 213002 (2015).
  * Haller  _et al._ [2015] E. Haller, J. Hudson, A. Kelly, D. A. Cotta, B. Peaudecerf, G. D. Bruce, and S. Kuhr, Single-atom imaging of fermions in a quantum-gas microscope, Nature Physics 11, 738 (2015).
  * Edge  _et al._ [2015] G. J. A. Edge, R. Anderson, D. Jervis, D. C. McKay, R. Day, S. Trotzky, and J. H. Thywissen, Imaging and addressing of individual fermionic atoms in an optical lattice, Physical Review A 92, 063406 (2015).
  * Omran  _et al._ [2015] A. Omran, M. Boll, T. A. Hilker, K. Kleinlein, G. Salomon, I. Bloch, and C. Gross, Microscopic Observation of Pauli Blocking in Degenerate Fermionic Lattice Gases, Physical Review Letters 115, 263001 (2015).
  * Cai  _et al._ [2022] Y. Cai, D. G. Allman, P. Sabharwal, and K. C. Wright, Persistent Currents in Rings of Ultracold Fermionic Atoms, Physical Review Letters 128, 150401 (2022).