# Evolution of Quantum Nonequilibrium for Coupled Harmonic Oscillators

Francisco Bento Lustosa1, Nelson Pinto-Neto1, Antony Valentini2 1Centro
Brasileiro de Pesquisas Físicas,
Rua Dr. Xavier Sigaud 150, Urca - CEP: 22290-180 - Rio de Janeiro-RJ, Brazil
2Augustus College, 14 Augustus Road, London SW19 6LN, UK.
Department of Physics and Astronomy,
Clemson University, Kinard Laboratory,
Clemson, SC 29634-0978, USA.

###### Abstract

In the context of de Broglie-Bohm pilot-wave theory, violations of the Born
rule are allowed and can be considered as describing nonequilibrium
distributions. We study the effects of interactions on quantum relaxation
towards equilibrium for a system of one-dimensional coupled harmonic
oscillators. We show by numerical simulations that interactions can delay or
even prevent complete relaxation for some initial states. We also discuss how
this effect might be relevant for cosmological scenarios and how
nonequilibrium could be detected in some models.

##  1 Introduction

The de Broglie-Bohm pilot-wave theory provides a description of quantum
systems in terms of trajectories (or configurations $q(t)$) that evolve
continuously through time guided by a wave function given by the usual
solution of the Schrödinger equation [1, 2, 3, 4, 5]. In that context, quantum
probability distributions can deviate from the usual Born rule and, through a
process of ‘quantum relaxation’, those distributions can relax to the state of
equilibrium $\rho(q,t)=|\psi(q,t)|^{2}$ [6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16]. This process is analogous to thermal relaxation and can be described in
terms of a ‘subquantum H-theorem’. Considering that all experimentally
observed quantum systems respect the Born rule and, therefore, are in quantum
equilibrium it is natural to assume that most systems in our universe have
already gone through the process of quantum relaxation. However, it is
possible that quantum nonequilibrium (or violations of the Born rule) have
existed in the very early universe [17, 7, 18, 8, 9, 19, 20].

According to standard inflationary cosmology , the temperature anisotropies
observed on the cosmic microwave background (CMB) are generated by quantum
vacuum fluctuations of the scalar field [21]. If some modes of this field were
out of equilibrium before the inflationary period it is possible that they
would transfer a nonequilibrium signature to the power spectrum of the CMB
[19, 22, 23, 24, 25, 26, 27, 28]. In another cosmological scenario,
nonequilibrium could be transferred from some modes of the inflaton to other
field modes during reheating, opening up the possibility that violations of
the Born rule could have an effect on other stages of the evolution of the
early universe (such as particle production or structure formation) or survive
in relic systems that could be detected today [22, 29, 30]. It has also been
proposed that black-hole formation and evaporation could lead to the
generation of quantum nonequilibrium in the resulting Hawking radiation [22,
31, 32, 33]. All of those scenarios, however, depend on a deeper understanding
of both the quantum relaxation process and the role of quantum mechanics in
extreme gravitational situations.

The pilot-wave theory has also been applied successfully to cosmology
providing a scenario where the initial singularity can be avoided due to
quantum effects in a bouncing model (see [34] for a recent review). The theory
has also been used to describe the transition of quantum to classical
perturbations during the early stages of the universe both in bouncing and
inflationary models [35, 36, 37, 38]. In those applications it was always
assumed that the wave function of the universe and all its subsystems were
initially in quantum equilibrium but it is also possible that violations of
the Born rule could produce detectable effects in these scenarios. In bouncing
models, one would have to study the evolution of quantum nonequilibrium of the
vacuum perturbations before and during the bounce to observe if they could
lead to detectable signatures on the CMB spectrum. It is also possible that,
even if quantum relaxation has taken place in the contracting phase,
deviations from equilibrium could be generated during the bounce by novel
quantum gravitational processes [33].

In de Broglie-Bohm theory a system is described by an evolving wave function
$\psi(q,t)$ that satisfies the usual Schrödinger equation
$\imath\partial\psi/\partial t=\mathit{H}\psi$ and also ‘guides’ an actual
configuration $q(t)$ of the system by de Broglie’s dynamical equation for the
trajectories:

$\dot{q}(t)=\frac{j(q,t)}{|\psi(q,t)|^{2}},$ (1)

where $j(q,t)$ is the usual Schrödinger current. From the Schrödinger equation
we can derive a continuity equation

$\frac{\partial|\psi(q,t)|^{2}}{\partial t}+\partial_{q}\cdot j(q,t)=0$ (2)

for the density $|\psi(q,t)|^{2}$, which has no intrinsic relation with
probabilities. A wave function guides the configuration of an _individual_
system. However, in actual quantum experiments we consider an ensemble of
systems guided by the same wave function but with different initial
configurations (or positions, in the case of a particle system). This ensemble
has a probability distribution that evolves according to

$\frac{\partial\rho(q,t)}{\partial
t}+\partial_{q}\cdot(\rho(q,t)\dot{q}(t))=0.$ (3)

Equations (1.2) and (1.3) take the same form, which means that if initially
$\rho(q,t)=|\psi(q,t)|^{2}$ this state of ‘quantum equilibrium’ will be
preserved for all times $t$. But in the context of de Broglie-Bohm theory we
can consider the possibility of nonequilibrium distributions that violate the
Born rule. If we explore that possibility, there must be a mechanism through
which those nonequilibrium states evolve to the usual distributions that we
observe the laboratory.

In classical statistical mechanics the process of relaxation to equilibrium
for an initial nonequilibrium probability distribution may be described by the
evolution of a coarse-grained H-function [39]. Using the fact that a
probability density $p$ and the volume element of phase space $d\Omega$ are
conserved along trajectories it can be shown that the classical H-function,
$H_{class}=\int d\Omega p\ln{p}$, which is essentially the negative of the
entropy, is constant in time. However, if we define a coarse-grained density
$\bar{p}$ averaged over small cells of the phase space and assume that
initially $\bar{p(0)}=p(0)$ it can be shown that the coarse-grained
$\bar{H}_{class}$ obeys

$\overline{H}_{class}(t)\leq\bar{H}_{class}(0),$ (4)

the classical H-theorem [40]. In the context of pilot-wave theory, the
continuity equations (1.2) and (1.3) imply that the ratio
$f=\rho(q,t)/|\psi(q,t)|^{2}$ is conserved along trajectories. Besides that,
analogous to the phase space volume element $d\Omega$, the element
$|\psi(q,t)|^{2}dq$ is also preserved along trajectories in configuration
space. It is then possible to define a subquantum H-function

$H=\int|\psi(q,t)|^{2}dqf\ln{f}=\int
dq\rho(q,t)\ln{\left(\frac{\rho(q,t)}{|\psi(q,t)|^{2}}\right)}.$ (5)

As the $H_{class}$, this quantity is constant in time. But dividing the
configuration space into cells and averaging $\rho$ and $|\psi|^{2}$ in each
cell, one can define the coarse-grained H-function

$\overline{H}(t)=\int
dq\overline{\rho(q,t)}\ln\left(\frac{\overline{\rho(q,t)}}{\overline{|\psi(q,t)|^{2}}}\right),$
(6)

that can be shown to obey a ‘subquantum H-theorem’

$\overline{H}(t)\leq\overline{H}(0)$ (7)

if the initial state has no fine-grained micro-structure
($\overline{\rho(q,0)}=\rho(q,0)$ and
$\overline{|\psi(q,0)|^{2}}=|\psi(q,0)|^{2}$) [6, 17, 7, 10, 9]. Quantum
equilibrium is reached when $\overline{H}(t)=0$. Whether or not equilibrium is
actually reached depends on the characteristics of the initial quantum state,
as does the relaxation timescale.

The expected quantum relaxation has been broadly confirmed by numerical
simulations for a number of systems, though there are some caveats. Two
classes of two-dimensional systems (the harmonic oscillator and a particle in
a box) have been extensively studied numerically for different superpositions
of energy modes of their quantum states. Rapid relaxation, with an
approximately-exponential decay of the coarse-grained H-function (1.6), has
been observed for most cases studied so far [10, 11, 41, 12, 15], with some
exceptions for certain states with a small number of superposed modes (for
which sometimes a small nonequilibrium ‘residue’ is found) [14]. Quantum
relaxation for a free scalar field on expanding space has also been studied
numerically, showing efficient relaxation at short wavelengths and a certain
suppression of relaxation at long (super-Hubble) wavelengths [23, 24]. More
recently, a system of two one-dimensional harmonic oscillators with a linear
time-dependent coupling was also studied and the results indicated an
approximate exponential decay to equilibrium but further simulations for
interacting systems are still needed [16].

Considering that the early universe is a prime candidate for the search for
nonequilibrium and its detectable effects, particular attention has been paid
in recent years to scalar fields initially out of equilibrium evolving on
expanding backgrounds, with a view to predicting possible anomalies in the
cosmic microwave background (CMB) caused by early violations of the Born rule
[24, 25, 27, 28]. A preliminary prediction for a large-scale power deficit was
obtained from a simplified model of a pre-inflationary era [25, 27], and a
detailed comparison with cosmological data has been made (with some intriguing
hints but inconclusive results) [28]. Studies of more complete cosmological
scenarios with early nonequilibrium, as well as more detailed predictions of
the effects on the CMB, are still needed.

In all these numerical simulations a number of factors can affect the
relaxation timescales obtained in each study. These include the number of
energy modes that are superposed, the relative phases in the initial wave
function, and the value of the coarse-graining length [12]. Underwood [42] has
also considered cases of ’extreme’ nonequilibrium, in which the density $\rho$
is initially concentrated in the tails of the wave function where
$|\psi(q,t)|^{2}$ is very small, for which the likelihood of relaxation seems
to depend on the total ’vorticity’ of $\psi$.

In a recent paper, one of the authors studied the effect of interactions for a
time-dependent coupling of two harmonic oscillators [16]. However, this study
had several numerical limitations due to the complexity of the wave function
that was given in terms of Bessel functions. Furthermore, as it was assumed
that the coupling increases linearly in time, it was not possible to separate
the effects of the coupling growth and the time-dependence. For these reasons,
in this work we take a step back, and focus exclusively on the possible
effects that the interaction has on the quantum relaxation process. We will
study a system composed of two one-dimensional harmonic oscillators with a
linear time-independent coupling. This example is of cosmological relevance,
since in the process of reheating in inflationary models, the energy transfer
of the inflaton field to other fields is usually modelled by this type of
coupling [43]. Our analysis will be based on the same numerical method used in
[16]. We will focus on the effects of the interaction on the relaxation
timescales and on the possible existence of nonequilibrium residues.

In Section II we describe our physical model and the solution of the
Schrödinger equation. In Section III we describe our numerical method and
setup, comparing its features with a different method used in previous work.
In Section IV we analyse the relaxation process for the system considered, by
calculating the evolution of nonequilibrium distributions and tracking the
behaviour of trajectories. We also obtain relaxation timescales and
nonequilibrium residues from the time evolution of the coarse-grained
H-function, and we study their dependence on the coupling and on the number of
modes superposed in the wave function.

##  2 Coupled Harmonic Oscillators

Harmonic oscillators are often used to model a wide range of systems. In
quantum field theory, a free scalar field mode on Minkowski space can be shown
to evolve in an analogous way to a standard harmonic oscillator, while on an
expanding background the corresponding oscillator has a time-dependent mass
and frequency [22]. In the inflationary models mentioned above, if the
inflaton was in a state of quantum nonequilibrium before inflation, it is
possible that some modes of the scalar field are still out of equilibrium
during and after inflation [19, 24, 25, 27]. Hence, in the process of
reheating, the nonequilibrium of the inflaton may be transferred to the other
fields when they interact [43]. It is also possible that nonequilibrium is
generated at the Planck scale [44, 33] during a cosmological bounce resulting
from the collapse of a previous contracting phase [45, 35, 36, 37, 38, 34]. In
any scenario where nonequilibrium survives the first stages of cosmological
expansion, it can result in observable effects on the CMB or in the initial
stages of particle creation.

In the case of inflation, both within the bouncing scenario or the usual
standard cosmological model, the decay of the scalar field happens through
interactions with other fields [21]. In models of multi-fluid inflation, the
effect of interaction between two scalar fields is crucial in the description
of particle creation [46]. In previous work, one of the authors considered the
effects of a time-dependent interaction on the relaxation process [16]. It was
in fact already known that a time-dependent expanding background can directly
affect the rate of quantum relaxation [24]. In this work we focus solely on
the effects of interaction in a simple model of coupled one-dimensional
oscillators.

This model is described by the following Hamiltonian:

$H(x_{a},x_{b},t)=\frac{p_{a}^{2}}{2m}+\frac{m\omega^{2}x_{a}^{2}}{2}+\frac{p_{b}^{2}}{2m}+\frac{m\omega^{2}x_{b}^{2}}{2}+\frac{mk}{2}x_{a}x_{b}.$
(8)

To decouple equation (2.1) we introduce the standard canonical transformation
of coordinates for the positions and momenta;

$\displaystyle\begin{pmatrix}x_{a}\\\
x_{b}\end{pmatrix}=\sqrt{\frac{1}{2m}}\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}\begin{pmatrix}x_{1}\\\ x_{2}\end{pmatrix},$ (9)
$\displaystyle\begin{pmatrix}p_{a}\\\
p_{b}\end{pmatrix}=\sqrt{\frac{m}{2}}\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}\begin{pmatrix}p_{1}\\\ p_{2}\end{pmatrix}.$ (10)

It can be shown that the transformed coordinates obey the same commutation
relations as the original ones. Writing the transformed Hamiltonian in terms
of the new coordinates $x_{r}$ ($r=1,2$) we have

$H(x_{1},x_{2},t)=\sum_{r=1}^{2}H_{r}=\sum_{r=1}^{2}\left(\frac{p_{r}^{2}}{2}+\frac{\Omega_{r}^{2}}{2}x_{r}^{2}\right),$
(11)

where the new frequencies are $\Omega_{1}=\sqrt{(\omega^{2}+k/2)}$ and
$\Omega_{2}=\sqrt{(\omega^{2}-k/2)}$. Because of the negative sign in
$\Omega_{2}$ the values of $k$ are restricted to the range $[0,2\omega^{2}]$.
After quantisation, we can write the Schrödinger equation as

$\imath\frac{\partial\Psi}{\partial
t}=\sum_{r=1}^{2}\left(-\frac{1}{2}\frac{\partial^{2}}{\partial
x_{r}^{2}}+\frac{1}{2}\Omega_{r}^{2}x_{r}^{2}\right)\Psi,$ (12)

where we take $\hbar=1$. The solution can be written as a superposition of
eingenstate solutions of the one-dimensional Schrödinger equations
$\imath\partial\psi_{r}/\partial t=\hat{H}_{r}\psi_{r}$. The full superposed
solution in terms of the original coordinates is then

$\Psi(x_{a},x_{b},t)=\sum_{n_{1},n_{2}}c_{n_{1},n_{2}}\psi_{n_{1}}(x_{1}(x_{a},x_{b}),t)\psi_{n_{2}}(x_{2}(x_{a},x_{b}),t),$
(13)

where the constants $c_{n_{1},n_{2}}$ are taken to be $(1/\sqrt{M})e^{\imath
2\pi\theta_{n_{1},n_{2}}}$ where $M$ is the number of superposed energy modes
and $\theta_{n_{1},n_{2}}$ are randomly chosen phases. For each value of $M$
we randomly choose a combination of quantum numbers $n_{1}$ and $n_{2}$ (using
the same method that was described in [16]). The functions
$\psi_{n_{r}}(x_{r},t)$ are the usual solutions for the 1D oscillator written
in terms of the Hermite polynomials.

To obtain actual trajectories $x_{a}(t)$ and $x_{b}(t)$ we construct guidance
equations in the transformed coordinates $x_{1}(t)$ and $x_{2}(t)$,
transforming the initial coordinates
$(x_{a}(0),x_{b}(0))\rightarrow(x_{1}(0),x_{2}(0))$ and solving the de Broglie
equations

$\dot{x_{r}}=\operatorname{Im}\left(\frac{\partial_{r}\Psi}{\Psi}\right),$
(14)

to obtain different mappings of
$(x_{1}(0),x_{1}(0))\rightarrow(x_{1}(t),x_{2}(t))$ for an ensemble of
trajectories. Before calculating the coarse-grained $\bar{\rho}$,
$\overline{|\psi|^{2}}$ and $\overline{H}(t)$ we transform back to the
original coordinates $(x_{a}(t),x_{b}(t))$, thus obtaining the evolution of a
chosen initial nonequilibrium distribution
$\rho(x_{a},x_{b},t)\neq|\Psi(x_{a},x_{b},t)|^{2}$.

##  3 Numerical Setup

###  3.1 Numerical Methods; backtracking and the forward trajectories method

The quantum relaxation process depends on a number of parameters affecting the
evolution of the quantum systems and of their probability distributions. In
order to observe how the approach to equilibrium changes with those
parameters, different methods of analysis have been developed.

Firstly, for a one-dimensional system of a particle in a box [8] an
approximate and partial approach to equilibrium was observed on a coarse-
grained level. Improved numerical simulations for a two-dimensional system
showed a good fit of the data to an H-function of the form
$\overline{H_{0}}e^{-t/\tau}$ [10, 12]. Furthermore, the effects of the number
of modes $M$ of the wave function and the coarse-graining length $\epsilon$ on
the relaxation timescales were demonstrated. The results showed a directly
correlation between the increase of $M$ and the decrease of $\tau$. A more
complete study of the evolution of nonequilibrium distributions for a 2D
harmonic oscillator for longer relaxation times was done in [14]. In that
work, the simulation results showed that for lower values of $M$ it was
possible that some residual nonequilibrium remained even after a long time.
Those simulations also provided a new best fit for the coarse-grained
H-function of the form $ae^{-t/\tau}+R$. The case of a free mode of a scalar
field on an expanding spatial background and on Minkowski space was also
studied [24]. Rapid relaxation was observed for the the scalar field (with
$M=25$) evolving in flat space-time but, for the time-dependent expanding
background, equilibrium was not reached for long-wavelength (super-Hubble)
modes of the scalar field.

All of those studies cited above used a similar _backtracking_ method (for a
detailed description see [10]). This method was successful in describing the
relaxation process and obtaining correlations between the relaxation
timescales $\tau$ and the values of $M$ and $\epsilon$. However, the
backtracking method has the disadvantage of taking a large amount of computer
time to calculate values of $\overline{H}(t)$ for different values of $t$. For
large number of modes, and for ensembles of trajectories, finding the
evolution of the coarse-grained H-function becomes computationally
intractable.

A way to avoid this problem is to study the character of the trajectories [14,
47, 11, 48, 15, 49] without calculating the evolution of actual complete
distributions. Due to the chaotic nature of many of the quantum trajectories,
it is possible to see indications of relaxation, or lack thereof, by analysing
selected sets of trajectories. If the trajectories explore large regions of
the configuration space this indicates that quantum relaxation can occur more
efficiently. On the other hand, if they are confined to small regions that is
an indication that quantum equilibrium cannot be reached. Although useful to
determine general properties of the trajectories, this method does not allow
for a more precise quantification of the relaxation process.

Another way to study the trajectories and quantum relaxation is related to the
_vorticity_ [50, 51] of the velocity fields and it was used to study the time
evolution of some extreme nonequilibrium distributions in [29]. The vorticity
is related to the number of modes $M$ as well as to the chosen random phases
$\theta_{n_{1},n_{2}}$ but it is not clear how we can use this method in a
more general manner to quantify relaxation timescales.

Finally, in recent work another method was used to study a system of two one-
dimensional oscillators with a time-dependent coupling (we shall refer to this
method as the Forward Trajectories Method, or FTM) [16]. Using this method,
one can calculate the trajectories from $t_{0}$ to a final time, recording the
distribution of trajectories at each time $t$ to obtain the H-function. One of
the reasons that the backtracking method was used in previous work was to
avoid the uneven sampling of distributions when some regions of configuration
space contain most trajectories and others have close to none. However, with a
large enough number of trajectories this problem does not appear relevant to
the accurate calculation of $\overline{H}(t)$. For each time $t$ that we
record the distribution, the coarse-grained value of $\rho$ is calculated by
summing the number of trajectories that are in each cell. It is the total
coarse-grained distribution that is relevant to the calculation of relaxation
timescales. The backtracking method is highly sensitive to changes in the
sampling within each cell whereas the FTM method gives approximately the same
values of $\overline{H}(t)$ for different coarse-graining lengths.

Both methods have advantages and limitations, however, the FTM is more
efficient and allowed for a wider analysis of the effect of the parameters on
the relaxation process for the system considered in [16]. In Figure 1 we show
a comparison of the two methods for a simple case with $M=9$ for the wave
function (2.6). The error bars are obtained from the standard deviation
calculated from three different coarse-graining lengths (which leads to
different number of points per cell) and the points shown on the graph are the
mean values of the H-functions. Both cases were calculated on a grid with a
total of 256 cells. The two H-functions, obtained in the interval $t_{0}=0$ to
$t_{f}=2\pi$, are in relatively good agreement with one another and the bigger
spread in some of the points for the backtracked case indicates that the FTM
is slightly more precise. The bigger error bars for the backtracking case are
a result of the sensitivity of the method to changes in the number of points
per cell (the error bars for the FTM are not visible because the values of
$\overline{H}(t)$ for different numbers of points per cell are almost the
same). Since those changes are practically irrelevant in the method we will
use from now on, we will calculate the error bars in a different manner in the
next sections.

Figure 1: $\overline{H}(t)$ functions calculated using the FTM method (red
stars) and the backtracking method (blue stars) The solid lines are the best
fit functions of the form $(\bar{H_{0}}-R)e^{-t/\tau}+R$ for each data set.

###  3.2 Our Setup

The simulations in the next section are all done using the FTM method.
Trajectories are calculated using a Runge-Kutta method of order 8 with an
initial absolute error tolerance of $10^{-9}$ and then again with a tolerance
of $10^{-10}$ [52] to solve the set of equations (2.7). If the results of the
two calculations have a difference of more then $\delta=5\times 10^{-3}$ in
either coordinate the calculation is repeated with an error tolerance one
order of magnitude smaller. This process is repeated until the absolute error
tolerance reaches the limit value of $10^{-16}$, if the difference between
calculations still exceeds $\delta$ the trajectory is excluded from the
distribution. In order to obtain simpler velocity fields we transform the
initial positions using equation (2.2) and we transform back to the original
coordinates at each time $t$ we wish to calculate the coarse-grained
H-function.

The calculation of distributions is done with $N=230400$ trajectories that
evolve in a two-dimensional box in the region $[x_{a}=-5..5,x_{b}=-5..5]$. To
obtain the coarse-grained functions $\overline{\rho(x_{a},x_{b},t)}$,
$\overline{|\psi(x_{a},x_{b},t)|^{2}}$ and $\overline{H}(t)$ we divide the
considered region in a grid of $16\times 16$ cells. The initial distribution
$\rho(x_{a},x_{b},0)$ is obtained using a random Gaussian generator for each
calculation. As the trajectories evolve some cells will have more points then
others and the coarse-grained value of the distribution in each cell is
obtained by summing the number of points that are in a given cell at each time
$t$.

The wave functions are defined up to factors $c_{n_{1},n_{2}}$. The phases
$\theta_{n_{1},n_{2}}$ and quantum numbers $[n_{1},n_{2}]$ are randomly
generated for each number of superposed modes $M$. In order to include the
random effect of the choice of phases and avoid any bias in the initial wave
functions [42, 16], we repeat the calculations for different sets of random
phases and quantum numbers. Error bars are obtained with the mean deviation
calculated from the different values of $\overline{H}(t)$ for different sets
of phases and quantum numbers. Finally, all simulations are done with
$m=\omega=1$ and for different values of the coupling constant $k$.

##  4 Interaction and the Evolution of Quantum Nonequilibrium

In this section we analyse how the variation of the value of the coupling
constant affects the evolution of an initial nonequilibrium distribution. We
will show how the interaction influences the behaviour of trajectories and
distributions, as well as affecting the relaxation timescales and
nonequilibrium residues.

###  4.1 Distributions and Confinement of Trajectories

The relaxation process can be analyzed by various methods but it is more
easily visualized by the direct observation of the evolution of the
nonequilibrium distribution. For a given initial ensemble with density $\rho$,
if the system evolves towards equilibrium we expect that the final coarse-
grained functions $\overline{\rho}$ and $\overline{|\psi|^{2}}$ will have
approximately the same form [10, 12, 14]. If they are exactly equal, from
equation (1.7), we find $\overline{H}(t)=0$ and the distribution has reached
the equilibrium state. However, previous simulations have shown that for a
sufficiently low number of superposed modes of the wave function, even after a
significantly large time, complete relaxation may not be achieved in some
cases [14]. Other factors may prevent complete relaxation, such as the form of
the initial distribution and the characteristics of the velocity field [42].
It has also been shown that a scalar field mode initially out of equilibrium
evolving on a time-dependent cosmological background does not relax completely
if the mode wavelength is larger than the Hubble radius [24], indicating that
time-dependent Hamiltonians can lead to distributions that avoid relaxation.
The time-dependent coupled system studied in [16] also indicated that
interactions may delay or even prevent complete relaxation.

|  |
---|---|---
(a) $\overline{|\psi(0)|^{2}}$ | (b) $\overline{|\psi(5\pi)|^{2}}$ | (c) $\overline{|\psi(10\pi)|^{2}}$
|  |
(d) $\overline{\rho(0)}$ | (e) $\overline{\rho(5\pi)}$ | (f) $\overline{\rho(10\pi)}$
Figure 2: Coarse-grained distribution and $|\psi|^{2}$ at times $t=0$ ((a) and (d)), $t=5\pi$ ((b) and (e)) and $t=10\pi$ ((c) and (f)) for the case of $24$ superposed modes of the wave function (2.6) and coupling constant $k=0.1$. The initial $\overline{\rho}$ is very close to $\overline{|\psi|^{2}}$ at $t=5\pi$ and is virtually indistinguishable from the equilibrium distribution at $t=10\pi.$ |  |
---|---|---
(a) $\overline{|\psi(0)|^{2}}$ | (b) $\overline{|\psi(5\pi)|^{2}}$ | (c) $\overline{|\psi(10\pi)|^{2}}$
|  |
(d) $\overline{\rho(0)}$ | (e) $\overline{\rho(5\pi)}$ | (f) $\overline{\rho(10\pi)}$
Figure 3: Coarse-grained distribution and $|\psi|^{2}$ at times $t=0$ ((a) and
(d)), $t=5\pi$ ((b) and (e)) and $t=10\pi$ ((c) and (f)) for the case of $24$
superposed modes of the wave function (2.6) and coupling constant $k=1.8$. The
initial coarse-grained nonequilibrium Gaussian distribution does not reach an
approximate equilibrium state at $t=10\pi$, indicating that the interaction
might delay or prevent the relaxation process.

In this work we focus on how a constant interaction, determined by the value
of the coupling constant $k$, affects the evolution of a nonequilibrium
distribution. As a preliminary illustration, Figure 2 shows the evolution of
an initial coarse-grained nonequilibrium Gaussian distribution and the coarse-
grained $|\psi|^{2}$ given by (2.6) from $t=0$ to $t=10\pi$ for the case with
$M=24$ and $k=0.1$. The distribution evolves efficiently to an approximate
equilibrium state and is virtually indistinguishable from the
$\overline{|\psi|^{2}}$ at the final time. In contrast, in Figure 3 we show
the evolution of $\overline{\rho}$ and $\overline{|\psi|^{2}}$ with the same
number of modes but with $k=1.8$. It is clear that in this case complete
relaxation is not achieved at time $t=10\pi$ and $\overline{\rho}$ and
$\overline{|\psi|^{2}}$ differ significantly.

The evolution of distributions shown in Figures 2 and 3 is a preliminary
indication of the effect of interactions on quantum relaxation. To further
analyse the relaxation process it is useful to observe how trajectories evolve
individually. Quantum trajectories are often chaotic [11, 48, 15] and tend to
explore most of the support region of the wave function in cases where
efficient relaxation takes place. The number of superposed modes of the wave
function directly influences the character of the trajectories since they
determine the number of nodes of the wave function [12, 42]. When trajectories
approach those nodes the velocity field becomes rapidly-varying and the
trajectories tend to be more chaotic for higher numbers of modes. Chaotic
trajectories usually explore a wide region of configuration space and
trajectories that are initially very close can have very distant final
positions. When most trajectories have those characteristics the approach to
equilibrium occurs efficiently. On the other hand, the confinement of
trajectories to certain regions is directly related to systems that avoid
complete relaxation [14].

To analyse the characteristics of the trajectories for our system of coupled
oscillators we make use of two tests that provide a visual representation of
how given sets of trajectories explore the configuration space [14]. The first
test consists of plotting the initial and final positions of five sets of
neighbouring trajectories to observe if they spread over the whole support
region of the wave function or if they stay confined and close to one another.
In the second test we observe the full evolution of 5 individual trajectories
to see how they explore the configuration space. In Figure 4 we plot the
results of those tests for the case where the wave function has a
superposition of 24 modes and for four values of the coupling constant $k$.
The right column shows 5 trajectories with randomly selected initial
positions111The initial positions are (2.39,-1.94), (1.27, 3.76),
(2.21,-0.01), (-0.89, 2.04) and (-2.30, -0.76). evolving from $t=0$ to
$t=10\pi$. On the left column we have 5 sets of 25 trajectories centered on
the initial positions used on the right column plots. For $k=0.1$ and $k=0.9$
most trajectories are not confined and spread through the bulk of the support
region of the wave function. For the case $k=1.1$ the blue trajectories are
strongly confined to a region close to their initial positions and most of the
green trajectories are also confined. The case with $k=1.8$ shows confinement
to regions close to the initial positions for the purple and blue sets of
trajectories and also shows confinement to a certain region for more than half
of the black, green and red sets.

These results indicate that for $k<1$ there is little or no confinement, while
for $k>1$ the confinement becomes increasingly stronger. Simulations with
different numbers of modes also corroborate the indication that there might be
two regimes ($k<1$ and $k>1$) for the system we are studying. As we will
discuss, this might be explained by the relation between the values of the
fundamental frequencies $\omega$ and $k$. Further analysis of this hypothesis
will be done in subsection 4(c) by analysing dependence of the relaxation
timescales and nonequilibrium residues on the coupling constant.

|
---|---
(a) | (b)
|
(c) | (d)
|
(e) | (f)
|
(g) | (h)
Figure 4: We display initial and final positions ((a), (c), (e) and (g)) of 5
sets of 25 trajectories with close initial positions. In ((b), (d), (f) and
(h)) we show 5 full trajectories. All calculations were done with
superposition of $20$ modes on the wave function and for increasing (top to
bottom) coupling constant $k=0.1,0.9,1.1$ and $1.8$.

###  4.2 Long Time Nonequilibrim Residues

In previous works, numerical simulations that calculate the coarse-grained
H-function were usually done for relatively short periods of time (except in
[14]). This can be mostly attributed to the computational difficulties related
to calculating a large enough number of trajectories to accurately model a
distribution, especially when considering wave functions with larger numbers
of modes. However, considering the approximate exponential decay of most
H-functions calculated so far, those short time simulations are very useful to
estimate relaxation timescales $\tau$ for a fitting function of the form
$\bar{H_{0}}e^{-t/\tau}$.

|
---|---
(a) $M=4$ and $k=0.5$ | (b) $M=4$ and $k=1.8$
|
(c) $M=12$ and $k=0.5$ | (d) $M=12$ and $k=1.8$
|
(e) $M=20$ and $k=0.5$ | (f) $M=20$ and $k=1.8$
Figure 5: $\overline{H}(t)$ functions averaged over three presets of random
phases $\theta_{n_{1},n_{2}}$ calculated in the time interval $[0,100\pi]$.
For $k=0.5$ there is complete relaxation at an efficient rate for the cases
$M=12$ and $M=20$. When $k=1.8$ relaxation timescales are considerably larger
and there are some nonequilibrium residues in all cases. The red line
represents the best fit function $(\bar{H_{0}}-R)e^{-t/\tau}+R$.

Longer time simulations are needed to verify if, after the initial decay of
the H-function, complete relaxation is actually achieved (within the precision
limits of the numerical simulations). If that is the case, the value of
$\overline{H}(t)$ becomes approximately zero and the system is considered to
be in complete equilibrium. Another possibility, which has been observed for
some cases with a small number of modes [14], is that a residual
nonequilibrium remains even after a long time, indicating that the Born rule
might cause observable effects for such systems. Systems that present this
feature can then plausibly model realistic situations where quantum
nonequilibrium might be detected.

In Figure 5 we show examples of the evolution of H-functions for a time
interval of $100\pi$.222This value is given in the units $\hbar=c=1$. The
calculations of $\overline{H}(t)$ are done for three different sets of
randomly-chosen initial phases $\theta_{n_{1},n_{2}}$ and the results are then
averaged. The error bars are given by the standard deviation of the three
values of the H-function (at each time $t$). The red lines represent the best
fit of the form $(\bar{H_{0}}-R)e^{-t/\tau}+R$. Values for $\tau$ and $R$ are
obtained by fitting each H-function result and then averaging over the three
relaxation timescales and residues obtained from the different sets of initial
phases.

The case with $M=4$ decays very slightly for both values of $k$ but the
average value of $\tau$ is considerably higher for the case with higher
coupling constant ($\tau\approx 56.43$ for $k=0.5$ and $\tau\approx 101.65$
for $k=1.8$). The averaged residues also represent a higher percentage of
$\overline{H}(0)$ for the case with stronger interaction ($R\approx 38.13\%$
for $k=0.5$ and $R\approx 72.48\%$ for $k=1.8$).

The cases with $M=12$ and $M=20$ for the coupling constant $k=0.5$ show
approximately complete relaxation. The average residues for both cases
represent approximately $0.5\%$ of the initial value of the coarse-grained
H-function. The relaxation timescale is smaller for the case with a higher
number of modes. The cases with $k=1.8$ show a considerable increase both in
the values of $\tau$ and $R$. However, for a sufficiently long time both cases
have a nonequilibrium residue of approximately $5\%$ of the initial value of
$\overline{H}(t)$, which cannot be considered negligible. These results show
that a larger coupling can significantly delay, and even somewhat suppress,
quantum relaxation (especially for a sufficiently small number of modes). This
might allow for nonequilibrium effects to be detectable in certain situations.
In the next subsection we will further explore the relation between $\tau$,
$R$ and $k$.

|
---|---
(a)$\tau(k)$ for $M=12$. | (b)$R(k)$ for $M=12$.
|
(c)$\tau(k)$ for $M=15$. | (d)$R(k)$ for $M=15$.
|
(e) $\tau(k)$ for $M=18$. | (f) $R(k)$ for $M=18$.
Figure 6: Relaxation timescales and nonequilibrium residues as functions of
the interaction coupling $k$ (varying from $0.00$ to $1.00$) obtained from the
fittings of $\overline{H}(t)$. The values of $\tau$ and $R$ were calculated
for each H-function and then averaged over ten sets of randomly-chosen initial
phases for each value of $k$. Error bars correspond to standard deviations
over these ten values for each point.

###  4.3 Relaxation Timescales, Residues and Interactions

In the last two subsections we have used the Forward Trajectories Method to
study the effect of the coupling constant on the evolution of distributions,
trajectories and H-functions showing how interactions might prevent
nonequilibrium initial distributions from reaching a final equilibrium state.
In order to properly quantify the effects indicated by those results we will
now study how the H-function evolves for different values of the coupling
constant and the number of superposed modes of the wave function. In this
subsection we study if and how interactions might be correlated to the
relaxation timescales and possible nonequilibrium residues.

We have studied the evolution of 37 different H-functions with increasing
coupling constants ranging from $k=0.0$ to $k=1.80$ during the time interval
$t=[0,12\pi]$. These H-functions were calculated for ten different sets of
random phases and then fitted to functions of the form
$(\bar{H_{0}}-R)e^{-t/\tau}+R$ which in turn provided ten different values for
the relaxation timescales $\tau$ and residues $R$. The values obtained from
the best fits of the coarse-grained H-functions were then averaged over and
errors bars were obtained from the standard deviation of the ten different
values for each $k$. The results of those simulations indicate the existence
of two regimes concerning the variation of the coupling constant. The analysis
of the evolution in those two regimes has led us to confirm the indications
from the trajectories shown in Figure 4.

|
---|---
(a)$\tau(k)$ for $M=12$. | (b)$R(k)$ for $M=12$.
|
(c)$\tau(k)$ for $M=15$. | (d)$R(k)$ for $M=15$.
|
(e) $\tau(k)$for $M=18$. | (f) $R(k)$ for $M=18$.
Figure 7: Relaxation timescales and nonequilibrium residues as functions of
the interaction coupling $k$ (varying from $1.05$ to $1.80$) obtained from the
fittings of $\overline{H}(t)$. The values of $\tau$ and $R$ were calculated
for each H-function and then averaged over the ten sets of random phases for
each value of $k$. Error bars correspond to standard deviations over these ten
values for each point.

This can be seen in Figures 6 and 7. In Figure 6 we observe a very slight
increase in the values of both $\tau$ and $R$ and an overall oscillating
evolution. This can be explained by the fact that the interaction does not
dominate over the fundamental frequency ($k<\omega^{2}$) and its increase does
not fundamentally affect the relaxation process. In Figure 7 the situation is
quite different, the increase in the coupling constant is directly related to
an increase in the relaxation timescale. The increase in the values of the
nonequilibrium residue is also noticeable, although less accentuated.

These results indicate that, when the oscillating frequency $\omega$ dominates
over the interaction, the effects of the coupling do not lead directly to an
increase in the relaxation timescale or the nonequilibrium residue. However,
when the interaction dominates it can significantly delay and perhaps even
prevent complete quantum relaxation. This conclusion suggests that quantum
nonequilibrium might survive situations where strong interactions dominate
over other characteristics of the system. Such effects might occur in the
interactions between scalar fields in the early universe or in the process of
parametric resonance characteristic of inflationary reheating scenarios.

##  5 Conclusions

In this paper we have studied the effects of interactions on quantum
relaxation in a system of coupled one-dimensional harmonic oscillators. This
study was done using a numerical Forward Trajectories Method that allowed us
to analyze the system in a number of different ways.

Our detailed study of the evolution of nonequilibrium distributions shows that
relaxation takes place efficiently for a significant number of superposed
modes of the wave function and for a small value of the coupling constant $k$.
Analysis of the trajectories shows that there is little or no confinement when
the coupling is smaller then the squared fundamental frequency of the
oscillators, indicating that in those cases relaxation is unlikelly to be
affected by the coupling. When $k>\omega^{2}$ a higher percentage of
trajectories stay confined to certain regions and complete relaxation is less
likely, although still possible after a long enough time and with a sufficient
number of modes.

The evolution of the coarse-grained H-function for a long time period
confirmed previous results which demonstrated that for a low number of modes
($M=4$) quantum equilibrium need not be reached exactly. For larger numbers of
modes ($M=12$ and $M=20$) the system approximately obeys the Born rule after a
sufficient time period. However, simulations with larger couplings show that
the relaxation timescale grows significantly and nonequilibrium might still be
large enough to produce detectable effects after a considerable time (for
example, the value of $\overline{H}(t)$ represents more then $20\%$ of the
initial value at $t=16\pi$ for the cases $M=12$ and $M=20$).

Finally, we analyzed the evolution of $\overline{H}(t)$ over a shorter time
period $t=[0,12\pi]$ for a variety of values of $k$, confirming the findings
cited above. Fitting $\overline{H}(t)$ with a function of the form
$(\bar{H_{0}}-R)e^{-t/\tau}+R$ we obtained the values of the timescale $\tau$
and of the nonequilibrium residues $R$ ranging from $k=0.00$ to $k=1.80$. The
growth of the coupling does not significantly affect the relaxation parameters
when $k<\omega^{2}$, indicating that when the interaction is not dominant it
has no clear effect on the evolution of $\overline{H}(t)$. However, when
$k>\omega^{2}$ the timescales increases significantly with the increase in the
coupling constant. The nonequilibrium residue also seems to be affected by the
growth of $k$, but it is generally quite small for the numbers of modes
considered $M=12,15,18$. In any case, the values of $R$ are generally higher
when the interaction is stronger and this further confirms the indications of
the rest of our analysis.

Our numerical results provide a picture of how a simple interaction may delay
or partially prevent quantum relaxation. Considering this as a toy model of
more complicated systems that existed in the very early universe, where
interactions might play an important role, the results suggest a way in which
quantum nonequilibrium could survive the early stages of cosmological
evolution. In standard inflationary scenarios, if some modes of the inflaton
were still out of equilibrium after inflation, this would affect the ‘decay’
of the inflaton field, which occurs via interactions with other fields and
which creates particles in the early universe [29]. If interactions play a
role in preventing the relaxation to the Born rule, this could provide a
feasible way through which nonequilibrium might survive the ‘reheating’ phase.

Another relevant scenario would be to consider multi-field inflation [53, 46].
In such models it is still an open question how quantum nonequilibrium would
evolve during inflation and after it. Nonetheless, in such models interactions
play an important role and our results indicate that this would also be a
situation where quantum relaxation could be delayed.

Finally, in bouncing cosmological models ([34] and references therein), it is
possible that quantum vacuum fluctuations, originating in the far past of a
contracting phase dominated by a pressureless matter field, generate almost
scale invariant power spectra of cosmological perturbations. As the initial
conditions are imposed in a distant past, it is plausible that relaxation had
already occurred. However, if during a quantum bounce quantum gravity effects
are significant, nonequilibrium could be generated during this phase [33]. In
the case of bouncing models with multiple fields where couplings are important
in the evolution of perturbations, nonequilibrium may then propagate to the
future, and the results presented here may be relevant to understand its
evolution in the universe. All these hypotheses can be studied in depth using
our numerical method if it is extended to more complex time-dependent wave
functions. We hope to return to this in future work.

DATA ACCESS:
All the data used in our analysis were produced by original simulations
developed by the authors. The data used for the plots in our paper as well as
the source codes of the simulations are accessible through the link:
(https://drive.google.com/file/d/1XsYLu1VfaNH5yvmrLi29mKPhj8KHzbmA/view?usp=sharing).
The fitting of the data was made using mathematical software as described in
the paper.

AUTOR CONTRIBUTIONS:
F.B.L.’s developed the code, produced the numerical results, worked on the
numerical and theoretical analysis of the results and drafted the manuscript.
N.P.N. worked on the theoretical analysis of the results and helped draft the
manuscript. A. V. also worked on the theoretical analysis of the results and
on the draft of the manuscript.

FUNDING:
F.B.L.’s work was funded by CNPq (PCI-DC grant 300769/2022-9 associated with
the PCI/MCTI/CBPF program). N.P.N. acknowledges the support of CNPq of Brazil
under grant PQ-IB 310121/2021-3.

AKNOWLEDGEMENTS:
F.B.L. would like to acknowledge the fundamental contribution of Samuel Colin
in the development of the numerical method used in our work.

## References

  * [1] L. de Broglie, La nouvelle dynamique des quanta, in: Électrons et Photons: Rapports et Discussions du Cinquième Conseil de Physique, Gauthier-Villars, Paris, 1928, [English translation in ref. [2]].
  * [2] G. Bacciagaluppi, A. Valentini, Quantum theory at the crossroads: Reconsidering the 1927 Solvay conference, Cambridge University Press, Cambridge, 2009.  arXiv:quant-ph/0609184.
  * [3] D. Bohm, A suggested interpretation of the quantum theory in terms of "hidden" variables. I, Phys. Rev. 85 (1952) 166–179.  doi:10.1103/PhysRev.85.166.
  * [4] D. Bohm, A suggested interpretation of the quantum theory in terms of hidden variables. II, Phys. Rev. 85 (1952) 180–193.  doi:10.1103/PhysRev.85.180.
  * [5] P. Holland, The Quantum Theory of Motion: An Account of the de Broglie-Bohm Causal Interpretation of Quantum Mechanics, Cambridge University Press, 1993.
  * [6] A. Valentini, Signal-locality, uncertainty, and the subquantum H-theorem. I, Phys. Lett. A156 (1991) 5–11.  doi:10.1016/0375-9601(91)90116-P.
  * [7] A. Valentini, On the pilot-wave theory of classical, quantum and subquantum physics, Ph.D. thesis, International School for Advanced Studies, Trieste, Italy, (1992).
URL http://hdl.handle.net/20.500.11767/4334

  * [8] A. Valentini, Hidden Variables, Statistical Mechanics and the Early Universe, in: Chance in Physics: Foundationns and Perspectives, eds. J Bricmont et al., Springer, Berlim, 2001.  arXiv:quant-ph/0104067.
  * [9] A. Valentini, Foundations of statistical mechanics and the status of the Born rule in de Broglie-Bohm pilot-wave theory, in: Statistical Mechanics and Scientific Explanation: Determinism, Indeterminism and Laws of Nature, ed. V. Allori, World Scientific, 2020.  arXiv:1906.10761.
  * [10] A. Valentini and H. Westman, Dynamical origin of quantum probabilities, Proc. R. Soc. A 461 (2005) 253–272.  arXiv:quant-ph/0403034.
  * [11] C. Efthymiopoulos and G. Contopoulos, Chaos in Bohmian quantum mechanics, J. Phys. A: Math. Gen 39 (8) (2006) 1819–1852.  doi:10.1088/0305-4470/39/8/004.
  * [12] M. D. Towler, N. J. Russell and A. Valentini, Time scales for dynamical relaxation to the Born rule, Proc. R. Soc. A 468 (2012) 990–1013.  arXiv:1103.1589.
  * [13] S. Colin, Relaxation to quantum equilibrium for dirac fermions in the de Broglie-Bohm pilot-wave theory, Proc. R. Soc. A 468 (2012) 1116–1135.  arXiv:1108.5496.
  * [14] E. Abraham, S. Colin and A. Valentini, Long-time relaxation in pilot-wave theory, J. Phys. A: Math. Theor. 47 (39) (2014) 395306.  arXiv:1310.1899.
  * [15] C. Efthymiopoulos, G. Contopoulos and A. C. Tzemos, Chaos in de Broglie-Bohm quantum mechanics and the dynamics of quantum relaxation, Ann. Fond. Louis de Broglie 42 (2017) 133–159.  arXiv:1703.09810.
  * [16] F. B. Lustosa, S. Colin and S. E. Perez Bergliaffa, Quantum relaxation in a system of harmonic oscillators with time-dependent coupling, Proc. R. Soc. A 477 (2248) (2021) 20200606.  arXiv:2007.02939.
  * [17] A. Valentini, Signal-locality, uncertainty, and the subquantum H-theorem. II, Phys. Lett. A 158 (1991) 1–8.  doi:10.1016/0375-9601(91)90330-B.
  * [18] A. Valentini, Pilot-wave theory of fields, gravitation and cosmology, in: Bohmian Mechanics and Quantum Theory: An Appraisal, eds. J. T. Cushing et al., Kluwer, Dordrecht, 1996.  doi:10.1007/978-94-015-8715-0_3.
URL https://doi.org/10.1007/978-94-015-8715-0_3

  * [19] A. Valentini, Inflationary cosmology as a probe of primordial Quantum Mechanics, Phys. Rev. D 82 (2010) 063513.  arXiv:0805.0163.
  * [20] A. Valentini, Beyond the quantum, Phys. World 22N11 (2009) 32–37.  arXiv:1001.2758.
  * [21] P. Peter, J.-P. Uzan, Primordial cosmology, Oxford Graduate Texts, Oxford Univ. Press, Oxford, 2009.
URL https://cds.cern.ch/record/1208401

  * [22] A. Valentini, Astrophysical and cosmological tests of quantum theory, J. Phys. A: Math. Theor. 40 (2007) 3285–3303.  arXiv:hep-th/0610032.
  * [23] A. Valentini, De Broglie-Bohm prediction of quantum violations for cosmological super-Hubble modes (2008).  arXiv:0804.4656.
  * [24] S. Colin, A. Valentini, Mechanism for the suppression of quantum noise at large scales on expanding space, Phys. Rev. D88 (2013) 103515.  arXiv:1306.1579.
  * [25] S. Colin, A. Valentini, Primordial quantum nonequilibrium and large-scale cosmic anomalies, Phys. Rev. D 92 (4) (2015) 043520.  arXiv:1407.8262.
  * [26] A. Valentini, Statistical anisotropy and cosmological quantum relaxation (2015).  arXiv:1510.02523.
  * [27] S. Colin, A. Valentini, Robust predictions for the large-scale cosmological power deficit from primordial quantum nonequilibrium, Int. J. Mod. Phys. D 25 (2016) 1650068.  arXiv:1510.03508.
  * [28] S. D. P. Vitenti, P. Peter and A. Valentini, Modeling the large-scale power deficit with smooth and discontinuous primordial spectra, Phys. Rev. D100 (4) (2019) 043506.  arXiv:1901.08885.
  * [29] N. G. Underwood and A. Valentini, Quantum field theory of relic nonequilibrium systems, Phys. Rev. D 92 (2015) 063531.  arXiv:1409.6817.
  * [30] N. G. Underwood and A. Valentini, Anomalous spectral lines and relic quantum nonequilibrium, Phys. Rev. D 101 (4) (2020) 043004.  arXiv:1609.04576.
  * [31] A. Valentini, Black holes, information loss, and hidden variables (2004).  arXiv:hep-th/0407032.
  * [32] A. Kandhadai and A. Valentini, Mechanism for nonlocal information flow from black holes, Int. J. Mod. Phys. A 35 (06) (2020) 2050031.  arXiv:1912.05374.
  * [33] A. Valentini, Quantum gravity and quantum probability (2021).  arXiv:2104.07966.
  * [34] N. Pinto-Neto, Bouncing quantum cosmology, Universe 7 (4) (2021).  doi:10.3390/universe7040110.
URL https://www.mdpi.com/2218-1997/7/4/110

  * [35] N. Pinto-Neto, G. Santos and W. Struyve, Quantum-to-classical transition of primordial cosmological perturbations in de broglie–bohm quantum theory, Phys. Rev. D 85 (2012) 083506.  arXiv:1309.2670.
  * [36] N. Pinto-Neto, J. C. Fabris, Quantum cosmology from the de Broglie-Bohm perspective, Class. Quant. Grav. 30 (2013) 143001.  arXiv:1306.0820.
  * [37] D. C. F. Celani, N. Pinto-Neto, S. D. P. Vitenti, Particle Creation in Bouncing Cosmologies, Phys. Rev. D 95 (2) (2017) 023523.  arXiv:1610.04933.
  * [38] N. Pinto-Neto, W. Struyve, Bohmian quantum gravity and cosmology, in: Applied Bohmian Mechanics: From Nanoscale Systems to Cosmology, eds. Xavier Oriols Pladevall and Jordi Mompart, Jenny Stanford Publishing, 2019\.  arXiv:1801.03353.
  * [39] P. Davies, The Physics of Time Asymmetry, University of California Press, Berkeley and Los Angeles, 1977.
  * [40] P. Ehrenfest, T. Ehrenfest, The Conceptual Foundations of the Statistical Approach in Mechanics, Dover Books on Physics, Dover Publications, 2014.
URL https://books.google.com.br/books?id=g0xoBQAAQBAJ

  * [41] S. Colin, W. Struyve, Quantum non-equilibrium and relaxation to equilibrium for a class of de Broglie-Bohm-type theories, New J. Phys. 12 (4) (2010) 043008.  arXiv:0911.2823.
  * [42] N. G. Underwood, Extreme quantum nonequilibrium, nodes, vorticity, drift and relaxation retarding states, J. Phys. A: Math. Theor. 51 (5) (2018) 055301.  arXiv:1705.06757.
  * [43] R. Allahverdi, R. Brandenberger, F.-Y. Cyr-Racine, A. Mazumdar, Reheating in inflationary cosmology: Theory and applications, Annu. Rev. Nucl. Part. Sci. 60 (1) (2010) 27–51.  arXiv:1001.2600.
  * [44] A. Valentini, Trans-Planckian fluctuations and the stability of quantum mechanics (2014).  arXiv:1409.7467.
  * [45] N. Pinto-Neto, The Bohm interpretation of Quantum Cosmology, Found. Phys. 35 (4) (2005) 577–603.  arXiv:0410117.
  * [46] P. Peter, N. Pinto-Neto and S. D. P. Vitenti, Quantum cosmological perturbations of multiple fluids, Phys. Rev. D 93 (2016) 023520.  arXiv:1510.06628.
  * [47] H. Frisk, Properties of the trajectories in Bohmian mechanics, Phys. Lett. A 227 (3) (1997) 139 – 142.  doi:10.1016/S0375-9601(97)00044-3.
  * [48] C. Efthymiopoulos, C. Kalapotharakos and G. Contopoulos, Nodal points and the transition from ordered to chaotic Bohmian trajectories, J. Phys. A: Math. Theor. 40 (43) (2007) 12945–12972.  arXiv:0709.2038.
  * [49] A. Kandhadai and A. Valentini, Perturbations and quantum relaxation, Found. Phys. 49 (1) (2019) 1–23.  arXiv:1609.04485.
  * [50] D. A. Wisniacki and E. R. Pujals, Motion of vortices implies chaos in Bohmian mechanics, Europhys. Lett. 71 (2) (2005) 159–165.  arXiv:quant-ph/0502108.
  * [51] D. A. Wisniacki, E. R. Pujals and F. Borondo, Vortex interaction, chaos and quantum probabilities, Europhys. Lett. 73 (5) (2006) 671–676.  arXiv:nlin/0507015.
  * [52] E. Hairer, S. P. Nørsett and G. Wanner, Solving Ordinary Differential Equations I (2nd Revised. Ed.): Nonstiff Problems, Springer-Verlag, Berlin, Heidelberg, 1993.
  * [53] D. Wands, Multiple field inflation, Lect. Notes Phys. 738 (2008) 275–304.  arXiv:astro-ph/0702187, doi:10.1007/978-3-540-74353-8_8.