email<EMAIL_ADDRESS>email<EMAIL_ADDRESS>
# Kinetic energy equipartion: a tool to characterize quantum thermalization

Carlos F. Destefani  [  Xavier Oriols  [  Department of Electronic
Engineering, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona,
Spain

###### Abstract

Weak values have gradually been transitioning from a theoretical curiosity to
a practical novel tool to characterize quantum systems in the laboratory. When
built by post-selecting the position, weak values of the momentum are linked
to the so-called hidden varaiables of the Bohmian and Stochastic quantum
mechanics. According to both theories, the Orthodox kinetic energy has, in
fact, two hidden-variable components: one linked to the current (or Bohmian)
velocity, and another linked to the osmotic velocity (or quantum potential),
and which are respectively identified with phase and amplitude of the
wavefunction. Inspired by such alternative formulations, we address what
happens to each of these two velocity components when the Orthodox kinetic
energy thermalizes in closed systems, and how the pertinent weak values yield
experimental information about them. We show that, after thermalization, the
expectation values of both the (squared) current and osmotic velocities
approach the same stationary value, that is, each of the Bohmian kinetic and
quantum potential energies approaches half of the Orthodox kinetic energy.
Such a ‘kinetic energy equipartition’ is a novel signature of quantum
thermalization that can empirically be tested in the laboratory, following a
well-defined operational protocol as given by the expectation values of
(squared) real and imaginary parts of the local-in-position weak value of the
momentum, which are respectively related to the current and osmotic
velocities. Thus, the kinetic energy equipartion presented here is independent
on any ontological status given to these hidden variables, and it could be
used as a novel element to characterize quantum thermalization in the
laboratory, beyond the traditional use of expectation values linked to
Hermitian operators. Numerical results for the nonequilibrium dynamics of a
few-particle harmonic trap under random disorder are presented as
illustration. And the advantages in using the center-of-mass frame of
reference for dealing with systems with many indistinguishable particles are
also discussed.

##  I Introduction

A renewed interest in statistical mechanics of closed quantum systems has
arisen [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15] as a consequence of
the successful experimental ability to isolate and manipulate bosonic [16, 17,
18, 19, 20] and fermionic [21, 22, 23, 24, 25] many-body systems built on
ultra-cold atomic gases subjected to optical lattices. The main question to be
addressed is when an initial nonequilibrium state thermalizes and, if so, in
which conditions. The Eigenstate Thermalization Hypothesis (ETH) [26, 27],
which has become a cornerstone in the study of quantum thermalization, claims
that all relevant energy eigenstates of a given Hamiltonian, in the
description of a quantum state, are thermal in the sense that they are similar
to an equilibrium state as long as one deals with macroscopic observables. In
recent years a large amount of numerical experiments has successfully tested
ETH by directly diagonalizing in physical space some sort of short-range many-
body lattice Hamiltonian, like Fermi- or Bose-Hubbard [17, 22, 23, 28, 29, 30]
and XXZ- or XYZ-Heisenberg [6, 15, 18, 31, 32, 33, 34], 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 the initial
nonequilibrium state. The true time-evolution of such a state is not as
widespread due to the inherently huge configuration space involved.

Our understanding about quantum thermalization has largely been based in terms
of expectation values of observables linked to Hermitian operators. It is
well-known that other explanations of quantum phenomena allow one to discuss
properties not directly linked to Hermitian operators. Such alternative
explanations are in general labeled as hidden-variable theories. The Bohmian
theory, formulated by de Broglie in 1927 [35] and further developed by Bohm in
1952 [36], is the most well-known example of a quantum theory with additional
microscopic variables: particles have always well-defined positions that
conform trajectories. Another example is the Stochastic quantum mechanics,
proposed by Nelson in 1966 [37]; although it also assumes particles with well-
defined trajectories, these are unknown and only their statistical behavior is
handled.

These two hidden-variable theories assume that the Orthodox kinetic energy in
fact has two components: in the Bohmian theory, it is computed as the square
of the so-called Bohmian velocity plus a quantum potential, while in the
Stochastic quantum mechanics, it is given by the square of a mean current
velocity plus the square of a so-called osmotic velocity. The central question
in this paper is whether these two hidden-variable components of the kinetic
energy also thermalize when the Orthodox kinetic energy thermalizes. We show
that, in fact, one can characterize the thermalization time $t_{eq}$ as the
time when the expectation value of the (squared) current and osmotic
velocities become equal or, similarly, when the expectation values of the
Bohmian kinetic and quantum potential energies become the same, with each
being equal to half of the Orthodox kinetic energy. Such a kinetic energy
equipartition is the hidden-variable signature of quantum thermalization.

Historically, hidden variables were thought to be hidden in the sense that
they were not measurable in the laboratory, so that they were useful only as a
complementary visualization tool of quantum dynamics. However, nowadays, it is
well-known that one may have empirical access to them through the so-called
weak values, which have a well-defined operational protocol, independent on
any ontology assigned to such hidden variables. Since the original single-
particle proposal [38] weak values have been attracting a lot of theoretical
[39, 40, 41, 42, 43] and experimental [44, 45, 46] interests in distinct
research fields, and a many-particle generalization has been presented by the
authors elsewhere [47]. One should clarify though that the fact that one may
have experimental access to hidden variables does not directly imply that they
are ontologically real. Each quantum theory defines its own ontological
elements. For example, the Bohmian velocity is part of the ontology of the
Bohmian mechanics (as the field which drives the particle trajectories), but
not of the Orthodox ontology. And, in fact, in the Stochastic quantum
mechanics the Bohmian velocity itself is not a direct ontological element, but
part of a mean current velocity. The discussion of such ontological status is
far from the scope of our paper. Our main focus is to emphasize that the
existence of the local-in-position weak values of the momentum opens a new
unexplored link between theoretical predictions and empirical data that allows
a novel characterization of quantum systems. One could use other types of weak
values, but working with the local-in-position weak values of the momentum
allows us to re-use all the mathematical machinery (without choosing their
ontology) of the Bohmian and Stochastic quantum mechanics. And for that, in
particular, we focus here on how the expectation values of the (squared) real
and imaginary parts of such new empirical data allow one to characterize the
process of quantum thermalization of a fermionic few-body closed system,
defined by an harmonic trap under random disorder. At the end of the paper, we
reformulate our proposal in the center-of-mass frame, so that our findings are
also extendable for larger systems with many indistinguishable particles.

The paper is organized as follows. Section II presents a brief summary of
needed theoretical background. Section III addresses the Ortodox kinetic
energy equipartition among its two hidden-variable components. Section IV
presents the theoretical model and the numerical results for a particular few-
body harmonic trap under random disorder, within three distinct scenarios for
the dynamics; the center-of-mass frame is then considered for approaching
larger systems. In Section V we conclude.

##  II Theoretical background

This section provides a summary of the needed theoretical background:
subsection II.1 presents the weak values from the polar form of the many-body
wave function, as well as the main equations of both Bohmian quantum mechanics
and stochastic quantum mechanics; subsection II.2 develops the non-hermitian
expectation values derived from such weak values. To simplify the notation,
atomic units are employed throughout the text, which assumes $N$
nonrelativistic spinless particles each in a $1$D physical space, so that the
position in configuration space is $\mathbf{x}=(x_{1},...,x_{N})$.

###  II.1 Weak values post-selected in position

A simple path to define the (complex) local-in-position weak value of the
momentum for particle $j$ at position $\mathbf{x}$, $p_{W,j}(\mathbf{x},t)$,
linked to the hermitian operator for the momentum,
$\hat{p}_{j}=-i\;\partial/\partial x_{j}$, comes after inserting the identity
$\int d\mathbf{x}|\mathbf{x}\rangle\langle\mathbf{x}|$ into the expectation
value $\langle\Psi(t)|\hat{p}_{j}|\Psi(t)\rangle$,

$\langle\hat{p}_{j}\rangle=\langle\Psi(t)|\hat{p}_{j}|\Psi(t)\rangle=\int
d\mathbf{x}|\Psi(\mathbf{x},t)|^{2}\frac{\langle\mathbf{x}|\hat{p}_{j}|\Psi(t)\rangle}{\langle\mathbf{x}|\Psi(t)\rangle}.$
(1)

From the polar form of the wave function,
$\Psi(\mathbf{x},t)=R(\mathbf{x},t)\text{exp}(iS(\mathbf{x},t))$ with
$|\Psi(\mathbf{x},t)|^{2}=R^{2}(\mathbf{x},t)$, the weak value as stated in
(1) decomposes like

$\displaystyle p_{W,j}(\mathbf{x},t)$ $\displaystyle\equiv$
$\displaystyle\frac{\langle\mathbf{x}|\hat{p}_{j}|\Psi(t)\rangle}{\langle\mathbf{x}|\Psi(t)\rangle}=\frac{\partial
S(\mathbf{x},t)}{\partial x_{j}}-i\;\frac{1}{R(\mathbf{x},t)}\frac{\partial
R(\mathbf{x},t)}{\partial x_{j}},$ (2) $\displaystyle=$ $\displaystyle
v_{j}(\mathbf{x},t)-i\;u_{j}(\mathbf{x},t),$

where the current $v_{j}(\mathbf{x},t)$ and osmotic $u_{j}(\mathbf{x},t)$
velocities for particle $j$ are identified in (2) as

$\displaystyle v_{j}(\mathbf{x},t)$ $\displaystyle=$
$\displaystyle\frac{\partial S(\mathbf{x},t)}{\partial
x_{j}}=\text{Im}\left[\frac{1}{\Psi(\mathbf{x},t)}\frac{\partial\Psi(\mathbf{x},t)}{\partial
x_{j}}\right],$ (3) $\displaystyle u_{j}(\mathbf{x},t)$ $\displaystyle=$
$\displaystyle\frac{1}{R(\mathbf{x},t)}\frac{\partial
R(\mathbf{x},t)}{\partial
x_{j}}=\text{Re}\left[\frac{1}{\Psi(\mathbf{x},t)}\frac{\partial\Psi(\mathbf{x},t)}{\partial
x_{j}}\right].$ (4)

Notice that $v_{j}(\mathbf{x},t)$ and $u_{j}(\mathbf{x},t)$ depend only,
respectively, on the phase and on the amplitude of $\Psi(\mathbf{x},t)$. The
strategy above in fact can be used for any operator; for the bilinear
momentum, $\hat{p}_{j}\hat{p}_{l}=-\partial/\partial x_{j}\partial/\partial
x_{l}$, one has

$\langle\hat{p}_{j}\hat{p}_{l}\rangle=\langle\Psi(t)|\hat{p}_{j}\hat{p}_{l}|\Psi(t)\rangle=\int
d\mathbf{x}|\Psi(\mathbf{x},t)|^{2}\frac{\langle\mathbf{x}|\hat{p}_{j}\hat{p}_{l}|\Psi(t)\rangle}{\langle\mathbf{x}|\Psi(t)\rangle},$
(5)

where the bilinear weak value $p_{W,jl}(\mathbf{x},t)$ is

$\displaystyle
p_{W,jl}(\mathbf{x},t)\equiv\frac{\langle\mathbf{x}|\hat{p}_{j}\hat{p}_{l}|\Psi(t)\rangle}{\langle\mathbf{x}|\Psi(t)\rangle}$
$\displaystyle=v_{j}(\mathbf{x},t)v_{l}(\mathbf{x},t)-u_{j}(\mathbf{x},t)u_{l}(\mathbf{x},t)-\frac{\partial
u_{l}(\mathbf{x},t)}{\partial x_{j}}$ (6)
$\displaystyle-i\;\left(u_{j}(\mathbf{x},t)v_{l}(\mathbf{x},t)+v_{j}(\mathbf{x},t)u_{l}(\mathbf{x},t)+\frac{\partial
v_{l}(\mathbf{x},t)}{\partial x_{j}}\right),$

which, for $j=l$, defines the local-in-position weak value of (twice) the
kinetic energy for particle $j$, $K_{W,j}(\mathbf{x},t)$,

$\displaystyle
2K_{W,j}(\mathbf{x},t)\equiv\frac{\langle\mathbf{x}|\hat{p}_{j}^{2}|\Psi(t)\rangle}{\langle\mathbf{x}|\Psi(t)\rangle}=v^{2}_{j}(\mathbf{x},t)-u^{2}_{j}(\mathbf{x},t)$
$\displaystyle-\frac{\partial u_{j}(\mathbf{x},t)}{\partial
x_{j}}-i\;\left(2u_{j}(\mathbf{x},t)v_{j}(\mathbf{x},t)+\frac{\partial
v_{j}(\mathbf{x},t)}{\partial x_{j}}\right).$ (7)

Subsecion II.2 will show that the imaginary part does not contribute for
ensemble values of (2), (6), or (7). One could also define from (2) and (6)
the (weak) correlation of the momentum post-selected in position,
$C_{p_{W,j};p_{W,l}}\equiv\langle p_{W,jl}(\mathbf{x},t)\rangle-\langle
p_{W,j}(\mathbf{x},t)\rangle\langle p_{W,l}(\mathbf{x},t)\rangle$, as

$C_{p_{W,j};p_{W,l}}=-\langle\frac{\partial u_{l}(\mathbf{x},t)}{\partial
x_{j}}\rangle-i\;\langle\frac{\partial v_{l}(\mathbf{x},t)}{\partial
x_{j}}\rangle,$ (8)

which could be used in situations where separated entanglements in phase and
amplitude of $\Psi(\mathbf{x},t)$ were accessible [48], since the real
(imaginary) part of $C_{p_{W,j};p_{W,l}}$ depends only on the amplitude
(phase) of $\Psi(\mathbf{x},t)$. Some given properties among current (3) and
osmotic (4) velocities will prove important in our derivation. For that, one
can use some basic elements of both Bohmian and Stochastic quantum mechanics.

####  II.1.1 Bohmian quantum mechanics

The Bohmian theory assumes that particles follow well-defined trajectories
guided by $\Psi(\mathbf{x},t)$, solution of the Schroedinger equation

$i\frac{\partial\Psi(\mathbf{x},t)}{\partial
t}=H(\mathbf{x})\Psi(\mathbf{x},t),$ (9)

whose Hamiltonian is

$H(\mathbf{x})=-\frac{1}{2}\sum_{j=1}^{N}\frac{\partial^{2}}{\partial
x_{j}^{2}}+V(\mathbf{x})=K(\mathbf{x})+V(\mathbf{x}),$ (10)

with $K(\mathbf{x})=\sum_{j=1}^{N}K_{j}(x_{j})$, and
$K_{j}(x_{j})\equiv-\partial^{2}/2\partial x_{j}^{2}$ the kinetic energy of
each of the $N$ particles, whose interaction is given by the potential energy
$V(\mathbf{x})$. From the polar form of $\Psi(\mathbf{x},t)$, equation (9)
gets rewritten as two coupled equations. On one hand, the imaginary part
yields the continuity equation,

$\frac{\partial R^{2}(\mathbf{x},t)}{\partial
t}+\sum_{j=1}^{N}R^{2}(\mathbf{x},t)\frac{\partial S(\mathbf{x},t)}{\partial
x_{j}}=0,$ (11)

in which from (3) one identifies the $j$-component of the current density as
$J_{j}(\mathbf{x},t)\equiv R^{2}(\mathbf{x},t)v_{j}(\mathbf{x},t)$. On the
other hand, the real part yields the quantum Hamilton-Jacobi equation,

$\frac{\partial S(\mathbf{x},t)}{\partial
t}+K_{B}(\mathbf{x},t)+Q_{B}(\mathbf{x},t)+V(\mathbf{x})=0,$ (12)

with $K_{B}(\mathbf{x},t)=\sum_{j=1}^{N}K_{B,j}(\mathbf{x},t)$ and
$Q_{B}(\mathbf{x},t)=\sum_{j=1}^{N}Q_{B,j}(\mathbf{x},t)$, and where the
$j$-component of the Bohmian kinetic $K_{B,j}(\mathbf{x},t)$ and quantum
potential $Q_{B,j}(\mathbf{x},t)$ energies are

$\displaystyle K_{B,j}(\mathbf{x},t)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\frac{\partial S(\mathbf{x},t)}{\partial
x_{j}}\right)^{2},$ (13) $\displaystyle Q_{B,j}(\mathbf{x},t)$
$\displaystyle=$
$\displaystyle-\frac{1}{2R(\mathbf{x},t)}\frac{\partial^{2}R(\mathbf{x},t)}{\partial
x^{2}_{j}}.$ (14)

The Bohmian trajectories, thanks to the current density $J_{j}(\mathbf{x},t)$
in (11), are obtained from the integration solely of $v_{j}(\mathbf{x},t)$
($u_{j}(\mathbf{x},t)$ plays no role); so, along the text
$v_{j}(\mathbf{x},t)$ is interchangeably identified either as Bohmian or
current velocity. In Bohmian theory, while the kinetic energy is determined by
the phase of the wavefunction, the quantum potential, determined by its
amplitude, has indeed the status of a potential energy since, from a time
derivative in (3) and by using (12), one gets the set of Newton-like
equations,

$\frac{d^{2}x_{j}(t)}{dt^{2}}=-\left[\frac{\partial}{\partial
x_{j}}\left(V(\mathbf{x},t)+Q_{B}(\mathbf{x},t)\right)\right]_{\mathbf{x}=\mathbf{x(t)}},$
(15)

which is another way for getting Bohmian trajectories.

####  II.1.2 Stochastic quantum mechanics

The Stochastic quantum mechanics may be understood as an attempt to give a
kinematic interpretation also to the quantum potential, although such a
potential is not explicit in its original derivation. Notice indeed that, from
(3)-(4), equations (13)-(14) can be cast into

$\displaystyle K_{B,j}(\mathbf{x},t)$ $\displaystyle=$
$\displaystyle\frac{1}{2}v_{j}^{2}(\mathbf{x},t),$ (16) $\displaystyle
Q_{B,j}(\mathbf{x},t)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}u_{j}^{2}(\mathbf{x},t)-\frac{1}{2}\frac{\partial
u_{j}(\mathbf{x},t)}{\partial x_{j}},$ (17)

so that $K_{B,j}(\mathbf{x},t)$ and $Q_{B,j}(\mathbf{x},t)$, respectively,
depend only on Bohmian and osmotic velocities. Notice that
$K_{B,j}(\mathbf{x},t)$ and $Q_{B,j}(\mathbf{x},t)$ are not weak values
computed from the hermitian operator of the kinetic energy, as it happens with
$K_{W,j}(\mathbf{x},t)$ in (7). One should keep in mind that, according to
this stochastic theory, $v_{j}(\mathbf{x},t)+u_{j}(\mathbf{x},t)$ has the
meaning of a mean velocity, while the true velocity is a random velocity
around it. Such a theory is defined in terms of a stochastic diffusion process
in the configuration space, which requires that the probability
$R^{2}(\mathbf{x},t)$ satisfies both forward ($+$) and backward ($-$) Fokker-
Plank equations for a parameter $\nu$,

$\displaystyle\frac{\partial R^{2}(\mathbf{x},t)}{\partial t}$
$\displaystyle=$ $\displaystyle-\sum_{j}\frac{\partial}{\partial
x_{j}}\left((v_{j}(\mathbf{x},t)\pm
u_{j}(\mathbf{x},t))R^{2}(\mathbf{x},t)\right)$ (18) $\displaystyle\pm$
$\displaystyle\nu\sum_{j}\frac{\partial^{2}R^{2}(\mathbf{x},t)}{\partial
x_{j}^{2}},$

whose sum yields the continuity equation in (11), irrespective of
$u_{j}(\mathbf{x},t)$ and $\nu$, and whose difference yields

$\nu\frac{\partial^{2}R^{2}(\mathbf{x},t)}{\partial
x_{j}^{2}}=\frac{\partial\left(u_{j}(\mathbf{x},t)R^{2}(\mathbf{x},t)\right)}{\partial
x_{j}},$ (19)

which is satisfied by the $u_{j}(\mathbf{x},t)$ in (4) and with $\nu=1/2$. In
other words, the Fokker-Planck equations in (18), with (4) and $\nu=1/2$,
pictures from (19) that the osmotic current
$u_{j}(\mathbf{x},t)R^{2}(\mathbf{x},t)$ is balanced by some diffusion current
$\nu\partial R^{2}(\mathbf{x},t)/\partial x_{j}$, so that the continuity
equation in (11) remains valid and determined solely by the current velocity
$v_{j}(\mathbf{x},t)$. As such, the kinematic interpretation of the quantum
potential, implicit in the derivation of the Stochastic quantum mechanics,
empirically reproduces both Bohmian and Orthodox quantum mechanics.

###  II.2 Expectation values from weak values

From (1)-(2) the expectation value for the momentum weak value post-selected
in position, $\langle p_{W,j}\rangle$, is

$\langle\hat{p}_{j}\rangle=\int
d\mathbf{x}\;p_{W,j}(\mathbf{x},t)|\Psi(\mathbf{x},t)|^{2}\equiv\langle
p_{W,j}\rangle.$ (20)

For this particular case one can directly compute the expectation values of
each of its real and imaginary parts, yielding from (1)-(2)

$\displaystyle\langle v_{j}\rangle$ $\displaystyle=$ $\displaystyle\int
d\mathbf{x}\;\text{Re}\left[p_{W,j}(\mathbf{x},t)\right]|\Psi(\mathbf{x},t)|^{2}=\langle\hat{p}_{j}\rangle,$
(21) $\displaystyle\langle u_{j}\rangle$ $\displaystyle=$ $\displaystyle\int
d\mathbf{x}\;\text{Im}\left[p_{W,j}(\mathbf{x},t)\right]|\Psi(\mathbf{x},t)|^{2}=0,$
(22)

so that, while $v_{j}(\mathbf{x},t)\neq p_{W,j}(\mathbf{x},t)$, $\langle
v_{j}\rangle=\langle p_{W,j}\rangle$ and the osmotic velocity has no role in
the expectation values of neither $\hat{p}_{j}$ nor $p_{W,j}(\mathbf{x},t)$
(for a wave function vanishing at the system borders). Strictly speaking,
$v_{j}(\mathbf{x},t)$ and $u_{j}(\mathbf{x},t)$ are not weak values themselves
but they are, respectively, post-processed real and imaginary parts of the
weak value $p_{W,j}(\mathbf{x},t)$. In other words, while
$p_{W,j}(\mathbf{x},t)$ is linked to the hermitian operator $\hat{p}_{j}$, no
hermitian operators can be linked to $v_{j}(\mathbf{x},t)$ and
$u_{j}(\mathbf{x},t)$. The connection among $p_{W,j}(\mathbf{x},t)$ and
$v_{j}(\mathbf{x},t)$ is approached elsewhere [49, 50, 51, 52, 53, 54, 55,
56], while more attention has recently been paid to the meaning of
$u_{j}(\mathbf{x},t)$ [57, 58, 59, 60, 61, 62, 63, 64, 65]. We distinguish in
the paper three types of expectation values: (i) $\langle\hat{a}\rangle$ with
“hat” for the operator $\hat{a}$; (ii) $\langle a_{W}\rangle$ with subindex
“W” for the weak value $a_{W}$; (iii) $\langle a\rangle$ for the value $a$
obtained by post-processing the weak value.

The expectation value of the bilinear weak value, $\langle p_{W,jl}\rangle$,
is from (5)-(6)

$\displaystyle\langle\hat{p}_{j}\hat{p}_{l}\rangle$ $\displaystyle=$
$\displaystyle\int
d\mathbf{x}\;p_{W,jl}(\mathbf{x},t)|\Psi(\mathbf{x},t)|^{2}\rangle,$ (23)
$\displaystyle=$ $\displaystyle\langle v_{j}v_{l}\rangle+\langle
u_{j}u_{l}\rangle\equiv\langle p_{W,jl}\rangle,$

where we have used

$\displaystyle\left\langle\frac{\partial u_{l}}{\partial x_{j}}\right\rangle$
$\displaystyle=$ $\displaystyle-2\langle u_{l}u_{j}\rangle,$ (24)
$\displaystyle\left\langle\frac{\partial v_{l}}{\partial x_{j}}\right\rangle$
$\displaystyle=$ $\displaystyle-\langle u_{l}v_{j}\rangle-\langle
v_{l}u_{j}\rangle;$ (25)

for an antisymmetric $\Psi(\mathbf{x},t)$ one has $\langle
u_{l}v_{j}\rangle=\langle v_{l}u_{j}\rangle$. Once more, the weak value
$p_{W,jl}(\mathbf{x},t)$ is linked to $\hat{p}_{j}\hat{p}_{l}$, but
$v_{j}(\mathbf{x},t)v_{l}(\mathbf{x},t)$ and
$u_{j}(\mathbf{x},t)u_{l}(\mathbf{x},t)$ are just post-processed data obtained
from the real and imaginary parts of the weak values $p_{W,j}(\mathbf{x},t)$
and $p_{W,l}(\mathbf{x},t)$. Interestingly, from (24) the expectation values
of the Bohmian kinetic $\langle K_{B,j}\rangle$ and quantum potential $\langle
Q_{B,j}\rangle$ energies in (16) and (17) become

$\displaystyle\langle K_{B,j}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2}\langle v_{j}^{2}\rangle,$ (26) $\displaystyle\langle
Q_{B,j}\rangle$ $\displaystyle=$ $\displaystyle\frac{1}{2}\langle
u_{j}^{2}\rangle,$ (27) $\displaystyle\langle\hat{K}_{j}\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{2}(\langle v^{2}_{j}\rangle+\langle
u^{2}_{j}\rangle),$ (28)

where the last equation is the expectation value of the Orthodox kinetic
energy, $\langle\hat{K}_{j}\rangle=\langle\hat{p}^{2}_{j}\rangle/2$, as
promptly obtained from (23). Neither $\langle v_{j}^{2}\rangle$ nor $\langle
u_{j}^{2}\rangle$ are expectation values of a weak value but instead,
respectively, of the (squared) real and imaginary parts of the weak value
$p_{W,j}(\mathbf{x},t)$. Notice that by integrating the weak values of the
kinetic energy in (7) with the use of (24)-(25), and comparing the result with
(28), one obtains

$\langle\hat{K}_{j}\rangle=\int
d\mathbf{x}\;K_{W,j}(\mathbf{x},t)|\Psi(\mathbf{x},t)|^{2}=\langle
K_{W,j}\rangle.$ (29)

Similar to the (weak) correlation in (8), one can define the correlations
$C_{\hat{a},\hat{b}}\equiv\langle\hat{a}\hat{b}\rangle-\langle\hat{a}\rangle\langle\hat{b}\rangle$
between operators and $C_{a,b}\equiv\langle ab\rangle-\langle a\rangle\langle
b\rangle$ between quantities post-processed from weak values. For the momentum
operator, from (23) and (21) one obtains

$\displaystyle C_{\hat{p}_{j},\hat{p}_{l}}=\langle v_{j}v_{l}\rangle+\langle
u_{j}u_{l}\rangle-\langle v_{j}\rangle\langle
v_{l}\rangle=C_{v_{j},v_{l}}+C_{u_{j},u_{l}},$ (30) $\displaystyle
C_{p_{W,j};p_{W,l}}=2\langle u_{j}u_{l}\rangle+i\left(\langle
u_{j}v_{l}\rangle+\langle v_{j}u_{l}\rangle\right),$ (31)

where (22) is used in (30), while (31) results from applying (24)-(25) in (8).
So, while the Bohmian velocity fully determines $\langle\hat{p}_{j}\rangle$ in
(21), the osmotic velocity, although satisfying (22), induces a deviation
between the quantum correlations of the momentum with respect to the curent
velocity in (30). Notice that from (30)-(31) one has
$\text{Re}[C_{p_{W,j};p_{W,l}}]=2C_{u_{j},u_{l}}$, while one should have
$\text{Im}[C_{p_{W,j};p_{W,l}}]\approx 0$.

##  III Equipartion at thermalization of the hidden-variable components

Since the expectation values $\langle K_{B,j}\rangle$ and $\langle
Q_{B,j}\rangle$ cannot be computed from hermitian operators, the typical
argumentation of ETH to discuss thermalization of observables cannot directly
be applied to these quantities, although its very same spirit can still be
used. For that, let us first summarize the standard interpretation of the
thermalization of expectation values.

From an initial nonequilibrium state $|\Psi(0)\rangle$, whose unitary
evolution is dictated by
$|\Psi(t)\rangle=\sum_{n}c_{n}e^{-i\omega_{n}t}|n\rangle$, being $|n\rangle$
an energy eigenstate with eigenvalue $\omega_{n}$, and $c_{n}=\langle
n|\Psi(0)\rangle$ defined by initial conditions, the expectation value of some
hermitian operator $\hat{a}$ reads

$\langle\hat{a}\rangle=\sum_{n}\rho_{n,n}(0)a_{n,n}+\sum_{n,m\neq
n}\rho_{n,m}(t)a_{m,n},$ (32)

with $a_{m,n}=\langle m|\hat{a}|n\rangle$, and
$\rho_{n,m}(t)=c_{m}^{*}c_{n}e^{i(\omega_{m}-\omega_{n})t}$ the density
matrix, composed by diagonal time-independent and nondiagonal time-dependent
terms; while the former can never be neglected (unless is zero by
construction), the latter needs to be negligible after some given time if one
expects $\langle\hat{a}\rangle$ to thermalize. A system is said to equilibrate
if, after some time $t_{eq}$ enough for full dephasing between different
energy eigenstates, the nondiagonal terms cancel out so that (32) can solely
be computed from the diagonal terms,
$\langle\hat{a}\rangle\approx\sum_{n}\rho_{n,n}(0)a_{n,n}$, for most times
$t>t_{eq}$. The system is then said to thermalize when $\langle\hat{a}\rangle$
becomes roughly equal to the expectation value as computed from its
microcanonical density matrix. The ETH states that such dephasing is more
typical to nondegenerated and chaotic many-body scenarios, where the
nondiagonals $a_{m,n}$ in (32) become exponentially smaller than the diagonals
$a_{n,n}$. In other words, a non-equilibrium state whose evolution involves a
large number of eigenstates (not necessarily of particles) is required. Such a
picture immediately applies to the thermalization of the orthodox kinetic
$K(\mathbf{x})$ and potential $V(\mathbf{x})$ operators in (10). For the
hidden-variable components of the Bohmian and Stochastic kinetic energies,
however, one needs an extra step, which renders us the main result of our work
in what follows.

The Bohmian velocity for particle $j$ is, from (3),

$v_{j}(\mathbf{x},t)=\text{Im}\left[\frac{1}{\Psi(\mathbf{x},t)}\frac{\partial\Psi(\mathbf{x},t)}{\partial
x_{j}}\right]=\frac{1}{2i}\left(\frac{\tilde{\Psi}_{j}}{\Psi}-\frac{\tilde{\Psi}^{*}_{j}}{\Psi^{*}}\right),$
(33)

where $\tilde{\Psi}_{j}\equiv\partial\Psi(\mathbf{x},t)/\partial x_{j}$ and
$\Psi\equiv\Psi(\mathbf{x},t)$. The ensemble value of the product $\langle
v_{j}v_{l}\rangle$ is

$\displaystyle\langle v_{j}v_{l}\rangle$ $\displaystyle=$ $\displaystyle\int
d\mathbf{x}|\Psi(\mathbf{x},t)|^{2}\>v_{j}(\mathbf{x},t)v_{l}(\mathbf{x},t)$
$\displaystyle=$ $\displaystyle-\frac{1}{4}\int
d\mathbf{x}\frac{1}{|\Psi|^{2}}\left(\tilde{\Psi}_{j}\Psi^{*}-\tilde{\Psi}_{j}^{*}\Psi\right)\left(\tilde{\Psi}_{l}\Psi^{*}-\tilde{\Psi}_{l}^{*}\Psi\right)$
$\displaystyle=$ $\displaystyle\frac{1}{2}\int
d\mathbf{x}\tilde{\Psi}_{j}\tilde{\Psi}_{l}^{*}$ $\displaystyle-$
$\displaystyle\frac{1}{4}\int
d\mathbf{x}\frac{1}{|\Psi|^{2}}\left(\tilde{\Psi}_{j}\Psi^{*}\tilde{\Psi}_{l}\Psi^{*}+\tilde{\Psi}_{j}^{*}\Psi\tilde{\Psi}_{l}^{*}\Psi\right).$

The first integral is just half of
$\langle\hat{p}_{j}\hat{p}_{l}\rangle=\langle(\hat{p}_{j}\Psi(t))^{*}|\hat{p}_{l}\Psi(t)\rangle=\int
d\mathbf{x}\tilde{\Psi}_{j}\tilde{\Psi}_{l}^{*}$, which then follows the same
thermalization process discussed after (32) for hermitian operators. The
second integral is the sum of one component plus its complex conjugate so that
(III) remains indeed real; from the decomposition
$\Psi(\mathbf{x},t)=\sum_{n}c_{n}e^{-i\omega_{n}t}\psi_{n}(\mathbf{x})$ one
can rewrite

$\displaystyle\int
d\mathbf{x}\frac{1}{|\Psi|^{2}}(\tilde{\Psi}_{j}\Psi^{*}\tilde{\Psi}_{l}\Psi^{*})=$
(35)
$\displaystyle\sum_{a,b,c,d}c_{a}c_{b}^{*}c_{c}c_{d}^{*}e^{-iw_{a,b,c,d}t}\int
d\mathbf{x}\;\frac{1}{|\Psi|^{2}}\tilde{\psi}_{a,j}\psi_{b}^{*}\tilde{\psi}_{c,l}\psi_{d}^{*},$

where $\tilde{\psi}_{s,q}\equiv\partial\psi_{s}(\mathbf{x})/\partial x_{q}$
and $\psi_{s}\equiv\psi_{s}(\mathbf{x})$, for $s=a,b,c,d$, $q=j,l$, and
$w_{a,b,c,d}=w_{a}-w_{b}+w_{c}-w_{d}$. Notice that contrarily to (32) and
(III), there is no time-independent term because the denominator
$|\Psi|^{2}\equiv|\Psi(\mathbf{x},t)|^{2}$ will be time-dependent even when
$w_{a,b,c,d}=0$. Due to the assumed chaotic nature of the eigenstates, one can
infer no spatial correlation between $\psi_{a,j}(\mathbf{x})$,
$\psi_{b}^{*}(\mathbf{x})$, $\tilde{\psi}_{c,l}(\mathbf{x})$,
$\psi_{d}^{*}(\mathbf{x})$; in addition, $|\Psi(\mathbf{x},t)|^{2}$ introduces
randomness in the configuration space as the system thermalizes in the
physical space [47]. With these in mind (35) simply yields a sum of random
numbers around zero with vanishing contribution, so that the thermalized value
of (III) is

$\langle
v_{j}v_{l}\rangle\approx\frac{\langle\hat{p}_{j}\hat{p}_{l}\rangle}{2}.$ (36)

From exactly the same procedure, since the osmotic velocity for particle $j$
is, from (4),

$u_{j}(\mathbf{x},t)=\text{Re}\left[\frac{1}{\Psi(\mathbf{x},t)}\frac{\partial\Psi(\mathbf{x},t)}{\partial
x_{j}}\right]=\frac{1}{2}\left(\frac{\tilde{\Psi}_{j}}{\Psi}+\frac{\tilde{\Psi}^{*}_{j}}{\Psi^{*}}\right),$
(37)

the ensemble value of the product $\langle u_{j}u_{l}\rangle$ will be exactly
the same as in (III), but with a positive sign in the second integral. Thus,
applying the same reasoning as above, one finds the thermalized value

$\langle
u_{j}u_{l}\rangle\approx\frac{\langle\hat{p}_{j}\hat{p}_{l}\rangle}{2},$ (38)

so that (23) remains satisfied after thermalization.

Applying (36) and (38) with $j=l$ into (26)-(28) yields

$\displaystyle\langle\hat{K}_{j}\rangle$ $\displaystyle=$
$\displaystyle\langle K_{B,j}\rangle+\langle Q_{B,j}\rangle,$ (39)
$\displaystyle\frac{\langle\hat{K}_{j}\rangle}{2}$ $\displaystyle\approx$
$\displaystyle\langle K_{B,j}\rangle\approx\langle Q_{B,j}\rangle.$ (40)

While (39) is already true in (26)-(28), and so it is valid at any scenario,
thermalized or not, the kinetic energy equipartition in (40), being only valid
at times $t>t_{eq}$, presents the hidden-variable signature of quantum
thermalization: Bohmian kinetic and quantum potential energies become equal,
with each being equal to half of the orthodox kinetic energy. Similarly, while
$\langle\hat{p}_{j}^{2}\rangle=\langle v_{j}^{2}\rangle+\langle
u_{j}^{2}\rangle$ applies at any time, at $t>t_{eq}$
$\langle\hat{p}_{j}^{2}\rangle/2\approx\langle v_{j}^{2}\rangle\approx\langle
u_{j}^{2}\rangle$ should also apply. That is, thermalization also implies that
(squared) Bohmian and osmotic velocities become equal, with each being equal
to half of the (squared) orthodox momentum; in other words, information from
phase and amplitude of the wave function become similar after $t>t_{eq}$,
which is a result of the initially localized wave function, in a
nonequilibrium dynamics, spreading almost homogeneously through the whole
configuration space after the onset of thermalization. In addition to
(39)-(40), and particularly for a randomly disordered harmonic trap, one
should also find that

$\displaystyle\langle\hat{H}\rangle$ $\displaystyle=$
$\displaystyle\langle\hat{K}\rangle+\langle\hat{V}\rangle,$ (41)
$\displaystyle\frac{\langle\hat{H}\rangle}{2}$ $\displaystyle\approx$
$\displaystyle\langle\hat{K}\rangle\approx\langle\hat{V}\rangle.$ (42)

Equation (41) expresses the conservation of total energy, from the Hamiltonian
in (10) in a unitary evolution, valid at any time. Equation (42) tells us that
the orthodox virial theorem, as one reaches some steady state at $t>t_{eq}$,
is expected to be restated; that is, potential and kinetic energies should
become equal, with each being equal to half of the total energy. It is assumed
in (41)-(42) that $\langle\hat{V}\rangle$ mostly comes from the confining
potential. At last, from (36) and (38) into (30)-(31), the thermalized
correlations should satisfy

$\frac{C_{\hat{p}_{j},\hat{p}_{l}}}{2}\approx C_{v_{j},v_{l}}\approx
C_{u_{j},u_{l}}\approx\frac{C_{p_{W,j};p_{W,l}}}{2}.$ (43)

##  IV Numerical results

We first describe the non-equilibrium initial state and the Hamiltonian of our
model; then we discuss our results according to three scenarios with distinct
initial conditions. At the end we reformulate our model in the center of mass
frame and address the respective results.

###  IV.1 Initial state and trap Hamiltonian

The initial $N$-electron non-equilibrium pure 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),$
(44)

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 (44) is

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

with spatial dispersion $\sigma_{j}$, central position $x_{0j}$, and central
velocity $p_{0j}$. Any nonzero $x_{0j}$ or $p_{0j}$ may activate the
nonequilibrium dynamics.

The Hamiltonian $H(\mathbf{x})$ propagating the many-body wavefunction
$\Psi(\mathbf{x},t)$ is given in (10), where the kinetic term $K(\mathbf{x})$
is already defined. The potential term $V(\mathbf{x})$, in our disordered
harmonic trap, is defined by

$V(\mathbf{x})\equiv V_{H}(\mathbf{x})+V_{I}(\mathbf{x})+V_{D}(\mathbf{x}),$
(46)

where the harmonic potential $V_{H}(\mathbf{x})$ with frequency $\omega$ is

$V_{H}(\mathbf{x})=\frac{1}{2}\omega^{2}\sum_{j=1}^{N}x_{j}^{2},$ (47)

and the electron-electron interaction potential $V_{I}(\mathbf{x})$ with
smooth parameter $\alpha$ is

$V_{I}(\mathbf{x})=\frac{1}{2}\sum_{j=1}^{N}\sum_{l\neq
j}^{N}\frac{1}{\sqrt{(x_{j}-x_{l})^{2}+\alpha^{2}}}.$ (48)

To ensure that the initial state in (44) is built as a superposition of a
large (and ‘chaotic’) number of eigenstates, we include the random disorder
potential $V_{D}(\mathbf{x})$,

$V_{D}(\mathbf{x})=\gamma_{D}\sum_{j=1}^{N}\sum_{l=1}^{M}b_{l}\exp[-\frac{4(x_{j}-g_{l})^{2}}{\sigma_{D}^{2}}],$
(49)

where $\gamma_{D}$ is its strength and $\sigma_{D}$ its spatial dispersion,
with $g_{l}$ running through $M$ grid points; the set of random numbers
$b_{l}$ satisfies $\langle b_{l}\rangle=0$ and $\langle b_{l}^{2}\rangle=1$,
and the disorder potential is normalized so that the integral of
$V^{2}_{D}(\mathbf{x})$ yields $\gamma_{D}^{2}$. Such a shape is typical of
speckle potentials [66, 67, 68, 69, 70, 71, 72] in fermionic traps.

###  IV.2 Expectations values in three scenarios

Our results initially focus on $N=2$. However our main conclusions, as the
kinetic energy equipartition at thermalization, are valid for any $N$ as the
discussion in the center of mass frame will latter settle. We focus on three
scenarios, D1, D2, D3, which are distinct by the initial values of $x_{0j}$
and $p_{0j}$, as summarized in Table 1, where each scenario considers both no-
disordered and disordered dynamics. Values of $x_{0j}$ in D1 and $p_{0j}$ in
D2 are chosen identical as to yield the same turning points in both dynamics;
the remaining simulation parameters are found in [73]. The no-disorder cases
employ smaller values of $x_{0j}$ and $p_{0j}$ to render the features more
visible.

The dynamics for each scenario D1, D2, D3 is respectively shown in figures 1,
2, 3: top panels (a)-(d) show a few initial cycles with no disorder, bottom
panels (e)-(h) show the full evolution with disorder. The structure of these
three figures, which focus on the time evolution of some pertinent expectation
values, is: panels (a), (e) show kinetic $\langle K\rangle$, harmonic $\langle
V_{H}\rangle$, interaction $\langle V_{I}\rangle$, and total $\langle
H\rangle=\langle K\rangle+\langle V_{H}\rangle+\langle V_{I}\rangle+\langle
V_{D}\rangle$ energies, with disorder energy $\langle V_{D}\rangle\approx 0$
at any $t$ not shown; panels (b), (f) show Bohmian velocity $\langle
v_{j}\rangle$, position $\langle x_{j}\rangle$, and osmotic velocity $\langle
u_{j}\rangle$, with momentum $\langle p_{j}\rangle=\langle v_{j}\rangle$ not
shown, and the label $j$ not indicated since it is redundant in our
antisymmetrized model; panels (c), (g) repeat the orthodox kinetic energy
$\langle K\rangle$, for its comparison with the Bohmian kinetic $\langle
K_{B}\rangle$ and quantum potential $\langle Q_{B}\rangle$ energies; panels
(d), (h) show the correlations for momentum $C_{\hat{p}_{1},\hat{p}_{2}}$,
Bohmian $C_{v_{1},v_{2}}$ and osmotic $C_{u_{1},u_{2}}$ velocities, and
position $C_{\hat{x}_{1},\hat{x}_{2}}$.

Scenario | D1 | D2 | D3
---|---|---|---
Disorder | No | Yes | No | Yes | No | Yes
($x_{01}$,$x_{02}$) | (-2,2) | (-20,20) | (-2,2) | (-2,2) | (-2,2) | (-20,20)
($p_{01}$,$p_{02}$) | (0,0) | (0,0) | (4,4) | (20,20) | (2,2) | (20,20)
($\langle x\rangle$, $\langle p\rangle$) |  |  |  |  |  |
($\langle K\rangle$,$\langle V_{H}\rangle$) |  |  |  |  |  |
($\langle K_{B}\rangle$,$\langle Q_{B}\rangle$) |  |  |  |  |  |
Table 1:  The three scenarios D1, D2, D3 for the dynamics with and without
disorder for $N=2$. Initial two lines: values of ($x_{01}$,$x_{02})$ and
($p_{01}$,$p_{02})$ for the initial non-equilibrium state. Last three lines:
features of expectation values of ($\langle x\rangle$, $\langle p\rangle$),
and of orthodox $(\langle K\rangle$, $\langle V_{H}\rangle$) and Bohmian
($\langle K_{B}\rangle$, $\langle Q_{B}\rangle$) energies, where symbols  and
respectively indicate whether or not such expectation values allow one to
identify $t_{eq}$.

####  IV.2.1 Dynamics from initial position

Figure 1: Expectation values from the dynamics in the D1 scenario, with
$(p_{01},p_{02})=(0,0)$. Top panels: few initial cycles under no disorder with
a smaller $(x_{01},x_{02})=(-2,2)$; Bottom panels: full dynamics under
disorder with a larger $(x_{01},x_{02})=(-20,20)$. Panels (a), (e) for
energies: Orthodox kinetic energy $\langle K\rangle$, harmonic confining
energy $\langle V_{H}\rangle$, Coulomb repulson energy $\langle V_{I}\rangle$,
total energy $\langle H\rangle$ ($\langle V_{I}\rangle$ is magnified by $100$
in (e)). Panels (b), (f) for velocities: Bohmian velocity $\langle v\rangle$,
position $\langle x\rangle$, osmotic velocity $\langle u\rangle$. Panels (c),
(g) for kinetic energies comparison: Orthodox $\langle K\rangle$, Bohmian
$\langle K_{B}\rangle$, quantum potential $\langle Q_{B}\rangle$ energies.
Panels (d), (h) for correlations: momentum $C_{\hat{p}_{1},\hat{p}_{2}}$,
Bohmian $C_{v_{1},v_{2}}$ and osmotic $C_{u_{1},u_{2}}$ velocities, position
$C_{\hat{x}_{1},\hat{x}_{2}}$.

Let us first focus on the D1 dynamics without disorder in Fig. 1(a)-(d), which
has $p_{01}=p_{02}=0$ and a small $x_{01}=-x_{02}=-2$. At $t=0$, 1(a) yields
$\langle V_{H}\rangle(0)=x_{01}^{2}/2+x_{02}^{2}/2+E_{V}/2=4.5$ and $\langle
K\rangle(0)=p_{01}^{2}/2+p_{02}^{2}/2+E_{V}/2=0.5$, where $E_{V}=\omega(n+1)$
is the ground state energy ($\omega=1$, $n=0$) which is equally shared among
potential and kinetic terms. That is, the virial theorem in (42) is not
satisfied (at any $t$) thanks to the non-equilibrium initial situation;
however, the total energy in (41), $\langle H\rangle(0)=5.27$ since $\langle
V_{I}\rangle(0)=0.27$, is conserved (at any $t$) thanks to the unitary
evolution. In 1(c), at $t=0$, one has $\langle Q_{B}\rangle(0)=0.5$ and
$\langle K_{B}\rangle(0)=0$ (since the initial velocities are zero, see 1(b)),
so that indeed $\langle K\rangle(0)=0.5$ and (39) is satisfied (at any $t$),
while (40) does not apply here. Each D1-cycle has a $\pi$-period and three
stages. For the first cycle: (i) at $t=0$ electrons are at
$(x_{1},x_{2})=(-2,2)$, with minimum $\langle K\rangle$ and maximum $\langle
V_{H}\rangle$; (ii) the dynamics pushes the electrons against each other until
that, at $t=\pi/2$, they try binding together at $(x_{1},x_{2})=(0,0)$, which
is avoided thanks to both exchange symmetry, Coulomb repulsion, and quantum
potential, as indicated by the peak in $\langle V_{I}\rangle$ and in $\langle
Q_{B}\rangle$, with maximum $\langle K\rangle$ and minimum $\langle
V_{H}\rangle$ and $\langle K_{B}\rangle$; (iii) at $t=\pi$ electrons are back
to $(x_{1},x_{2})=(-2,2)$ and a new cycle starts. Notice that $\langle
K_{B}\rangle$ has a double peak around the peak of $\langle Q_{B}\rangle$
because the velocity acquires a first maximum from $t=0$ to $t=\pi/2$, which
is the time electrons ‘stop’ to reverse their movements, and a second maximum
from $t=\pi/2$ to $t=\pi$. One sees how the quantum potential, acting ‘in-
phase’ with the Coulomb repulsion, carries the quantumness of the two-body
entanglement. The correlations in 1(d) are a mirror of the above discussions;
they are all negative since, as one variable increases, the other decreases in
the D1-dynamics. In terms of moduli, at $t=0$, $C_{\hat{x}1,\hat{x}2}$ is at
its maxima due to the maximum electron separation, and reaches its vanishing
minima at $t=\pi/2$ when electrons are closest to each other. The three
kinetic correlations are zero at $t=0$ thanks to the same initial zero
velocity for the electrons; at $t=\pi/2$, since at that time it is only
$\langle Q_{B}\rangle$ that contributes to $\langle K\rangle$, $C_{u1,u2}$
reaches its maximum while $C_{v1,v2}$ vanishes. Notice that (30) is satisfied
at any $t$, while (43) does not apply here. In 1(b), $\langle v\rangle=\langle
x\rangle=\langle u\rangle=0$ at any $t$ is a trivial consequence of the D1
dynamics being antidiagonal in the configuration space $x_{1}$$x_{2}$, with
electrons initially equidistant, what anticipates that such terms reflect
center of mass properties.

Figure 2: Expectation values from the dynamics in the D2 scenario, with
$(x_{01},x_{02})=(-2,2)$. Top panels: few initial cycles under no disorder
with a smaller $(p_{01},p_{02})=(4,4)$; Bottom panels: full dynamics under
disorder with a larger $(p_{01},p_{02})=(20,20)$. Panels (a), (e) for
energies: Orthodox kinetic energy $\langle K\rangle$, harmonic confining
energy $\langle V_{H}\rangle$, Coulomb repulson energy $\langle V_{I}\rangle$,
total energy $\langle H\rangle$ ($\langle V_{I}\rangle$ is magnified by $100$
in (e)). Panels (b), (f) for velocities: Bohmian velocity $\langle v\rangle$,
position $\langle x\rangle$, osmotic velocity $\langle u\rangle$. Panels (c),
(g) for kinetic energies comparison: Orthodox $\langle K\rangle$, Bohmian
$\langle K_{B}\rangle$, quantum potential $\langle Q_{B}\rangle$ energies.
Panels (d), (h) for correlations: momentum $C_{\hat{p}_{1},\hat{p}_{2}}$,
Bohmian $C_{v_{1},v_{2}}$ and osmotic $C_{u_{1},u_{2}}$ velocities, position
$C_{\hat{x}_{1},\hat{x}_{2}}$.

The full D1 dynamics with disorder is shown in Fig. 1(e)-(h). Thanks to the
larger $x_{01}=-x_{02}=-20$ the oscillation amplitudes are larger but with the
same $\pi$-period. The magnitude of $\langle V_{I}\rangle$ though remain
similar to the no-disorder case, such that we magnify it by $100$ here, while
the initial peaks in $\langle V_{I}\rangle$ and in $\langle Q_{B}\rangle$ as
well as the minima in $\langle K_{B}\rangle$ become steeper in the disordered
scenario. As the cycles succeed all expectation values overall reach quasi-
stationary magnitudes once thermalization is set at $t_{eq}\approx 55$. In
1(e), while maintaining $\langle H\rangle$ constant as in (41), kinetic and
potential energies interchange their magnitudes until (42) becomes valid, and
the virial theorem looks restated after thermalization, seeming an indication
that a steady state is reached; the Coulomb repulsion $\langle V_{I}\rangle$
smears out and its peaks disappear after thermalization. In 1(g), while (39)
remains true at any $t$, the peaks in $\langle Q_{B}\rangle$ and minima in
$\langle K_{B}\rangle$ also smear out and disappear after thermalization; most
interestingly, this happens in a way as to satisfy (40), allowing one to
visualize the central result of our paper: in addition to the ‘restatement’ of
the virial theorem, thermalization also implies a kinetic energy
equipartition, the hidden-variable signature of quantum thermalization. The
trivial values of $\langle x\rangle=\langle v\rangle=\langle u\rangle\approx
0$ remain at any $t$ in 1(f), and so presents no feature to identify
thermalization. The correlations in 1(h) are once more a mirror of the above
discussions, so the maxima and minima of the initial cycles smear out after
thermalization is set, when one also finds that, in addition to (30), equation
(43) is also satisfied, as well as $C_{\hat{p}_{1},\hat{p}_{2}}\approx
C_{\hat{x}_{1},\hat{x}_{2}}$; the correlations remain negative, reaching small
magnitudes after thermalization.

####  IV.2.2 Dynamics from initial velocity

Figure 3: Expectation values from the dynamics in the D3 scenario, with mixed
$x_{0}$, $p_{0}$ influence. Top panels: few initial cycles under no disorder
with small $(x_{01},x_{02})=(-2,2)$ and small $(p_{01},p_{02})=(2,2)$; Bottom
panels: full dynamics under disorder with a large $(x_{01},x_{02})=(-20,20)$
and a large $(p_{01},p_{02})=(20,20)$. Panels (a), (e) for energies: Orthodox
kinetic energy $\langle K\rangle$, harmonic confining energy $\langle
V_{H}\rangle$, Coulomb repulson energy $\langle V_{I}\rangle$, total energy
$\langle H\rangle$ ($\langle V_{I}\rangle$ is magnified by $100$ in (e)).
Panels (b), (f) for velocities: Bohmian velocity $\langle v\rangle$, position
$\langle x\rangle$, osmotic velocity $\langle u\rangle$. Panels (c), (g) for
kinetic energies comparison: Orthodox $\langle K\rangle$, Bohmian $\langle
K_{B}\rangle$, quantum potential $\langle Q_{B}\rangle$ energies. Panels (d),
(h) for correlations: momentum $C_{\hat{p}_{1},\hat{p}_{2}}$, Bohmian
$C_{v_{1},v_{2}}$ and osmotic $C_{u_{1},u_{2}}$ velocities, position
$C_{\hat{x}_{1},\hat{x}_{2}}$.

Let us start the discussion of the D2 dynamics once more from the no-disorder
case in Fig. 2(a)-(d), which has $x_{01}=-x_{02}=-2$ (as $(0,0)$ cannot be
used in a fermionic trap) and a small $p_{01}=p_{02}=4$ (such that both
electrons start moving in the same direction). As such, at $t=0$, one still
has $\langle V_{H}\rangle(0)=4.5$ and $\langle V_{I}\rangle(0)=0.27$, but a
larger $\langle K\rangle(0)=16.5$, built from $\langle K_{B}\rangle(0)=16$ and
$\langle Q_{B}\rangle(0)=0.5$, so yielding $\langle H\rangle(0)=21.27$. Each
D2-cycle has a $2\pi$-period and five stages. The first cycle is: (i) at $t=0$
electrons at $(x_{1},x_{2})=(-2,2)$ with minimum $\langle V_{H}\rangle$ and
maximum $\langle K\rangle$; (ii) at $t=\pi/2$ electrons reach the positive
turning point and try to collide at $(x_{1},x_{2})=(4,4)$, yielding peaks in
$\langle V_{I}\rangle$ and in $\langle Q_{B}\rangle$, while $\langle
K_{B}\rangle=0$ as the electrons stop at that time (the fact that the peak in
$\langle Q_{B}\rangle$ is now at the minimum of $\langle K\rangle$ is the
reason why a double peak is no longer seen in $\langle K_{B}\rangle$
contrarily to the D1-scenario); (iii) at $t=\pi$ electrons pass back at
$(x_{1},x_{2})=(-2,2)$; (iv) at $t=3\pi/2$ the electrons reach the negative
turning point and try to collide at $(x_{1},x_{2})=(-4,-4)$, with new peaks in
$\langle V_{I}\rangle$ and in $\langle Q_{B}\rangle$ and again with $\langle
K_{B}\rangle=0$; (v) at $t=2\pi$ the electrons are back to
$(x_{1},x_{2})=(-2,2)$ and a new cycle starts. The cycle shows that, although
the initial velocity is the same for both electrons, they acquire different
velocities in the dynamics as one moves in favour and the other against the
harmonic potential. Results in 2(b) are trivially understood from the previous
analysis, confirming the center of mass character of those quantities: the D2
dynamics, being diagonal in the configuration space $x_{1}$$x_{2}$, has that
$\langle x\rangle$ ranges from $\langle x\rangle=0$ when electrons are at
$(2,-2)$ to $\langle x\rangle=4$ ($\langle x\rangle=-4$) when electrons are at
the positive $(4,4)$ (negative $(-4,-4)$) turning point, while $\langle
v\rangle$ is just out-of-phase and $\langle u\rangle=0$ from (22).
Interestingly the correlations in 2(d) do not change from the D1 dynamics. So
let us explain the origin of the double-peak feature in $C_{v_{1},v_{2}}$: at
$t=0$, $C_{v_{1},v_{2}}=0$ since both electrons have same initial velocity to
the right; as time evolves, the left (right) electron gains (loses) velocity,
so negatively increasing $C_{v_{1},v_{2}}$, which reaches its first peak as
the left electron crosses the origin; then the left electron also starts to
decrease its velocity until both electrons reach the positive turning point,
causing $C_{v_{1},v_{2}}=0$ at $t=\pi/2$; electrons then start moving to the
left, again with a higher velocity for the left electron, inducing a new
negative increase in $C_{v_{1},v_{2}}$, until the left electron passes by the
origin at its highest velocity, causing the second peak in $C_{v_{1},v_{2}}$.

The full D2 dynamics with disorder is seen in Fig. 2(e)-(h), in which the
larger $p_{01}=p_{02}=20$ induces larger oscillation amplitudes with the same
$2\pi$-period, except for $\langle V_{I}\rangle$ which is once more magnified
by $100$. Overall, the smearing out of the oscillations in $\langle K\rangle$
and $\langle V_{H}\rangle$, of the peaks in $\langle V_{I}\rangle$ and
$\langle Q_{B}\rangle$, and of the minima in $\langle K_{B}\rangle$, happens
as the system thermalizes, with such quantities seeming to reach stationary
values at $t_{eq}\approx 220$. Thermalization in the D2 scenario also induces
the feature $\langle x\rangle\approx\langle v\rangle\approx 0$ in 2(f)
(compare with 1(f)). All four equations (39)-(42) are also satisfied in the D2
scenario at $t>t_{eq}$, such that the virial theorem seems restated in 2(e),
and the hidden-variable signature of thermalization is once more seen in 2(g).
While (30) and (43) are also valid in the D2 scenario, as well as
$C_{\hat{p}_{1},\hat{p}_{2}}\approx C_{\hat{x}_{1},\hat{x}_{2}}$, the
correlations in 2(h) become positive (since both electrons move in the same
direction) and stabilize at larger magnitudes (from the larger initial
velocity).

####  IV.2.3 Dynamics from initial position and initial velocity

After detailing the scenarios D1 and D2 we will only describe here the
distinct features of the D3 dynamics. Once more starting from the no-disorder
case in Fig. 3(a)-(d), which has both a small $p_{01}=p_{02}=2$ and a small
$x_{01}=-x_{02}=-2$. Such identical initial values obviously yield $\langle
V_{H}\rangle(0)=\langle K\rangle(0)=4.5$, and since $\langle
V_{I}\rangle(0)=0.27$ one has $\langle H\rangle(0)=9.27$, while $\langle
K_{B}\rangle(0)=4$ and $\langle Q_{B}\rangle(0)=0.5$ are the values building
$\langle K\rangle(0)$. The dynamics though is more involved since it is a
superposition of the two previous scenarios, with cycles in a $2\pi$-period.
At the initial cycles, $\langle K\rangle$ looks almost constant since $\langle
Q_{B}\rangle$ reaches magnitudes comparable to $\langle K_{B}\rangle$ (compare
3(c) with 2(c) and 1(c)), as $\langle K_{B}\rangle$ may still vanish whenever
$\langle v\rangle=0$ in 3(b). The correlation plots in 3(d), interestingly,
are identical to 2(d) and 1(d), while $\langle v\rangle$ and $\langle
x\rangle$ in 3(b) oscillate in between the turning points at $(2,2)$ and
$(-2,-2)$, with center-of-mass-like values.

The full D3 dynamics with disorder is shown in Fig. 3(e)-(h), which considers
both a large $p_{01}=p_{02}=20$ and a large $x_{01}=-x_{02}=-20$. With such
values the constancy of the kinetic and potential energies in 3(e) becomes
evident such that, besides the trivial energy conservation in (41), the virial
theorem in (42) looks ‘as if’ is always satisfied at any $t$, even before
thermalization is set at $t_{eq}\approx 110$, mistakenly telling that some
stationary state could be present since $t=0$. In other words, the onset of
thermalization is no longer identifiable in the energy expectation values of
3(e) (where $\langle V_{I}\rangle$ is magnified again by $100$), although it
remains identifiable in 3(f). The central result of our paper, however, is
once more found in 3(g), where the kinetic energy equipartition in (40) is
verified although following a different thermalizing path if compared to 2(g)
and 1(g). The correlations in 3(h) are about the same as in the D1 scenario,
where (43), as well as $C_{\hat{p}_{1},\hat{p}_{2}}\approx
C_{\hat{x}_{1},\hat{x}_{2}}$, apply once more.

###  IV.3 Revisiting the three scenarios within the center of mass framework

####  IV.3.1 The issue of empirical indistinguishibility

All expectation values presented in the previous figures refer to a single
particle but, thanks to the antisymmetry of $\Psi(\mathbf{x},t)$, there is no
need to identify ‘which’ particle. This is true no matter they are linked to
Hermitian operators like $\hat{K}$ or $\hat{p}$ or to post-processings of weak
values like $p_{W}(\mathbf{x},t)$, as the exchange symmetry of both phase and
amplitude of $\Psi(\mathbf{x},t)$ is directly translated, respectively, to the
current and osmotic velocities in (3)-(4) (and to the post-processing data).
However, if one is interested on how expectation values like $\langle
v\rangle$ or $\langle u\rangle$ are obtained in a laboratory, one needs to
differentiate between the ones computed from hermitian operators to the ones
computed from post-processings of weak values, as indicated in Table 2.

Expectation | Eq. | Operator | From weak values
---|---|---|---
$\langle p\rangle$ | (20) | $\hat{p}$ | $p_{W}(\mathbf{x},t)$
$\langle K\rangle$ | (28) | $\hat{K}$ | Re$[p_{W}(\mathbf{x},t)]^{2}$+Im$[p_{W}(\mathbf{x},t)]^{2}$
$\langle v\rangle$ | (21) |  | Re$[p_{W}(\mathbf{x},t)]$
$\langle u\rangle$ | (22) |  | Im$[p_{W}(\mathbf{x},t)]$
$\langle K_{B}\rangle$ | (26) |  | Re$[p_{W}(\mathbf{x},t)]^{2}$
$\langle Q_{B}\rangle$ | (27) |  | Im$[p_{W}(\mathbf{x},t)]^{2}$
Table 2:  The main expectation values relevant in the discussion of the
kinetic energy equipartition, addressed in the previous figures. Expectation
values marked with  cannot directly be computed from an Hermitian operator,
but they can be accessed by post-processing the local-in-position weak value
of the momentum $p_{W}(\mathbf{x},t)$.

Figure 4: Center of mass and relative energy expectation values for scenarios
D1, D2, and D3, respectively, in left ((a),(d),(g)), middle ((b),(e),(h)), and
right ((c),(f),(i)) panels. Upper panels: kinetic energy $\langle
K_{c}\rangle$ and $\langle K_{r}\rangle$, confining potential energy $\langle
V_{H,c}\rangle$ and $\langle V_{H,r}\rangle$), respective total energy
$E_{c}=\langle K_{c}\rangle+\langle V_{H,c}\rangle$ and $E_{r}=\langle
K_{r}\rangle+\langle V_{H,r}\rangle$ (black lines); $\langle V_{I,r}\rangle$
is not included. Middle panels: Bohmian kinetic $\langle K_{B,c}\rangle$ and
quantum potential $\langle Q_{B,c}\rangle$ energies in comparison with
orthodox $\langle K_{c}\rangle$ center of mass energies. Lower panels: Bohmian
kinetic $\langle K_{B,r}\rangle$ and quantum potential $\langle
Q_{B,r}\rangle$ energies in comparison with orthodox $\langle K_{r}\rangle$
relative energies. Legend in (c) applies to (a),(b); legend in (f) applies to
(d),(e); legend in (i) applies to (g),(h).

Let us ilustrate the problem with $N=2$. The experimental weak value
$p_{W,1}(x_{1},x_{2},t)$ requires two measurements: a first weak measurement
linked to $\hat{p}_{1}$, plus a second strong measurement linked to
$\hat{x}_{1}$ and $\hat{x}_{2}$. Such a weak value is obtained after a pre-
selection (repeating the same experiment for a large number of identically
prepared initial wave functions) plus a post-selection (to obtain an ensemble
value of all weakly measured momentum whose subsequent strongly measured
position yields a specific location); a similar procedure is required for
getting $p_{W,2}(x_{1},x_{2},t)$ by interchanging the measurements on
particles $1$ and $2$. Although $\langle p_{W,1}\rangle=\int
dx_{1}\;dx_{2}p_{W,1}(x_{1},x_{2},t)|\Psi(x_{1},x_{2},t)|^{2}$ is the same as
$\langle p_{W,2}\rangle=\int
dx_{2}\;dx_{1}p_{W,2}(x_{2},x_{1},t)|\Psi(x_{2},x_{1},t)|^{2}$, their
empirical evaluation requires one to identify which are the particles $1$ and
$2$ being weakly or strongly measured, which is impossible for systems with
indistinguishable particles. One should keep in mind though that while Bohmian
particles are ontologically distinguishable, the Bohmian dynamical laws ensure
that expectations values are empirically indistinguishable.

To handle this issue we consider an effective single-particle weak value,
$\tilde{p}_{W,j}(x_{k},t)$, by defining
$\mathbf{X_{k}}=(x_{1},..,x_{k-1},x_{k+1},..,x_{N})$ such that
$\mathbf{x}=(x_{k},\mathbf{X_{k}})$, as

$\displaystyle\tilde{p}_{W,j}(x_{k},t)$ $\displaystyle=$
$\displaystyle\frac{\int
d\mathbf{X_{k}}\;p_{W,j}(x_{k},\mathbf{X_{k}},t)|\Psi(x_{k},\mathbf{X_{k}},t)|^{2}}{P(x_{k},t)},$
(50) $\displaystyle P(x_{k},t)$ $\displaystyle=$ $\displaystyle\int
d\mathbf{X_{k}}\;|\Psi(x_{k},\mathbf{X_{k}},t)|^{2}.$ (51)

It is straightforward to show that

$\displaystyle\langle\tilde{p}_{W,j}\rangle$ $\displaystyle=$
$\displaystyle\int dx_{k}\;\tilde{p}_{W,j}(x_{k},t)P(x_{k},t)\rangle$ (52)
$\displaystyle=$ $\displaystyle\int
dx_{k}d\mathbf{X_{k}}\;p_{W,j}(x_{k},\mathbf{X_{k}},t)|\Psi(x_{k},\mathbf{X_{k}},t)|^{2}$
$\displaystyle=$ $\displaystyle\langle p_{W,j}\rangle.$

Here we only need the case $j=k$; the general case $j\neq k$ is found
elsewhere [47]. We remind that $\langle\tilde{p}_{W,k}\rangle=\langle
p_{W,k}\rangle$ does not imply that $\tilde{p}_{W,k}(x_{k},t)$ can promptly be
used in place of $p_{W,k}(x_{k},\mathbf{X_{k}},t)$ to obtain expectation
values like $\langle v\rangle$ or $\langle u\rangle$. This is only the case
when the degree $x_{k}$ is decoupled from the rest $N-1$ degrees in
$\mathbf{X_{k}}$, that is, when one has
$\Psi(x_{k},\mathbf{X_{k}},t)=\psi_{k}(x_{k},t)\psi_{\mathbf{k}}(\mathbf{X_{k}},t)$
. Then, it is straightforward to show that
$p_{W,k}(x_{k},\mathbf{X_{k}},t)=\tilde{p}_{W,k}(x_{k},t)$ with
$P(x_{k},t)=|\psi_{k}(x_{k},t)|^{2}$, so that the evaluation of the weak value
linked to the weak momentum and position of the particle $x_{k}$ does not
require the measurements of the others positions $\mathbf{X_{k}}$, but, still
requires to identify the particle $x_{k}$. This last difficulty can be solved
by dealing with the center of mass frame.

It is well-know [74, 75, 76, 77, 78, 79] that our trap Hamiltonian in the
absence of the disorder potential $V_{D}(\mathbf{x})$, that is, when
$H(\mathbf{x})\equiv K(\mathbf{x})+V_{H}(\mathbf{x})+V_{I}(\mathbf{x})$, can
be written in terms of center of mass $x_{c}$ and relative
$\mathbf{X_{r}}=(x_{r,1},x_{r,2},..,x_{r,N-1})$ coordinates, where
$x_{c}=\sum_{i=1}^{N}x_{i}/N$ and $x_{r,j}=x_{j}-x_{j+1}$. In such a situation
(top panels in all previous figures), the discussion above promptly applies by
taking $x_{c}\equiv x_{k}$ and $\mathbf{X_{r}}\equiv\mathbf{X_{k}}$. Then one
indeed could consider $\tilde{p}_{W,c}(x_{c},t)$, which one has experimental
access, instead of the intricate $p_{W,c}(x_{c},\mathbf{X_{r}},t)$, to compute
quantities from Table 2 like
$\tilde{v}_{c}^{2}(x_{c},t)=\text{Re}(\tilde{p}_{W,c}(x_{c},t))^{2}$ and
$\tilde{u}_{c}^{2}(x_{c},t)=\text{Im}(\tilde{p}_{W,c}(x_{c},t))^{2}$, instead
of $v_{c}^{2}(x_{c},\mathbf{X_{r}},t)$ or $u_{c}^{2}(x_{c},\mathbf{X_{r}},t)$.
One would also have $H=H_{c}+H_{\mathbf{r}}$ and
$\Psi(x_{c},\mathbf{X_{r}},t)=\psi_{c}(x_{c},t)\psi_{\mathbf{r}}(\mathbf{X_{r}},t)$,
such that $\tilde{p}_{W,c}(x_{c},t)=p_{W,c}(x_{c},\mathbf{X_{r}},t)$ in (50)
and $P(x_{c},t)=|\psi_{c}(x_{c},t)|^{2}$ in (51). That is, a weak measurement
of the momentum of the center of mass followed by a strong measurement of its
position, without measuring the other $N-1$ degrees, should yield information
about the $N$-body trap. The kinetic energy equipartition discussed so far
should remain valid in any frame of reference. So the question one needs to
pose now is whether or not the inclusion of disorder, which drives the system
towards thermalization (bottom panels in all previous figures), induces a
coupling between center of mass and relative coordinates.

####  IV.3.2 Results in the center of mass frame

To answer that question we once more make use of $N=2$, with
$\mathbf{X_{r}}=(x_{r,1})\equiv x_{r}$. The disordered trap Hamiltonian
becomes simply

$H(x_{c},x_{r})=H_{c}(x_{c})+H_{r}(x_{r})+V_{D}(x_{c},x_{r}),$ (53)

with

$\displaystyle H_{c}(x_{c})$ $\displaystyle=$
$\displaystyle-\frac{1}{4}\frac{\partial^{2}}{\partial
x_{c}^{2}}+\omega^{2}x_{c}^{2},$ (54) $\displaystyle H_{r}(x_{r})$
$\displaystyle=$ $\displaystyle-\frac{\partial^{2}}{\partial
x_{r}^{2}}+\frac{1}{4}\omega^{2}x_{r}^{2}+\frac{1}{|x_{r}+\alpha|};$ (55)

With $V_{D}(x_{c},x_{r})=0$, one would obtain
$\Psi(x_{c},x_{r},t)=\psi_{c}(x_{c},t)\psi_{r}(x_{r},t)$ with
$i\partial\psi_{c}(x_{c},t)/\partial t=H_{c}\psi_{c}(x_{c},t)$ and
$i\partial\psi_{r}(x_{r},t)/\partial t=H_{r}\psi_{r}(x_{r},t)$. Figure 4 shows
the dynamics for the scenarios D1, D2, D3, respectively, at left, middle,
right panels, in which $\langle V_{H,c}\rangle$ and $\langle K_{c}\rangle$
($\langle V_{H,r}\rangle$ and $\langle K_{r}\rangle$) label the confining
potential and kinetic components of (54) (of (55)); the component related to
the Coulomb term $\langle V_{I,r}\rangle$ in (55) is not shown since, as we
have seen in all previous figures, it remains small. The labels $\langle
K_{B,c}\rangle=\langle(v_{1}+v_{2})^{2}\rangle/4$ and $\langle
K_{B,r}\rangle=\langle(v_{1}-v_{2})^{2}\rangle/4$, as well as $\langle
Q_{B,c}\rangle=\langle(u_{1}+u_{2})^{2}\rangle/4$ and $\langle
Q_{B,r}\rangle=\langle(u_{1}-u_{2})^{2}\rangle/4$, respectively for the center
of mass and relative components of both Bohmian kinetic and quantum potential
energies, are also employed.

The upper panels (Fig. 4-(a),(b),(c)) show that, indeed, both center of mass
and relative coordinates thermalize, that is, $\langle V_{H,r}\rangle$ and
$\langle K_{r}\rangle$ on one hand, while $\langle V_{H,c}\rangle$ and
$\langle K_{c}\rangle$ on the other hand, reach about same stationary values
at the same $t_{eq}$, with relative energies being higher in scenarios D1 and
D3; the opposite happens in the D2 dynamics, where the relative energies also
present a much smaller $t_{eq}$. The black lines, which present the total
energies $E_{c}=\langle K_{c}\rangle+\langle V_{H,c}\rangle$ and
$E_{r}=\langle K_{r}\rangle+\langle V_{H,r}\rangle$, show that, in fact,
center of mass and relative coordinates are not decoupled at $t<t_{eq}$, since
their respective energies do not remain constant in the time evolution. So, in
principle, the picture made in the previous subsection should not apply.
However, such energies do stabilize after $t_{eq}$, seemingly indicating that
thermalization implies a negligible dependence of $x_{c}$ on $x_{r}$. This
should not come as a surprise since, first, the physics should not be
different by changing a frame of reference; second, as the disorder potential
is random and does not privilege any degree of freedom, it seems natural to
expect that, while coupling is to be found at the beginning, the time
evolution should homogeneously spread the probabilities over the configuration
space, no matter if among $x_{1}$ and $x_{2}$ or $x_{c}$ and $x_{r}$. So, the
empirical measurements as pictured in the previous subsection do apply at
thermalization and, more importantly, the negligible dependence of $x_{c}$ on
the relative coordinates becomes an even more robust result for thermalized
systems with large $N$ [80]. The middle (Fig. 4-(d),(e),(f)) and lower (Fig.
4-(g),(h),(i)) panels, respectively, confirm that the kinetic energy
equipartition remains verified for the center of mass and relative
coordinates, for each of the three scenarios, and it should also be verified
for any $N$ in this new reference frame. One also has to remind, though, that
thermalization has been studied in small systems with as little as 6 [18], 5
[6], or 2-4 bosons [81, 82], and even single-particle systems [83, 84, 85].

##  V Conclusions

Weak values have gradually been transitioning from a theoretical curiosity to
a practical tool in the laboratory allowing novel characterizations of quantum
systems, as they can provide information beyond the traditional expectation
values linked to Hermitian operators. In particular, we have shown that weak
values of the momentum post-selected in positions, without linking the
discussion to any specific ontology but re-using the mathematical machinery of
the Bohmian and Stochastic quantum theories, can be used as a relevant tool to
characterize quantum thermalization in closed systems. As example, we have
addressed the monopole oscillations in the configuration space of a two-
electron harmonic trap under random disorder, under different initial
conditions employed to initiate three distinct nonequilibrium dynamics.

The expectation values from the orthodox operators not always can be employed
to access the onset of thermalization. For example, $\langle v\rangle$ and
$\langle x\rangle$ in Fig. 1(f) for the D1 scenario, or $\langle K\rangle$ and
$\langle V_{H}\rangle$ in Fig. 3(e) for the D3 scenario; on the other hand,
for the D2 scenario, both Figs. 2(e)-(f) identify such onset. The differences
among these dynamics and their path to thermalization result from the
different initial conditions, in which in the configuration space either the
dynamics is antidiagonal (D1), diagonal (D2), or both (D3).

On the other hand, the onset of thermalization is always accessible from its
hidden-variable signature in every scenario, irrespective to initial
conditions. Not only (39) is obviously satisfied at any time, but the kinetic
energy equipartition in (40), stating that Bohmian kinetic and quantum
potential energies should become equal, with each equalizing half of the
orthodox kinetic energy, is always true after thermalization is set. This is
the same of saying that the squared values of osmotic and current velocities
become equal, with each equalizing half of the squared momentum, which also
implies that the correlations obey (43) in addition to (30); the validity of
(42) in addition to (41) can be taken as a restatement of the virial theorem
when reaching some steady state.

These hidden variables, linked separately to amplitude (osmotic) and phase
(current) components of the many-body wave function, are not linked to
orthodox operators, but are accessible in the laboratory through a post-
proccessing of local-in-position weak values protocols for momentum and
kinetic energy; the real and imaginary parts of the momentum weak value are
respectively tied to the current and osmotic hidden variables. Thermalization
is, so to speak, a manifestation of both real (amplitude) and imaginary
(phase) parts of the wavefunction becoming completely and basically
homogeneously spread through the whole configuration space, making it hard for
one to differentiate among one or another.

In order to properly understand the merits of our work it is essential to
notice that all hidden-variable results in Secs. II and III are not just
visualizing tools, but they open a new unexplored link between theoretical
predictions and empirical data. Notice that the use of term ‘hidden variables’
has no ontological implication (only historical reasons). All quantum theories
with empirical agreement with experiments do exactly predict the same weak
values. In simpler words, the link between theoretical predictions and
empirical data do not need to choose any particular ontology. To emphasize the
accessibility of weak values in the laboratory, we have also addressed, by
moving to the center of mass frame of reference, how the weak values of such
center of mass can be employed to approach larger systems with $N$ identical
particles, where the novel kinetic energy equipartion here presented should
also be verified.

###### Acknowledgements.

This research was funded by Spain’s Ministerio de Ciencia, Innovación y
Universidades under Grant PID2021-127840NB-I00 (MICINN/AEI/FEDER, UE),
Generalitat de Catalunya and FEDER for project 001-P-001644 (QUANTUMCAT),
European Union’s Horizon 2020 research and innovation programme under Grant
881603 GrapheneCore3.

## References

  * [1] M. Ueda, Nat. Rev. Phys. 2, 669 (2020).
  * [2] C. Gogolin and J. Eisert, Rep. Prog. Phys. 79, 056001 (2016).
  * [3] J. Eisert, M. Friesdorf, and C. Gogolin, Nat. Phys. 11, 124 (2015).
  * [4] V. I. Yukalov, Laser Phys. Lett. 8, 485 (2011).
  * [5] P. Reimann, New J. Phys. 21, 053014 (2019).
  * [6] M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 (2008).
  * [7] R. Nandkishore and D. A. Huse, Annu. Rev. Condens. Matter Phys. 6, 15 (2015).
  * [8] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Adv. Phys. 65, 239 (2016).
  * [9] J. M. Deutsch, Rep. Prog. Phys. 81, 082001 (2018).
  * [10] T. Langen, R. Geiger, and J. Schmiedmayer, Annu. Rev. Condens. Matter Phys. 6, 201 (2015).
  * [11] M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. Sen(De), and U. Sen, Adv. Phys. 56, 243 (2007).
  * [12] L. Sanchez-Palencia, D. Clément, P. Lugan, P. Bouyer, and A. Aspect, New J. Phys. 10, 045019 (2008).
  * [13] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
  * [14] P. Reimann, New J. Phys. 17, 055025 (2015).
  * [15] T. R. Oliveira, C. Charalambous, D. Jonathan, M. Lewenstein, and A. Riera, New J. Phys. 20, 033032 (2018).
  * [16] T. Kinoshita, T. Wenger, and D. Weiss, Nature 440, 900 (2006).
  * [17] S. Trotzky, Y-A. Chen, A. Flesch, I. P. McCulloch, U. Schollwöck, J. Eisert, and I. Bloch, Nat. Phys. 8, 325 (2012).
  * [18] A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss, and M. Greiner, Science 353, 794 (2016).
  * [19] 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).
  * [20] D. Dries, S. E. Pollack, J. M. Hitchcock, and R. G. Hulet, Phys. Rev. A 82, 033603 (2010).
  * [21] S. Will, D. Iyer, and M. Rigol, Nat. Commun. 6, 6009 (2015).
  * [22] J. Kajala, F. Massel, and P. Torma, Phys. Rev. Lett. 106, 206401 (2011).
  * [23] 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).
  * [24] B. Nagler, K. Jagering, A. Sheikhan, S. Barbosa, J. Koch, S. Eggert, I. Schneider, and A. Widera, Phys. Rev. A 101, 053633 (2020).
  * [25] 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).
  * [26] M. Srednicki, Phys. Rev. E 50, 888 (1994).
  * [27] J. M. Deutsch, Phys. Rev. A 43, 2046 (1991).
  * [28] C. Nation and D. Porras, Phys. Rev. E 102, 042115 (2020).
  * [29] M. Rigol, Phys. Rev. A 80, 053607 (2009).
  * [30] L. F. Santos and M. Rigol, Phys. Rev. E 81, 036206 (2010).
  * [31] 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).
  * [32] C. Gogolin, M. P. Muller, and J. Eisert, Phys. Rev. Lett. 106, 040401 (2011).
  * [33] M. Brenes, T. LeBlond, J. Goold, and M. Rigol, Phys. Rev. Lett. 125, 070605 (2020).
  * [34] M. Brenes, J. Goold, and M. Rigol, Phys. Rev. B 102, 075127 (2020).
  * [35] L. de Broglie, J. Phys. Radium 8, 225 (1927).
  * [36] D. Bohm, Phys. Rev. 85, 166 (1952).
  * [37] E. Nelson, Phys. Rev. B 150, 1079 (1966).
  * [38] Y. Aharonov, D. Z. Albert, and L. Vaidman, Phys. Rev. Lett. 60, 1351 (1988).
  * [39] D. Pandey, R. Sampaio, T. Ala-Nissila, G. Albareda, and X. Oriols, Phys. Rev. A 103, 052219 (2021).
  * [40] H. M. Wiseman, New J. Physics 9, 165 (2007).
  * [41] D. Dürr, S. Goldstein, and N. Zanghì, J. Stat. Phys. 134, 1023 (2009).
  * [42] D. Marian, N. Zanghi, and X. Oriols, Phys. Rev. Lett. 116, 110404 (2016).
  * [43] F. L. Traversa, G. Albareda, M. Di Ventra, and X. Oriols, Phys. Rev. A 87, 052124 (2013).
  * [44] S. Kocsis, B. Braverman, S. Ravets, M. J. Stevens, R. P. Mirin, L. K. Shalm, and A. M. Steinberg, Science 332, 1170 (2011).
  * [45] A. Hariri, D. Curic, L. Giner, and J. S. Lundeen, Phys. Rev. A 100, 032119 (2019).
  * [46] R. Ramos, D. Spierings, I. Racicot, and A. M. Steinberg, Nature 583, 529 (2020).
  * [47] C. F. Destefani and X. Oriols, Phys. Rev. A 107, 012213 (2023).
  * [48] A. Valdes-Hernandez, L. de la Peña, and A. M. Cetto, Phys. Lett. A 383, 838 (2019).
  * [49] A.G. Kofman, S. Ashhab, and F. Nori, Physics Reports 520, 43 (2012).
  * [50] J. Dressel, M. Malik, F. M. Miatto, A. N. Jordan, and R. W. Boyd, Rev. Mod. Phys. 86, 307 (2014).
  * [51] K. Renziehausen and I. Barth, Found. Phys. 50, 772 (2020).
  * [52] O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, Phys. Rev. X 6, 041065 (2016).
  * [53] V. Alba, B. Bertini, M. Fagotti, L. Piroli, and P. Ruggiero, J. Stat. Mech., 114004 (2021).
  * [54] I. Bouchole and J. Dubail, J. Stat. Mech., 014003 (2022).
  * [55] X. Oriols and J. Mompart, Applied Bohmian Mechanics: From Nanoscale Systems to Cosmology (Pan Stanford, Singapore, 2012).
  * [56] 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).
  * [57] J. S. Lundeen, B. Sutherland, A. Patel, C. Stewart, and C. Bamber, Nature 474, 188 (2011).
  * [58] J. Dressel and A. N. Jordan, Phys. Rev. A 85, 012107 (2012).
  * [59] E. Cohen and E. Pollak, Phys. Rev. A 98, 042112 (2018).
  * [60] G. Dennis, M. A. de Gosson, and B. J. Hiley, Phys. Lett. A 379, 1224 (2015).
  * [61] E. Heifetz and E. Cohen, Found. Phys. 45, 1514 (2015).
  * [62] G. Grossing, Phys. A 388, 811 (2009).
  * [63] A. S. Sanz, Front. Phys. 14(1), 11301 (2019).
  * [64] B. J. Hilley, J. Phys.: Conf. Ser. 361, 012014 (2012).
  * [65] R. Flack and B. J. Hiley, Entropy 20, 367 (2018).
  * [66] Y.-W. Hsueh, C.-H. Hsueh, and W.-C. Wu, Entropy 22, 855 (2020).
  * [67] L. Pezze, B. Hambrecht, and L. Sanchez-Palencia, EPL 88, 30009 (2009).
  * [68] L. Pezze, and L. Sanchez-Palencia, Phys. Rev. Lett. 106, 040601 (2011).
  * [69] C.-H. Hsueh, R. Ong, J.-F. Tseng, M. Tsubota, and W.-C. Wu, Phys. Rev. A 98, 063613 (2018).
  * [70] R. Ong, C.-H. Hsueh, and W.-C. Wu, Phys. Rev. A 100, 053619 (2019).
  * [71] D. Clément, A. F. Varón, J. A. Retter, L. Sanchez-Palencia, A. Aspect, and P. Bouyer, New J. Phys. 8, 165 (2006).
  * [72] 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).
  * [73] All figures in the paper consider $\omega=1.0$, $\alpha=0.01$, $\sigma_{1}=\sigma_{2}=1/\sqrt{\omega}=1.0$, $\gamma_{D}=25$, $\sigma_{D}\approx\Delta x$, with $\Delta x$ being the position resolution of the grid, defined in a way that its momentum resolution is $\Delta p\approx\Delta x$. The initial wave packet is antisymmetrized, such that the label $j$ in quantities like $\langle x_{1,2}\rangle$ is redundant. The final simulation time is $t_{f}=100$ for the D1 scenario, $t_{f}=250$ for the D2 scenario, and $t_{f}=150$ for the D3 scenario, with time step $\Delta t=0.001$. Simulation boxes are $[-60,60]$ with $2048$ points for each degree of freedom in the configuration space $x_{1}$$x_{2}$, so that for $N=2$ its size is $\approx 4\;10^{6}$.
  * [74] J. W. Abraham, K. Balzer, D. Hochstuhl, and M. Bonitz, Phys. Rev. B 86, 125112 (2012).
  * [75] C. R. McDonald, G. Orlando, J. W. Abraham, D. Hochstuhl, M. Bonitz, and T. Brabec, Phys. Rev. Lett. 111, 256801 (2013).
  * [76] J. W. Abraham, M. Bonitz, C. McDonald, G. Orlando, and T. Brabec, New J. Phys. 16, 013001 (2014).
  * [77] S. Bauch, K. Balzer, C. Henning, and M. Bonitz, Phys. Rev. B 80, 054515 (2009).
  * [78] M. Taut, Phys. Rev. A 48, 3561 (1993).
  * [79] S. Mandal, P. K. Mukherjee, and G. H. F. Diercksen, J. Phys. B: Mol. Opt. Phys. 36, 4483 (2003).
  * [80] X. Oriols and A. Benseny, New Journal of Physics, 19 063031 (2017).
  * [81] T. Fogarty, M. Á. García-March, L. F. Santos, and N. L. Harshman, Quantum 5, 486 (2021).
  * [82] G. Zisling, L. F. Santos, and Y. B. Lev, SciPost. Phys. 10, 088 (2021).
  * [83] Md. M. Ali, W.-M. Huang, and W.-M. Zhang, Sci. Rep. 10, 13500 (2020).
  * [84] J.-Q. Liao, H. Dong, and C. P. Sun, Phys. Rev. A 81, 052121 (2010).
  * [85] P. Lydzba, Y. Zhang, M. Rigol, and L. Vidmar, Phys. Rev. B 104, 214203 (2021).