Dynamical Thermalization of Interacting Fermionic Atoms in a Sinai Oscillator Trap

Abstract

We study numerically the problem of dynamical thermalization of interacting cold fermionic atoms placed in an isolated Sinai oscillator trap. This system is characterized by a quantum chaos regime for one-particle dynamics. We show that, for a many-body system of cold atoms, the interactions, with a strength above a certain quantum chaos border given by the Åberg criterion, lead to the Fermi–Dirac distribution and relaxation of many-body initial states to the thermalized state in the absence of any contact with a thermostate. We discuss the properties of this dynamical thermalization and its links with the Loschmidt–Boltzmann dispute.

1. Introduction

The problem of emergence of thermalization in dynamical systems started from the Loschmidt–Boltzmann dispute about time reversibility and thermalization in an isolated system of moving and colliding classical atoms [1,2] (see the modern overview in [3,4]). The modern resolution of this dispute is related to the phenomenon of dynamical chaos where an exponential instability of motion breaks the time reversibility at infinitely small perturbation (see e.g., [5,6,7,8]). The well known example of such a chaotic system is the Sinai billiard in which a particle moves inside a square box with an internal circle colliding elastically with all boundaries [9].

The properties of one-particle quantum systems, which are chaotic in the classical limit, have been extensively studied in the field of quantum chaos during the last few decades and their properties have been mainly understood (see, e.g., [10,11,12]). Thus, it was shown that the level spacing statistics in the regime of quantum chaos [13] is the same as for Random Matrix Theory (RMT) invented by Wigner for a description of spectra of complex nuclei [14,15]. This result became known as the Bohigas–Giannoni–Schmit conjecture [13,16]. Thus, classically chaotic systems ( e.g., Sinai billiard) are usually characterized by Wigner–Dyson (RMT) statistics with level repulsion [13,14,15] while the classically integrable systems usually show Poisson statistics for level spacing distribution [11,12,16]. In this way, the level spacing statistics gives a direct indication for ergodicity (Wigner–Dyson statistics) or non-ergodicity (Poisson statistics) of quantum eigenstates. It was also established that the classical chaotic diffusion can be suppressed by quantum interference effects leading to an exponential localization of eigenstates [17,18,19,20] being similar to the Anderson localization in disordered solid-state systems [21]. The localized phase is characterized by Poisson statistics and the delocalized or metallic phase has RMT statistics. For billiard systems, the localized (nonergodic) and delocalized (ergodic) regimes appear in the case of rough billiards as described in [22,23].

It was also shown that, in the regime of quantum chaos, the Bohr correspondence principle [24] and the fully correct semiclassical description of quantum evolution remain valid only for a logarithmically short Ehrenfest time scale $t_E\sim \ln (1/\hslash )/h$ [17,19]. Here, ℏ is an effective dimensionless Planck constant and h in the Kolmogorov–Sinai entropy characterizing the exponential divergence of classical trajectories. This result is in agreement with the Ehrenfest theorem, which states that the classical-quantum correspondence works on a time scale during which the wave packet remains compact [25]. However, for the classically chaotic systems, the Ehrenfest time scale is rather short due to an exponential instability of classical trajectories. After the Ehrenfest time scale $t_E$, the quantum out-of-time correlations (or OTOC as it is used to say now) stop decaying exponentially in contrast to exponentially decaying classical correlators [26,27]. For $t>t_E$, the decay of quantum correlations stops and they remain on the level of quantum fluctuations being proportional to ℏ [26,27,28]. Since the level of quantum fluctuations is proportional to ℏ, the classical diffusive spreading over the momentum is affected by quantum corrections only on a significantly larger diffusive time scale $t_D\propto 1/{\hslash }^2\gg t_E\propto \ln (1/\hslash )$ [17,19,26,27,28].

The problem of the emergence of RMT statistics and quantum ergodicity in many-body quantum systems is more complex and intricate as compared to one-particle quantum chaos. Indeed, it is well known that, in many-body quantum systems, the level spacing between nearest energy levels drops exponentially with the increase of number of particles or with energy excitation $\delta E$ above the Fermi level in finite size Fermi systems, e.g., in nuclei [29]. Thus, at first glance, it seems that an exponentially small interaction between fermions should mix many-body quantum levels leading to RMT level spacing statistics (see, e.g., [30]).

Furthermore, the size of the Hamiltonian matrix of a many-body system grows exponentially with the number of particles, but, since all interactions have a two-body nature, the number of nonzero interaction elements in this matrix does not grow faster than the number of particles in the fourth power. Thus, we have a very sparse matrix being rather far from the RMT type. A two-body random interaction model (TBRIM) was proposed in [31,32] to consider the case of generic random two-body interactions of fermions in the limiting case of strong interactions when one-particle orbital energies are neglected. Even if the TBRIM matrix is very sparse, it was shown that the level spacing statistics $p\left(s\right)$ is described by the Wigner–Dyson or RMT distribution [33,34].

However, it is also important to analyze another limiting case when the two-body interaction matrix elements of strength U are weak or comparable with one-particle energies with an average level spacing ${\mathsf{\Delta }}_1$. In metallic quantum dots, this case with $U/{\mathsf{\Delta }}_1\approx 1/g$ corresponds to a large conductance of a dot $g=E_{Th}/{\mathsf{\Delta }}_1\gg 1$, where $E_{Th}=\hslash /t_D$ is the Thouless energy with $t_D$ being a diffusion spread time over the dot [35,36,37]. In this case, the main question is about critical interaction strength U or excitation energy $\delta E$ above the Fermi level of the dot at which the RMT statistics becomes valid. First, numerical results and simple estimates for a critical interaction strength in a model similar to TBRIM were obtained by Sven Åberg in [38,39]. The estimate of a critical interaction $U_c$, called the Åberg criterion [40], compares the typical two-body matrix elements with the energy level spacing ${\mathsf{\Delta }}_c$ between quantum states directly coupled by two-body interactions. Thus, the Åberg criterion tells that the Poisson statistics is valid for many-body energy levels for $U<U_c\sim {\mathsf{\Delta }}_c$ and the RMT statistics sets in for $U>U_c\sim {\mathsf{\Delta }}_c$. In [41], this criterion, proposed independently of [38,39], was applied to the TBRIM of weakly interacting fermions in a metallic quantum dot being confirmed by extensive numerical simulations. It was also argued that the dynamical thermalization sets in an isolated finite fermionic system for energy excitations $\delta E$ above the critical border $\delta E_{ch}$ determined from the above criterion [41]:

$$\delta E>\delta E_{ch}\approx g^{2/3}{\mathsf{\Delta }}_1,g={\mathsf{\Delta }}_1/U.$$

(1)

The emergence of thermalization in an isolated many-body system induced by interactions between particles without any contact with an external thermostat represents the Dynamical Thermalization Conjecture (DTC) proposed in [41]. The validity of the Åberg criterion was numerically confirmed for various physical models (see [40] and references therein). An additional confirmation was given by the analytical derivation presented in [42] showing that, for three interacting particles, in a metallic dot the RMT sets in when the two-body matrix elements U become larger than the two-particle level spacing ${\mathsf{\Delta }}_2\sim {\mathsf{\Delta }}_c$ being parametrically larger than the three-particle level spacing ${\mathsf{\Delta }}_3\ll {\mathsf{\Delta }}_2$. The advanced theoretical arguments developed in [43,44] confirm the relation (1) for interacting fermions in a metallic quantum dot.

The test for the transition from Poisson to RMT statistics is rather direct and needs only the knowledge of energies’ eigenvalues. However, the verification of DTC is much more involved since it requires the computation of system eigenstates. Thus, it is much more difficult to check numerically the relation (1) for DTC. However, it is possible to show that there is a transition from non-thermalized eigenstates at weak interactions (presumably for $\delta E<\delta E_{ch}$) to dynamically thermalized individual eigenstates at relatively strong interactions (presumably for $\delta E>\delta E_{ch}$). Thus, for the TBRIM with fermions, the validity of DTC for individual eigenstates at $U>U_c\sim {\mathsf{\Delta }}_c$ has been demonstrated in [45,46] by the computation of energy E and entropy S of each eigenstate and its comparison with the theoretical result given by the Fermi–Dirac thermal distribution [47].

Even if the TBRIM represents a useful system for DTC tests, it is not so easy to realize it in real experiments. Thus, in this work, we investigate the DTC features in a system of cold fermionic atoms placed in the Sinai oscillator trap created by a harmonic two-dimensional potential with a repulsive circular potential created by a laser beam in a vicinity of the trap center. In such a case, the repulsive potential in the center is modeled as an elastic circle as in the case of Sinai billiard [9]. For one particle, it has been shown in [48] that the Sinai oscillator has an almost fully chaotic phase space and that, in the quantum case, the level spacing statistics is described by the RMT distribution. Due to one-particle quantum chaos in the Sinai oscillator, we expect that this system will have properties similar to the TBRIM. On the other side, the Sinai oscillator trap has been already experimentally realized with Bose–Einstein condensate of cold bosonic atoms [49,50,51]. At present, cold atom techniques allow for investigating various properties of cold interacting fermionic atoms [52,53] and we argue that the investigation of dynamical thermalization of such fermionic atoms, e.g., ${}^{6}Li$, in a Sinai oscillator trap is now experimentally possible. Thus, in this work, we study properties of DTC of interacting fermionic atoms in a Sinai oscillator trap. Here, we consider the two-dimensional (2D) case of such a system assuming that the trap frequency in the third direction is small and that the 2D dynamics are not significantly affected by the adiabatically slow motion in the third dimension.

Finally, we note that, at present, the TBRIM model in the limit of strong interactions attracts a high interest in the context of field theory since, in this limit, it can be mapped on a black hole model of quantum gravity in $1+1$ dimensions known as the Sachdev–Ye–Kitaev (SYK) model linked also to a strange metal [54,55,56,57,58,59]. In fact, the SYK model, in its fermionic formulation [56], corresponds to the TBRIM considered with a conductance close to zero $g\to 0$. In these lines, the dynamical thermalization in TBRIM and SYK systems has been discussed in [45,46]. Furthermore, there is also a growing interest in dynamical thermalization for various many-body systems known also as the eigenstate thermalization hypothesis (ETH) and many-body localization (MBL) (see, e.g., [60,61,62,63]). We think that the system of interacting fermionic atoms in a Sinai oscillator trap captures certain features of TBRIM and SYK models and thus represents an interesting test ground to investigate nontrivial physics of these systems in real cold atom experiments.

This paper is composed as follows: in Section 2, we describe the properties of the one-particle dynamics in a Sinai oscillator; numerical results for dynamical thermalization on interacting atoms in this oscillator are presented in Section 3; the conditions of thermalization for fermionic cold atoms in realistic experiments are given in Section 4; the discussion of the results is presented in Section 5.

2. Quantum Chaos in Sinai Oscillator

The model of one particle in the 2D Sinai oscillator is described in detail in [48] with the Hamiltonian:

$$H_1=\frac{1}{2m}(p_x^2+p_y^2)+\frac{m}{2}({\omega }_x^2{x}^2+{\omega }_y^2{y}^2)+V_d(x,y).$$

(2)

Here, the first two terms describe the 2D oscillator with frequencies ${\omega }_x,{\omega }_y$ and the last term gives the potential wall of elastic disk of radius $r_d$. We choose the dimensionless units with mass $m=1$, frequencies ${\omega }_x=1$, ${\omega }_y=\sqrt{2}$ and disk radius $r_d=1$. The disk center is located at $(x_d,y_d)=(-1/2,-1/2)$ so that the disk bungs a hole in the center as it was the case in the experiments [49]. The Poincare sections at different energies are presented in [48] showing that the phase space is almost fully chaotic (see Figure 1 there). The quantum evolution is described by the Schrödinger equation with the quantized Hamiltonian (2), where the conjugate momentum and coordinate variables become operators with the commutation relation $[x,p_x]=[y,p_y]=i\hslash$ [48]. For the quantum problem, we use the value of the dimensionless Planck constant $\hslash =1$ so that the ground state energy is $E_g=1.685$. In the following, the energies are expressed in atomic like units of energy $E_u=\hslash {\omega }_x$ (for our choice of Sinai oscillator parameters, we also have $E_u=\hslash {\omega }_x={\hslash }^2/\left(m{r_d}^2\right)$) [48] with the typical size of oscillator ground state being equal to the disk radius: $a_0=\mathsf{\Delta }{x}_{osc}={(\hslash /m{\omega }_x)}^{1/2}=r_d$.

In [48], it is shown that the classical dynamics of this system are almost fully chaotic. In the quantum case, the level spacing statistics is well described by the RMT distribution. The average dependence of energy level number k is well described by the theoretical dependence $k\left(\epsilon \right)={\epsilon }^2/\left(2\sqrt{2}\right)-\epsilon /2$ [48]. Thus, the one-particle density of states ${\rho }_1\left(\epsilon \right)$ and corresponding level spacing ${\mathsf{\Delta }}_1$ are:

$${\rho }_1\left(\epsilon \right)=\frac{dk}{d\epsilon }=\frac{\epsilon }{\sqrt{2}}-\frac{1}{2},{\mathsf{\Delta }}_1=\frac{1}{\rho }\approx \frac{\sqrt{2}}{\epsilon }\approx \frac{0.84}{\sqrt{k}}.$$

(3)

Examples of several eigenstates, computed on a numerical grid of 28,341 spatial points, are shown in Figure 1. More details on the numerical diagonalization of (2) and other example eigenstates can be found in [48].

3. Sinai Oscillator with Interacting Fermionic Atoms

3.1. Two-Body Interactions of Fermionic Atoms

The two-body interaction of atoms appears usually due to van der Waals forces which drop rapidly with the distance between two atoms and the short ranged interaction can be described in the framework of the scattering length approach (see, e.g., [64,65]). Therefore, we assume that the finite effective interaction range $r_c$ is significantly smaller than the disk radius $r_d$ and the typical size of the wave function, i.e., $r_c\ll r_d$. Such a short range interaction is indeed used to modelize atomic interactions in harmonic traps (see, e.g., [66]). For example, in a typical experimental situation, the disk radius is of the order of micron $r_d\sim 1\mathsf{\mu }\mathrm{m}={10}^{-4}$ cm, while, for Li and other alkali atoms, we have $r_c\sim 3\times {10}^{-7}$ cm [64,65]. Of course, in the limit of small $r_c\ll r_d$, the interaction between two atoms takes place mainly in the s-wave scattering, so effectively the interaction operates between fermions with different quantum numbers of the Sinai oscillator. This feature is of course taken into account in our numerical simulations.

In the following, we use a simple interaction function having a constant amplitude U for $r\le r_c$ and being zero for $r>r_c$, where we simply choose $r_c=0.2r_d$, which corresponds well to the short range interaction regime. The precise value of $r_c$ is not very important since a slight modification $r_c\to {\overline{r}}_c$ can be absorbed in a modified amplitude according to $U\to \overline{U}=U{({\overline{r}}_c/r_c)}^2$, a relation we verified numerically for various values of ${\overline{r}}_c<r_d$. We numerically checked that, at a used $r_c=0.2r_d$ value, we are in the regime when the interaction matrix elements are proportional to ${r_c}^2$ so that we are in the regime of short-range interactions. We mention that, in experiments, the strength of the interaction amplitude can be changed by a variation of the magnetic field via the Feshbach resonance [67].

3.2. Reduction to TBRIM Like Case and Its Analysis

Using the methods described in [48], we numerically compute a certain number of one-particle or orbital energy eigenvalues ${\epsilon }_k$ and corresponding eigenstates ${\phi }_k\left(\mathbf{r}\right)$ of the Sinai oscillator (2). As repulsive interaction potential $v\left(\mathbf{r}\right),$ we choose the short ranged box function $v\left(\mathbf{r}\right)=U$ if $\left|\mathbf{r}\right|\le r_c=0.2$ (since $r_d=1$) and $v\left(\mathbf{r}\right)=0$, otherwise. Here, the parameter $U>0$ gives the overall scale of the interaction strength depending on the charge of the particles and eventually other physical parameters.

Therefore, the corresponding many-body Hamiltonian with M one-particle orbitals and $0\le L\le M$ spinless fermions takes the form:

$$H=\sum _{k=1}^M{\epsilon }_k{c}_k^{†}{c}_k^{\phantom{†}}+\sum _{i<j,k<l}{V}_{ij,kl}{c}_i^{†}{c}_j^{†}{c}_l^{\phantom{†}}{c}_k^{\phantom{†}},$$

(4)

where, for $i<j$ and $k<l$, we have the interaction matrix elements:

$$V_{ij,kl}={\overline{V}}_{ij,kl}-{\overline{V}}_{ij,lk},{\overline{V}}_{ij,kl}=\int d{\mathbf{r}}_{\mathbf{1}}\int d{\mathbf{r}}_{\mathbf{2}}{\phi }_i^{*}\left({\mathbf{r}}_{\mathbf{1}}\right){\phi }_j^{*}\left({\mathbf{r}}_{\mathbf{2}}\right)v({\mathbf{r}}_{\mathbf{1}}-{\mathbf{r}}_{\mathbf{2}}){\phi }_k\left({\mathbf{r}}_{\mathbf{1}}\right){\phi }_l\left({\mathbf{r}}_{\mathbf{2}}\right)$$

(5)

and $c_k^{†}$, $c_k^{\phantom{†}}$ are fermion operators for the M orbitals satisfying the usual anticommutation relations. We note that, in the literature, when expressing a two-body interaction potential in second quantization, one usually uses the raw matrix elements ${\overline{V}}_{ij,kl}$ with an additional prefactor of $1/2$ and full independent sums for the four indices i, j, k and l. Using the particle exchange symmetry: ${\overline{V}}_{ij,kl}={\overline{V}}_{ji,lk}$, one can reduce the $i,j$ sums to $i<j$, which removes the prefactor $1/2$ (after exchanging the index names l and k for the $i>j$ contributions and exploiting that contributions at $i=j$ or $l=k$ obviously vanish). The definition of the anti-symmetrized interaction matrix elements $V_{ij,kl}$ according to (5) allows for also reducing the $k,l$ sums to $k<l$. Furthermore, the ordering of the two fermion operators $c_l^{\phantom{†}}{c}_k^{\phantom{†}}$ in (4) is also important and necessary to obtain positive expectation values if the interaction is repulsive. The anti-symmetrized matrix elements $V_{ij,kl}$ correspond to a $M_2\times M_2$ matrix with $M_2=M(M-1)/2)$. In order to avoid a global shift of the non-interacting eigenvalue spectrum due to the interaction, we also apply a diagonal shift $V_{ij,ij}\to V_{ij,ij}-(1/M_2){\sum }_{k<l}{V}_{kl,kl}$ to ensure that this matrix has a vanishing trace ( One can easily show that the trace of the $M_2\times M_2$ anti-symmetrized interaction matrix is proportional to the trace of the interaction operator in the many-body Hilbert space with a factor depending on M and L). Of course, such a global energy shift does not affect the issues of thermalization, interaction induced eigenfunction mixing or the quantum time evolution with respect to the Hamiltonian H, etc.

We note that the transition from the classical Hamiltonian (2) to the quantum one (4) is done by the standard procedure of second quantization (see, e.g., [68]).

3.3. Aberg Parameter

In absence of interaction, the energy eigenvalues of (4) are given as the sum of occupied orbital energies:

$$E\left(\left\{n_k\right\}\right)=\sum _{k=1}^M{\epsilon }_k{n}_k,$$

(6)

where $\left\{n_k\right\}$ represents a configuration such that $n_k\in \{0,1\}$ and ${\sum }_k{n}_k=L$. The associated eigenstates are the basis states where each orbital is either occupied (if $n_k=1$) or unoccupied (if $n_k=0$) and, in this work, we will denote these states in the usual occupation number representation: $|n_M\cdots n_2{n}_1>$, where, for convenience, we write the lower index orbitals starting from the right side.

The distribution of the total one-particle energies (6) is numerically rather close to a Gaussian (since $n_k$ act as quasi-random numbers) with mean and variance (see also Equation (A.4) of Ref. [46]):

$$E_{\mathrm{mean}}=L\overline{\epsilon },{\sigma }_0^2=\frac{L(M-L)}{M-1}\left(\overline{{\epsilon }^2}-{\overline{\epsilon }}^2\right),\overline{{\epsilon }^n}=\frac{1}{M}\sum _{k=1}^M{\epsilon }_k^n,n=1,2.$$

(7)

Therefore, the many-body level spacing ${\mathsf{\Delta }}_{\mathrm{MB}}$ or inverse Heisenberg time at the band center $E=E_{\mathrm{mean}}$ is given by ${\mathsf{\Delta }}_{\mathrm{MB}}=1/t_{\mathrm{H}}=\sqrt{2\pi }({\sigma }_0/d),$ where $d=M!/(L!(M-L)!)$ is the dimension of the fermion Hilbert space in the sector of M orbitals and L particles. In our numerical computations, we simply evaluated the quantities $\overline{{\epsilon }^n}$ of (7) using the exact one-particle energy eigenvalues obtained from the numerical diagonalization of the one-particle Sinai Hamiltonian $H_1$ given in (2). However, to get some analytical simplification for large M, one may use the one-particle density of states (3), which gives, after replacing the sums by integrals and neglecting the constant term, $\overline{{\epsilon }^n}\approx 2{\epsilon }_M^n/(n+2)$ and $\overline{{\epsilon }^2}-{\overline{\epsilon }}^2\approx {\epsilon }_M^2/18\approx \sqrt{2}M/9$.

For the question of whether the interaction strength is sufficiently strong to mix the non-interacting basis states, the important quantity is the effective level spacing of states coupled directly by the interaction ${\mathsf{\Delta }}_c=\sqrt{2\pi }[{\sigma }_0(L=2)/K],$ where $K=1+L(M-L)+L(L-1)(M-L)(M-L-1)/4$ is the number of nonzero elements for a column (or row) of H [41,69] and we need to use the variance for only two particles:

$${\sigma }_0^2(L=2)=\frac{2(M-2)}{M-1}\left(\overline{{\epsilon }^2}-{\overline{\epsilon }}^2\right)\Rightarrow \frac{{\sigma }_0^2(L=2)}{{\sigma }_0^2}=\frac{2(M-2)}{L(M-L)}$$

(8)

because the interaction only couples states where (at least) $L-2$ particles are on the same orbital such that (at most) only the partial sum of two one-particle energies is different between two coupled states. Even though for two particles the hypothesis of a Gaussian distribution is theoretically not justified, the distribution is still sufficiently similar to a Gaussian and it turns out that the value of $1/{\mathsf{\Delta }}_c=K/\left[\sqrt{2\pi }{\sigma }_0(L=2)\right]$ as the coupled two-particle density of states in the band center is numerically quite accurate with an error below 10% (for $M=16$ and our choice of ${\epsilon }_k$ values).

According to the Åberg criterion [38,39,41], the onset of chaotic mixing happens for typical interaction matrix elements U comparable to ${\mathsf{\Delta }}_c$. Therefore, we compute the quantity $V_{\mathrm{mean}}=\sqrt{\langle |V_{ij,kl}{|}^2\rangle }$ (which is proportional to the interaction amplitude U) where the average is done with respect to all $M_2^2$ matrix elements of the interaction matrix. This quantity might be problematic and not correspond to a typical interaction matrix element in the case of a long tail distribution. However, in our case, it turns out that $V_{\mathrm{mean}}\approx 2exp(\langle \ln |V_{ij,kl}|\rangle )$, which excludes this scenario. Using this quantity, we introduce the dimensionless Åberg parameter and the critical interaction amplitude $U_c$ by $A=V_{\mathrm{mean}}/{\mathsf{\Delta }}_c=U/U_c$ such that $A=1$ if $U=U_c$. We expect [38,39,41] to be the onset of strong/chaotic mixing at $A\gg 1$ and a perturbative regime for $A\ll 1$, while, at $A=1$, we have the critical interaction strength $U=U_c$. The value of $U_c$ depends on the parameters L, M, ${\sigma }_0$ and the overlap of the one-particle eigenstates according to (5). To obtain some useful analytical expression of $U_c$, we note that the quantity $V_{\mathrm{mean}}$, numerically computed for $4\le M\le 30$, can be quite accurately fitted by $V_{\mathrm{mean}}\approx 3\times {10}^{-4}U/{\epsilon }_M$. Furthermore, we remind readers of the expression ${\mathsf{\Delta }}_c=(1/K)\sqrt{4\pi (M-2)(\overline{{\epsilon }^2}-{\overline{\epsilon }}^2)/(M-1)},$ which can be simplified in the limit $M\gg 1$ and $L\gg 1$, such that $K\approx {(M-L)}^2{L}^2/4$, resulting in: ${\mathsf{\Delta }}_c=4/3\sqrt{2\pi }{\epsilon }_M/[{(M-L)}^2{L}^2]$. Here, we also used the found expression above $\overline{{\epsilon }^2}-{\overline{\epsilon }}^2\approx {\epsilon }_M^2/18$. From this, we find that $U_c={\mathsf{\Delta }}_{c}U/V_{\mathrm{mean}}\approx CM/[{(M-L)}^2{L}^2]$ with a numerical constant $C\approx 16\times {10}^4\sqrt{\pi }/9\approx 3.15\times {10}^4$, where we also used ${\epsilon }_M^2\approx 2\sqrt{2}M$ according to (3). Below, we will give more accurate numerical values of $V_{\mathrm{mean}}$, ${\mathsf{\Delta }}_c$ and $U_c$ for the parameter choice of M and L numerically relevant in this work.

We note that this estimate for $A=U/U_c$ applies to energies close to the many body band center of H and that, for energies away from the band center, the value of ${\mathsf{\Delta }}_c$ is enhanced, thus reducing the effective value of A. Furthermore, we computed $V_{\mathrm{mean}}$ by a simplified average over all interacting matrix elements not taking into account an eventual energy dependence according to the index values of $i,j,k,l$ in (5).

3.4. Density of States

In this work, we present numerical results for the case of $M=16$ orbitals and $L=7$ particles corresponding to a many-body Hilbert space of dimension $d=M!/(L!(M-L)!)=11440$ and the number $K=820$ of directly coupled states of a given initial state by non-vanishing interaction matrix elements in (4). Thus, in our studies, the whole Hilbert space is built only on these $M=16$ orbitals. We diagonalize numerically the many-body Hamiltonian (4) for various values of A in the range $0.025\le A\le 200$. We have also verified that the results and their physical interpretation are similar for smaller cases such as $M=12$, $L=5$ (with $d=792$, $K=246$) or $M=14$, $L=6$ ($d=3003$, $K=469$).

We mention that, for $M=16$ and $L=7$, we find numerically that $V_{\mathrm{mean}}=3.865\times {10}^{-5}U$ and, from (8), that ${\mathsf{\Delta }}_c=\sqrt{2\pi }[{\sigma }_0(L=2)/K]=6.1706\times {10}^{-3}$, where the quantities $\overline{{\epsilon }^n}$ were exactly computed from the numerical orbitals energies ${\epsilon }_k$. From this, we find that $U_c={\mathsf{\Delta }}_{c}U/V_{\mathrm{mean}}\approx 159.65$. This expression is more accurate than the more general analytical estimate for arbitrary $M\gg 1$ and $L\gg 1$ given in the last section (which would provide $U_c\approx 127$ for $M=16$ and $L=7$).

Our first observation is that, even in the presence of interactions, the density of states has approximately a Gaussian form with the same center $E_{\mathrm{mean}}$ given in (7) for the case $A=0$. This is simply due the fact that the interaction matrix has, by choice, a vanishing trace and does not provide a global shift of the spectrum. We determine the variance ${\sigma }^2\left(A\right)$ of the Gaussian density of states by a fit of the integrated density of states $P\left(E\right)$ using

$$P\left(E\right)=(1+\mathrm{erf}\left[q\left(E\right)\right])/2,q\left(E\right)=(E-E_{\mathrm{mean}})/\left[\sqrt{2}\sigma \left(A\right)\right],$$

(9)

such that $P^{\prime }\left(E\right)$ is a Gaussian of width $\sigma \left(A\right)$ and center $E_{\mathrm{mean}}$ (see Appendix A of Ref. [46] for more details). From this, we find the behavior:

$${\sigma }^2\left(A\right)={\sigma }_0^2(1+\alpha A^2),$$

(10)

where $\alpha$ is a constant depending on M and L; for $M=16$, $L=7$ the fit values of ${\sigma }_0$ and $\alpha$ are ${\sigma }_0=3.013\pm 0.009$ and $\alpha =0.00877\pm 0.00010$. It is also possible to determine $\sigma \left(A\right)$ using the expression ${\sigma }^2\left(A\right)={Tr}_{\mathrm{Fock}}[{(H-E_{\mathrm{mean}}\mathbf{1})}^2]/d$ where the trace in Fock space can be evaluated either by using the matrix H before diagonalizing it or using its exact energy eigenvalues $E_m$. This provides the same behavior as (10) with the very similar numerical values ${\sigma }_0=3.013\pm 0.007$ and $\alpha =0.00858\pm 0.00008$ (for $M=16$, $L=7$). We mention that the integrated Gaussian density of states (9) is not absolutely exact but quite accurate for values $A\le 10$. For larger values of A, the deviations increase, but the overall form is still correct. As described in [46], the quality of the fit can be considerably improved if we replace in (9) the linear function $q\left(E\right)$ by a polynomial of degree 5. In this case, the precision of the fit is highly accurate for the full range of A values we consider. In particular, we use this improved fit to perform the spectral unfolding when computing the nearest level spacing distribution (shown below).

To obtain some theoretical understanding of (10), one can consider a model where the initial interaction matrix elements (5) are replaced by independent Gaussian variables with identical variance $V_{\mathrm{mean}}^2$. In this case, one can show theoretically [46] that ${\sigma }^2\left(A\right)={\sigma }_0^2+K_2{V}_{\mathrm{mean}}^2$, where $K_2=L(L-1)[1+M-L+(M-L)(M-L-1)/4]$ is a number somewhat larger than K taking into account that certain non-vanishing interaction matrix elements in Fock space are given as a sum of several initial interaction matrix elements (5) (see Appendix A of [46] for details). The parameter $K_2$ takes for $M=16$, $L=7$ ($M=14$, $L=6$ or $M=12$, $L=5$) the value $K_2=1176$ ($K_2=690$ or $K_2=370$, respectively). Since $V_{\mathrm{mean}}=A{\mathsf{\Delta }}_c=A\sqrt{2\pi }{\sigma }_0(L=2)/K,$ we indeed obtain (10) with $\alpha ={\alpha }_{\mathrm{th}}=4\pi (M-2)K_2/[K^2\left(L(M-L)\right]$. For $M=16$, $L=7$, we find ${\sigma }_0=3.0279$ (see (7)) and ${\alpha }_{\mathrm{th}}=0.00488$. The latter is roughly by a factor of 2 smaller than the numerical value. We attribute this to the fact that the real initial interaction matrix elements (5) are quite correlated, and not independent uniform Gaussian variables, leading therefore to an effective increase of the number $K_2$ due to hidden correlations. The important point is that theoretically at very large values values of M and L, e.g., $M\approx 2L\gg 1$, we have $K_2\approx K\approx L^4/4$ and ${\alpha }_{\mathrm{th}}\approx 32\pi /L^5$, which is parametrically small for very large L. Therefore, there is a considerable range of values $1<A<1/\sqrt{\alpha }$ where the interaction strongly mixes the non-interacting many-body eigenstates but where the density of states is only weakly affected by the interaction. This regime is also known as the Breit–Wigner regime (see, e.g., [40], for the case of interacting Fermi systems).

3.5. Thermalization and Entropy of Eigenstates

In the following, we mostly concentrate on values $A\le 10$ such that the effect of the increase of the spectral width $\sigma \left(A\right)$ is still small or at least quite moderate. The question arises if a given many-body state, either an exact eigenstate of H or a state obtained from a time evolution with respect to H, is thermalized according to the Fermi–Dirac distribution [47]. As in [45,46], we determine the occupation numbers $n_k=\langle c_k^{†}{c}_k^{\phantom{†}}\rangle$ for such a state, as well as the corresponding fermion entropy S [47] and the effective total one-particle energy $E_{1p}$ by:

$$S=-\sum _{k=1}^M\left(n_k\ln n_k+(1-n_k)\ln (1-n_k)\right),E_{1p}=\sum _{k=1}^M{\epsilon }_k{n}_k$$

(11)

based on the assumption of weakly interacting fermions. In the regime of modest interaction $A\lesssim 5$ (for $M=16$, $L=7$), corresponding to a constant spectral width $\sigma \left(A\right)\approx {\sigma }_0$, we have typically $E_{1p}\approx E_{\mathrm{ex}}$ (for exact eigenstates of H) or $E_{1p}\approx \langle H\rangle$ (for other states). If the given state is thermalized, its occupation numbers $n_k$ should be close to the theoretical Fermi–Dirac filling factor $n\left({\epsilon }_k\right)$ with $n\left(\epsilon \right)=1/(1+exp[\beta (\epsilon -\mu )\left]\right)$, where inverse temperature $\beta =1/T$ and chemical potential $\mu$ are determined by the conditions:

$$L=\sum _{k=1}^{M}n\left({\epsilon }_k\right),E=\sum _{k=1}^M{\epsilon }_{k}n\left({\epsilon }_k\right).$$

(12)

Here, E is normally given by $E_{1p}$, but one may also consider the value $E_{\mathrm{ex}}$ (or $\langle H\rangle$) provided the latter is in the energy interval where the conditions (12) allow for a unique solution. Furthermore, for a given energy E, we can also determine the theoretical (or thermalized) entropy $S_{\mathrm{th}}\left(E\right)$ using (11) with $n_k$ being replaced by $n\left({\epsilon }_k\right)$ (where $\beta$, $\mu$ are determined from (12) for the energy E).

The many-body states with energies above $E_{\mathrm{mean}}$ are artificial since they correspond to negative temperatures due to the finite number of orbitals considered in our model. Therefore, we limit our studies to the lower half of the energy spectrum $29\le E\le 39\approx E_{\mathrm{mean}}$ (for $M=16$, $L=7$). In Figure 2, we compare the thermalized Fermi–Dirac occupation number $n\left(\epsilon \right)$ with the occupation numbers $n_k$ for two eigenstates at level numbers $m=123$ (1354, with $m=1$ corresponding to the ground state) with approximate energy eigenvalue $E\approx 32$ ($E\approx 35$) for three different Åberg parameters $A=0.35$, $A=3.5$ and $A=10$. These states are not too close to the ground state but still quite far below the band center.

We note that the regime of negative temperatures is natural for the TBRIM where the energy spectrum is inside a finite energy band (this regime has been discussed in [45,46]). However, for the Sinai oscillator, the energy spectrum is unbounded and, due to that, the regime of negative temperatures, appearing in the numerical simulations due to a finite number of one-particle orbital, is artificial.

At weak interaction, $A=0.35$, both states are not at all thermalized with occupation numbers being either close to 1 or 0. Apparently, these states result from weak perturbations of the non-interacting eigenstates $|0000011000110111>$ or $|1000100011001011>$, where the $n_k$ values are rounded to 1 (or 0) if $n_k>0.5$ ($n_k<0.5$). For $m=1354$, the values of $n_k$ are a little bit farther away from the ideal values 0 or 1 as compared to $m=123$ but still sufficiently close to be considered as perturbative. Apparently, the state $m=123$, which is lower in the spectrum (with larger effective two-body level spacing), is less affected by the interaction than the state 1354. In both cases, the entropy S is quite below the thermalized entropy $S_{\mathrm{th}}$ (see

Table 1. Parameters of the eigenstates corresponding to Figure 2. S is the entropy, $E_{1p}$ the effective total one-particle energy, both given by (11), and $E_{\mathrm{ex}}$ is the exact energy eigenvalue. Inverse temperature $\beta$, chemical potential $\mu$, theoretical entropy $S_{\mathrm{th}}$ are determined by (12) or (11) (with $n_k$ replaced by $n\left({\epsilon }_k\right)$) for both energies $E_{1p}$, $E_{\mathrm{ex}}$.

A	m	S	${\mathit{S}}_{\mathbf{th}}\left({\mathit{E}}_{1\mathbf{p}}\right)$	${\mathit{E}}_{1\mathbf{p}}$	$\mathit{\mu }\left({\mathit{E}}_{1\mathbf{p}}\right)$	$\mathit{\beta }\left({\mathit{E}}_{1\mathbf{p}}\right)$	${\mathit{S}}_{\mathbf{th}}\left({\mathit{E}}_{\mathbf{ex}}\right)$	${\mathit{E}}_{\mathbf{ex}}$	$\mathit{\mu }\left({\mathit{E}}_{\mathbf{ex}}\right)$	$\mathit{\beta }\left({\mathit{E}}_{\mathbf{ex}}\right)$
0.35	122	0.95	7.91	32.15	5.31	1.05	7.89	32.13	5.31	1.05
0.35	1353	4.91	10.16	35.29	4.98	0.45	10.16	35.30	4.98	0.45
3.5	122	6.99	8.28	32.52	5.28	0.95	7.54	31.81	5.34	1.15
3.5	1353	10.16	10.23	35.45	4.95	0.43	10.10	35.15	5.00	0.47
10	122	8.91	8.98	33.33	5.22	0.77	4.96	30.10	5.46	2.02
10	1353	10.52	10.54	36.28	4.75	0.32	9.53	34.12	5.14	0.63

for numerical values of entropies, energies, inverse temperature and chemical potential for the states shown in Figure 2).

At intermediate interaction, $A=3.5$, the occupation numbers are closer to the theoretical Fermi–Dirac values but still with considerable deviations. Here, both entropy values S are rather close to $S_{\mathrm{th}}$. The state 1354 seems to be better thermalized than the state $m=123$, the latter having a slightly larger deviation between both entropy values. At stronger interaction, $A=10$, both states are very well thermalized with a good matching of both entropy values (again with the state 1354 being a bit better thermalized than the state $m=123$) provided we use $E_{1\mathrm{p}}$ as reference energy to compute temperature and chemical potential. The temperature obtained from $E_{\mathrm{ex}}$ is too small because here the increase of $\sigma \left(A\right)$ is already quite strong and $E_{\mathrm{ex}}$ rather strongly deviates from $E_{1\mathrm{p}}$. In addition, the value of $S_{\mathrm{th}}$ using $E_{\mathrm{ex}}$ does not match S. Obviously, at stronger interaction values, it is necessary to use $E_{1\mathrm{p}}$ to test the thermalization hypothesis of a given state.

Figure 3 shows the mutual dependence between the three parameters $\beta$, $\mu$ on E when solving the conditions (12). The chemical potential as a function of $\beta =1/T$ is rather constant except for smallest values of $\beta$ where $\mu \sim 1/\beta$ with a negative prefactor. One can actually easily show from (12) that in the limit $\beta \to 0$ the chemical potential does not depend on ${\epsilon }_k$ and is given by $\mu =-\ln [1+(M-2L)/L]/\beta$ providing a singularity if $L\ne M/2$ with negative (positive) prefactor for $L<M/2$ ($L>M/2$) and $\mu =0$ for $L=M/2$. The temperature (${\beta }^{-1}$) vanishes for E close to the lower energy border and diverges for E close to the band center $E_{\mathrm{mean}}$.

In Figure 4, we present the nearest level spacing distribution $p\left(s\right)$ for different values of the Åberg parameter. To compute $p\left(s\right)$, we have used only the “physical” levels in the lower half of the energy spectrum and the unfolding has been done with the integrated density of states (9), where $q\left(E\right)$ is replaced by a fit polynomial of degree 5. For the smallest value $A=0.2$, the distribution $p\left(s\right)$ is very close to the Poisson distribution with some residual level repulsion at very small spacings. This is a quite well known effect because typically the transition from Wigner–Dyson to Poisson statistics (when tuning some suitable parameter such as the Åberg parameter from strong to weak coupling) is non-uniform in energy and happens first at larger spacings (energy differences) and then at smaller spacings. The reason is simply that two levels which by chance are initially very close are easily repelled by a small residual coupling matrix element (when slightly changing a disorder realization or similar). For $A=0.5$, there is somewhat more level repulsion at small spacings, but the distribution is still rather close to the Poisson distribution with some modest deviations for $s\le 1.2$. For the larger Åberg values $A=3.5$ and $A=10$, we clearly obtain Wigner–Dyson statistics (taking into account the quite limited number of only $d/2-1=5719$ level spacing values for the histograms). These results clearly confirm that the transition from $A<1$ to $A>1$ corresponds indeed to a transition from a perturbative regime to a regime of chaotic mixing with Wigner–Dyson level statistics [14].

A further confirmation that $A=1$ is critical can be seen in Figure 5, which compares the dependence of the entropy S of exact eigenstates (lower half of the spectrum) on $E_{1\mathrm{p}}$ or $E_{\mathrm{ex}}$ with the theoretical thermalized entropy $S_{\mathrm{th}}\left(E\right)$. For the Åberg values $A=0.2$ (and $A=0.5$), the entropy S of all (most) states is significantly below its theoretical value $S_{\mathrm{th}}$. Actually, the distribution of data points is considerably concentrated at smaller entropy values which is not so clearly visible in Figure. In particular, the average of the ratio of $S/S_{\mathrm{th}}\left(E_{1\mathrm{p}}\right)$ is $0.178$ for $A=0.2$ and $0.522$ for $A=0.5$. For the Åberg values $A=3.5$ and $A=10$, most or nearly all entropy values (for $E_{1\mathrm{p}}$) are very close to the theoretical line with the average ratio $S/S_{\mathrm{th}}\left(E_{1\mathrm{p}}\right)$ being $0.990$ for $A=3.5$ and $0.998$ for $A=10$. For $A=3.5$, the states with lowest energies are not yet perfectly thermalized and the data points for $E_{\mathrm{ex}}$ and $E_{1\mathrm{p}}$ are still rather close. For $A=10$, all states are well thermalized (when using the energy $E_{1\mathrm{p}}$) while the data points for $E_{\mathrm{ex}}$ are quite outside the theoretical curve simply due to the overall increase of the width of the energy spectrum. This observation is also in agreement with the discussion of Figure 2. For smaller values $A<0.2$ (not shown in Figure 5), we find that the data points are still closer to the E-axis while, for larger values $A>10$, the data points are clearly on the theoretical curve for $E_{1\mathrm{p}}$ (but more concentrated on energy values closer to the center with larger entropy values and larger temperatures), while, for $E_{\mathrm{ex}}$, according to (10), the overall width of the exact eigenvalue spectrum increases strongly and the data points are clearly outside the theoretical curve (except for a few states close to the band center).

We note that the data points in Figure 5a,b significantly deviate from the theoretical thermalization curve since the Åberg criterion is not satisfied $(A<1)$. For $A=3.5;10$, the deviations are significantly reduced, but they are still more pronounced in a vicinity of the ground state in agreement with the relation (1).

In Figure 6, the occupation numbers $n_k$ (averaged over several energy eigenvalues inside a given energy cell) are shown in the plane of energy E and orbital index k as color density plot for the Åberg parameter $A=3.5$. The comparison with the theoretical occupation numbers $n\left({\epsilon }_k\right)$ (shown in the same way) provides further confirmation that, at $A=3.5$, there is indeed already a quite strong thermalization of most eigenstates.

3.6. Thermalization of Quantum Time Evolution

The question arises how or if a time dependent state $\left|\psi \right(t)>=exp(-iHt\left)\right|\psi \left(0\right)>$, obeying the quantum time evolution with the Hamiltonian H and an initial state $\left|\psi \right(0)>$ being a non-interacting eigenstate $|n_M\cdots n_2{n}_1)>$ (with all $n_k\in \{0,1\}$ and ${\sum }_k{n}_k=L$), evolves eventually into a thermalized state. We have computed such time dependent states using the exact eigenvalues and eigenvectors of H to evaluate the time evolution operator. As initial states, we have chosen four states (for $M=16$, $L=7$): (i) $|{\varphi }_1>=|0000100000111111>$ where a particle at orbital 7 is excited from the non-interacting ground state (with all orbitals from 1 to 7 occupied) to the orbital 12, (ii) $|{\varphi }_2>=|0010100000011111>$ where two particles at orbitals 6 and 7 are excited from the non-interacting ground state to the orbitals 12 and 14, (iii) $|{\varphi }_3>=|0000011000110111>$ and (iv) $|{\varphi }_4>=|1000100011001011>$. The states $|{\varphi }_3>$ and $|{\varphi }_4>$ are obtained from the exact eigenstate of H for $A=0.35$ at level number $m=123$ and 1354, respectively, by rounding the occupation numbers $n_k$ to 1 (or 0) if $n_k>0.5$ ($n_k<0.5$) (states of top panels in Figure 2). The approximate energies (6) of these four states are $E\approx 30$ ($|{\varphi }_1>$), $E\approx 32$ ($|{\varphi }_2>$ and $|{\varphi }_3>$) and $E\approx 35$ ($|{\varphi }_4>$).

It is useful to express the time in multiples of the elementary quantum time step defined as:

$$\mathsf{\Delta }t=\frac{t_{\mathrm{H}}}{d}=\frac{1}{\sqrt{2\pi {\sigma }^2\left(A\right)}},$$

(13)

where $t_{\mathrm{H}}$ is the Heisenberg time (at the given value of A), d the dimension of the Hilbert space and $\sigma \left(A\right)$ the width of the Gaussian density of states given in (10). The quantity $\mathsf{\Delta }t$ is the shortest physical time scale of the system (inverse of the largest energy scale) and obviously for $t\ll \mathsf{\Delta }t$ the unitary evolution operator is close to the unit matrix multiplied by a uniform phase factor: $exp(-iHt)\approx exp(-i{E}_{\mathrm{mean}}t)\mathbf{1}$ since the eigenvalues $E_m$ of H satisfy $|E_m-E_{\mathrm{mean}}|\lesssim \sigma \left(A\right)$. We expect that any signification deviation of $\left|\psi \right(t)>$ with respect to the initial condition $\left|\psi \right(0)>$ happens at $t\ge \mathsf{\Delta }t$ (or later in case of very weak interaction). Furthermore, by analyzing the time evolution in terms of the ratio $t/\mathsf{\Delta }t$ the results do not depend on the global energy scale of the spectral width. The longest time scale is the Heisenberg time $t_{\mathrm{H}}\approx {10}^4\mathsf{\Delta }t$ (since $d=11440$ for $M=16$, $L=7$). Later, we also discuss intermediate time scales such as the inverse decay rate obtained from the Fermi golden rule.

To show graphically the time evolution, we compute the time dependent occupation numbers $n_k\left(t\right)=<\psi \left(t\right)\left|c_k^{†}{c}_k^{\phantom{†}}\right|\psi \left(t\right)>$ and present them in a color density plot in the plane $(k,t/\mathsf{\Delta }t)$. In addition, at the time value used last, we compute the effective total one-particle energy $E_{1\mathrm{p}}$ using the relation (11) (note that $E_{1\mathrm{p}}$ is not conserved with respect to the time evolution except for very weak interaction) and use this value to determine from (12) the inverse temperature $\beta$, chemical potential $\mu$ and the thermalized Fermi–Dirac filling factor $n\left({\epsilon }_k\right)$ at each k value for the orbital index. These values of ideally thermalized occupation numbers will be shown in an additional vertical bar ( Note that this additional bar is not related to the usual color bar that provides the translation of colors to $n_k$ values). The latter is shown in Figure 6 and applies also to all subsequent figures with color density plots for $n_k$ values right behind the data for the last time values separated by a vertical white line. This presentation allows for an easy verification if the occupation numbers at the last time values are indeed thermalized or not.

In Figure 7, we show the time evolution for the initial state $|{\varphi }_1>$ and the two Åberg parameter values $A=1$ and $A=3.5$ using a linear time scale with integer multiples of $\mathsf{\Delta }t$ and for $t\le 2000\mathsf{\Delta }t\approx t_{\mathrm{H}}/6$. At $A=1$, the occupation number $n_{12}$ (of the excited particle) shows, at the beginning, a periodic structure, with an approximate period $400\mathsf{\Delta }t$ for $t<1000\mathsf{\Delta }t$, and a modest decay for $t>1000\mathsf{\Delta }t$. At the same time, the first orbitals above the Fermi sea are slightly excited. At final $t=2000\mathsf{\Delta }t$, the state is clearly not thermalized. For $A=3.5$, we see a very rapid partial decay of $n_6$ and $n_{12}$ together with an increase of $n_7$. Furthermore, for $n_k$ with $8\le k\le 11$, there are later and more modest excitations with a periodic time structure. Here, the final state at $t=2000\mathsf{\Delta }t$ is also not thermalized, but it is closer to thermalization as for the case $A=1$.

The linear time scale used in Figure 7 is not very convenient since it cannot capture a rapid decay/increase of $n_k$ well at small times and its maximal time value is also significantly limited below the Heisenberg time. Therefore, we use in Figure 8 and Figure 9 a logarithmic time scale with $0.1\mathsf{\Delta }t\le t\le {10}^6\mathsf{\Delta }t\approx {10}^2{t}_{\mathrm{H}}$. Note that, in these figures, the different $n_k$ values for each cell are not time averaged but represent the precise values for certain, exponentially increasing, discrete time values (see caption of Figure 8 for the precise values). Therefore, in case of periodic oscillations of $n_k$, there will be, for larger time values, a quasi random selection of different time positions with respect to the period.

In Figure 8, the time evolution for the initial states $|{\varphi }_1>$ and $|{\varphi }_2>$ is shown for the Åberg values $A=1,3.5,10$. For $|{\varphi }_1>$ at $A=1$ and $A=3.5$, the observations of Figure 7 are confirmed with the further information that the absence of thermalization in these cases is also valid for time scales larger than $2000\mathsf{\Delta }t$ up to ${10}^6\mathsf{\Delta }t$ and for $A=3.5$ the initial decay of $n_6$ and $n_{12}$ happens at $t\approx 10\mathsf{\Delta }t$. For $|{\varphi }_1>$ at $A=10$, the decay starts at $t\approx 3\mathsf{\Delta }t$ and an approximate thermalization happens at $t>40\mathsf{\Delta }t$. However, here there is still some time periodic structure and it would be necessary to do some time average to have perfect thermalization. For $|{\varphi }_2>$ at $A=1$, the decay of excited orbitals 12 and 14 starts at $t\approx 100\mathsf{\Delta }t$ and saturates at $t\approx 1000\mathsf{\Delta }t$ at which time also orbitals 6 and 7 are excited. After this, there are very small excitations of orbitals 8, 9, 10 and maybe 13, 15. There is also some very modest decay of the Fermi sea orbitals 2, 4 and 5 at $t>1000\mathsf{\Delta }t$. The final state at $t={10}^6\mathsf{\Delta }t$ is not thermalized even though some orbitals have $n_k$ values close to thermalization. For $|{\varphi }_2>$ at $A=3.5$, the decay of excited orbitals 12 and 14 starts at $t\approx 10\mathsf{\Delta }t$ and, for $t>300\mathsf{\Delta }t$, there is thermalization (but requiring some time average as for $|{\varphi }_1>$ at $A=10$). Interestingly, at intermediate times $10\mathsf{\Delta }t<t<100\mathsf{\Delta }t$, the high orbitals 13 and 16 are temporarily slightly excited and decay afterwards rather quickly to their thermalized values. For $|{\varphi }_2>$ at $A=10$, the decay of excited orbitals 12 and 14 starts even at $t\approx 3\mathsf{\Delta }t$ and thermalization seems to set in at $t>30\mathsf{\Delta }t$.

Figure 9 is similar to Figure 8 except for the initial states $|{\varphi }_3>$ and $|{\varphi }_4>$ which have occupation numbers $n_k\in \{0,1\}$ obtained by rounding the $n_k$ values of the two eigenstates visible in the two top panels of Figure 2. Here, the initial decay of excited orbitals starts roughly at $t\approx 300\mathsf{\Delta }t$ ($t\approx (10-20)\mathsf{\Delta }t$ or $t\approx (2-3)\mathsf{\Delta }t$) for $A=1$ ($A=3.5$ or $A=10$, respectively). There is no thermalization for both states at $A=1$ (but some $n_k$ values are close to thermalized values), approximate thermalization for $A=3.5$ and $|{\varphi }_3>$ and good thermalization for $A=3.5$ and $|{\varphi }_4>$ as well as $A=10$ (both states).

Using the time dependent values $n_k\left(t\right)$, one can immediately determine the corresponding entropy $S\left(t\right)$ using (11). At $t=0$, we have obviously $S\left(0\right)=0$, since, for all four initially considered states, we have perfect occupation number values of either $n_k=0$ or $n_k=1$. Naturally, one would expect that the entropy increases with a certain rate and saturates then at some maximal value which may correspond (or be lower) to the thermalized entropy $S_{\mathrm{th}}\left(E_{1\mathrm{p}}\right)$ (with $E_{1\mathrm{p}}$ determined for the state $\left|\psi \right(t)>$ at large times) depending if there is presence (or absence) of thermalization according to the different cases visible in Figure 8 and Figure 9. However, in the absence of thermalization, we see that there may also exist periodic oscillations with a finite amplitude at very long time scales.

In Figure 10, we show the time dependent entropy $S\left(t\right)$ for the two initial states $|{\varphi }_1>$, $|{\varphi }_4>$ and the three values $A=1$, $A=3.5$ and $A=10$ of the Åberg parameter. For $A=10$, there is indeed a rather rapid saturation of the entropy of both states at a maximal value that is indeed close to the thermalized entropy $S_{\mathrm{th}}\left(E_{1\mathrm{p}}\right)$. We note that $E_{1\mathrm{p}}$ is not conserved at strong interactions and that its initial value $E_{1\mathrm{p}}\approx 30$ ($E_{1\mathrm{p}}\approx 35$) at $t=0$ evolves to $E_{1\mathrm{p}}\approx 33.5$ ($E_{1\mathrm{p}}\approx 37$) at large times for $|{\varphi }_1>$ ($|{\varphi }_4>$) corresponding roughly to $S\approx S_{\mathrm{th}}\left(E_{1\mathrm{p}}\right)\approx 9.2$ ($10.8$) visible as thin blue horizontal lines in Figure 10. For $A=3.5$ (or $A=1$), the thermalized entropy values, visible as thin green (red) lines, are lower as compared to the case $A=10$ due to different final $E_{1\mathrm{p}}$ values. For $A=3.5$ and $|{\varphi }_4>$, there is also saturation of S to its thermalized value. For $A=3.5$ and $|{\varphi }_1>$, there seems to be an approximate saturation at a quite low value $S\approx 6$ but with periodic fluctuations in the range $6\pm 0.3$. For $A=1$ and $|{\varphi }_4>$, there is a quite late and approximate saturation with some fluctuations that are visible for $t>{10}^4\mathsf{\Delta }t$ and with $S\approx 10\pm 0.2$. For $A=1$ and $|{\varphi }_1>$, there is a late periodic regime for $t>{10}^3\mathsf{\Delta }$ with a quite large amplitude $S\approx 3\pm 1$ and with $S_{\max }\approx 4$ significantly below the thermalized entropy $S_{\mathrm{th}}\left(E_{1\mathrm{p}}\right)\approx 5.5$. The panels using a normal (instead of logarithmic) time scale with $t\le 200\mathsf{\Delta }t$ miss completely the long time limits for $A=1$ and might incorrectly suggest that there is an early saturation at quite low values of S.

The periodic (or quasi-periodic) time dependence of $n_k\left(t\right)$ or $S\left(t\right)$, for the cases with lower values of A and/or an initial state with lower energy, indicates that, for such states, only a small number ($2,3,\dots$) of exact eigenstates of H contribute mostly in the expansion of $\left|\psi \right(t)>$ in terms of these eigenstates.

Figure 10 also shows that the initial increase of $S\left(t\right)$ is rather comparable between the two states for identical values of A even though the long time limit might be very different. Furthermore, a closer inspection of the data indicates that typically $S\left(t\right)$ is close to a quadratic behavior for $t\lesssim \mathsf{\Delta }t$ but which immediately becomes linear for $t\gtrsim \mathsf{\Delta }t$ similarly as the transition probabilities between states in the context of time dependent perturbation theory. To study the approximate slope in the linear regime, we define (For practical reasons, we decide to incorporate the quantum time step $\mathsf{\Delta }t$ in the definition of ${\mathrm{\Gamma }}_c$, i.e., ${\mathrm{\Gamma }}_c$ is defined as the ratio of the initial slope $S^{\prime }\left(t\right)$ over the global spectral bandwidth $\sim \sigma \left(A\right)\sim 1/\mathsf{\Delta }t$) the quantity ${\mathrm{\Gamma }}_c=dS\left(t\right)/d(t/\mathsf{\Delta }t)=\mathsf{\Delta }t{S}^{\prime }\left(t\right)$ for $t=\xi \mathsf{\Delta }t$, where $\xi \gtrsim 1$ is a numerical constant of order one. To determine ${\mathrm{\Gamma }}_c$ practically, we perform first the fit $S\left(t\right)={\overline{S}}_{\infty }(1-exp[-{\overline{\gamma }}_1(t/\mathsf{\Delta }t)])$ for $0\le t/\mathsf{\Delta }t\le 100$ and use the exponential decay rate ${\overline{\gamma }}_1$ to perform a refined fit $S\left(t\right)=S_{\infty }(1-exp[-{\gamma }_1(t/\mathsf{\Delta }t)-{\gamma }_2{(t/\mathsf{\Delta }t)}^2])$ for the interval $0\le t/\mathsf{\Delta }t\le 5/{\overline{\gamma }}_1$. From this, we determine ${\mathrm{\Gamma }}_c=S_{\infty }{\gamma }_1$, which is rather close to ${\overline{S}}_{\infty }{\overline{\gamma }}_1$ for $A\le 2$ but not for larger values of A where the decay time is reduced and not sufficiently large in comparison to the initial quadratic regime. Therefore, the quadratic term in the exponential is indeed necessary to obtain a reasonable fit quality. This procedure corresponds to an effective average of the value of $\xi$ between 1 and roughly $1/{\gamma }_1$, which is indeed useful to smear out some oscillations in the initial increase of $S\left(t\right)$ for smaller values $A\le 1$.

We note that, for many-body quantum systems, the exponential growth of entropy with time had been also discussed and numerically illustrated in [40] (see also related publications in References 25 and 26 there). Recently, such an exponential growth of entropy has been discussed in [62,70].

Figure 11 shows the dependence of these values of ${\mathrm{\Gamma }}_c$ on the parameter A for our four initial states. At first sight, one observes a behavior ${\mathrm{\Gamma }}_c\propto A^2$ for $A\lesssim 2$ and a saturation for larger values of A. However, a more careful analysis shows that there are modest but clearly visible deviations with respect to the quadratic behavior in A (power law fits ${\mathrm{\Gamma }}_c\propto A^p$ for $A\le 2$ provide exponents close to $p\approx 1.75-1.85$) and it turns out that these deviations correspond to a logarithmic correction: ${\mathrm{\Gamma }}_c=f\left(A\right)=(C_1-C_2\ln \left[g\left(A\right)\right])g\left(A\right)$ with $g\left(A\right)=A^2$ (for fits with $A\le 1$) or with $g\left(A\right)=A^2/(1+C_3{A}^2)$ (for fits with all A values).

To understand this behavior, we write for sufficiently small times $n_k\left(t\right)\approx 1-\delta n_k\left(t\right)$ (if $n_k\left(0\right)=1$) or $n_k\left(t\right)\approx \delta n_k\left(t\right)$ (if $n_k\left(0\right)=0$), where $\delta n_k\left(t\right)$ is the small modification of $n_k\left(t\right)$. Time dependent perturbation theory suggests that $\delta n_k\left(t\right)\sim {(t/\mathsf{\Delta }t)}^2$ for $t\lesssim \mathsf{\Delta }t$ and $\delta n_k\left(t\right)\approx a_k{A}^{2}t/\mathsf{\Delta }t$ for $t\gtrsim \mathsf{\Delta }t$ such that still $\delta n_k\left(t\right)\ll 1$ with coefficients $a_k$ dependent on k (and also on M, L) and satisfying a linear relation to ensure the conservation of particle number. Using (11) and neglecting corrections of order $\delta n_k^2$, we obtain: $S\approx -{\sum }_k\left(\delta n_k\ln \delta n_k-\delta n_k\right)$ and ${\mathrm{\Gamma }}_c=\mathsf{\Delta }t{S}^{\prime }\left(t\right)\approx -\mathsf{\Delta }t{\sum }_k\delta n_k^{\prime }\left(t\right)\ln \delta n_k\left(t\right)$ with $t=\xi \mathsf{\Delta }t$. Since $\delta n_k^{\prime }\left(t\right)\approx a_k{A}^2/\mathsf{\Delta }t$, we find indeed the behavior:

$${\mathrm{\Gamma }}_c=[C_1-C_2\ln \left(A^2\right)]A^2,C_1=-\sum _k{a}_k\langle \ln \left(a_k\xi \right)\rangle ,C_2=\sum _k{a}_k,$$

(14)

where $\langle \cdots \rangle$ indicates an average over some modest values of $\xi \gtrsim 1$. The precise values of $a_k$ may depend rather strongly on the orbital index k and the initial state (see also Figure 8 and Figure 9), but the coefficients $C_1$, $C_2$ depend only slightly on the initial state (see

Table 2. Values of the fit parameters $C_1,\dots ,C_5$ for the initial states $|{\varphi }_1>,\dots ,|{\varphi }_4>$ used for the analytical fits of ${\mathrm{\Gamma }}_c$ (${\mathrm{\Gamma }}_F$) in Figure 11 (see also Figure 12 and Figure 13).

Initial State	${\mathit{C}}_1$	${\mathit{C}}_2$	${\mathit{C}}_3$	${\mathit{C}}_4$	${\mathit{C}}_5$
$|{\varphi }_1>=|0000100000111111>$	$0.107\pm 0.009$	$0.0081\pm 0.0023$	$0.092\pm 0.017$	$0.0048\pm 0.0001$	$0.0054\pm 0.0003$
$|{\varphi }_2>=|0010100000011111>$	$0.100\pm 0.003$	$0.0103\pm 0.0011$	$0.069\pm 0.007$	$0.0068\pm 0.0001$	$0.0099\pm 0.0003$
$|{\varphi }_3>=|0000011000110111>$	$0.110\pm 0.005$	$0.0130\pm 0.0016$	$0.080\pm 0.012$	$0.0076\pm 0.0001$	$0.0098\pm 0.0003$
$|{\varphi }_4>=|1000100011001011>$	$0.103\pm 0.003$	$0.0210\pm 0.0007$	$0.018\pm 0.005$	$0.0094\pm 0.0001$	$0.0140\pm 0.0002$

). Furthermore, by replacing $A^2\to g\left(A\right)=A^2/(1+C_3{A}^2)$ to allow for a saturation at large A and with a further fit parameter $C_3$, it is possible to describe the numerical data by the more general fit ${\mathrm{\Gamma }}_c=f\left(a\right)$ for the full range of A values.

3.7. Survival Probability and Fermi’s Golden Rule

The knowledge of the time dependent states $\left|\psi \right(t)>$ allows us also to compute the decay function $p_{\mathrm{dec}}\left(t\right)={|<\psi \left(0\right)|\psi \left(t\right)>|}^2$ which represents the survival probability of the initial non-interacting eigenstate due to the influence of interactions. Again, for the very short time window $t\lesssim \mathsf{\Delta }t$, we expect a quadratic decay: $1-p_{\mathrm{dec}}\left(t\right)\approx \langle {(H-E_{\mathrm{mean}})}^2\rangle t^2\approx \mathrm{const}.{(t/\mathsf{\Delta }t)}^2$ with $\langle \cdots \rangle$ being the quantum expectation value with respect to $\left|\psi \right(0)>$ and a numerical constant $\lesssim 1$ since $1/\mathsf{\Delta }t$ represents roughly the spectral width of H. For $t\gtrsim \mathsf{\Delta }t$, but such that $1-p_{\mathrm{dec}}\left(t\right)\ll 1$, we have, according to Fermi’s golden rule: $1-p_{\mathrm{dec}}\left(t\right)={\mathrm{\Gamma }}_F(t/\mathsf{\Delta }t),$ where ${\mathrm{\Gamma }}_F$ is the decay rate (Again, for practical reasons and similarly to ${\mathrm{\Gamma }}_c$, we incorporate in the definition of ${\mathrm{\Gamma }}_F$ the time scale $\mathsf{\Delta }t$, i.e., ${\mathrm{\Gamma }}_F=\mathsf{\Delta }t\times$ usual decay rate found in the literature and meaning that ${\mathrm{\Gamma }}_F$ is defined as the ratio of the usual decay rate over the global spectral bandwidth) of the state.

To determine ${\mathrm{\Gamma }}_F$ numerically, we apply the fit: $p_{\mathrm{dec}}\left(t\right)=Cexp(-{\mathrm{\Gamma }}_{\mathrm{F}}t/\mathsf{\Delta }t)$ in two steps. First, we use the interval $1\le t/\mathsf{\Delta }t\le 50$ and, if $5/{\mathrm{\Gamma }}_F<50$, corresponding to a rapid decay (which happens for larger values of A), we repeat the fit for the reduced interval $1\le t/\mathsf{\Delta }t\le 5/{\mathrm{\Gamma }}_F$. The choice of the Amplitude $C\ne 1$ and the condition $t\ge \mathsf{\Delta }t$ for the fit range allow for taking into account the effects due to the small initial window of quadratic decay. In Figure 12, we show two examples for the initial state $|{\varphi }_1>$ and the Åberg values $A=3.5$ and $A=10$. In both cases, the shown maximal time value $t_{\max }=50\mathsf{\Delta }t$ (if $A=3.5$) or $t_{\max }\approx 13.5$ (if $A=10$) defines the maximal time value for the fit range. For $A=10$, the fit nicely captures the decay for $1\le t/\mathsf{\Delta }t\le 6$, while, for $A=3.5$, there are also some oscillations in the decay function for which the fit procedure is equivalent to some suitable average in the range $1\le t/\mathsf{\Delta }t\le 30$. For very small values of A, the fit procedure also works correctly since it captures only the initial decay that is important if $p_{\mathrm{dec}}\left(t\right)$ does not decay completely at large times and which typically happens in the perturbative regime $A\lesssim 1$.

Figure 13 shows the dependence of ${\mathrm{\Gamma }}_F$ on A for the usual four initial states together with the fit ${\mathrm{\Gamma }}_F=f\left(A\right)=C_4{A}^2/(1+C_5{A}^2)$ for the data with initial state $|{\varphi }_1>$. The values of the parameters $C_4$, $C_5$ for this and the other initial states are given in Table 2. Here, the initial quadratic dependence ${\mathrm{\Gamma }}_F\propto A^2$ is highly accurate (with no logarithmic correction). Similarly to ${\mathrm{\Gamma }}_c$, there is only a slight dependence of the values of ${\mathrm{\Gamma }}_F$ and the fit values on the choice of initial state.

Theoretically, we expect according to Fermi’s Golden rule that: ${\mathrm{\Gamma }}_F\approx \left(\mathsf{\Delta }t\right)2\pi V_{\mathrm{Fock}}^2{\rho }_c\left(E\right),$ where $V_{\mathrm{Fock}}^2={Tr}_{\mathrm{Fock}}\left(V^2\right)/\left(Kd\right)={\sigma }_0^2\alpha A^2/K$ according to the discussion below (10) and ${\rho }_c\left(E\right)$ is the effective two-body density of states for states directly coupled by the interaction such that ${\rho }_c\left(E_{\mathrm{mean}}\right)=1/{\mathsf{\Delta }}_c$ (see discussion below (7)). We note that $V_{\mathrm{Fock}}$ is the typical interaction matrix element in Fock space which is slightly larger than $V_{\mathrm{mean}}$ (see the theoretical discussion above for the computation of the coefficient $\alpha$ used in (10) and Appendix A of [46]). The factor $\mathsf{\Delta }t$ is due to our particular definition (Again, for practical reasons and similarly to ${\mathrm{\Gamma }}_c$, we incorporate in the definition of ${\mathrm{\Gamma }}_F$ the time scale $\mathsf{\Delta }t$, i.e., ${\mathrm{\Gamma }}_F=\mathsf{\Delta }t\times$ usual decay rate found in the literature and meaning that ${\mathrm{\Gamma }}_F$ is defined as the ratio of the usual decay rate over the global spectral bandwidth) of decay rates. The expression of ${\mathrm{\Gamma }}_F$ is actually also valid for larger values of A provided we use the density of states ${\rho }_c\left(E\right)$ in the presence of interactions that provides an additional factor $1/\sqrt{1+\alpha A^2}$ according to (10). Therefore, at the band center, we have: $2\pi \mathsf{\Delta }t{\rho }_c=K/\left[{\sigma }_0\left(L\right){\sigma }_0(L=2)(1+\alpha A^2)\right],$ which gives, together with (8):

$${\mathrm{\Gamma }}_F=\frac{{\sigma }_0\left(L\right)}{{\sigma }_0(L=2)}\left(\frac{\alpha A^2}{1+\alpha A^2}\right)=\sqrt{\frac{L(M-L)}{2(M-2)}}\left(\frac{\alpha A^2}{1+\alpha A^2}\right).$$

(15)

For $M=16$ and $L=7$, the square root factor is 1.5 and we have to compare $1.5\alpha \approx 0.0132$ with the values of $C_4$ in Table 1, which are somewhat smaller, probably due to a reduction factor for the energy dependent density of states since the energies of the initial states have a certain distance to the band center. Furthermore, according to (15), we have to compare $C_5$ with $\alpha \approx 0.00877$ which is not perfect but gives the correct order of magnitude. For both parameters, the numerical matching is quite satisfactory taking into account the very simple argument using the same typical value of the interaction matrix elements for all cases of initial states.

Finally, we mention that, for the three Åberg parameter values $A=1$, $A=3.5$, $A=10$ used in Figure 8 and Figure 9, we have typical decay times in units of $\mathsf{\Delta }t$ being $1/{\mathrm{\Gamma }}_F\approx 300$, 30, 3, respectively (with some modest fluctuations depending on initial states). These values match quite well the observed time values at which the initially occupied orbitals start to decay (see the above discussion of Figure 8 and Figure 9).

3.8. Time Evolution of Density Matrix and Spatial Density

We now turn to the effects of the many-body time evolution in position space (see for example Figure 1). For this, we compute the spatial density

$$\rho (x,y,t)=<\psi \left(t\right)|{\mathrm{\Psi }}^{†}(x,y)\mathrm{\Psi }(x,y)|\psi \left(t\right)>,{\mathrm{\Psi }}^{(†)}(x,y)=\sum _k{\phi }_k^{(*)}(x,y)c_k^{(†)},$$

(16)

where ${\phi }_k(x,y)$ is the one-particle eigenstate of orbital k, with some examples shown in Figure 1. Here, ${\mathrm{\Psi }}^{(†)}(x,y)$ denotes the usual fermion field operators (in the case of continuous x, y variables) or standard fermion operators for discrete position basis states (when using a discrete grid for x and y positions as we did for the numerical solution of the non-interacting Sinai oscillator model in Section 2). The sum over orbital index k in (16) requires in principle a sum over a full complete basis set of orbitals with infinite number (case of continuous x, y values) or a very large number (case of discrete x-y grid) significantly larger than the very modest number of orbitals M we used for the numerical solution of the many-body Sinai oscillator.

However, we can simply state that, in our model, by construction, all orbitals with $k>M$ are never occupied such that, in the expectation value for $\rho (x,y,t),$ only the values $k\le M$ are necessary. Taking this into account together with the fact that the one-particle eigenstates are real valued, we obtain the more explicit expression:

$$\rho (x,y;t)=\sum _{k,l=1}^M{\phi }_k(x,y){\phi }_l(x,y)n_{kl}\left(t\right),n_{kl}\left(t\right)=<\psi \left(t\right)\left|c_k^{†}{c}_l^{\phantom{†}}\right|\psi \left(t\right)>,$$

(17)

where $n_{kl}\left(t\right)$ is the density matrix in orbital representation generalizing the occupation numbers $n_k\left(t\right)$ which are its diagonal elements. Due to the complex phases of $\left|\psi \right(t)>$ (when expanded in the usual basis of non-interacting many-body states), the density matrix is complex valued but hermitian: $n_{kl}^{*}\left(t\right)=n_{lk}\left(t\right)$. Therefore, its anti-symmetric imaginary part does not contribute in $\rho (x,y;t)$. We have numerically evaluated (17) and we present in Figure 14 color plots of the density matrix and the spatial density $\rho (x,y;t)$ for $A=3.5$, the initial state $|\psi \left(0\right)>=|{\varphi }_2>$ and four time values $t/\mathsf{\Delta }t=0.1,30,100,1000$. Since the density $\rho (x,y;t)$ does not provide a lot of spatial structure, we also show in Figure 14 the density difference with respect to the initial condition $\mathsf{\Delta }\rho (x,y;t)=\rho (x,y;t)-\rho (x,y;0)$ which reveals more of its structure (figures and videos for the time evolution of this and other cases are available for download at the web page [71]).

At the time $t/\mathsf{\Delta }t=0.1$, density matrix and spatial density are essentially identical to the initial condition at $t=0$. For $\mathsf{\Delta }\rho$, we see a non-trivial structure since there is a small difference with the initial condition and the color plot simply amplifies small maximal amplitudes to maximal color values (red/yellow for strongest positive/negative values even if the latter are small in an absolute scale). The density matrix is diagonal and its diagonal values are either 1 (for initially occupied orbitals) or 0 (for initially empty orbitals) and the spatial density simply gives the sum of densities due to the occupied eigenstates.

At $t/\mathsf{\Delta }t=30$, we see a non-trivial structure in the density matrix with a lot of non-vanishing values in certain off-diagonal elements. Furthermore, the orbitals 13 and 16 are also slightly excited (see also discussion of Figure 8) and there is a significant change of the spatial density.

Later, at $t/\mathsf{\Delta }t=100$, the number/values of off-diagonal elements in the density matrix is somewhat reduced, but they are still visible. Especially between orbitals 12 and 13 as well as 14 and 15, there is a rather strong coupling. Orbital 13 is now more strongly excited than the initially excited orbital 12. In addition, orbitals 14 and 15 are quite strong. The spatial density has become smoother and the structure of $\mathsf{\Delta }\rho$ is roughly close to the case at $t/\mathsf{\Delta }t=30$ but with some significant differences.

Finally, at $t/\mathsf{\Delta }t=1000$, the density matrix seems be diagonal with values close to the thermalized values. There is a further increase of the density smoothness and $\mathsf{\Delta }\rho$ has a similar but different structure as for $t/\mathsf{\Delta }t=100$ or $t/\mathsf{\Delta }t=30$.

Apparently, at intermediate times $20\le t/\mathsf{\Delta }t\le 100$, there are some quantum correlations between certain orbitals, visible as off-diagonal elements in the density matrix which disappear for later times. This kind of decoherence is similar to the exponential decay observed in [46] for the off-diagonal element of the $2\times 2$ density matrix for a qubit coupled to a chaotic quantum dot or the SYK black hole. However, to study this kind of decoherence more carefully in the context here, it would be necessary to use as initial state a non-trivial linear combination of two non-interacting eigenstates and not to rely on the creation of modest off-diagonal elements for intermediate time scales as we see here.

The spatial density is globally rather smooth and typically quite well given by the “classical” relation $\rho (x,y;t)\approx {\sum }_k{\phi }_k^2(x,y)n_k\left(t\right)$ in terms of the time dependent occupation numbers. Only for intermediate time scales with more visible quantum coherence (more off-diagonal elements $n_{kl}\left(t\right)\ne 0$), this relation is less accurate. However, at $A=3.5$, the density still exhibits small but regular fluctuations in its detail structure as can be seen in the structure of $\mathsf{\Delta }\rho$ for later time scales. A closer inspection of the data (for time values not shown in Figure 14) also shows that, even at long time scales, there are significant fluctuations of $\rho$ when t is slightly changed by a few multiples of $\mathsf{\Delta }t$.

In Figure 14, we also show for comparison the theoretical thermalized quantities where, in (17), the density matrix is replaced by a diagonal matrix with entries being the thermalized occupations numbers $n_{kk}=n\left({\epsilon }_k\right)$ at energy $E=32.9$, which is the typical total one-particle average energy of $\left|\psi \right(t)>$ for the long time limit $t/\mathsf{\Delta }t\ge 1000$ showing that, at $t/\mathsf{\Delta }t=1000$, the state is very close to thermalization but still with small significant differences (see also discussion of Figure 7 for this case).

We may also generalize the spatial density (16) to a spatial density correlator which we define as:

$${\rho }_{\mathrm{corr}}(x,y;x_0,y_0;t)=<\psi \left(t\right)\left|{\mathrm{\Psi }}^{†}(x,y)\mathrm{\Psi }(x_0,y_0)\right|\psi \left(t\right)>=\sum _{k,l=1}^M{\phi }_k(x,y){\phi }_l(x_0,y_0)n_{kl}\left(t\right)$$

(18)

depending on initial $(x_0,y_0)$ and final position $(x,y)$. As an illustration, we choose the fixed value $(x_0,y_0)=(1.22,0.15)$ which is very close to the maximal position (center of the red area) of the one-particle ground state ${\phi }_1(x,y)$ visible in panel (a) of Figure 1. The spatial density correlator is potentially complex with a non-vanishing imaginary part in case of non-vanishing off-diagonal matrix elements of $n_{kl}\left(t\right)\ne 0$ for $k\ne l$. In Figure 15, we present density plots of absolute value, real and imaginary part of ${\rho }_{\mathrm{corr}}(x,y;x_0,y_0;t)$ in $(x,y)$ plane and with the given value $(x_0,y_0)=(1.22,0.15)$ for the same parameters of Figure 14 (concerning initial state, Åberg parameter, time values and also thermalized case). However, for the thermalized case, the density matrix is diagonal by construction and the imaginary part of ${\rho }_{\mathrm{corr},\mathrm{th}}(x,y;x_0,y_0)$ vanishes (giving a blue panel due to zero values).

There are significant time dependent fluctuations of ${\rho }_{\mathrm{corr}}(x,y;x_0,y_0;t)$ for all time scales with real part and absolute value being dominated by rather strong maximal values for positions close to the initial position. However, the imaginary part (which vanishes at $t=0$ and is typically smaller than the values in maximum domain of a real part) shows a more interesting structure since the color plot amplifies small amplitudes (in absolute scale). Apart from this, the absolute and real part values for positions outside the maximum domain (far away from the initial position) seem to decay for long timescales, which is also confirmed by the thermalized case. Even though the case for $t/\mathsf{\Delta }t=1000$ seems to be rather close to the thermalized case (for absolute value and real part), there are still differences that are more significant here as in Figure 14.

Globally, the obtained results show that the dynamical thermalization takes place well, leading to the usual Fermi–Dirac thermal distribution when the Åberg criterion is satisfied and interactions are sufficiently strong to drive the system into the thermal state.

4. Estimates for Cold Atom Experiments

We discuss here typical parameters for cold atoms in a trap following [72]. Thus, for sodium atoms, we have $\omega \approx {\omega }_{x,y,z}\approx 2\pi 10Hz$ with $a_0=\sqrt{\hslash /\left(m\omega \right)}=6.5\mu m$ and oscillator level spacing $E_u=\hslash \omega \approx 0.5nK$ (nanoKelvin). The typical scattering length is $a_s\approx 3nm$ being small compared to $a_0$. The atomic density is ${\rho }_0=1/{\left(a_0\right)}^3\approx 4\times {10}^{9}c{m}^{-3}$. Since $a_s\ll a_0$, the two-body interaction is of $\delta$-function type with $v\left({\mathbf{r}}_{\mathbf{1}}-{\mathbf{r}}_{\mathbf{2}}\right)=(4\pi {\hslash }^2{a}_s/m)\delta \left({\mathbf{r}}_{\mathbf{1}}-{\mathbf{r}}_{\mathbf{2}}\right)$ [66,72]. In our numerical simulations, $\delta \left(\mathbf{r}\right)$ is replaced by a function $H(r_c-\left|\mathbf{r}\right|)/\left(C{2}^d{r}_c^d\right)$ with a small $r_c$, volume C of the unit sphere in d dimensions and with the Heaviside function $H\left(x\right)=1$ (or 0) for $x\ge 0$ and $H=0$ for $x<0$. Hence, our parameter U introduced in Section 3 corresponds to $U=\left(C{2}^d{r}_c^d\right)(4\pi {\hslash }^2{a}_s/m)$.

Below, we present the estimates for the dynamical thermalization border for excitations of fermionic atoms in a vicinity of their Fermi energy in a 3D Sinai oscillator following the lines of Equation (1). In such a case, the two-body interaction energy scale between atoms is $U_s=4\pi {\hslash }^2{\rho }_0{a}_s/m=4\pi (a_s/a_0)\hslash \omega$ [66,72] so that $U_s/\hslash \omega \sim 6\times {10}^{-3}$. Compared to sodium, the mass of Li atoms is approximately three times smaller so that, for the same $\omega$, we have $a_0\approx 10\mathsf{\mu }\mathrm{m}$ and $U_s/\hslash \omega \approx 4\times {10}^{-3}$. We think that the scattering length can be significantly increased via the Feshbach resonance, allowing for reaching effective interaction values $U_s/\hslash \omega \sim 1$ being similar to the value $A\sim 3$ used in our numerical studies with the onset of dynamical thermalization.

Usually, a 3D trap with fermionic atoms can capture about $N_a\sim {10}^5$ atoms with $\omega \approx {\omega }_x\sim {\omega }_y\sim {\omega }_z\sim 2\pi 10Hz$. Following the result (1), it is interesting to determine the DTC border dependence on $N_a\gg 1$ for a Sinai oscillator with $r_d\sim 1\mu m\sim a_0/5$. We assume that, similar to a 2D case, the scattering on an elastic ball in the trap center leads to quantum chaos and chaotic eigenstates with $\ell \le N_a$ components (e.g., in the basis of oscillator eigenfunctions). The Fermi energy of the trap is then $E_F=\hslash {\left(N_a{\omega }_x{\omega }_y{\omega }_z\right)}^{1/3}\approx \hslash \omega {N_a}^{1/3}$ [52,53]. Assuming that all these components have random amplitudes of a typical size $1/\sqrt{\ell }$, we then obtain an estimate for a typical matrix element of two-body interaction between one-particle eigenstates

$$U_2\approx {\alpha }_s\hslash \omega /{\ell }^{3/2},{\alpha }_s=4\pi (a_s/a_0),a_0=\sqrt{\hslash /m\omega }.$$

(19)

The derivation of this estimate is very similar to the case of two interacting particles in a disordered potential with localized eigenstates [73]. At the same time, in the vicinity of the Fermi energy $E_F$, we have the one-particle level spacing ${\mathsf{\Delta }}_1=d{E}_F/d{N}_a\approx \hslash \omega /\left(3{N_a}^{2/3}\right)$. Hence, the effective conductance appears in in (1) is $g={\mathsf{\Delta }}_1/U_2\approx {\ell }^{3/2}/\left(3{\alpha }_s{N_a}^{2/3}\right)$. Thus, from (1), we obtain the dynamical thermalization border for excitation energy $\delta E$ in a 3D Sinai oscillator trap with $N_a$ fermionic atoms:

$$\delta E>\delta E_{ch}\approx {\mathsf{\Delta }}_1{g}^{2/3}\approx 2\ell {\mathsf{\Delta }}_1/({{\alpha }_s}^{2/3}{N_a}^{4/9})\sim {N_a}^{5/9}{\mathsf{\Delta }}_1/{{\alpha }_s}^{2/3}\sim \hslash \omega /({{\alpha }_s}^{2/3}{N_a}^{1/9})\sim E_F/({{\alpha }_s}^{2/3}{N_a}^{4/9}).$$

(20)

It is assumed that $\delta E\ll E_F$. Here, the last three relations are written in an assumption that $\ell \sim N_a$. Thus, at large $N_a$ values and not too small ${\alpha }_s$, the critical energy border $\delta E_{ch}$ for dynamical thermalization is rather small compared to $E_F$. However, still $\delta E_{ch}\gg {\mathsf{\Delta }}_1$. Here, we used the maximal value for the number of components $\ell \sim N_a$. It is possible that, in reality, ℓ can be significantly smaller than $N_a$. However, the determination of the dependence $\ell \left(N_a\right)$ requires separate studies taking into account the properties of chaotic eigenstates and their spreading over the energy surface. This spreading can have rather nontrivial properties (see, e.g., [23]). This is confirmed by the results presented in Appendix A for the 2D case of Sinai oscillator showing the numerically obtained dependence of two-body matrix elements on energy for transitions in a vicinity of Fermi energy $E_F$ (see Figure A1 there).

5. Conclusions

In this work, we demonstrated the existence of interaction induced dynamical thermalization of fermionic atoms in a Sinai oscillator trap if the interaction strength between atoms exceeds a critical border determined by the Åberg criterion [38,39,41]. This thermalization takes place in a completely isolated system in the absence of any contact with an external thermostat. In the context of the Loschmidt–Boltzmann dispute [1,2], we should say that formally this thermalization is reversible in time since the Schrodinger equation of the system has symmetry $t\to -t$. The classically chaotic dynamics of atoms in the Sinai oscillator trap breaks in practice this reversibility due to exponential growth of errors induced by chaos. In the regime of quantum chaos, there is no exponential growth of errors due to the fact that the Ehrenfest time scale of chaos is logarithmically short [17,19,26,28]. An example for the stability of time reversibility is given in [27,28]. In fact, the experimental reversal of atom waves in the regime of quantum chaos has been even observed with cold atoms in [74]. In view of this and the fact that the spectrum of atoms in the Sinai oscillator trap is discrete, we can say that dynamical thermalization will have obligatory revivals in time returning from the thermalized state (e.g., bottom panels in Figure 8) to the initial state (top panels in Figure 8). This is the direct consequence of the Poincare recurrence theorem [75]. However, the time for such a recurrence grows exponentially with the number of components contributing to the initial state (which is also exponentially large in the regime of dynamical thermalization) and thus, during such a long time scale, external perturbations (coming from outside of our isolated system, e.g., not perfect isolation) will break in practice this time reversibility.

We hope that our results will initiate experimental studies of dynamical thermalization with cold fermionic atoms in systems such as the Sinai oscillator trap.