email<EMAIL_ADDRESS>email<EMAIL_ADDRESS>
# Assessing quantum thermalization in physical and configuration spaces
via many-body weak values

Carlos F. Destefani1 [  Xavier Oriols1 [  1 Department of Electronic
Engineering, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona,
Spain

###### Abstract

We explore the origin of the arrow of time in an isolated quantum system
described by the Schrödinger equation. We provide an explanation from weak
values in the configuration space, which are understood as operational
properties obtained in the laboratory following a well-defined protocol. We
show that quantum systems satisfying the eigenstate thermalization hypothesis
can simultaneously provide thermalized ensemble expectation values and
nonthermalized weak values of the momentum, both from the same operational
probability distribution. The reason why weak values of the momentum may
escape from the eigenstate thermalization hypothesis is because they are
linked only to off-diagonal elements of the density matrix in the energy
representation. For indistinguishable particles, however, operational
properties can not be defined in the configuration space. Therefore, we state
that the origin of the arrow of time in isolated quantum systems described by
the Schrödinger equation comes from dealing with properties obtained by
averaging (tracing out) some degrees of freedom of the configuration space. We
then argue that thermalization does not occur in the properties defined in the
configuration space, and our argument is compatible with defending that
thermalization is a real phenomenon in the properties defined in the physical
space. All of these conclusions are testable in the laboratory through many-
body weak values.

##  I Introduction

The arrow of time has always been a topic of lively debate [1, 2, 3, 4, 5, 6,
7, 8]. It appears in many disciplines as, for example, a cosmological arrow of
time pointing in the direction of the Universe expansion [2]. A related arrow
of time appears in the time evolution of time-irreversible macroscopic systems
governed by the second law of thermodynamics, where entropy always increases
with time [1, 7]. The puzzle implicit in such irreversibility is that most
fundamental microscopic laws have no arrow of time. They are time-reversible
laws, in the sense that time appears as a variable, just like position,
without any privileged direction. But, if macroscopic laws emerge from
microscopic laws, what can make them so different? A possible explanation is
that such fundamental laws are in fact not time-reversal. For example, it has
been argued that the time-reversible Schrödinger equation is not the true law
at a fundamental level, and that it should be substituted by laws from
spontaneous collapse theories which, by construction, are time-irreversible
[3]. Another explanation argues that real systems are never perfectly
isolated, so that time-reversible fundamental laws, despite being the true
laws, are not directly applicable [4]. And yet another argumentation claims
that the evolution of a real system depends, apart from the true time-
reversible microscopic laws, on the initial conditions which provide an
irreversible time evolution [6].

In this paper we explore a different path, via weak values in the
configuration space, to understand the physical origins of the arrow of time
in perfectly isolated nonrelativistic systems described by the Schrödinger
equation, assumed as a true reversible law where initial conditions are not
relevant to explain the observed time-irreversibility in our results. Our goal
is to link the evolution of microscopic (or macroscopic) properties of a model
system with the configuration (or physical) space where such properties are
defined. The wave function solution of the Schrödinger equation is defined in
a $3N$-dimensional configuration space, while typical properties where time-
irreversibility is observed, are defined in smaller spaces where some degrees
of freedom of the configuration space are integrated. Thus, the question that
motivates us is whether the presence or absence of an arrow of time in the
time evolution of the properties of a quantum system is a consequence of
defining them in the full configuration space or in a smaller space.

For our goal, the renewed interest in statistical mechanics of closed quantum
systems [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]
provides the perfect scenario. Such interest has been generated by the
successful experimental ability to isolate and manipulate bosonic [25, 26, 27,
28, 29] and fermionic [30, 31, 32, 33, 34] many-body systems built on ultra-
cold atomic gases subjected to optical lattices. Such quantum systems can be
described by the many-particle Schrödinger equation; an arrow of time appears
because, despite these systems being expected to present unitary evolution,
some of their initial nonequilibrium nonthermalized properties may later
thermalize.

A preliminary consideration is that it is not at all obvious whether the
configuration space is more or less fundamental than the ordinary physical
space. The nonrelativistic Schrödinger equation can be seen as a sort of
approximation to the relativistic quantum field theory [35]. For
indistinguishable particles, quantum field theory does not require knowledge
of the exact position of each particle in the configuration space, but only of
how many particles are present in a position of the physical space and, as
such, configuration space seems less fundamental than physical space.

Another consideration is that thermalization is typically reported in
properties mensurable in a laboratory, so that one needs to discuss, within
similar empirical protocols, nonthermalized properties also mensurable in a
laboratory. But making conclusions testable in a laboratory opens new
difficulties since, strictly speaking, a measured closed system is no longer a
closed system because of its interaction with the measuring apparatus [36].
Such measurement produces a collapse of the quantum state of the isolated
system, which is at the origin of quantum randomness and backaction. Within
the orthodox theory, such a collapse requires a time-asymmetric law, different
from the time-symmetric Schrödinger equation, so that an “orthodox quantum-
mechanical arrow of time” seems to enter into play [37]. We will see how the
weak values protocol provides operational properties of the quantum system
without backaction or quantum randomness (avoiding the role of the collapse in
our discussion).

It is important to differentiate between operational (empirical) and
ontological (real) properties of a system. Operational properties are those
whose definition comes exclusively from operations done on the laboratory over
the system (with or without ontological meaning for that property), which are
defined independently on any quantum theory. Ontological properties, on the
other hand, are those that a specific quantum theory postulates to exist.
Therefore, it is possible that a given operational property coincides with an
ontological property in one theory, but not in another. A typical example is
the velocity of a particle, known to be an operational property computed from
weak values, independently on any quantum theory; such operational property
coincides with an ontological property in the Bohmian theory (the velocity
itself), but it is not an ontological property in the Orthodox theory. It is
far from our scope to discuss whether or not thermalization is an ontological
phenomenon or not occurring in the configuration space or in the physical
space, since this depends on which quantum theory is invoked. Our less
controversial focus here is to approach thermalization in closed systems from
operational properties testable in laboratory.

The structure of the paper is the following. Sect. II presents the many-body
generalization of the single-particle weak values, stating them as operational
properties in configuration space without explicit dependence on quantum
randomness and backaction, for both distinguishable and indistinguishable
particles. Sect. III discusses thermalization and equilibration concepts as
found in the literature for closed quantum systems, and address the eigenstate
thermalization hypothesis [38, 39]. Sect. IV defines our model system and its
nonequilibrium dynamics, and summarizes our results for both expectation
values and weak values from the Schrödinger equation dynamics. In Sect. V we
conclude.

##  II Operational properties in the configuration space

To simplify notation, along the paper we use natural units and consider a
$1$-dimensional physical space with degree of freedom $x$, so that
$\mathbf{x}=\\{x_{1},...,x_{N}\\}$ is the position in the $N$-dimensional
configuration space; the extension to a $3N$-dimensional space should be
straightforward. As already mentioned, a measured closed system is no longer a
closed system, and a strong measurement, for example, of the momentum operator
$\hat{p}$ yields the eigenvalue $p$, and produces the initial state to be
converted into the momentum eigenstate. On the other hand, expectation values
and weak values yield operational information of the system without backaction
and quantum randomness.

###  II.1 Expectation values

We define the expectation value of the momentum $\langle p(t)\rangle$ as an
operational property of the system in the sense that it is linked to a well
defined protocol in the laboratory,

$\langle p(t)\rangle=\int dp\;p\;\mathbb{P}(p,t),$ (1)

where the probability distribution $\mathbb{P}(p,t)$ can be obtained as
follows: i) many identical initial states $|\Psi(t)\rangle$ are prepared at
time $t$; ii) for each initial state, a (weak or strong) measurement of
momentum is done, yielding the value $p$ at time $t$; iii) $\mathbb{P}(p,t)$
is constructed by counting how many $p$ occurs when repeating ii) on ensemble
i).

When then applying Born law to predict the value of $\mathbb{P}(p,t)$ one can
easily identifies

$\langle p(t)\rangle=\int
dp\;p\;\mathbb{P}(p,t)=\langle\Psi(t)|\hat{p}|\Psi(t)\rangle.$ (2)

Notice that the right hand side of (2) depends on the state of the system
$|\Psi(t)\rangle$ before a measurement is done, without neither randomness nor
backaction. In fact, $\langle p(t)\rangle$ is a typical property used to
analyze when an isolated quantum system thermalizes. We are here interested in
discussing thermalization from operational properties of the isolated quantum
system requiring the measurement of both momentum and position simultaneously.
Let us then start by discussing weak values in physical space.

###  II.2 Weak values in the physical space

It has recently been shown that weak values [40] are able to yield dynamic
information on two noncommuting operators at a single time avoiding backaction
[41, 42] and quantum randomness. Weak values have attracted a lot of
theoretical [41, 42, 43, 44, 45, 46] and experimental [47, 48, 49] interests
in many research fields.

At the laboratory the single-particle weak value of momentum $p_{W}(x,t)$ is
given by

$p_{W}(x,t)=\frac{\int dp\;p\;\mathbb{P}(p,x,t)}{\int dp\;\mathbb{P}(p,x,t)},$
(3)

computed from the probability distribution $\mathbb{P}(p,x,t)$ via the
following procedure: i) many identical initial states $|\Psi(t)\rangle$ are
prepared at time $t$; ii) for each initial state, a weak measurement of the
momentum is done, yielding the value $p$ at time $t$; iii) subsequently, a
strong measurement of the position is done, yielding the value $x$ at time
$t$; iv) $\mathbb{P}(p,x,t)$ is constructed by counting how many $p$ and $x$
occurs when repeating ii) and iii) on ensemble i).

Since $x$ is post-selected in (3), but not integrated out, one gets
information on how the expectation value of the momentum is distributed in the
physical space. Again, when applying Born law to predict the value of
$\mathbb{P}(p,x,t)$ one can easily identifies

$p_{W}(x,t)=\frac{\int dp\;p\;\mathbb{P}(p,x,t)}{\int
dp\;\mathbb{P}(p,x,t)}=\text{Real}\left(\frac{\langle
x|\hat{p}|\Psi(t)\rangle}{\langle x|\Psi(t)\rangle}\right).$ (4)

Similarly to (2), the ensemble-over-identical-experiments in the right hand
side in (4) eliminates the undesired backaction and quantum randomness induced
by the first measuring apparatus [40, 42, 46, 50, 51], so that $p_{W}(x,t)$ in
(4) depends only on the initial state before the measurements took place.
Expression (4) allows us to give the weak value a very simple interpretation.
In a single-particle system, from $\mathbb{I}=\int dx|x\rangle\langle x|$, one
can rewrite the expectation value of the momentum in (2) (see Appendix A) as

$\langle p\rangle(t)=\int dx\langle\Psi(t)|x\rangle\langle
x|\hat{p}|\Psi(t)\rangle=\int dx|\Psi(x,t)|^{2}p_{W}(x,t),$ (5)

so that the same probability distribution $\mathbb{P}(p,x,t)$ used to compute
$p_{W}(x,t)$ in (4) can be employed to obtain $\mathbb{P}(p,t)$ used to
compute $\langle p(t)\rangle$ in (2), since

$\mathbb{P}(p,t)=\int dx\;\mathbb{P}(p,x,t).$ (6)

Using the mathematics (and not necessarily the ontology) of Bohmian mechanics,
one can also re-interpret $p_{W}(x,t)$ as the (operational) velocity of the
particle at position $x$ and time $t$ [43, 44, 45, 46, 52, 53, 54, 55, 56,
57],

$p_{W}(x,t)=\text{Imag}\left(\frac{1}{\Psi(x,t)}\frac{\partial\Psi(x,t)}{\partial
x}\right)=\frac{J(x,t)}{|\Psi(x,t)|^{2}},$ (7)

with $J(x,t)=\text{Imag}(\Psi^{*}(x,t)\partial\Psi(x,t)/\partial x)$ the
current density (see Appendix B).

###  II.3 Weak values in the configuration space for distinguishable particles

Notice that $p_{W}(x,t)$ is an operational property in the ordinary physical
space, but we now need to define an operational property in the configuration
space for dealing with an isolated quantum system with $N$ particles.
Therefore, in this paper, we extend the original single-particle weak values
in (4) to $N$-particle scenarios for both distinguishable and
indistinguishable cases. For the former case, the weak values for the
$j$-particle is

$p_{W}^{j}(\mathbf{x},t)=\frac{\int
dp_{j}\;p_{j}\;\mathbb{P}(p_{j},\mathbf{x},t)}{\int
dp_{j}\;\mathbb{P}(p_{j},\mathbf{x},t)},$ (8)

where now the probability distribution $\mathbb{P}(p_{j},\mathbf{x},t)$ is
obtained as follows: i) many identical initial states $|\Psi(t)\rangle$ are
prepared at time $t$; ii) for each initial state, a weak measurement of the
momentum of the $j$-particle is done, yielding the value $p_{j}$ at time $t$;
iii) subsequently, a strong measurement of the positions of particles
$1$,$2$,…,$N$ yielding respectively the values $x_{1}$,$x_{2}$,…,$x_{N}$ is
done; iv) $\mathbb{P}(p_{j},\mathbf{x},t)$ is constructed by counting how many
$p_{j}$,$x_{1}$,$x_{2}$,…,$x_{N}$ occurs when repeating ii) and iii) on
ensemble i).

Our definition of distinguishable particles above is operational in the sense
that the measuring apparatus is somehow able to distinguish particles, for
example, by measuring their masses but, to avoid unnecessary notation, we have
not indicated this extra measurement in the above protocol. Again Born law
allows us to rewrite (8) as

$p_{W}^{j}(\mathbf{x},t)=\frac{\int
dp_{j}\;p_{j}\;\mathbb{P}(p_{j},\mathbf{x},t)}{\int
dp_{j}\;\mathbb{P}(p_{j},\mathbf{x},t)}=\text{Real}\left(\frac{\langle\mathbf{x}|\hat{p}_{j}|\Psi(t)\rangle}{\langle\mathbf{x}|\Psi(t)\rangle}\right),$
(9)

where once more from the mathematics (and not necessarily from the ontology)
of Bohmian mechanics, one can also re-interpret $p_{W}^{j}(\mathbf{x},t)$ as
the (operational) velocity of the $j$-particle at the position $\mathbf{x}$ in
the configuration space and time $t$ [43, 44, 45, 46, 52, 53, 54, 55, 56, 57],

$p_{W}^{j}(\mathbf{x},t)=\frac{J^{j}(\mathbf{x},t)}{|\Psi(\mathbf{x},t)|^{2}},$
(10)

with
$J^{j}(\mathbf{x},t)=\text{Imag}(\Psi^{*}(\mathbf{x},t)\partial\Psi(\mathbf{x},t)/\partial
x_{j})$ the current density in the $x_{j}$ direction.

###  II.4 Weak values for indistinguishable particles

The most common situation in the laboratory however relates to identical
particles, for which a proper many-body wave function should implicitly
include the exchange symmetry among the particles. Equation (9) then becomes
inaccessible in a laboratory because it is no longer possible to know,
operationally, which position belongs to each particle. To deal with
indistinguishable particles, one needs to construct a many-body weak value
defined in physical space coordinate $x$ by averaging (integrating) all
degrees of freedom (see Appendix A). By doing so one obtains

$\tilde{p}_{W}(x,t)=\frac{\int dp\;p\;\mathbb{\tilde{P}}(p,x,t)}{\int
dp\;\mathbb{\tilde{P}}(p,x,t)},$ (11)

where the probability distribution $\mathbb{\tilde{P}}(p,x,t)$ is now obtained
as follows: i) many identical initial states $|\Psi(t)\rangle$ are prepared at
time $t$; ii) for each initial state, a weak measurement of the momentum of
one nonidentified particle is done, yielding the value $p$ at time $t$; iii)
subsequently, a strong measurement of the position of the same or another
nonidentified particle is done, yielding the value $x$ at time $t$; iv)
$\mathbb{\tilde{P}}(p,x,t)$ is constructed by counting how many $p$ and $x$
occurs when repeating ii) and iii) on ensemble i).

Born law again allows us to rewrite (11) as (see Appendix A)

$\tilde{p}_{W}(x,t)=\frac{\int dp\;p\;\mathbb{\tilde{P}}(p,x,t)}{\int
dp\;\mathbb{\tilde{P}}(p,x,t)}=\frac{1}{N^{2}}\sum_{j=1}^{N}\sum_{k=1}^{N}p_{W}^{j,k}(x,t),$
(12)

with

$p_{W}^{j,k}(x,t)=\frac{\int dx_{1}...\int dx_{k-1}\int dx_{k+1}...\int
dx_{N}\;p_{W}^{j}(..,x_{k-1},x,x_{k+1},..,t)|\Psi(..,x_{k-1},x,x_{k+1},..,t)|^{2}}{\int
dx_{1}...\int dx_{k-1}\int dx_{k+1}...\int
dx_{N}|\Psi(..,x_{k-1},x,x_{k+1},..,t)|^{2}}.$ (13)

Notice that $\tilde{p}_{W}(x,t)$ is, in fact, the local velocity as used in
quantum hydrodynamic models [52, 53, 54, 55], being empirically accessible in
both distinguishable and indistinguishable scenarios. The operational protocol
for computing $\mathbb{\tilde{P}}(p,x,t)$ for indistinguishable particles is
related to $\mathbb{P}(p_{k},\mathbf{x},t)$ for distinguishable particles as

$\mathbb{\tilde{P}}(p,x,t)=\frac{1}{N^{2}}\sum_{j=1}^{N}\sum_{k=1}^{N}\int
dx_{1}...\int dx_{k-1}\int dx_{k+1}...\int
dx_{N}\;\mathbb{P}(p_{j},..,x_{k-1},x,x_{k+1},..,t).$ (14)

##  III Thermalization in isolated systems from expectation values

Our main contribution in the study of quantum thermalization, as detailed in
next section, is the inclusion of many-body weak values in the configuration
space as operational properties. However, such a study has usually been done
in the literature in terms of expectation values as in (2), and as such we
summarize in this section the role of expectation values to characterize
thermalization.

For an initial nonequilibrium pure state $|\Psi(0)\rangle$, the Schrödinger
equation provides its unitary evolution as
$|\Psi(t)\rangle=\sum_{n}c_{n}e^{-iE_{n}t}|n\rangle$, where $|n\rangle$ is an
energy eigenstate with eigenvalue $E_{n}$, and $c_{n}=\langle
n|\Psi(0)\rangle$ is defined by the initial conditions. The expectation value
of some observable $\hat{A}$ is given by

$\langle A\rangle(t)=\sum_{n}\rho_{n,n}A_{n,n}+\sum_{n,m\neq
n}\rho_{n,m}(t)A_{m,n},$ (15)

with the time-dependent off-diagonal elements of the density matrix in the
energy representation defined as

$\rho_{n,m}(t)=c_{m}^{*}c_{n}e^{i(E_{m}-E_{n})t},$ (16)

and the time-independent diagonal elements as

$\rho_{n,n}(t)=c_{n}^{*}c_{n}=\rho_{n,n}(0),$ (17)

with the operator $\hat{A}$ in the energy representation being

$A_{m,n}=\langle m|\hat{A}|n\rangle.$ (18)

A system is said to equilibrate if, after some time $t_{eq}$ enough for full
dephasing between different energy eigenstates, the off-diagonal terms
(coherences) cancel out so that (15) can be computed solely from the time-
independent diagonal terms (populations), that is, $\langle
A\rangle(t)\approx\sum_{n}\rho_{n,n}A_{n,n}$ for most times $t>t_{eq}$ (except
for some recurrence times). The properties of the system after equilibration
are fully determined by the initial conditions $\rho_{n,n}(0)=|c_{n}|^{2}$,
since the density matrix populations are time-independent. The closed system
is then said to thermalize when $\langle A\rangle(t)$ becomes roughly equal to
the expectation value as computed from the classical density matrix in the
microcanonical ensemble, $\rho_{cl}$, in which equal probabilities are
attached to each microstate within an energy window defined by the initial
conditions; that is, a subset $N_{act}$ of the energy eigenstates are
initially activated at $t=0$, and they will remain as the only states
dictating the system dynamics at any $t$. It is important to notice that
$N_{act}$ refers to a given number of relevant elements in a basis set, but it
has no relation with the number of particles $N$ of the system. In other
words, one can also expect thermalization even in few-particle systems, as
detailed in next section; in fact, thermalization has also been studied in
laboratories in small systems with as little as 6 [27], 5 [14], or 2-4 bosons
[58, 59], 3 qubits [60], and even single-particle systems [61, 62, 63].

The eigenstate thermalization hypothesis (ETH) [38, 39] has become the
standard theory dealing with quantum thermalization in closed systems. It
states that the dephasing above mentioned is typical to nondegenerate and
chaotic many-body nonintegrable systems, where the off-diagonals $A_{m,n}$ in
(18) become exponentially smaller than $A_{n,n}$. In recent years a large
amount of numerical experiments has successfully tested such a hypothesis by
directly diagonalizing some sort of short range many-body lattice Hamiltonian,
like Fermi- or Bose-Hubbard [26, 31, 32, 64, 65, 66], and XXZ- or XYZ-
Heisenberg [14, 23, 27, 67, 68, 69, 70], in the search of chaotic signatures
in the statistics of their spectra, as in general induced by local impurities,
without the need to explicitly evolve $|\Psi(0)\rangle$. The ETH states that
nonintegrable systems (where total energy may be the only conserved quantity),
after a quench (which may create a nonequilibrium initial state by activating
a subset $N_{act}$ of excited eigenstates), can present a ‘chaotic’ spectrum
ruled by the Wigner-Dyson statistics (which contains level repulsion), and
that the long time average of the expectation value of some observable roughly
equals its thermal equilibrium value in an (microcanonical) ensemble. Such a
hypothesis claims that thermalization is indeed hidden in the chaotic initial
nature of the Hamiltonian eigenstates themselves. Lattice models typically
handle $\approx 24$ sites with $\approx 1/3$ filling, and the above local
impurities are added to break their otherwise integrable character, as the ETH
overall claims that integrable systems are not expected to thermalize.

On the other hand, our time evolution dwells in true configuration space with
an antisymmetrized wave function and full long range electron-electron
interaction; due to the tensorial nature of our problem, and since we need to
employ a grid with $M\approx 10^{3}$ points per degree of freedom for decent
position and momentum resolutions, we can realistically only deal with
$N\lesssim 4$ particles ($N=3$ already implies $M\approx 10^{9}$ grid points
at each time step).

##  IV Numerical results for expectation values and weak values

We now apply the many-body weak values machinery for the analysis of quantum
thermalization in a model with $N$ spinless electrons [20, 71, 72, 73, 74, 75,
76, 77], typical of condensates in harmonic oscillator traps under a speckle
field. Such a field translates to a ‘chaotic’ random disorder potential, which
can yield a chaotic energy spectrum as requested by the ETH and so induce
thermalization even in systems with small $N$. Our model can attach a disorder
to each grid point, and an initial velocity at $t=0$ is given to the electrons
as to simulate the initial quench of the confining potential, for each of the
considered $N=1,2,3$ systems. In fact, most of the thermalization literature
dealing with identical particles employs some lattice-based model, since such
models avoid the need of explicit knowledge of the exact position of each
particle in a point of the configuration space; instead, they only need to
know how many particles are present in each site of the physical space. From a
computational point of view lattice models have unquestionable advantages but
are inappropriate for our goal of discussing whether or not thermalization
occurs in configuration space, a goal that forces us to directly solve the
time-evolution of the few-body Schrödinger equation in configuration space,
and to analyze thermalization by monitoring the time-evolution of both
expectation values and weak values.

###  IV.1 Initial state

The pure initial $N$-electron nonequilibrium antisymmetric state is

$\langle\mathbf{x}|\Psi(0)\rangle=\frac{1}{\mathcal{C}}\sum_{n=1}^{N!}\text{sign}(\vec{p}(n))\prod_{j=1}^{N}\psi_{j}(x_{p(n)_{j}},0),$
(19)

with $\mathcal{C}$ a normalization constant and $\text{sign}(\vec{p}(n))$ the
sign of the permutation $\vec{p}(n)=\\{p(n)_{1},..,P(n)_{N}\\}$. Each initial
Gaussian state in (19) is

$\psi_{j}(x)=\exp\left[-\frac{(x-x_{0j})^{2}}{2\sigma_{j}^{2}}\right]\exp\left[ip_{0j}(x-x_{0j})\right],$
(20)

with spatial dispersion $\sigma_{j}$, central position $x_{0j}$, and central
velocity $p_{0j}$. The dynamical evolution of $|\Psi(t)\rangle$ is determined
by the Schrödinger equation, $i\partial|\Psi(t)\rangle/\partial
t=\hat{H}|\Psi(t)\rangle$, where the Hamiltonian $\hat{H}$ is described in
next section.

###  IV.2 Full Hamiltonian

The $N$-electron Hamiltonian of our model system is
$\hat{H}=\hat{H}_{0}+\hat{D}$, with

$\hat{H}_{0}=\sum_{j=1}^{N}\left[\hat{k}_{j}+\hat{v}_{j}+\sum_{k<j}^{N}\hat{e}_{k,j}\right],\;\;\;\;\;\;\hat{D}=\sum_{j=1}^{N}\hat{d}_{j}.$
(21)

In $\hat{H}_{0}$, $\langle
x_{j}|\hat{e}_{k,j}|x_{k}\rangle=1/\sqrt{(x_{j}-x_{k})^{2}+\alpha^{2}}$ takes
care of the Coulomb repulsion with a smooth parameter $\alpha$, $\langle
x_{j}|\hat{k}_{j}|x_{j}\rangle=-\partial^{2}/(2\partial x_{j}^{2})$ stands for
the kinetic energy, and $\langle
x_{j}|\hat{v}_{j}|x_{j}\rangle=\omega^{2}x_{j}^{2}/2$ is the harmonic trap
potential. On the other hand, $\hat{D}$ introduces random disorder at every
grid point, where $\langle
x_{j}|\hat{d}_{j}|x_{j}\rangle=\gamma_{D}\sum_{k=1}^{M}a_{k}\exp[-4(x_{j}-g_{k})^{2}/\sigma_{D}^{2}]$,
with $\gamma_{D}$ its strength and $\sigma_{D}$ its spatial dispersion, while
$g_{k}$ runs through $M$ grid points; the set of random numbers $a_{k}$
satisfies $\langle a_{k}\rangle=0$ and $\langle a_{k}^{2}\rangle=1$, and the
disorder potential is normalized via $\int\langle
x_{j}|\hat{d}_{j}|x_{j}\rangle^{2}dx_{j}=\gamma_{D}^{2}$. Such random disorder
can be mapped onto speckle field potentials typical in some optical lattice
experiments. All simulation parameters are found in [78].

Figure 1: (a): Zoom of a disordered harmonic potential for $N=1$, $\langle
x_{1}|\hat{v}_{1}+\hat{d}_{1}|x_{1}\rangle$, from a larger simulation box.
(b): Successive energy splittings $\Delta E=E_{n}-E_{n-1}$ from a grid
diagonalization of the respective Hamiltonian $\hat{H}$. (d): zoom around the
peak of the corresponding density matrix modulus $|\rho_{n,m}(t=0)|$, which
remains the same at any time $t$ in an unitary evolution. (c): initial
positions $x_{0j}$ for $N=1,2,3$, and initial velocity $p_{0j}$ which is the
same in every case.

Figure 1(a) exemplifies, for $N=1$, a typical shape of the disordered harmonic
potential, $\langle x_{1}|\hat{v}_{1}+\hat{d}_{1}|x_{1}\rangle$, while figure
1(b) shows the corresponding successive energy splittings $\Delta
E=E_{n}-E_{n-1}$, as obtained from a direct diagonalization of $\hat{H}$,
which oscillate around the pure harmonic oscillator value of $1$. Figure 1(d)
shows the related shape of the density matrix modulus $|\rho_{n,m}(t=0)|$,
from where one identifies $N_{act}\approx 70$ (counting every level above
$10\%$ of peak value), while figure 1(c) shows initial positions $x_{0j}$,
$j=1...N$, and the initial velocity $p_{0}$, the same for any $j$ and $N$, for
the $N=1,2,3$ systems. Notice that $p_{0}$ not only takes the role of
simulating the initial quench of the trap and so to initiate the
nonequilibrium dynamics, but it is also responsible for determining the size
of the energy window and so the value of $N_{act}$. The initial energy $E_{0}$
of the wave packet defines, in the language of the ETH, the center of the
microcanonical energy window; since $E_{0}\approx\omega+p_{0}^{2}/2=201$, the
peak is at $n=m=201$. As mentioned in (15), the diagonal terms are time-
independent and, although real and imaginary parts of the off-diagonals terms
oscillate in time, their modulus remain constant so that, in an unitary
evolution, the modulus seen in figure 1(d) remains the same at any $t$ (see
Appendix B).

Both initial wave packet and Hamiltonian in (19)-(21) have many adjustable
parameters in [78] upon which the value of $t_{eq}$ depends on: (i) Coulomb
correlation by varying $\omega$ or $\alpha$; (ii) initial Gaussian by varying
$\sigma_{j}$, $x_{0j}$, or $p_{0j}$; (iii) disorder potential by varying
$\gamma_{D}$ or $\sigma_{D}$. All of them play a role in determining the
$\Delta E$ splittings of the involved $N_{act}$ activated eigenstates, which
is what drives the nonequilibrium dynamics of both expectation values and weak
values. It is not our goal to fully characterize the equilibration process as
a function of all those parameters. Neither to address the topic of many-body
localization [79, 80, 81], which should work against thermalization, nor to go
deeper in the issue of quantum-to-classical transition at $t\gg t_{eq}$; these
two latter issues are beyond the scope of our work and are focus of extensive
studies elsewhere. The main goal of our paper is to present a distinct
perspective in the understanding of quantum thermalization by looking at the
many-body weak values of the momentum in the configuration space.

Each numerical experiment corresponds to a static realization of a disorder
pattern; $t_{eq}$ hardly changes among different runs, given that all other
model parameters remain unchanged. The disorder amplitude should not be too
strong, to avoid localization due to all small wells created on the top of the
trap, neither too weak, to avoid too long simulation times until reaching
$t_{eq}$. We emphasize that even in the absence of electron-electron
collision, disorder collision is able to create correlation among distinct
degrees of freedom in configuration space when $N>1$.

###  IV.3 Time evolution

Figure 2: Wave packet dynamics. Upper panels for $N=3$: initial
$|\Psi(x_{1},x_{2},[x_{3}],t=0)|^{2}$ (a) and final
$|\Psi(x_{1},x_{2},[x_{3}],t=150)|^{2}$ (b) shapes. Middle panels for $N=2$:
initial $|\Psi(x_{1},x_{2},t=0)|^{2}$ (c) and final
$|\Psi(x_{1},x_{2},t=150)|^{2}$ (d) shapes. Lower panels: (e) for $N=1$:
initial (dotted) $|\Psi(x_{1},t=0)|^{2}$ and final (solid)
$|\Psi(x_{1},t=150)|^{2}$ shapes; (f): $1$D-view for the three systems:
$|\Psi(x_{1},[x_{2}],[x_{3}],t)|^{2}$ for $N=3$ (red),
$|\Psi(x_{1},[x_{2}],t)|^{2}$ for $N=2$ (blue), and $|\Psi(x_{1},t)|^{2}$ for
$N=1$ (green), with solid (dotted) lines for $t=150$ ($t=0$); also shown the
classical microcanonical harmonic oscillator distribution $\rho_{cl}(x_{1})$
(black dashed line). All plots are in log-scale, and the horizontal axis in
(a),(c) ((b),(d)) is the same as in (e) ((f)).

Figure 2 plots in configuration space the time evolution of the $N$-electron
wave function for the three $N=1,2,3$ systems from the initial nonequilibrium
state. We show $|\Psi(x_{1},x_{2},[x_{3}],t)|^{2}$ for $N=3$ and
$|\Psi(x_{1},x_{2},t)|^{2}$ for $N=2$ at initial ($t=0$, panels (a),(c)) and
final ($t=150$, panels (b),(d)) simulation times, while panel (e) shows
$|\Psi(x_{1},t)|^{2}$ for $N=1$ for both initial ($t=0$, dotted) and final
($t=150$, solid) simulation times; the notation $[x_{j}]$ means that the
degree $j$ is integrated out, with results independent on chosen $j$ due to
the antisymmetry of the problem. The initially localized wave packets fully
spread out due to random scattering generated by both disorder potential and
Coulomb repulsion, with such spreading becoming more homogeneous as $N$
increases. In panel (f) we show the respective $1$D plots of
$|\Psi(x_{1},[x_{2}],[x_{3}],t)|^{2}$ ($N=3$, red),
$|\Psi(x_{1},[x_{2}],t)|^{2}$ ($N=2$, blue), and $|\Psi(x_{1},t)|^{2}$ ($N=1$,
green), with initial (final) results at dotted (solid) lines; the
probabilities at large $t$, as one expects to have crossed $t_{eq}$, approach
the microcanonical distribution of a classical harmonic oscillator (dashed
black line), $\rho_{cl}(x)=[\pi l_{0}\sqrt{p_{0}^{2}-x^{2}/l_{0}^{2}}]^{-1}$,
so that the dynamics develops in between the classical turning points at
$x_{TP}=\pm 20$ ($=p_{TP}$ since $\omega=1/l_{0}^{2}=1$); the deviation from
$\rho_{cl}(x)$ decreases as $N$ increases. In Appendix C we analyze the same
nonequilibrium dynamics in momentum representation.

###  IV.4 Thermalized expectation values

Figures 3 and 4 show time evolution and thermalization of some typical
expectation values $\langle
A\rangle(t)=\langle\Psi(t)|\hat{A}|\Psi(t)\rangle$, for $N=1,2,3$ respectively
in panels (a), (b), (c). On one hand, figure 3 focus on energy terms
normalized by $N$: kinetic $\langle K\rangle(t)$, with
$\hat{K}=\sum_{j}\hat{k}_{j}$, potential $\langle V\rangle(t)=\langle
V_{HO}\rangle(t)+\langle V_{D}\rangle(t)+\langle V_{Cou}\rangle(t)$, with
$\hat{V}_{HO}=\sum_{j}\hat{v}_{j}$, $\hat{V}_{D}=\sum_{j}\hat{d}_{j}$, and
$\hat{V}_{Cou}=\sum_{j,k<j}\hat{e}_{k,j}$, and half of total energy $\langle
E\rangle(t)=\langle K\rangle(t)+\langle V\rangle(t)$. On other hand, figure 4
focus on position $\langle x_{j}\rangle(t)$ and momentum $\langle
p_{j}\rangle(t)$ for the $j$-electron, and on their RMS values
$z_{j,RMS}(t)=\sqrt{\langle z^{2}_{j}\rangle(t)-\langle z_{j}\rangle^{2}(t)}$,
with $z=x,p$, and results independing on chosen $j$. _Without_ random
disorder, such expectation values would only exhibit harmonic oscillations
with period $2\pi/\omega$: while $\langle x_{j}\rangle(t)$ and $\langle
p_{j}\rangle(t)$ would respectively oscillate within position $\pm x_{TP}=\pm
20$ and momentum $\pm p_{TP}=\pm 20$ turning points, $\langle V\rangle(t)/N$
and $\langle K\rangle(t)/N$ would respectively oscillate within $0$ and
$x_{TP}^{2}/2=200$ and within 0 and $p_{TP}^{2}/2=200$; these latter values
increase a little upon $N$ due to Coulomb repulsion, whose isolated
contribution is shown ($100\times$-magnified) in figure 3.

Figure 3: Energy expectation values from the dynamics in figure 2. Panels (a),
(b), (c) respectively for $N=1,2,3$. Kinetic $\langle K\rangle(t)$, potential
$\langle V\rangle(t)$, half of total energy $\langle E\rangle(t)=\langle
K\rangle(t)+\langle V\rangle(t)$, and isolated contribution of $\langle
V_{Cou}\rangle(t)$ ($100\times$-magnified) are all normalized by $N$. Inset
for $N=2$ shows a longer propagation time, $t=[150-300]$. Legend in (a) and
horizontal axis in (c) apply to all panels.

The _presence_ of random disorder, even though $\langle V_{D}\rangle(t)\approx
0$ at any $t$, brings the initial nonequilibrium state into a final
_equilibrium_ state after a relaxation time $t_{eq}$. We know from the
discussions in (15) and in figure 1(d) that thermalization of an observable
$\hat{A}$ is determined by the diagonal populations of the density matrix,
while its off-diagonal coherences should dephase and only yield small
fluctuations around the relaxed value (see Appendix B). So one may estimate
the value of $t_{eq}$ either from figure 3, when the virial theorem $\langle
K\rangle\approx\langle V\rangle\approx\langle E\rangle/2$ is roughly satisfied
(since $\langle V_{Cou}\rangle\ll\langle V_{HO}\rangle$) [74], or from figure
4, when $\langle p_{j}\rangle\approx\langle x_{j}\rangle\approx 0$ seemingly
indicating a frozen dynamics after thermalization. From this latter result the
RMS values become $z_{j,RMS}(t)\approx\sqrt{\langle z^{2}_{j}\rangle(t)}$,
becoming also constant in figure 4 at $t>t_{eq}$ (from $\approx 14.2$ for
$N=1$ to $\approx 14.4$ for $N=3$), so that the values of
$p^{2}_{j,RMS}/2=\langle p^{2}_{j}\rangle/2$ and $x^{2}_{j,RMS}/2=\langle
x^{2}_{j}\rangle/2$ roughly yield the respective values of $\langle
K\rangle/N$ and $\langle V\rangle/N$ at $t>t_{eq}$ in figure 3. As expected
for an unitary evolution, $\langle E\rangle(t)/N$ remains conserved at any
$t$, from $\approx 201$ for $N=1$ to $\approx 207$ for $N=3$.

Figure 4: Position $\langle x_{j}\rangle(t)$ and momentum $\langle
p_{j}\rangle(t)$ expectation values from the dynamics in figure 2. Panels (a),
(b), (c) respectively for $N=1,2,3$. Their respective RMS values, as defined
in the text, are also shown, where results do not depend on chosen $j$. Inset
for $N=2$ shows a longer propagation time, $t=[150-300]$. Legend in (a) and
horizontal axis in (c) apply to all panels.

The value of $t_{eq}$ depends on all parameters [78] in (19)-(21), e.g., the
smaller is $p_{0j}$ or the higher is $\gamma_{D}$ the smaller is $t_{eq}$. By
increasing the influence of $\langle V_{Cou}\rangle$ in comparison to $\langle
K\rangle$ (by decreasing $\omega$ or $\alpha$), $t_{eq}$ increases since the
oscillation period and so $x_{TP}$ increases. The plots of $\langle
V_{Cou}\rangle(t)$ in figure 3 show that the Coulomb repulsion is more
effective at the turning points for $t\ll t_{eq}$, where electrons spend more
time reversing their movements, while the disorder potential overall acts
through a whole oscillation, although one may take it as more effective at the
origin. Coulomb correlation has a striking influence on the thermalization
process: in configuration space the only scattering mechanism for $N=1$ is due
to disorder, while for $N>1$ Coulomb scattering also makes more difficult for
the system to relax. This is seen by the slightly increasing values of
$t_{eq}$ as one moves in figure 3 from (a) ($t_{eq}\approx 70$) to (b)
($t_{eq}\approx 80$) to (c) ($t_{eq}\approx 90$). The vanishing of $\langle
p\rangle(t)$ in figure 4 seems more effective as one moves from (a) to (b) to
(c) but, however, it does not necessarily imply that electrons have achieved a
stationary-state null velocity at $t\gg t_{eq}$, as our following analysis on
weak values of the momentum will clarify (see also ‘phase-space’ in Appendix
C).

###  IV.5 Nonthermalized weak values

We can at last elaborate on how the many-body weak values of the momentum may
improve our understanding on thermalization. In table 1 we summarize the five
types of operational properties accessible for the three different types of
quantum systems considered in this section: single-particle, distinguishable
particles, indistinguishable particles.

$N=1$ |  $\langle p(t)\rangle$ (2) |  $p_{W}(x,t)$ (4) |  |  |
---|---|---|---|---|---
$N>1$ (Dis) |  $\langle p(t)\rangle$ (2) |  |  $p_{W}^{j}(\mathbf{x},t)$ (9) |  $\tilde{p}_{W}(x,t)$ (12) |  $p_{W}^{j,k}(x,t)$ (13)
$N>1$ (Ind) |  $\langle p(t)\rangle$ (2) |  |  |  $\tilde{p}_{W}(x,t)$ (12) |
Table 1: Five operational properties accessible in the laboratory for each of
the three quantum systems considered in our work: (i) with $N=1$ particles,
(ii) with $N>1$ distinguishable particles, and (iii) with $N>1$
indistinguishable particles.

The plot in figure 5(a) corresponds to the quantum system $N=1$ in table 1.
Although the expectation value $\langle p_{j}\rangle(t)$ seems to indicate
that the quantum behavior at $t\gg t_{eq}$ roughly equals the behavior of a
diagonal density matrix in (17), the weak values $p_{W}^{1,1}(0,t)$ (which
obviously corresponds to $p_{W}(0,t)$ in (4)) certifies that the off-diagonal
terms in (16) do not vanish after thermalization. For $N=1$, the fact that
expectation values thermalize is just a result that positive and negative off-
diagonals elements can roughly compensate each other, but they certainly do
not disappear as indicated by $p_{W}^{1,1}(x,t)$. We remind that $\langle
p\rangle(t)=\int dx\int dp\;p\;\mathbb{P}(p,x,t)$ and $p_{W}(x,t)$ can, both,
be computed from the same empirical probability $\mathbb{P}(p,x,t)$. In other
words, $\mathbb{P}(p,x,t)$ provides simultaneous thermalized and
nonthermalized results depending on how it is treated.

Figure 5: Local-in-position many-body weak values of the momentum from the
dynamics in figure 2 for $N=1$ in (a), and for distinguishable $N=2$ particles
in (b),(c) where Coulomb/exchange terms are disconsidered. Panels (a), (b)
show $p_{W}^{j,j}(x,t)$ from (13), which does not depend on $j$. Panel (c)
shows both $p_{W}^{j,k}(x,t)$ from (13) and $\tilde{p}_{W}(x,t)$ from (12).
Values of $x$ are the respective initial $x_{0j}$ values. The expectation
value of the momentum $\langle p_{j}\rangle(t)$ is also shown. Horizontal axis
in (a),(b) the same as in (c).

The plots in figures 5(b) and 5(c) correspond to the quantum system
$N>1(\text{Dis})$ in table 1, because neither exchange nor Coulomb interaction
among the particles are included; that is, the many-body wave function here
could be written as $\Psi(\mathbf{x},t)=\psi_{1}(x_{1},t)\psi_{2}(x_{2},t)$.
The weak values $p_{W}^{1,1}(0,t)$ and $p_{W}^{1,1}(-4,t)$ in panel (b)
confirm the nonthermalized operational properties even at $t\gg t_{eq}$. In
panel (c), one notices that $p_{W}^{1,2}(0,t)$, which corresponds to the weak
measurement of the momentum of particle 1 and strong measurement of the
position of particle 2, shows a thermalized behavior as it overlaps $\langle
p_{j}\rangle(t)$; but this is because here, as the particles have no
correlations among them, one gets $p_{W}^{j,k}(x,t)=\langle p_{j}\rangle(t)$
when $j\neq k$, as discussed in (32). Also in panel (c) $\tilde{p}_{W}(0,t)$,
from (12), which here is just a particle-average of the oscillating term
$p_{W}^{1,1}(0,t)$ and the non-oscillating term $p_{W}^{1,2}(0,t)$, presents
less oscillations but still clearly showing a non-thermalized behavior of
these operational properties.

Figure 6: Local-in-position many-body weak values of the momentum from the
dynamics in figure 2 for indistinguishable particles with $N=2$ in (a),(b) and
$N=3$ in (c). Panels (a), (c) show $p_{W}^{j,j}(x,t)$ from (13), which does
not depend on $j$, for the respective initial $x_{0j}$ values. Panel (b) shows
both $p_{W}^{j,k}(x,t)$ from (13) and $\tilde{p}_{W}(x,t)$ from (12). The
expectation value of the momentum $\langle p_{j}\rangle(t)$ is also shown.
Inset in (a) shows a longer propagation time, $t=[150-300]$. Horizontal axis
in (a),(b) the same as in (c).

Figure 6 corresponds to the (most common) quantum system $N>1(\text{Ind})$ in
table 1, in which one only has operational access in the laboratory to the
expectation value $\langle p\rangle(t)$ in (2) and to the weak value
$\tilde{p}_{W}(x,t)$ in (12). Panels (a) and (b) for $N=2$; panel (c) for
$N=3$. Since we also have mathematical (not operational) access in our
simulations (directly from configuration space) to the weak value
$p_{W}^{j,k}(x,t)$, we have also plotted it to help us to understand the
behavior of the operational weak value $\tilde{p}_{W}(x,t)$, which is a
particle average over different $p_{W}^{j,k}(x,t)$. We see in Fig. 6(a) that
$p_{W}^{1,1}(x,t)$ has smaller oscillations than in Fig. 5(b), but it still
does not thermalize, while $\langle p_{j}\rangle(t)$ effectivelly thermalize
(this later result is independent on $j$); the larger time window in the inset
so confirms. In figure 6(b) one notices that both $p_{W}^{1,2}(0,t)$ and
$p_{W}^{1,1}(0,t)$ have a similar non-thermalizing behavior. From these two
latter values we obtain
$\tilde{p}_{W}(0,t)=(p_{W}^{1,1}(0,t)+p_{W}^{1,2}(0,t))/2$ (thanks to the
properties in (31)), whose plot in panel (b) shows that this operational
parameter for identical particles does not fully thermalize. The $N=3$ case in
panel (c) starts to show the trend that, at higher $N$, $p_{W}^{1,1}(x,t)$
will approach the expectation value $\langle p_{1}\rangle(t)$ and so will also
thermalize, which is seem for all used $x_{0j}$ values. This is understood
from the fact that, for identical particles, $p_{W}^{j,k}(x,t)$ (and
$p_{W}^{j,j}(x,t))$ contains $N-1$ spatial integrals, while $\langle
p_{j}\rangle(t)$ contains $N$ and so, as one increases $N$, one gets
$p_{W}^{j,k}(x,t)\approx\langle p_{j}\rangle(t)$ because $N-1\approx N$.

The fact that the weak value of the momentum does not thermalize in the
configuration space can be understood when $p_{W}^{j}(\mathbf{x},t)$ in (9) is
mathematically interpreted as a Bohmian velocity, satisfying all of its
mathematical properties without any ontologic implications. Without external
perturbation, the initial state $|\Psi(0)\rangle$ of a closed system cannot
change with time its condition of being or not an energy eigenstate [36]. In
other words, the only energy eigenstates are the ones that are so at all
times. We know that for a closed system with $\Psi(\mathbf{x},t)$ vanishing at
the boundaries, the only Bohmian velocities in (9) that are zero are those
linked to an energy eigenstate [56]. From (10) we see that the weak value of
the momentum depends on the current density, and a time-dependent
$J^{j}(\mathbf{x},t)$ (strictly different from zero) is required for a time-
dependent probability presence in such a space, due to the continuity equation
in the configuration space implicit in the Schrödinger equation. Thus, since
our initial nonequilibrium state is not an energy eigenstate, from all
previous arguments, we conclude that $p_{W}^{j}(\mathbf{x},t)$ will never
vanish in the configuration space no matter the value of $N$, independently on
the thermalization or not of the expectation value of the momentum $\langle
p_{j}\rangle(t)$; that is, $p_{W}^{j}(\mathbf{x},t)$ will never thermalize
when evaluated at a point $\mathbf{x}$ of the configuration space (the same
conclusions apply to the presence probability in such space). Most of the
developments done in this paper come from the fact that, in most cases, the
weak value $p_{W}^{j}(\mathbf{x},t)$ in the configuration space is not
empirically accessible in the laboratory, and as such it is not an operational
property for indistinguishable particles; in such cases, though, the weak
value $\tilde{p}_{W}(x,t)$ in the physical space tends to thermalize as $N$
increases.

##  V Conclusions

We have seen how, from the empirical knowledge of a unique distribution
probability $\mathbb{P}(p,\textbf{x},t)$, it is possible to get,
simultaneously, both thermalized (expectation values) and nonthermalized (weak
values) operational properties of a closed quantum system described by the
Schrödinger equation. A possible origin of the arrow of time in such systems
comes from dealing with operational properties obtained by averaging (tracing
out) some of the degrees of freedom defined in configuration space. As such
there is no contradiction in that there are properties defined in the
configuration space that do not thermalize, while some properties defined in
the physical space (by averaging or tracing out some degrees of freedom) have
their own time irreversible equations of motion. Such conclusion can be tested
in the laboratory through the many-body weak values of the momentum for
distinguishable particles where, at least conceptually, it is possible to get
information of the position of all particles in the laboratory after their
strong measurement.

However, for indistinguishable particles, a position measurement at $x$ cannot
be linked to a specific particle since there is no experimental protocol to
tag identical particles and, as such, instead of $\mathbb{P}(p,\textbf{x},t)$
in configuration space, one has operational access to $\mathbb{P}(p,x,t)$ in
physical space. But we have shown that $\mathbb{P}(p,x,t)$ can be understood
as an averaging of $\mathbb{P}(p,\textbf{x},t)$ over degrees of freedom in the
configuration space. For the simpler single-particle case, where obviously the
distinction between distinguishable/indistinguishable particles and between
configuration/physical spaces makes no sense, we have also seen that
expectation values can thermalize while weak values may not. The simple
explanation is that the expectation value has one integration over the single
valid coordinate while the weak value has not. Even for a system with $N=2$
identical particles, we see different behaviors for the expectation value
$\langle p\rangle(t)$ and the weak value $\tilde{p}_{W}(x,t)$ in the physical
space. In general, for identical particles, $\tilde{p}_{W}(x,t)$ contains
$N-1$ spatial integrals, while $\langle p_{j}\rangle(t)$ contains $N$ and so,
as one increases $N$, one gets $p_{W}^{j,k}(x,t)\approx\langle
p_{j}\rangle(t)$ because $N-1\approx N$; that is, the former thermalizes when
the latter also does. But such thermalization is seen in the physical space,
not in the configuration space.

So does thermalization occur in configuration space? No thermalization occurs
in the most common operational properties of a quantum system when evaluated
in a point $\mathbf{x}$ of the configuration space, even when other properties
defined in simpler spaces are thermalized. Notice that we have avoided along
the paper to discuss ontologic properties. Instead, our operational approach
has the advantage that the conclusions do not depend on which quantum theory
is invoked, but it also forbids the discussion whether or not thermalization
is a real (ontologic) phenomenon occurring in physical space. This latter
discussion depends on the ontology of each quantum theory; there are quantum
theories where the configuration space is not the fundamental space.

In summary, as mentioned in our Intro [1, 2, 3, 4, 5, 6, 7, 8], there are many
arrows of time which could require different explanations. We here only
discuss the origin of the arrow of time in the non-relativistic many-body
Schrödinger equation in closed quantum systems: a non-thermalized property in
the configuration space, when some of its degrees of freedom are averaged,
leads to a thermalized property in the physical space. We have also developed
the operational many-body weak values of the momentum to make such a
conclusion testable in the laboratory.

###### Acknowledgements.

This research was funded by Spain’s Ministerio de Ciencia, Innovación y
Universidades under Grant No. RTI2018-097876-B-C21 (MCIU/AEI/FEDER, UE), Grant
PID2021-127840NB-I00 (MICINN/AEI/FEDER, UE), the “Generalitat de Catalunya”
and FEDER for the project 001-P-001644 (QUANTUMCAT), the European Union’s
Horizon 2020 research and innovation programme under Grant No. 881603
GrapheneCore3 and under the Marie Skłodowska-Curie Grant No. 765426 TeraApps.

##  Appendix A Weak values equations in many-body systems

In this appendix we develop the main weak values equations in the paper.

###  A.1 Development of equations (4) and (9) in the paper

For a single-particle system described by the wave function $\psi(x,t)=\langle
x|\psi(t)\rangle$, the expression for the weak values of the momentum,
$p_{W}(x,t)$, can be obtained from the position distribution of the mean
momentum $\langle p\rangle(t)$ of the single-particle operator,
$\hat{p}=p|p\rangle\langle p|$, as

$\langle p\rangle(t)=\langle\psi(t)|\left(\int dx|x\rangle\langle
x\right)|\hat{p}|\psi(t)\rangle=\int dx|\psi(x,t)|^{2}\frac{\langle
x|\hat{p}|\psi(t)\rangle}{\langle x|\psi(t)\rangle}=\int
dx|\psi(x,t)|^{2}p_{W}(x,t),$ (22)

where we have defined the weak values as
$p_{W}(x,t)=\text{Real}\left(\frac{\langle x|\hat{p}|\psi(t)\rangle}{\langle
x|\psi(t)\rangle}\right)$. Since $\langle p\rangle(t)$ is real, as $\hat{p}$
is an hermitian operator, we have $\int
dx|\psi(x,t)|^{2}\text{Imag}\left(\frac{\langle
x|\hat{p}|\psi(t)\rangle}{\langle x|\psi(t)\rangle}\right)=0$. Thus, only the
real part is considered for defining the weak values and so we reproduce
equation (4) in the paper.

From a similar development we can find the weak values for an $N$-particle
system, with $\Psi(\mathbf{x},t)=\langle\mathbf{x}|\psi(t)\rangle$. The mean
momentum $\langle p_{j}\rangle(t)$ of degree of freedom $j$ belonging to
operator
$\hat{P}_{j}\equiv\hat{1}\otimes...\otimes\hat{p}_{j}\otimes...\otimes\hat{1}$
with $\hat{p}_{j}=p_{j}|p_{j}\rangle\langle p_{j}|$ is

$\displaystyle\langle p_{j}\rangle(t)$ $\displaystyle=$
$\displaystyle\langle\Psi(t)|\hat{P}_{j}|\Psi(t)\rangle=\int
d\mathbf{x}\langle\Psi(t)|\mathbf{x}\rangle\langle\mathbf{x}|\hat{P}_{j}|\Psi(t)\rangle=\int
d\mathbf{x}\langle\Psi(t)|\mathbf{x}\rangle\langle
x_{1}|\otimes...\otimes\left(\langle
x_{j}|\hat{p}_{j}\right)\otimes...\otimes\langle x_{N}|\Psi(t)\rangle$ (23)
$\displaystyle=$ $\displaystyle\int
d\mathbf{x}\Psi^{*}(\mathbf{x},t)(-i)\frac{\partial\Psi(\mathbf{x},t)}{\partial
x_{j}}=\int
d\mathbf{x}|\Psi(\mathbf{x},t)|^{2}\text{Imag}\left(\frac{\frac{\partial\Psi(\mathbf{x},t)}{\partial
x_{j}}}{\Psi(\mathbf{x},t)}\right)=\int
d\mathbf{x}|\Psi(\mathbf{x},t)|^{2}p_{W}^{j}(\mathbf{x},t),$

where we have used $\int d\mathbf{x}=\int dx_{1}...\int dx_{N}$,
$|\mathbf{x}\rangle=|x_{1}\rangle\otimes...\otimes|x_{N}\rangle$, and
$\left(\langle x_{j}|\hat{p}_{j}\right)=\int dx_{j}^{\prime}\langle
x_{j}|\hat{p}_{j}|x_{j}^{\prime}\rangle\langle x_{j}^{\prime}|=\langle
x_{j}|(-i)\frac{\partial}{\partial x_{j}}$. This result shows that $\langle
p_{j}\rangle(t)$ can be decomposed into different components along the
positions $\mathbf{x}$ on the configuration space. Each component
$p_{W}^{j}(\mathbf{x},t)$ is the many-body weak values of the momentum of the
$j$-th particle,

$p_{W}^{j}(\mathbf{x},t)=\text{Real}\left(\frac{\langle\mathbf{x}|\hat{P}_{j}|\Psi(t)\rangle}{\langle\mathbf{x}|\Psi(t)\rangle}\right)=\text{Imag}\left(\frac{\frac{\partial\Psi(\mathbf{x},t)}{\partial
x_{j}}}{\Psi(\mathbf{x},t)}\right)=\frac{J^{j}(\mathbf{x},t)}{|\Psi(\mathbf{x},t)|^{2}},$
(24)

with $|\Psi(\mathbf{x},t)|^{2}$ the probability of finding a particle at the
given configuration position $\mathbf{x}$. Expression (24) so reproduces
equation (9) in the paper. The last identity in (24) shows that the weak
values of the momentum is just the Bohmian velocity of the $j$-th particle in
the configuration position $\mathbf{x}$. If one deals with few particles well-
separated in the physical space, then the measurement of the many-body weak
values of each particle in the laboratory is unproblematic.

###  A.2 Development of equation (13) in the paper

The problem with the weak values in (24) appears in the laboratory when we
consider $N$ particles in the same region of the physical space, since
$p_{W}^{j}(\mathbf{x},t)$ depends on all positions of the $N$ particles. Then,
it seems impossible for practical purposes to develop a measurement protocol
identifying the $N$ positions of the particles simultaneously. Thus, we want
to rewrite (23) in a way that it only depends on one position of one of the
$N$ particles. We are interested in an expression for computing $\langle
p_{j}\rangle(t)$ as a product of a probability in the physical space,
$\mathbb{P}^{k}(x,t)$, by a weak values which is also local in the physical
space, $p_{W}^{j,k}(x,t)$, which goes like

$\langle p_{j}\rangle(t)=\int dx\mathbb{P}^{k}(x,t)\;p_{W}^{j,k}(x,t),$ (25)

where

$\mathbb{P}^{k}(x,t)=\int dx_{1}...\int dx_{k-1}\int dx_{k+1}...\int
dx_{N}\;|\Psi(\mathbf{x},t)|^{2}.$ (26)

From (23), (24), and (26), one exactly gets equation (13) in the paper,

$p_{W}^{j,k}(x,t)=\frac{\int dx_{1}...\int dx_{k-1}\int dx_{k+1}...\int
dx_{N}\;p_{W}^{j}(..,x_{k-1},x,x_{k+1},..,t)|\Psi(..,x_{k-1},x,x_{k+1},..,t)|^{2}}{\int
dx_{1}...\int dx_{k-1}\int dx_{k+1}...\int
dx_{N}|\Psi(..,x_{k-1},x,x_{k+1},..,t)|^{2}}.$ (27)

It is straightforward to show that $p_{W}^{j,k}(x,t)$ in (27) pondered by
$\mathbb{P}^{k}(x,t)$ in (26) exactly gives $\langle p_{j}\rangle(t)$.

###  A.3 Development of equation (12) in the paper

The problem with the weak values in (27) is that it seems very difficult to
identify on which $j$-th particle the momentum is (weakly) measured, and on
which $k$-particle the position is (strongly) measured. Even worst, it seems
not possible to repeat the experiment and get the position and momentum of the
same two particles as in the previous measurement. In fact, the evaluation of
$\langle p_{j}\rangle(t)$ is independent on which $k$-th particle one does the
position measurement. If there are $N$ particles in the system we can write

$\langle p_{j}\rangle(t)=\frac{1}{N}\sum_{k=1}^{N}\int
dx\;\mathbb{P}^{k}(x,t)p_{W}^{j,k}(x,t)=\int
dx\;\mathbb{P}(x,t)\frac{1}{N}\sum_{k=1}^{N}p_{W}^{j,k}(x,t),$ (28)

since $\mathbb{P}^{k}(x,t)=\mathbb{P}(x,t)$ for identical particles. Finally,
if we assume that we will not identify the $j$-th particle for the momentum
measurement, then instead of trying to get $\langle p_{j}\rangle(t)$ we will
get an average over the $N$ particles,

$\displaystyle\langle p\rangle(t)\equiv\frac{1}{N}\sum_{j=1}^{N}\langle
p_{j}\rangle(t)=\frac{1}{N}\frac{1}{N}\sum_{j=1}^{N}\sum_{k=1}^{N}\int
dx\;\mathbb{P}(x,t)p_{W}^{j,k}(x,t)=\int
dx\;\mathbb{P}(x,t)\frac{1}{N^{2}}\sum_{j=1}^{N}\sum_{k=1}^{N}p_{W}^{j,k}(x,t).$
(29)

As such, we arrive to the weak values of the momentum as in equation (12) in
the paper,

$\tilde{p}_{W}(x,t)=\frac{1}{N^{2}}\sum_{j=1}^{N}\sum_{k=1}^{N}p_{W}^{j,k}(x,t),$
(30)

which satisfies $\langle p\rangle(t)=\int
dx\;\mathbb{P}(x,t)\;\tilde{p}_{W}(x,t)$.

For identical particles, either fermions or bosons, we have
$\mathbb{P}^{k}(x,t)=\mathbb{P}^{j}(x,t)\equiv\mathbb{P}(x,t)$ for all $j,k$
since $|\Psi(..,x_{k},..,x_{j},..,t)|^{2}=|\Psi(..,x_{j},..,x_{k},..,t)|^{2}$.
We also have
$p_{W}^{j}(..,x_{l},..,x_{j},..,t)=p_{W}^{l}(..,x_{j},..,x_{l},..,t)$ when
$j,l\neq k$ because
$J^{j}(..,x_{l},..,x_{j},..,t)=J^{l}(..,x_{j},..,x_{l},..,t)$ when $j,l\neq
k$. As such, we obtain

$\displaystyle p_{W}^{j,k}(x,t)$ $\displaystyle=$ $\displaystyle
p_{W}^{l,k}(x,t)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\text{for
all }j,l\neq k,$ $\displaystyle p_{W}^{j,k}(x,t)$ $\displaystyle=$
$\displaystyle
p_{W}^{j,l}(x,t)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\text{for
all }k,l\neq j,$ $\displaystyle p_{W}^{j,k}(x,t)$ $\displaystyle=$
$\displaystyle
p_{W}^{k,j}(x,t)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\text{for
all }j\neq k,$ $\displaystyle p_{W}^{j,j}(x,t)$ $\displaystyle=$
$\displaystyle
p_{W}^{k,k}(x,t)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\text{for
all }j,k.$ (31)

For a separable wave function,
$\Psi(\mathbf{x},t)=\psi_{1}(x_{1},t)...\psi_{N}(x_{N},t)$, we obtain

$\displaystyle p_{W}^{j,k}(x,t)$ $\displaystyle=$ $\displaystyle\frac{\int
dx\;p_{W}^{j}(x,t)|\psi_{j}(x,t)|^{2}}{\int dx|\psi_{j}(x,t)|^{2}}=\int
dx\;J^{j}(x,t)=\langle p_{j}\rangle(t)\;\;\;\;\;\;\;\text{for all }j\neq k,$
(32) $\displaystyle p_{W}^{j,j}(x,t)$ $\displaystyle=$ $\displaystyle
p_{W}^{j}(x,t)=\frac{J^{j}(x,t)}{|\psi_{j}(x,t)|^{2}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\text{for
all }j,$ (33)

where $J^{j}(x,t)$ and $p_{W}^{j}(x,t)$ are the current density and the weak
values, respectively, linked to the single-particle $\psi_{j}(x_{j},t)$. So,
for separable systems, the terms $p_{W}^{j,k}(x,t)$ are spatially uniform,
while $p^{j}_{W}(x,t)$ depend strongly on the position. For nonseparable
systems such differences are not true. Thus, $p_{W}^{j,k}(x,t)$ also provides
a procedure to quantify the interaction between distinct particles.

###  A.4 Time-averaging of weak values

Figure 7: Time-averaged many-body weak values of the momentum,
$\bar{f}_{W}(\\{x_{0j}\\},t,T)=\bar{p}_{W}^{j,k}(\\{x_{0j}\\},t,T)$, from the
weak values presented in figure 6 in the paper. Only $j=k=1$ is plotted, which
is the same as $j=k=2$ for $N=2$ and $j=k=2,3$ for $N=3$. The values of
$\\{x_{0j}\\}$ are taken as the set of initial values of the respective wave
packets for $N=1$ in (a),(b), $N=2$ in (c),(d), and $N=3$ in (e),(f). Left and
right panels consider respectively a time $t=20<t_{eq}$ and $t=120>t_{eq}$,
while the integration period $T/2$ in all panels runs from $0$ to $20$.

We have seen in figure 6 in the paper that, contrary to the momentum
expectation value in figure 4 in the paper, the weak values of the momentum
only approach $0$ at higher $N$. For $N=2$ they remain within a finite range
of values after $t_{eq}$, while for $N=1$ such range is much larger. We have
also discussed how the weak values (as well as the density matrix coherences)
have a random nature for $t>t_{eq}$. Thus we compute the time-average of such
weak values on a period of time $T$,

$\bar{f}_{W}(x,t,T)=\frac{1}{T}\int_{t-T/2}^{t+T/2}dt^{\prime}\;f_{W}(x,t^{\prime}),$
(34)

where $f_{W}(x,t)$ can be any of the $4$ different types of weak values
expressed in the paper, that is, $p_{W}(x,t)$, $p_{W}^{j}(x,t)$,
$p_{W}^{j,k}(x,t)$, and $\tilde{p}_{W}(x,t)$. We consider in figure 7 the case
$f_{W}(x,t)=p_{W}^{j,k}(x,t)$, and then plot
$\bar{f}_{W}(x,t,T)=\bar{p}_{W}^{j,k}(x,t,T)$ from the weak values presented
in figure 6 in the paper, as a function of $T$. Upper, middle, lower panels
respectively for $N=1$, $N=2$, $N=3$; panels in left and right column
respectively related to a time before ($t=20<t_{eq}$) and after
($t=120>t_{eq}$) thermalization. For $N>1$, the time-average of the weak
values always vanish at $t>t_{eq}$, even at small $T$ ((d),(f)); for
$t<t_{eq}$ this may only happen at much higher $T$ ((c),(e)). The $N=1$ case
is, once more, much less smooth due to the many nodes of the wave function.

##  Appendix B Diagonal and off-diagonal elements of the density matrix in the
energy representation and its connection with weak values of the momentum

The dynamics presented in our paper is based on an initial nonequilibrium
_pure_ state evolving in a closed system, while being affected by some
‘chaotic’ disorder potential. Without disorder, such a pure state simply
evolves periodically in the underneath harmonic oscillator potential. With
disorder, we have seen its role in bringing the nonequilibrium state into an
equilibrium regime characterized e.g. by the behaviour of the expectation
values shown in figures 3 and 4 in the paper. But what is the role of disorder
on the density matrix in the energy representation?

###  B.1 Density matrix populations and coherences

We consider a pure state as a superposition of (single-particle or many-
particle) energy eigenstates,

$\langle x_{1}|\otimes...\otimes\langle
x_{N}|\Psi(t)\rangle=\langle\mathbf{x}|\Psi(t)\rangle=\sum_{n}c_{n}\langle\mathbf{x}|n\rangle
e^{-iE_{n}t}=\sum_{n}\;c_{n}R_{n}(\mathbf{x})e^{-iE_{n}t},$ (35)

being $|n\rangle$ an energy eigenstate with eigenvalue $E_{n}$, $c_{n}=\langle
n|\Psi(0)\rangle$, and $\langle\mathbf{x}|n\rangle=R_{n}(\mathbf{x})$ a real
function. The $n$-sum runs within a set $N_{act}$ of activated states, which
is defined by some quench responsible for creating the nonequilibrium initial
state; in our model, the quench is a sudden shift of the harmonic trap
translated to an initial velocity for the electrons at $t=0$. The pure state
density matrix, $\hat{\rho}(t)=|\Psi(t)\rangle\langle\Psi(t)|$, in the energy
basis becomes

$\rho_{n,m}(t)=\langle m|\hat{\rho}(t)|n\rangle=\langle
m|\Psi(t)\rangle\langle\Psi(t)|n\rangle=c_{n}\;c_{m}^{*}\;e^{i(E_{m}-E_{n})t},$
(36)

where the diagonal elements, $\rho_{n,n}$ (_populations_), are clearly time-
independent, while the off-diagonal terms (_coherences_) are not. When time-
averaged, within a time interval $T\to\infty$, the off-diagonal elements
vanish and the density matrix becomes diagonal,

$\lim_{T\to\infty}\frac{1}{T}\int_{-T/2}^{T/2}\;\rho_{n,m}(t)\;dt=c_{n}\;c_{m}^{*}\;\lim_{T\to\infty}\frac{1}{T}\int_{-T/2}^{T/2}e^{i(E_{m}-E_{n})t}dt=c_{n}\;c_{m}^{*}2\pi\delta_{n,m},$
(37)

with $\delta_{n,m}$ the Kronecker delta. This does not imply that the density
matrix in (36) (without time averaging) becomes diagonal as $t\to\infty$, but
only that the off-diagonal elements oscillate around zero.

In figure 8 we show the evolution of real and imaginary parts of a few off-
diagonal elements $\rho_{n,m}(t)$ ($N=1$), which can be written from (36) by
employing $c_{l}=|c_{l}|e^{i\theta_{l}}$ as

$\rho_{n,m}(t)=|c_{n}|\;|c_{m}|\text{cos}(\theta_{n}-\theta_{m}+(E_{m}-E_{n})t)+i|c_{n}|\;|c_{m}|\text{sin}(\theta_{n}-\theta_{m}+(E_{m}-E_{n})t),$
(38)

such that the modulus $|\rho_{n,m}(t)|=|c_{n}||c_{m}|$ is also time-
independent, as seen in panel (a). The oscillation period depends inversely on
their level splittings, $\Delta E=E_{m}-E_{n}$, as shown in panels (b)-(d),
which relate to the same level $n_{p}=201$ at the peak of $\rho_{n,m}(t)$: as
one increases the splitting, by considering a matrix element with $1$ (b), $3$
(c), and $6$ (d) levels apart, the oscillation period decreases. Also, if
considering both $n,m$ values away from the peak at $n_{p}=201$, the
respective matrix elements have exponentially smaller amplitudes. In a pure
harmonic oscillator all level splittings are simply proportional to the
difference in number of levels, but the disorder potential creates _random_
splittings. The full density matrix is given in figure 1(d) in the paper,
which shows $N_{act}\approx 70$ activated states at $t=0$, which will remain
as the main states determining the system unitary evolution at any $t$. From
these results we conclude that the thermalized system keeps memory of its
initial state through the propagation of its off-diagonal coherences (in
particular by keeping the information $c_{l}=|c_{l}|e^{i\theta_{l}}$), which
never disappear. Most importantly for the current operator as discussed in the
paper, since its diagonal populations are zero by construction.

Figure 8: Density matrix coherences for $N=1$. From the grid-diagonalization
as depicted in the inset (b) of figure 1 in the paper, the time-evolution of a
few of the off-diagonal coherences $\rho_{n,m}(t)$ is shown with respect to
the peak at $n_{p}=201$, with real (imaginary) part in blue (red). As the
energy splitting, $\Delta E=E_{m}-E_{n}$, increases from (b) to (c) to (d),
the respective oscillation period decreases. In (a) one verifies that the
respective moduli $|\rho_{n,m}(t)|^{2}$ remain constant.

###  B.2 Connection between density matrix coherences and weak values

First, we need the presence probability density $|\Psi(\mathbf{x},t)|^{2}$,
given by

$\displaystyle|\Psi(\mathbf{x},t)|^{2}$ $\displaystyle=$
$\displaystyle\text{trace}\left[\hat{\rho}(|x_{1}\rangle\otimes...\otimes|x_{N}\rangle\langle
x_{1}|\otimes...\otimes\langle x_{N}|)\right]$ (39) $\displaystyle=$
$\displaystyle\sum_{n}\sum_{m}c_{n}c_{m}^{*}\;e^{i(E_{m}-E_{n})\;t}R_{n}(\mathbf{x})R^{*}_{m}(\mathbf{x}),$

and that can be decomposed into time-independent populations and time-
dependent coherences as

$|\Psi(\mathbf{x},t)|^{2}=|\Psi_{dia}(\mathbf{x},t)|^{2}+|\Psi_{off}(\mathbf{x},t)|^{2},$
(40)

with

$\displaystyle|\Psi_{dia}(\mathbf{x},t)|^{2}$ $\displaystyle=$
$\displaystyle\sum_{n}\rho_{n,n}|R_{n}(\mathbf{x})|^{2},$
$\displaystyle|\Psi_{off}(\mathbf{x},t)|^{2}$ $\displaystyle=$
$\displaystyle\sum_{n}\sum_{m\neq
n}\rho_{n,m}(t)R_{n}(\mathbf{x})R^{*}_{m}(\mathbf{x}).$ (41)

We also need the current density $J^{j}(\mathbf{x},t)$ due to the $j$-th
particle at one given position in the configuration space,

$\displaystyle J^{j}(\mathbf{x},t)$ $\displaystyle=$
$\displaystyle\text{trace}\left[\hat{\rho}\left(\hat{p}_{j}|x_{1}\rangle\otimes...\otimes|x_{N}\rangle\langle
x_{1}|\otimes...\otimes\langle
x_{N}|+|x_{1}\rangle\otimes...\otimes|x_{N}\rangle\langle
x_{1}|\otimes...\otimes\langle x_{N}|\hat{p}_{j}\right)\right]$ (42)
$\displaystyle=$
$\displaystyle\frac{i}{2}\sum_{n}\sum_{m}c_{n}c_{m}^{*}\;e^{i(E_{m}-E_{n})\;t}\left[R_{n}(\mathbf{x},t)\;\frac{\partial
R_{m}(\mathbf{x},t)}{\partial x_{j}}-R_{m}(\mathbf{x},t)\;\frac{\partial
R_{n}(\mathbf{x},t)}{\partial x_{j}}\right],$

that can also be decomposed in terms of diagonal and off-diagonal components
as

$J^{j}(\mathbf{x},t)=J^{j}_{dia}(\mathbf{x},t)+J^{j}_{off}(\mathbf{x},t),$
(43)

with

$\displaystyle J^{j}_{dia}(\mathbf{x},t)$ $\displaystyle=$
$\displaystyle\frac{i}{2}\sum_{n}\rho_{n,n}\left[R_{n}(\mathbf{x})\;\frac{\partial
R_{n}(\mathbf{x})}{\partial x_{j}}-R_{n}(\mathbf{x})\;\frac{\partial
R_{n}(\mathbf{x})}{\partial
x_{j}}\right]=\sum_{n}\rho_{n,n}(t)\;J^{j}_{n,n}(\mathbf{x}),$ $\displaystyle
J^{j}_{off}(\mathbf{x},t)$ $\displaystyle=$
$\displaystyle\frac{i}{2}\sum_{n}\sum_{m\neq
n}\rho_{n,m}(t)\;\left[R_{n}(\mathbf{x})\;\frac{\partial
R_{m}(\mathbf{x})}{\partial x_{j}}-R_{m}(\mathbf{x})\;\frac{\partial
R_{n}(\mathbf{x})}{\partial x_{j}}\right]=\sum_{n}\sum_{m\neq
n}\rho_{n,m}(t)\;J^{j}_{m,n}(\mathbf{x}),$ (44)

where
$J^{j}_{m,n}(\mathbf{x})=\frac{i}{2}\left[R_{n}(\mathbf{x})\;\frac{\partial
R_{m}(\mathbf{x})}{\partial x_{j}}-R_{m}(\mathbf{x})\;\frac{\partial
R_{n}(\mathbf{x})}{\partial x_{j}}\right]$. The relevant point is that only
off-diagonal elements $\rho_{n,m}(t)\;J^{j}_{m,n}(\mathbf{x})$ for $n\neq m$
provide current densities, since the contribution of diagonal elements vanish,
$J^{j}_{dia}(\mathbf{x},t)=0$, in closed systems with wave function vanishing
at the boundaries. The diagonal terms $\rho_{n,n}(t)\;J^{j}_{n,n}(\mathbf{x})$
do not contribute to the total current because energy eigenstates are pure
real (or pure imaginary), so that their current $J^{j}_{n,n}(\mathbf{x})=0$ in
first equation in (44).

Both results from (40) and (43) are in agreement with the well-known
continuity equation,

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\frac{\partial|\Psi(\mathbf{x},t)|^{2}}{\partial
t}+\sum_{j=1}^{N}\frac{\partial J^{j}(\mathbf{x},t)}{\partial x_{j}}$ (45)
$\displaystyle=$
$\displaystyle\frac{\partial|\Psi_{dia}(\mathbf{x},t)|^{2}}{\partial
t}+\sum_{j=1}^{N}\frac{\partial J^{j}_{dia}(\mathbf{x},t)}{\partial
x_{j}}+\frac{\partial|\Psi_{off}(\mathbf{x},t)|^{2}}{\partial
t}+\sum_{j=1}^{N}\frac{\partial J^{j}_{off}(\mathbf{x},t)}{\partial x_{j}}$
$\displaystyle=$
$\displaystyle\frac{\partial|\Psi_{off}(\mathbf{x},t)|^{2}}{\partial
t}+\sum_{j=1}^{N}\frac{\partial J^{j}_{off}(\mathbf{x},t)}{\partial x_{j}},$

where we have used the trivial results
$\frac{\partial|\Psi_{dia}(\mathbf{x},t)|^{2}}{\partial t}=0$ (because
$|\Psi_{dia}(\mathbf{x},t)|$ is time-independent) and
$\sum_{j=1}^{N}\frac{\partial J^{j}_{dia}(\mathbf{x},t)}{\partial x_{j}}=0$
(because $J^{j}_{dia}(\mathbf{x},t)=0$).

Since we know from (36) and (38) that the coherences never vanish
($\rho_{n,m}(t)\neq 0$ for $n\neq m$) and always oscillate, we conclude that
$\frac{\partial|\Psi_{off}(\mathbf{x},t)|^{2}}{\partial t}\neq 0$, so that
(45) means that the off-diagonal probability presence
$|\Psi_{off}(\mathbf{x},t)|^{2}$ and the off-diagonal current density
$J^{j}_{off}(\mathbf{x},t)$ are dynamically changing during the whole
simulation, before and after $t_{eq}$. We notice now that
$|\Psi(\mathbf{x},t)|^{2}$ and $J^{j}_{off}(\mathbf{x},t)$ are the elements
that define the weak values in equation (3) in the paper,

$p_{W}^{j}(\mathbf{x},t)=\frac{J^{j}(\mathbf{x},t)}{|\Psi(\mathbf{x},t)|^{2}}=\frac{J^{j}_{off}(\mathbf{x},t)}{|\Psi(\mathbf{x},t)|^{2}}=\frac{1}{|\Psi(\mathbf{x},t)|^{2}}\left(\sum_{n}\sum_{m\neq
n}\rho_{n,m}(t)\;J^{j}_{m,n}(\mathbf{x})\right).$ (46)

This provides the required connection between weak values and coherences.

##  Appendix C Momentum representation and ‘phase-space’

In this appendix we provide the dynamics evolution in momentum representation
and construct a pseudo phase-space of the system.

Figure 9: Wave packet dynamics in momentum representation. This figure is
similar to figure 2 in the paper, which employed position representation, but
instead it considers momentum representation. Upper panels: wave packet
$|\Psi(p_{1},p_{2},[p_{3}],t)|^{2}$ for $N=3$ at $t=0$ in (a) and $t=150$ in
(b). Middle panels: $|\Psi(p_{1},p_{2},t)|^{2}$ for $N=2$ at $t=0$ in (c) and
$t=150$ in (d). Panel (e): wave packet $|\Psi(p_{1},t)|^{2}$ for $N=1$ at
$t=0$ (dotted) and at $t=150$ (solid). Panel (f): $1$D-view of
$|\Psi(p_{1},[p_{2}],[p_{3}],t)|^{2}$ for $N=3$ (red),
$|\Psi(p_{1},[p_{2}],t)|^{2}$ for $N=2$ (blue), and $|\Psi(p_{1},t)|^{2}$ for
$N=1$ (green), with solid (dotted) lines for the wave packet at $t=150$
($t=0$), where all $t=0$ plots overlap at $p_{0j}=20$. All plots are in log-
scale. Horizontal axis in (a),(c) ((b),(d)) is the same as in (e) ((f)).
Figure 10: ‘Phase-space’ analysis. From position $\langle x_{1}\rangle$(t) and
momentum $\langle p_{1}\rangle(t)$ expectation values in figure 4 in the
paper, we compile a respective ‘phase-space’ for $N=1$ in (a), $N=2$ in (c),
and $N=3$ in (e). Panels (b),(d),(f) are a zoom at the origin for the
respective $N$. The initial value of $\langle x_{1}\rangle(t=0)$ for each $N$
is the average of the set of respective initial positions $\\{x_{0j}\\}$,
while $\langle p_{1}\rangle(t=0)=20$ for any $N$. As one approaches the final
simulation time, $t=150>t_{eq}$, the ‘phase-space’ looks noisier for $N=1,2$
as a consequence of the recurrences as seen in the respective expectation
values in the paper. Horizontal axis in (a),(c) ((b),(d)) is the same as in
(e) ((f)).

###  C.1 Dynamics in momentum representation

Figure 2 in the paper summarizes the dynamics for the systems with $N=1,2,3$
by showing the initial and final snapshots of the respective wave functions in
_position_ representation. Since our algorithm for propagating the Schrödinger
equation uses a split-operator method where the kinetic energy is handled in
momentum representation, where it is diagonal, and then Fourier-transformed
back to position representation (where the potential terms are diagonal), for
consistency, we show in figure 9 exactly the same plots as in figure 2 in the
paper, but in _momentum_ representation. That is, we plot in the upper panels
$|\Psi(p_{1},p_{2},[p_{3}],t)|^{2}$, with $[p_{3}]$ integrated out, for $N=3$
at $t=0$ in (a) and $t=150$ in (b) (results would be the same if integrating
on $[p_{1}]$ or $[p_{2}]$). The middle panels show $|\Psi(p_{1},p_{2},t)|^{2}$
for $N=2$ at $t=0$ in (c) and $t=150$ in (d). Panel (e) shows
$|\Psi(p_{1},t)|^{2}$ for $N=1$ at $t=0$ ($t=150$) in dotted (solid) lines,
while panel (f) compiles all respective $1$D plots of
$|\Psi(p_{1},[p_{2}],[p_{3}],t)|^{2}$ for $N=3$, $|\Psi(p_{1},[p_{2}],t)|^{2}$
for $N=2$, and $|\Psi(p_{1},t)|^{2}$ for $N=1$, with solid (dotted) lines for
$t=150$ ($t=0$). Since the initial velocity is the same and centered at
$p_{0j}=20$ for each degree of freedom and for every $N$, all $1$D plots at
$t=0$ overlap. One also notices in momentum representation the same full
spread of the wave function after thermalization, at $t\gg t_{eq}$, which then
remains in between the ’turning points’ $p_{TP}=\pm 20$ (since $\omega=1$). As
in position representation, the higher the $N$ the more homogeneous is the
wave packet spread.

###  C.2 Pseudo phase-space

We have shown in the paper that thermalization provides time-independent
expectation values, so that another useful plot is a pseudo ‘phase-space’ as
if that could be simply compiled from the expectation values of $\langle
p_{j}\rangle(t)$ and $\langle x_{j}\rangle(t)$ in figure 4 in the paper.
Figure 10 shows the ‘phase space’ for $N=1,2,3$ respectively in upper, middle,
lower panels; left column the full range, right column a zoom at the origin
$(0,0)$. The figures look the same no matter which $j$-electron is considered
at each $N$. For any $N$ the curves start at $\langle p_{1}\rangle(0)=20$ and
$\langle x_{1}\rangle(0)=0$, only exception being $N=2$ which starts at
$\langle x_{1}\rangle(0)=-2$. The curves after thermalization (right column)
are about the same at any $N$, with $\langle p_{1}\rangle(t)\approx\langle
x_{1}\rangle(t)\approx 0$; they look visually different since some recurrences
may happen causing a noisier ‘phase space’ as time keeps rolling after
$t_{eq}$; however, once more, the behaviour is more smooth as $N$ increases.

## References

  * [1] R. Penrose, The Emperor’s New Mind: Concerning computers, Brains and the Laws of Physics (Oxford Univ. Press, Oxford 1989).
  * [2] S. W. Hawking, Phys. Rev. D 32, 2489 (1985).
  * [3] G. C. Ghirardi, A. Rimini, and T. Weber, Phys. Rev. D 34, 470 (1986).
  * [4] L. Sklar, Physics and Change: Philosophical Issues in the Foundations of Statistical Mechanics (Cambridge Univ. Press, Cambridge 1993).
  * [5] D. Z. Albert, Time and Change (MA: Harvard Univ. Press, Cambridge 2000).
  * [6] J. Bricmont, D. Durr, M. Galavotti, F. Petruccione, and N. Zanghi, Change in Physics: Foundations and Perspectives (Springer, Berlin 2001).
  * [7] R. Penrose, J. Stat. Phys. 77, 217 (1994).
  * [8] H. D. Zeh, The Physical Basis of the Direction of Time (5th ed., Springer, Berlin 2007).
  * [9] M. Ueda, Nat. Rev. Phys. 2, 669 (2020).
  * [10] C. Gogolin and J. Eisert, Rep. Prog. Phys. 79, 056001 (2016).
  * [11] J. Eisert, M. Friesdorf, and C. Gogolin, Nat. Phys. 11, 124 (2015).
  * [12] V. I. Yukalov, Laser Phys. Lett. 8, 485 (2011).
  * [13] P. Reimann, New J. Phys. 21, 053014 (2019).
  * [14] M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 (2008).
  * [15] R. Nandkishore and D. A. Huse, Annu. Rev. Condens. Matter Phys. 6, 15 (2015).
  * [16] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Adv. Phys. 65, 239 (2016).
  * [17] J. M. Deutsch, Rep. Prog. Phys. 81, 082001 (2018).
  * [18] T. Langen, R. Geiger, and J. Schmiedmayer, Annu. Rev. Condens. Matter Phys. 6, 201 (2015).
  * [19] M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. Sen(De), and U. Sen, Adv. Phys. 56, 243 (2007).
  * [20] L. Sanchez-Palencia, D. Clément, P. Lugan, P. Bouyer, and A. Aspect, New J. Phys. 10, 045019 (2008).
  * [21] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
  * [22] P. Reimann, New J. Phys. 17, 055025 (2015).
  * [23] T. R. Oliveira, C. Charalambous, D. Jonathan, M. Lewenstein, and A. Riera, New J. Phys. 20, 033032 (2018).
  * [24] H. Kim, T. N. Ikeda, and D. A. Huse, Phys. Rev. E 90, 052105 (2014).
  * [25] T. Kinoshita, T. Wenger, and D. Weiss, Nature 440, 900 (2006).
  * [26] S. Trotzky, Y-A. Chen, A. Flesch, I. P. McCulloch, U. Schollwöck, J. Eisert, and I. Bloch, Nat. Phys. 8, 325 (2012).
  * [27] A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss, and M. Greiner, Science 353, 794 (2016).
  * [28] J. P. Ronzheimer, M. Schreiber, S. Braun, S. S. Hodgman, S. Langer, I. P. McCulloch, F. Heidrich-Meisner, I. Bloch, and U. Schneider, Phys. Rev. Lett. 110, 205301 (2013).
  * [29] D. Dries, S. E. Pollack, J. M. Hitchcock, and R. G. Hulet, Phys. Rev. A 82, 033603 (2010).
  * [30] S. Will, D. Iyer, and M. Rigol, Nat. Commun. 6, 6009 (2015).
  * [31] J. Kajala, F. Massel, and P. Torma, Phys. Rev. Lett. 106, 206401 (2011).
  * [32] U. Schneider, L. Hackermuller, S. Will, T. Best, I. Bloch, T. A. Costi, R. W. Helmes, D. Rasch, and A. Rosch, Science 322, 1520 (2008).
  * [33] B. Nagler, K. Jagering, A. Sheikhan, S. Barbosa, J. Koch, S. Eggert, I. Schneider, and A. Widera, Phys. Rev. A 101, 053633 (2020).
  * [34] U. Schneider, L. Hackermuller, J. P. Ronzheimer, S. Will, S. Braun, T. Best, I. Bloch, E. Demler, S. Mandt, D. Rasch, and A. Rosch, Nat. Phys. 8, 213 (2012).
  * [35] W. C. Myrvold, Synthese 192, 3247 (2015).
  * [36] H-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford Univ. Press, Oxford, 2002).
  * [37] In other quantum theories, like Bohmian mechanics or many-worlds, the collapse is just an ‘effective’ phenomenon whose consequences can still be treated from a bigger unitary Schrödinger equation that includes the measuring apparatus.
  * [38] M. Srednicki, Phys. Rev. E 50, 888 (1994).
  * [39] J. M. Deutsch, Phys. Rev. A 43, 2046 (1991).
  * [40] Y. Aharonov, D. Z. Albert, and L. Vaidman, Phys. Rev. Lett. 60, 1351 (1988).
  * [41] B. E. Y. Svensson, Quanta 2, 18 (2013).
  * [42] D. Pandey, R. Sampaio, T. Ala-Nissila, G. Albareda, and X. Oriols, Phys. Rev. A 103, 052219 (2021).
  * [43] H. M. Wiseman, New J. Physics 9, 165 (2007).
  * [44] D. Dürr, S. Goldstein, and N. Zanghì, J. Stat. Phys. 134, 1023 (2009).
  * [45] D. Marian, N. Zanghi, and X. Oriols, Phys. Rev. Lett. 116, 110404 (2016).
  * [46] F. L. Traversa, G. Albareda, M. Di Ventra, and X. Oriols, Phys. Rev. A 87, 052124 (2013).
  * [47] S. Kocsis, B. Braverman, S. Ravets, M. J. Stevens, R. P. Mirin, L. K. Shalm, and A. M. Steinberg, Science 332, 1170 (2011).
  * [48] A. Hariri, D. Curic, L. Giner, and J. S. Lundeen, Phys. Rev. A 100, 032119 (2019).
  * [49] R. Ramos, D. Spierings, I. Racicot, and A. M. Steinberg, Nature 583, 529 (2020).
  * [50] A.G. Kofman, S. Ashhab, and F. Nori, Physics Reports 520, 43 (2012).
  * [51] J. Dressel, M. Malik, F. M. Miatto, A. N. Jordan, and R. W. Boyd, Rev. Mod. Phys. 86, 307 (2014).
  * [52] K. Renziehausen and I. Barth, Found. Phys. 50, 772 (2020).
  * [53] O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, Phys. Rev. X 6, 041065 (2016).
  * [54] V. Alba, B. Bertini, M. Fagotti, L. Piroli, and P. Ruggiero, J. Stat. Mech., 114004 (2021).
  * [55] I. Bouchole and J. Dubail, J. Stat. Mech., 014003 (2022).
  * [56] X. Oriols and J. Mompart, Applied Bohmian Mechanics: From Nanoscale Systems to Cosmology (Pan Stanford, Singapore, 2012).
  * [57] A. O. T. Pang, H. Ferretti, N. Lupu-Gladstein, W.-K. Tham, A. Brodutch, K. Bonsma-Fisher, J. E. Sipe, and A. M. Steinberg, Quantum 4, 365 (2020).
  * [58] T. Fogarty, M. Á. García-March, L. F. Santos, and N. L. Harshman, Quantum 5, 486 (2021).
  * [59] G. Zisling, L. F. Santos, and Y. B. Lev, SciPost. Phys. 10, 088 (2021).
  * [60] C. Neill, P. Roushan, M. Fang, Y. Chen, M. Kolodrubetz, Z. Chen, A. Megrant, R. Barends, B. Campbell, B. Chiaro, A. Dunsworth, E. Jeffrey, J. Kelly, J. Mutus, P. J. J. O’Malley, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, A. Polkovnikov, and J. M. Martinis, Nat. Phys. 12, 1037 (2016).
  * [61] Md. M. Ali, W.-M. Huang, and W.-M. Zhang, Sci. Rep. 10, 13500 (2020).
  * [62] J.-Q. Liao, H. Dong, and C. P. Sun, Phys. Rev. A 81, 052121 (2010).
  * [63] P. Lydzba, Y. Zhang, M. Rigol, and L. Vidmar, Phys. Rev. B 104, 214203 (2021).
  * [64] C. Nation and D. Porras, Phys. Rev. E 102, 042115 (2020).
  * [65] M. Rigol, Phys. Rev. A 80, 053607 (2009).
  * [66] L. F. Santos and M. Rigol, Phys. Rev. E 81, 036206 (2010).
  * [67] Y. Tang, W. Kao, K.-Y. Li, S. Seo, K. Mallayya, M. Rigol, S. Gopalakrishnan, and B. L. Lev, Phys. Rev. X 8, 021030 (2018).
  * [68] C. Gogolin, M. P. Muller, and J. Eisert, Phys. Rev. Lett. 106, 040401 (2011).
  * [69] M. Brenes, T. LeBlond, J. Goold, and M. Rigol, Phys. Rev. Lett. 125, 070605 (2020).
  * [70] M. Brenes, J. Goold, and M. Rigol, Phys. Rev. B 102, 075127 (2020).
  * [71] Y.-W. Hsueh, C.-H. Hsueh, and W.-C. Wu, Entropy 22, 855 (2020).
  * [72] L. Pezze, B. Hambrecht, and L. Sanchez-Palencia, EPL 88, 30009 (2009).
  * [73] L. Pezze, and L. Sanchez-Palencia, Phys. Rev. Lett. 106, 040601 (2011).
  * [74] C.-H. Hsueh, R. Ong, J.-F. Tseng, M. Tsubota, and W.-C. Wu, Phys. Rev. A 98, 063613 (2018).
  * [75] R. Ong, C.-H. Hsueh, and W.-C. Wu, Phys. Rev. A 100, 053619 (2019).
  * [76] D. Clément, A. F. Varón, J. A. Retter, L. Sanchez-Palencia, A. Aspect, and P. Bouyer, New J. Phys. 8, 165 (2006).
  * [77] J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht, P. Lugan, D. Clément, L. Sanchez-Palencia, P. Bouyer, and Alain Aspect, Nature 453, 891 (2008).
  * [78] All figures in the paper consider $\alpha=0.1$, $\omega=1.0$, $\sigma_{j}=l_{0}=1/\sqrt{\omega}=1.0$; initial velocities are $p_{0j}=20$ for any $j$ and $N$; initial positions are $x_{01}=0$ for $N=1$, $(x_{01},x_{02})=(-4,0)$ for $N=2$, $(x_{01},x_{02},x_{03})=(-4,0,4)$ for $N=3$, having in mind that the initial wave packet is antisymmetrized such that the label $j$ in $x_{0j}$ is redundant; for the disorder, $\gamma_{D}=25$ and $\sigma_{D}\approx\Delta x$, with $\Delta x$ being the position resolution of the respective grid, defined in a way that its momentum resolution is $\Delta p\approx\Delta x$, since our propagator uses both position and momentum representations. Initial time is $t_{0}=0$, final time is $t_{f}=150$ (for $N=2$ $t_{f}=300$), with time step $\Delta t=0.001$. Simulation boxes are roughly twice the extension of the turning points, and $\approx 10^{3}$ points per dimension are used; so e.g for $N=3$ our configuration space has $\approx 10^{9}$ points.
  * [79] A. Maksymov, P. Sierant, and J. Zakrzewski, Phys. Rev. B 102, 134205 (2020).
  * [80] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, Nat. Phys. 12, 907 (2016).
  * [81] L. F. Santos and M. Rigol, Phys. Rev. E 82, 031130 (2010).