Dynamical thermalization of interacting fermionic atoms in a Sinai-oscillator trap

TL;DR

This paper numerically investigates how interacting fermionic atoms in a Sinai-oscillator trap undergo dynamical thermalization, resulting in Fermi-Dirac distribution without external heat baths, linked to quantum chaos.

Contribution

It demonstrates that interactions above a certain threshold induce thermalization in a quantum chaotic trap, connecting many-body dynamics with statistical mechanics.

Findings

Interactions above the Aberg criterion lead to Fermi-Dirac distribution.

The system relaxes to thermal equilibrium without contact with a thermostat.

Links between dynamical thermalization and quantum chaos are discussed.

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 Aberg criterion, lead to the Fermi-Dirac distribution and relaxation of many-body initial states to the thermalized state in absence of any contact with a thermostate. We discuss the properties of this dynamical thermalization and its links with the Loschmidt-Boltzmann dispute.

Tables

Table 1. Table 1: Parameters of the eigenstates corresponding to Figure  2 . S 𝑆 S is the entropy, E 1 ​ p subscript 𝐸 1 𝑝 E_{1p} the effective total one-particle energy, both given by ( 11 ), and E ex subscript 𝐸 ex E_{\rm ex} is the exact energy eigenvalue. Inverse temperature β 𝛽 \beta , chemical potential μ 𝜇 \mu , theoretical entropy S th subscript 𝑆 th S_{\rm th} are determined by ( 12 ) or ( 11 ) (with n k subscript 𝑛 𝑘 n_{k} replaced by n ​ ( ε k ) 𝑛 subscript 𝜀 𝑘 n(\varepsilon_{k}) ) for both energies E 1 ​ p subscript 𝐸 1 𝑝 E_{1p} , E ex subscript 𝐸 ex E_{\rm ex} .

$A$	$m$	$S$	$S_{th}\left(E_{1p}\right)$	$E_{1p}$	$\mu \left(E_{1p}\right)$	$\beta \left(E_{1p}\right)$	$S_{\mathrm{th}}\left(E_{\mathrm{ex}}\right)$	$E_{\mathrm{ex}}$	$\mu \left(E_{\mathrm{ex}}\right)$	$\beta \left(E_{\mathrm{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

Table 2. Table 2: Values of the fit parameters C 1 , … , C 5 subscript 𝐶 1 … subscript 𝐶 5 C_{1},\ldots,C_{5} for the initial states | ϕ 1 > , … , | ϕ 4 > ket subscript italic-ϕ 1 … ket subscript italic-ϕ 4 |\phi_{1}\!\!>,\ldots,|\phi_{4}\!\!> used for the analytical fits of Γ c subscript Γ 𝑐 \Gamma_{c} ( Γ F subscript Γ 𝐹 \Gamma_{F} ) in Figure  11 (Figure  13 ).

initial state	$C_1$	$C_2$	$C_3$	$C_4$	$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$

Equations

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

$$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).$$

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

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

$$V_{ij,kl}={\bar{V}}_{ij,kl}-{\bar{V}}_{ij,lk},{\bar{V}}_{ij,kl}=\int d{\mathbf{r}}_1\int d{\mathbf{r}}_2{\phi }_i^{\ast }({\mathbf{r}}_1){\phi }_j^{\ast }({\mathbf{r}}_2)v({\mathbf{r}}_1-{\mathbf{r}}_2){\phi }_k({\mathbf{r}}_1){\phi }_l({\mathbf{r}}_2)$$

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

$$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.$$

$${\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)}$$

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

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

$$S=-\sum_{k=1}^{M}\Bigl{(}n_{k}\,\ln n_{k}+(1-n_{k})\,\ln(1-n_{k})\Bigr{)}\quad,\quad E_{1p}=\sum_{k=1}^{M}\,\varepsilon_{k}\,n_{k}$$

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

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

$${\Gamma }_c=[C_1-C_2\ln (A^2)]A^2,C_1=-\sum _k{a}_k⟨\ln (a_k\xi )⟩,C_2=\sum _k{a}_k$$

$${\Gamma }_F=\frac{{\sigma }_0(L)}{{\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).$$

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

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

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

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

$$\delta E>\delta E_{ch}\approx {\Delta }_1{g}^{2/3}\approx 2\ell {\Delta }_1/({{\alpha }_s}^{2/3}{N_a}^{4/9})\sim {N_a}^{5/9}{\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}).$$

Full text

.1 I 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 loschmidt ; boltzmann (see the modern overview in mayer ; jalabert ). 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. arnold ; sinaibook ; chirikov ; lichtenberg ). 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 sinaibil .

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 decades and their properties have been mainly understood (see e.g. gutzwiller ; haake ; stockmann ). Thus it was shown that the level spacing statistics in the regime of quantum chaos bohigas1984 is the same as for Random Matrix Theory (RMT) invented by Wigner for a description of spectra of complex nuclei wigner ; mehta . This result became known as the Bohigas-Giannoni-Schmit conjecture bohigas1984 ; ullmo . Thus classically chaotic systems (e.g. Sinai billiard) are usually characterized by Wigner-Dyson (RMT) statistics with level repulsion wigner ; mehta ; bohigas1984 while the classically integrable systems usually show Poisson statistics for level spacing distribution haake ; stockmann ; ullmo . 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 chirikov1981 ; fishman ; chirikov1988 ; fishmanschol being similar to the Anderson localization in disordered solid-state systems anderson . 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 rough1 ; rough2 .

It was also shown that in the regime of quantum chaos the Bohr correspondence principle bohrcorres and the fully correct semiclassical description of quantum evolution remain valid only for a logarithmically short Ehrenfest time scale $t_{E}\sim\ln(1/\hbar)/h$ chirikov1981 ; chirikov1988 . Here $\hbar$ 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 ehrenfest . 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 to decay exponentially in contrast to exponentially decaying classical correlators dls1981 ; dls1983 . For $t>t_{E}$ the decay of quantum correlations stops and they remain on the level of quantum fluctuations being proportional to $\hbar$ dls1981 ; dls1983 ; stmapscholar . Since the level of quantum fluctuations is proportional to $\hbar$ the classical diffusive spreading over the momentum is affected by quantum corrections only on a significantly larger diffusive time scale $t_{D}\propto 1/\hbar^{2}\gg t_{E}\propto\ln(1/\hbar)$ chirikov1981 ; chirikov1988 ; dls1981 ; dls1983 ; stmapscholar .

The problem of 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 bohr . Thus on a 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. guhr ).

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 grows not faster than the number of particles in 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 french1 ; bohigas1 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, is was shown that the level spacing statistics $p(s)$ is described by the Wigner-Dyson or RMT distribution french2 ; bohigas2 .

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 $\Delta_{1}$. In metallic quantum dots this case with $U/\Delta_{1}\approx 1/g$ corresponds to a large conductance of a dot $g=E_{Th}/\Delta_{1}\gg 1$ where $E_{Th}=\hbar/t_{D}$ is the Thouless energy with $t_{D}$ being a diffusion spread time over the dot thouless ; imry ; akkermans . 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 aberg1 ; aberg2 . The estimate of a critical interaction $U_{c}$, called the Åberg criterion dlsnobel , compares the typical two-body matrix elements with the energy level spacing $\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\Delta_{c}$ and the RMT statistics sets in for $U>U_{c}\sim\Delta_{c}$. In jacquod this criterion, proposed independently of aberg1 ; aberg2 , 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 jacquod :

[TABLE]

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 jacquod . The validity of the Åberg criterion was numerically confirmed for various physical models (see dlsnobel and Refs. therein). An additional confirmation was given by the analytical derivation presented in sushkov showing that for 3 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 $\Delta_{2}\sim\Delta_{c}$ being parametrically larger than the three-particle level spacing $\Delta_{3}\ll\Delta_{2}$. The advanced theoretical arguments developed in mirlin1 ; mirlin2 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 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\Delta_{c}$ has been demonstrated in kolovsky2017 ; frahmtbrim 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 landau .

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 sinaibil . For one particle it has been shown in ermannsinai 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 of the TBRIM. On the other side the Sinai-oscillator trap has been already experimentally realized with Bose-Einstein condensate of cold bosonic atoms ketterle1 ; ketterle2 ; ketterle3 . At present cold atom techniques allow to investigate various properties of cold interacting fermionic atoms roati1 ; roati2 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 is 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 sachdevprl ; kitaev ; sachdevprx ; rosenhaus ; maldacena ; garcia1 . In fact, the SYK model, in its fermionic formulation sachdevprx , corresponds to the TBRIM considered with a conductance close to zero $g\rightarrow 0$. In these lines the dynamical thermalization in TBRIM and SYK systems has been discussed in kolovsky2017 ; frahmtbrim . 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. huse ; polkovnikov ; borgonovi ; alet ). 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.

I II Quantum chaos in Sinai oscillator

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

[TABLE]

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 ketterle1 . The Poincare sections at different energies are presented in ermannsinai 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\hbar$ ermannsinai . For the quantum problem we use the value of the dimensionless Planck constant $\hbar=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}=\hbar\omega_{x}$ (for our choice of Sinai oscillator parameters we also have $E_{u}=\hbar\omega_{x}=\hbar^{2}/(m{r_{d}}^{2})$) ermannsinai with the typical size of oscillator ground state being equal to the disk radius: $a_{0}=\Delta x_{osc}=(\hbar/m\omega_{x})^{1/2}=r_{d}$.

In ermannsinai it is shown that the classical dynamics of this system is 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(\varepsilon)=\varepsilon^{2}/(2\sqrt{2})-\varepsilon/2$ ermannsinai . Thus the one-particle density of states $\rho_{1}(\varepsilon)$ and corresponding level spacing $\Delta_{1}$ are:

[TABLE]

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

II III Sinai oscillator with interacting fermionic atoms

II.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 frame work of the scattering length approach (see e.g. gribakin1 ; gribakin2 ). Therfore 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. wilkens ). For example in a typical experimental situation the disk radius is of the order of micron $r_{d}\sim 1\mu m=10^{-4}cm$ while for Li and other alkali atoms we have $r_{c}\sim 3\times 10^{-7}cm$ gribakin1 ; gribakin2 .

In the following, we use a simple interaction function having a constant amplitude $U$ for $r\leq 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\bar{r}_{c}$ can be absorbed in a modified amplitude according to $U\to\bar{U}=U(\bar{r}_{c}/r_{c})^{2}$, a relation we verified numerically for various values of $\bar{r}_{c}<r_{d}$. 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 kohler .

II.2 Reduction to TBRIM like case and its analysis

Using the methods described in ermannsinai we numerically compute a certain number of one-particle or orbital energy eigenvalues $\varepsilon_{k}$ and corresponding eigenstates $\varphi_{k}({\bf r})$ of the Sinai oscillator (2). As repulsive interaction potential $v({\bf r})$ we choose the short ranged box function $v({\bf r})=U$ if $|{\bf r}|\leq r_{c}=0.2$ (since $r_{d}=1$) and $v({\bf r})=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\leq L\leq M$ spinless fermions takes the form:

[TABLE]

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

[TABLE]

and $c^{\dagger}_{k}$, $c^{\phantom{\dagger}}_{k}$ 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 $\bar{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: $\bar{V}_{ij,kl}=\bar{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 to reduce also the $k,l$ sums to $k<l$. Furthermore, the ordering of the two fermion operators $c^{\phantom{\dagger}}_{l}c^{\phantom{\dagger}}_{k}$ 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 trace111One 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 and does not affect the issues of thermalization, interaction induced eigenfunction mixing or the quantum time evolution with respect to the Hamiltonian $H$ etc.

II.3 Åberg parameter

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

[TABLE]

where $\{n_{k}\}$ 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 Eq. (A.4) of Ref. frahmtbrim ):

[TABLE]

Therefore the many-body level spacing $\Delta_{\rm MB}$ or inverse Heisenberg time at the band center $E=E_{\rm mean}$ is given by $\Delta_{\rm MB}=1/t_{\rm 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{\varepsilon^{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{\varepsilon^{n}}\approx 2\,\varepsilon_{M}^{n}/(n+2)$ and $\overline{\varepsilon^{2}}-\overline{\varepsilon}^{2}\approx\varepsilon_{M}^{2}/18\approx\sqrt{2}\,M/9$.

For the question if 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 $\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$ jacquod ; flambaum and we need to use the variance for only two particles:

[TABLE]

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/\Delta_{c}=K/[\sqrt{2\pi}\,\sigma_{0}(L=2)]$ 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 $\varepsilon_{k}$ values).

According to the Åberg criterion aberg1 ; aberg2 ; jacquod the onset of chaotic mixing happens for typical interaction matrix elements $U$ comparable to $\Delta_{c}$. Therefore we compute the quantity $V_{\rm 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 case of a long tail distribution. However, in our case it turns out that $V_{\rm mean}\approx 2\exp(\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_{\rm mean}/\Delta_{c}=U/U_{c}$ such that $A=1$ if $U=U_{c}$. We expect aberg1 ; aberg2 ; jacquod 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_{\rm mean}$, numerically computed for $4\leq M\leq 30$, can be quite accurately fitted by $V_{\rm mean}\approx 3\times 10^{-4}\,U/\varepsilon_{M}$. Furthermore, we remind the expression $\Delta_{c}=(1/K)\,\sqrt{4\pi(M-2)(\overline{\varepsilon^{2}}-\overline{\varepsilon}^{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: $\Delta_{c}=4/3\,\sqrt{2\pi}\,\varepsilon_{M}/[(M-L)^{2}\,L^{2}]$. Here we also used the above found expression $\overline{\varepsilon^{2}}-\overline{\varepsilon}^{2}\approx\varepsilon_{M}^{2}/18$. From this we find that $U_{c}=\Delta_{c}U/V_{\rm mean}\approx C\,M/[(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 $\varepsilon_{M}^{2}\approx 2\sqrt{2}\,M$ according to (3). Below we will give more accurate numerical values of $V_{\rm mean}$, $\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 $\Delta_{c}$ is enhanced thus reducing the effective value of $A$. Furthermore, we computed $V_{\rm 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).

II.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\leq A\leq 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_{\rm mean}=3.865\times 10^{-5}\,U$ and from (8) that $\Delta_{c}=\sqrt{2\pi}\,[\sigma_{0}(L=2)/K]=6.1706\times 10^{-3}$ where the quantities $\overline{\varepsilon^{n}}$ where exactly computed from the numerical orbitals energies $\varepsilon_{k}$. From this we find that $U_{c}=\Delta_{c}\,U/V_{\rm 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 presence of interactions, the density of states has approximately a Gaussian form with the same center $E_{\rm 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}(A)$ of the Gaussian density of states by a fit of the integrated density of states $P(E)$ using

[TABLE]

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

[TABLE]

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(A)$ using the expression $\sigma^{2}(A)=\mbox{Tr}_{\rm Fock}\left[(H-E_{\rm mean}{\bf 1})^{2}\right]/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\leq 10$. For larger values of $A$ the deviations increase but the overall form is still correct. As described in frahmtbrim the quality of the fit can be considerably improved if we replace in (9) the linear function $q(E)$ 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_{\rm mean}^{2}$. In this case one can show theoretically frahmtbrim that $\sigma^{2}(A)=\sigma_{0}^{2}+K_{2}V_{\rm 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 frahmtbrim 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_{\rm mean}=A\Delta_{c}=A\sqrt{2\pi}\,\sigma_{0}(L=2)/K$ we indeed obtain (10) with $\alpha=\alpha_{\rm th}=4\pi(M-2)K_{2}/[K^{2}(L(M-L)]$. For $M=16$, $L=7$ we find $\sigma_{0}=3.0279$ (see (7)) and $\alpha_{\rm 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_{\rm 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. dlsnobel for the case of interacting Fermi systems).

II.5 Thermalization and entropy of eigenstates

In the following, we mostly concentrate on values $A\leq 10$ such that the effect of the increase of the spectral width $\sigma(A)$ 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 landau . As in kolovsky2017 ; frahmtbrim we determine the occupation numbers $n_{k}=\langle c^{\dagger}_{k}c^{\phantom{\dagger}}_{k}\rangle$ for such a state, as well as the corresponding fermion entropy $S$ landau and the effective total one-particle energy $E_{1p}$ by :

[TABLE]

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(A)\approx\sigma_{0}$, we have typically $E_{1p}\approx E_{\rm 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(\varepsilon_{k})$ with $n(\varepsilon)=1/(1+\exp[\beta(\varepsilon-\mu)])$ where inverse temperature $\beta=1/T$ and chemical potential $\mu$ are determined by the conditions:

[TABLE]

Here $E$ is normally given by $E_{1p}$ but one may also consider the value $E_{\rm 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_{\rm th}(E)$ using (11) with $n_{k}$ being replaced by $n(\varepsilon_{k})$ (where $\beta$, $\mu$ are determined from (12) for the energy $E$).

The many-body states with energies above $E_{\rm 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\leq E\leq 39\approx E_{\rm mean}$ (for $M=16$, $L=7$). In Figure 2 we compare the thermalized Fermi-Dirac occupation number $n(\varepsilon)$ with the 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.

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 [math] 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_{\rm th}$ (see Table 1 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_{\rm 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_{\rm 1p}$ as reference energy to compute temperature and chemical potential. The temperature obtained from $E_{\rm ex}$ is too small because here the increase of $\sigma(A)$ is already quite strong and $E_{\rm ex}$ rather strongly deviates from $E_{\rm 1p}$. Also the value of $S_{\rm th}$ using $E_{\rm ex}$ does not match $S$. Obviously at stronger interaction values it is necessary to use $E_{\rm 1p}$ 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 $\varepsilon_{k}$ and is given by $\mu=-\,\ln[1+(M-2L)/L]/\beta$ providing a singularity if $L\neq 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_{\rm mean}$.

In Figure 4 we present the nearest level spacing distribution $p(s)$ for different values of the Åberg parameter. To compute $p(s)$ 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(E)$ is replaced by a fit polynomial of degree 5. For the smallest value $A=0.2$ the distribution $p(s)$ 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\leq 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 wigner .

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_{\rm 1p}$ or $E_{\rm ex}$ with the theoretical thermalized entropy $S_{\rm th}(E)$. 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_{\rm th}$. Actually the distribution of data points is considerably concentrated at smaller entropy values which is not so clearly visible in the Figure. In particular the average of the ratio of $S/S_{\rm th}(E_{\rm 1p})$ 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_{\rm 1p}$) are very close to the theoretical line with the average ratio $S/S_{\rm th}(E_{\rm 1p})$ 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_{\rm ex}$ and $E_{\rm 1p}$ are still rather close. For $A=10$ all states are well thermalized (when using the energy $E_{\rm 1p}$) while the data points for $E_{\rm 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_{\rm 1p}$ (but more concentrated on energy values closer to the center with larger entropy values and larger temperatures) while for $E_{\rm 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).

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(\varepsilon_{k})$ (shown in the same way) provides further confirmation that at $A=3.5$ there is indeed already a quite strong thermalization of most eigenstates.

II.6 Thermalization of quantum time evolution

The question arises how or if a time dependent state $|\psi(t)\!\!>=\exp(-iHt)|\psi(0)\!\!>$, obeying the quantum time evolution with the Hamiltonian $H$ and an initial state $|\psi(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) $|\phi_{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) $|\phi_{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) $|\phi_{3}\!\!>=|0000011000110111\!\!>$ and (iv) $|\phi_{4}\!\!>=|1000100011001011\!\!>$. The states $|\phi_{3}\!\!>$ and $|\phi_{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$ ($|\phi_{1}\!\!>$), $E\approx 32$ ($|\phi_{2}\!\!>$ and $|\phi_{3}\!\!>$) and $E\approx 35$ ($|\phi_{4}\!\!>$).

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

[TABLE]

where $t_{\rm H}$ is the Heisenberg time (at the given value of $A$), $d$ the dimension of the Hilbert space and $\sigma(A)$ the width of the Gaussian density of states given in (10). The quantity $\Delta t$ is the shortest physical time scale of the system (inverse of the largest energy scale) and obviously for $t\ll\Delta t$ the unitary evolution operator is close to the unit matrix multiplied by a uniform phase factor: $\exp(-iHt)\approx\exp(-iE_{\rm mean}t)\,{\bf 1}$ since the eigenvalues $E_{m}$ of $H$ satisfy $|E_{m}-E_{\rm mean}|\lesssim\sigma(A)$. We expect that any signification deviation of $|\psi(t)\!\!>$ with respect to the initial condition $|\psi(0)\!\!>$ happens at $t\geq\Delta t$ (or later in case of very weak interaction). Furthermore, by analyzing the time evolution in terms of the ratio $t/\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_{\rm H}\approx 10^{4}\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}(t)=<\!\!\psi(t)|c^{\dagger}_{k}c^{\phantom{\dagger}}_{k}|\psi(t)\!\!>$ and present them in a color density plot in the plane $(k,\;t/\Delta t)$. Also, at the last used time value we compute the effective total one-particle energy $E_{\rm 1p}$ using the relation (11) (note that $E_{\rm 1p}$ 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(\varepsilon_{k})$ at each $k$ value for the orbital index. These values of ideally thermalized occupation numbers will be shown in an additional vertical bar 222Note 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 $|\phi_{1}\!\!>$ and the two Åberg parameter values $A=1$ and $A=3.5$ using a linear time scale with integer multiples of $\Delta t$ and for $t\leq 2000\,\Delta t\approx t_{\rm H}/6$. At $A=1$ the occupation number $n_{\rm 12}$ (of the excited particle) shows at the beginning a periodic structure, with an approximate period $400\,\Delta t$ for $t<1000\,\Delta t$, and a modest decay for $t>1000\,\Delta t$. At the same time the first orbitals above the Fermi sea are slightly excited. At final $t=2000\,\Delta t$ the state is clearly not thermalized. For $A=3.5$ we see a very rapid partial decay of $n_{6}$ and $n_{\rm 12}$ together with an increase of $n_{\rm 7}$. Furthermore for $n_{k}$ with $8\leq k\leq 11$ there are later and more modest excitations with a periodic time structure. Here the final state at $t=2000\,\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 well capture a rapid decay/increase of $n_{k}$ at small times and its maximal time value is also significantly limited below the Heisenberg time. Therefore we use in Figures 8 and 9 a logarithmic time scale with $0.1\,\Delta t\leq t\leq 10^{6}\,\Delta t\approx 10^{2}\,t_{\rm 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 $|\phi_{1}\!\!>$ and $|\phi_{2}\!\!>$ is shown for the Åberg values $A=1,\,3.5,\,10$. For $|\phi_{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\,\Delta t$ up to $10^{6}\,\Delta t$ and for $A=3.5$ the initial decay of $n_{6}$ and $n_{12}$ happens at $t\approx 10\,\Delta t$. For $|\phi_{1}\!\!>$ at $A=10$ the decay starts at $t\approx 3\,\Delta t$ and an approximate thermalization happens at $t>40\,\Delta t$. But here there is still some time periodic structure and it would be necessary to do some time average to have perfect thermalization. For $|\phi_{2}\!\!>$ at $A=1$ the decay of excited orbitals 12 and 14 starts at $t\approx 100\,\Delta t$ and saturates at $t\approx 1000\,\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\,\Delta t$. The final state at $t=10^{6}\,\Delta t$ is not thermalized even though some orbitals have $n_{k}$ values close to thermalization. For $|\phi_{2}\!\!>$ at $A=3.5$ the decay of excited orbitals 12 and 14 starts at $t\approx 10\,\Delta t$ and for $t>300\,\Delta t$ there is thermalization (but requiring some time average as for $|\phi_{1}\!\!>$ at $A=10$). Interestingly at intermediate times $10\,\Delta t<t<100\,\Delta t$ the high orbitals 13 and 16 are temporarily slightly excited and decay afterwards rather quickly to their thermalized values. For $|\phi_{2}\!\!>$ at $A=10$ the decay of excited orbitals 12 and 14 starts even at $t\approx 3\,\Delta t$ and thermalization seems to set in at $t>30\,\Delta t$.

Figure 9 is similar to Figure 8 but for the initial states $|\phi_{3}\!\!>$ and $|\phi_{4}\!\!>$ which have occupations 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\,\Delta t$ ($t\approx(10-20)\,\Delta t$ or $t\approx(2-3)\,\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 $|\phi_{3}\!\!>$ and good thermalization for $A=3.5$ and $|\phi_{4}\!\!>$ as well as $A=10$ (both states).

Using the time dependent values $n_{k}(t)$ one can immediately determine the corresponding entropy $S(t)$ using (11). At $t=0$ we have obviously $S(0)=0$ since 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_{\rm th}(E_{\rm 1p})$ (with $E_{\rm 1p}$ determined for the state $|\psi(t)\!\!>$ at large times) depending if there is presence (or absence) of thermalization according to the different cases visible in Figures 8 and 9. However, in 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(t)$ for the two initial states $|\phi_{1}\!\!>$, $|\phi_{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 which is indeed close to the thermalized entropy $S_{\rm th}(E_{\rm 1p})$. We note that $E_{\rm 1p}$ is not conserved at strong interactions and that its initial value $E_{\rm 1p}\approx 30$ ($E_{\rm 1p}\approx 35$) at $t=0$ evolves to $E_{\rm 1p}\approx 33.5$ ($E_{\rm 1p}\approx 37$) at large times for $|\phi_{1}\!\!>$ ($|\phi_{4}\!\!>$) corresponding roughly to $S\approx S_{\rm th}(E_{\rm 1p})\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_{\rm 1p}$ values. For $A=3.5$ and $|\phi_{4}\!\!>$ there is also saturation of $S$ to its thermalized value. For $A=3.5$ and $|\phi_{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 $|\phi_{4}\!\!>$ there is a quite late and approximate saturation with some fluctuations which is visible for $t>10^{4}\,\Delta t$ and with $S\approx 10\pm 0.2$. For $A=1$ and $|\phi_{1}\!\!>$ there is a late periodic regime for $t>10^{3}\,\Delta$ with a quite large amplitude $S\approx 3\pm 1$ and with $S_{\rm max}\approx 4$ significantly below the thermalized entropy $S_{\rm th}(E_{\rm 1p})\approx 5.5$. The panels using a normal (instead of logarithmic) time scale with $t\leq 200\,\Delta t$ miss completely the long time limits for $A=1$ and might wrongly suggest that there is an early saturation at quite low values of $S$.

The periodic (or quasi-periodic) time dependence of $n_{k}(t)$ or $S(t)$, 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,\,\ldots$) of exact eigenstates of $H$ contribute mostly in the expansion of $|\psi(t)\!\!>$ in terms of these eigenstates.

Figure 10 also shows that the initial increase of $S(t)$ 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(t)$ is close to a quadratic behavior for $t\lesssim\Delta t$ but which immediately becomes linear for $t\gtrsim\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 333For practical reasons we decide to incorporate the quantum time step $\Delta t$ in the definition of $\Gamma_{c}$, i.e. $\Gamma_{c}$ is defined as the ratio of the initial slope $S^{\prime}(t)$ over the global spectral bandwidth $\sim\sigma(A)\sim 1/\Delta t$. the quantity $\Gamma_{c}=dS(t)/d(t/\Delta t)=\Delta tS^{\prime}(t)$ for $t=\xi\Delta t$ where $\xi\gtrsim 1$ is a numerical constant of order one. To determine $\Gamma_{c}$ practically we perform first the fit $S(t)=\bar{S}_{\infty}\,(1-\exp[-\bar{\gamma}_{1}(t/\Delta t)])$ for $0\leq t/\Delta t\leq 100$ and use the exponential decay rate $\bar{\gamma}_{1}$ to perform a refined fit $S(t)=S_{\infty}\,(1-\exp[-\gamma_{1}(t/\Delta t)-\gamma_{2}(t/\Delta t)^{2}])$ for the interval $0\leq t/\Delta t\leq 5/\bar{\gamma}_{1}$. From this we determine $\Gamma_{c}=S_{\infty}\,\gamma_{1}$ which is rather close to $\bar{S}_{\infty}\,\bar{\gamma}_{1}$ for $A\leq 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(t)$ for smaller values $A\leq 1$.

Figure 11 shows the dependence of these values of $\Gamma_{c}$ on the parameter $A$ for our four initial states. At first sight on observes a behavior $\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 $\Gamma_{c}\propto A^{p}$ for $A\leq 2$ provide exponents close to $p\approx 1.75-1.85$) and it turns out that these deviations correspond to a logarithmic correction: $\Gamma_{c}=f(A)=(C_{1}-C_{2}\,\ln[g(A)])\,g(A)$ with $g(A)=A^{2}$ (for fits with $A\leq 1$) or with $g(A)=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}(t)\approx 1-\delta n_{k}(t)$ (if $n_{k}(0)=1$) or $n_{k}(t)\approx\delta n_{k}(t)$ (if $n_{k}(0)=0$) where $\delta n_{k}(t)$ is the small modification of $n_{k}(t)$. Time dependent perturbation theory suggests that $\delta n_{k}(t)\sim(t/\Delta t)^{2}$ for $t\lesssim\Delta t$ and $\delta n_{k}(t)\approx a_{k}\,A^{2}\,t/\Delta t$ for $t\gtrsim\Delta t$ such that still $\delta n_{k}(t)\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 $\Gamma_{c}=\Delta t\,S^{\prime}(t)\approx-\Delta t\sum_{k}\delta n_{k}^{\prime}(t)\ln\delta n_{k}(t)$ with $t=\xi\,\Delta t$. Since $\delta n_{k}^{\prime}(t)\approx a_{k}A^{2}/\Delta t$ we find indeed the behavior :

[TABLE]

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 Figures 8 and 9) but the coefficients $C_{1}$, $C_{2}$ depend only slightly on the initial state (see Table 2). Furthermore, by replacing $A^{2}\to g(A)=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 $\Gamma_{c}=f(a)$ for the full range of $A$ values.

II.7 Survival probability and Fermi’s golden rule

The knowledge of the time dependent states $|\psi(t)\!\!>$ allows us also to compute the decay function $p_{\rm dec}(t)=|<\!\!\psi(0)|\psi(t)\!\!>|^{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\Delta t$ we expect a quadratic decay: $1-p_{\rm dec}(t)\approx\langle(H-E_{\rm mean})^{2}\rangle t^{2}\approx{\rm const.}\,(t/\Delta t)^{2}$ with $\langle\cdots\rangle$ being the quantum expectation value with respect to $|\psi(0)\!\!>$ and a numerical constant $\lesssim 1$ since $1/\Delta t$ represents roughly the spectral width of $H$. For $t\gtrsim\Delta t$ but such that $1-p_{\rm dec}(t)\ll 1$ we have according to Fermi’s golden rule: $1-p_{\rm dec}(t)=\Gamma_{F}\,(t/\Delta t)$ where $\Gamma_{F}$ is the decay rate 444 Again for practical reasons and similarly to $\Gamma_{c}$ we incorporate in the definition of $\Gamma_{F}$ the time scale $\Delta t$, i.e. $\Gamma_{F}=\Delta t\ \times$ usual decay rate found in the literature and meaning that $\Gamma_{F}$ is defined as the ratio of the usual decay rate over the global spectral bandwidth. of the state.

To determine numerically $\Gamma_{F}$ we apply the fit: $p_{\rm dec}(t)=C\,\exp(-\Gamma_{\rm F}\,t/\Delta t)$ in two steps. First we use the interval $1\leq t/\Delta t\leq 50$ and if $5/\Gamma_{F}<50$, corresponding to a rapid decay (which happens for larger values of $A$) we repeat the fit for the reduced interval $1\leq t/\Delta t\leq 5/\Gamma_{F}$. The choice of the Amplitude $C\neq 1$ and the condition $t\geq\Delta t$ for the fit range allow to take into account the effects due to the small initial window of quadratic decay. In Figure 12 we show two examples for the initial state $|\phi_{1}\!\!>$ and the Åberg values $A=3.5$ and $A=10$. In both cases the shown maximal time value $t_{\rm max}=50\,\Delta t$ (if $A=3.5$) or $t_{\rm 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\leq t/\Delta t\leq 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\leq t/\Delta t\leq 30$. For very small values of $A$ the fit procedure works also correctly since it captures only the initial decay which is important if $p_{\rm dec}(t)$ does not decay completely at large times and which typically happens in the perturbative regime $A\lesssim 1$.

Figure 13 shows the dependence of $\Gamma_{F}$ on $A$ for the usual four initial states together with the fit $\Gamma_{F}=f(A)=C_{4}A^{2}/(1+C_{5}A^{2})$ for the data with initial state $|\phi_{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 $\Gamma_{F}\propto A^{2}$ is highly accurate (with no logarithmic correction). Similarly to $\Gamma_{c}$ there is only a slight dependence of the values of $\Gamma_{F}$ and the fit values on the choice of initial state.

Theoretically, we expect according to Fermi’s Golden rule that: $\Gamma_{F}\approx(\Delta t)\,2\pi V_{\rm Fock}^{2}\,\rho_{c}(E)$ where $V_{\rm Fock}^{2}=\mbox{Tr}_{\rm Fock}(V^{2})/(Kd)=\sigma_{0}^{2}\,\alpha A^{2}/K$ according to the discussion below (10) and $\rho_{c}(E)$ is the effective two-body density of states for states directly coupled by the interaction such that $\rho_{c}(E_{\rm mean})=1/\Delta_{c}$ (see discussion below (7)). We note that $V_{\rm Fock}$ is the typical interaction matrix element in Fock space which is slightly larger than $V_{\rm mean}$ (see the theoretical discussion above for the computation of the coefficient $\alpha$ used in (10) and Appendix A of frahmtbrim ). The factor $\Delta t$ is due to our particular definition ${}^{\ref{foot4}}$ of decay rates. The expression of $\Gamma_{F}$ is actually also valid for larger values of $A$ provided we use the density of states $\rho_{c}(E)$ in the presence of interactions which provides an additional factor $1/\sqrt{1+\alpha A^{2}}$ according to (10). Therefore, at the band center we have: $2\pi\,\Delta t\,\rho_{c}=K/[\sigma_{0}(L)\,\sigma_{0}(L=2)(1+\alpha A^{2})]$ which gives together with (8):

[TABLE]

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 Figures 8 and 9 we have typical decay times in units of $\Delta t$ being $1/\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 above discussion of Figures 8 and 9).

II.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

[TABLE]

where $\varphi_{k}(x,y)$ is the one-particle eigenstate of orbital $k$, with some examples shown in Figure 1. Here $\Psi^{(\dagger)}(x,y)$ denote the usual fermion field operators (in 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 I). 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\leq 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:

[TABLE]

where $n_{kl}(t)$ is the density matrix in orbital representation generalizing the occupation numbers $n_{k}(t)$ which are its diagonal elements. Due to the complex phases of $|\psi(t)\!\!>$ (when expanded in the usual basis of non-interacting many-body states) the density matrix is complex valued but hermitian: $n_{kl}^{*}(t)=n_{lk}(t)$. 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(0)\!\!>=|\phi_{2}\!\!>$ and four time values $t/\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 $\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 quantware ).

At the time $t/\Delta t=0.1$ density matrix and spatial density are essentially identical to the initial condition at $t=0$. For $\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 [math] (for initially empty orbitals) and the spatial density simply gives the sum of densities due to the occupied eigenstates.

At $t/\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/\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 stronger excited than the initially excited orbital 12. Also orbital 14 and 15 are quite strong. The spatial density has become smoother and the structure of $\Delta\rho$ is roughly close to the case at $t/\Delta t=30$ but with some significant differences.

Finally at $t/\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 $\Delta\rho$ has a similar but different structure as for $t/\Delta t=100$ or $t/\Delta t=30$.

Apparently at intermediate times $20\leq t/\Delta t\leq 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 frahmtbrim 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}\varphi_{k}^{2}(x,y)\,n_{k}(t)$ in terms of the time dependent occupation numbers. Only for intermediate time scales with more visible quantum coherence (more off-diagonal elements $n_{kl}(t)\neq 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 $\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 $\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(\varepsilon_{k})$ at energy $E=32.9$ which is the typical total one-particle average energy of $|\psi(t)\!\!>$ for the long time limit $t/\Delta t\geq 1000$ showing that at $t/\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:

[TABLE]

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 $\varphi_{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}(t)\neq 0$ for $k\neq l$. In Figure 15 we present density plots of absolute value, real and imaginary part of $\rho_{\rm 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_{\rm corr,th}(x,y;x_{0},y_{0})$ vanishes (giving a blue panel due to zero values).

There are significant time dependent fluctuations of $\rho_{\rm 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 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/\Delta t=1000$ seems to be rather close to the thermalized case (for absolute value and real part) there are still differences which are more significant here as in Figure 14.

In global the obtained results show that the dynamical thermalization well takes place 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.

III IV Estimates for cold atom experiments

We discuss here typical parameters for cold atoms in a trap following ketterlefermi . Thus for sodium atoms we have $\omega\approx\omega_{x,y,z}\approx 2\pi 10Hz$ with $a_{0}=\sqrt{\hbar/(m\omega)}=6.5\mu m$ and oscillator level spacing $E_{u}=\hbar\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/(a_{0})^{3}\approx 4\times 10^{9}cm^{-3}$. Since $a_{s}\ll a_{0}$ the two-body interaction is of $\delta$-function type with $v({{\bf r_{1}-r_{2}}})=(4\pi\hbar^{2}a_{s}/m)\delta({\bf r_{1}-r_{2}})$ wilkens ; ketterlefermi . In our numerical simulations $\delta({\bf r})$ is replaced by a function $H(r_{c}-|{\bf r}|)/(C\,2^{d}r_{c}^{d})$ with a small $r_{c}$, volume $C$ of the unit sphere in $d$ dimensions and with the Heaviside function $H(x)=1$ (or [math]) for $x\geq 0$ and $H=0$ for $x<0$. Hence, our parameter $U$ introduced in Section II corresponds to $U=(C\,2^{d}r_{c}^{d})(4\pi\hbar^{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 3D Sinai oscillator following the lines of Eq.(1). In such a case the two-body interaction energy scale between atoms is $U_{s}=4\pi\hbar^{2}\rho_{0}a_{s}/m=4\pi(a_{s}/a_{0})\hbar\omega$ wilkens ; ketterlefermi so that $U_{s}/\hbar\omega\sim 6\times 10^{-3}$. Comparing to sodium the mass of Li atoms is approximately 3 times smaller so that for the same $\omega$ we have $a_{0}\approx 10\mu m$ and $U_{s}/\hbar\omega\approx 4\times 10^{-3}$. We think that the scattering length can be significantly increased via the Feshbach resonance allowing to reach effective interaction values $U_{s}/\hbar\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 2D case the scattering on an elastic ball in the trap center leads to quantum chaos and chaotic eigenstates with $\ell\leq N_{a}$ components (e.g. in the basis of oscillator eigenfunctions). The Fermi energy of the trap is then $E_{F}=\hbar(N_{a}\omega_{x}\omega_{y}\omega_{z})^{1/3}\approx\hbar\omega{N_{a}}^{1/3}$ roati1 ; roati2 . 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

[TABLE]

The derivation of this estimate is very similar to the case of two interacting particles in a disordered potential with localized eigenstates tip . At the same time in a vicinity of the Fermi energy $E_{F}$ we have the one-particle level spacing $\Delta_{1}=dE_{F}/dN_{a}\approx\hbar\omega/(3{N_{a}}^{2/3})$. Hence the effective conductance appearing in (1) is $g=\Delta_{1}/U_{2}\approx\ell^{3/2}/(3\alpha_{s}{N_{a}}^{2/3})$. 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:

[TABLE]

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\Delta_{1}$. Here we used the maximal value for the number of components $\ell\sim N_{a}$. It is possible that in a reality $\ell$ can be significantly smaller than $N_{a}$. However, the determination of the dependence $\ell(N_{a})$ requires separate studies taking into account the properties of chaotic eigenstates and their spreading over the energy surface. These spreading can have rather nontrivial properties (see e.g. rough2 ). This is confirmed by the results presented in Appendix for the 2D case of Sinai oscillator showing the numerically obtained dependence of two-body matrix element on energy for transitions in a vicinity of Fermi energy $E_{F}$ (see Fig. 16 there).

IV V Discussion

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 aberg1 ; aberg2 ; jacquod . This thermalization takes place in a completely isolated system in absence of any contact with an external thermostat. In the context of the Loschmidt-Boltzmann dispute loschmidt ; boltzmann we should say that formally this thermalization is reversible in time since the Schrodinger equation of the system has symmetry $t\rightarrow-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 chirikov1981 ; chirikov1988 ; dls1981 ; stmapscholar . An example for the stability of time reversibility is given in dls1983 ; stmapscholar . In fact the experimental reversal of atom waves in the regime of quantum chaos has been even observed with cold atoms in hoogerland . 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 poincare . 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.

Acknowledgments

We are thankful to Shmuel Fishman for deep discussions of quantum chaos problems and related scientific topics during many years.

NOTE: This work represents the contribution to the Special Issue in memory of Shmuel Fishman to appear at MDPI Condenced Matter journal.

This work was supported in part by the Programme Investissements d’Avenir ANR-11-IDEX-0002-02, reference ANR-10-LABX-0037-NEXT (project THETRACOM). This work was granted access to the HPC GPU resources of CALMIP (Toulouse) under the allocation 2018-P0110.

A1. Two-body matrix element near the Fermi energy

In this Appendix we present numerical results for the dependence of the quantity $V_{\rm mean}(\varepsilon)=\sqrt{\langle V_{ij,kl}^{2}\rangle}$ as a function of $\varepsilon$ where the average is done only for orbitals with energies $\varepsilon_{n}$ close to $\varepsilon$ (for $n\in\{i,j,k,l\}$), i.e.: $|\varepsilon_{n}-\varepsilon|\leq\Delta\varepsilon$ with $\Delta\varepsilon=2$. We note that this is different from the quantity $V_{\rm mean}$ used in Section II where the average was done over all orbitals (up to a maximal number being $M$). The reason for the special average with orbital energies close to $\varepsilon$ (which will be identified with the Fermi energy $E_{F}$) is that these transitions are dominant in presence of the Pauli blockade near the Fermi level.

We remind that according to the discussion of Sections I and II the matrix elements $V_{ij,kl}$ were computed for an interaction potential of amplitude $U$ for $|{\bf r_{1}-r_{2}}|<r_{c}$ (with the radius $r_{c}=0.2r_{d}=0.2$) and being zero for $|{\bf r_{1}-r_{2}}|\geq r_{c}$. Furthermore, they have been anti-symmetrized and a diagonal shift $V_{ij,ij}\to V_{ij,ij}-(1/M_{2})\sum_{k<l}V_{kl,kl}$ was applied to ensure that the interaction matrix has a vanishing trace.

Due to this shift and the precise average procedure there is a slight (purely theoretical) dependence on the maximal orbital number $M$ for this average (there is a cut-off effect for $\varepsilon$ close to the maximal orbital energy $\varepsilon_{M}$). Due to this we considered two values of $M=30$ and $M=60$.

The numerically obtained dependence is shown in Fig. 16 and is well described by the fit $V_{\rm mean}/U=a/\varepsilon^{b}$ with $a=1.56\times 10^{-4}$ and $b=0.78$. The small value of $a$ is due to antisymmetry of two-particle fermionic states and due to a small value of $r_{c}=0.2r_{d}$ which leads to a decrease of the effective interaction strength being proportional to ${r_{c}}^{2}$.

We note that the Fermi energy is determined by the number of fermionic atoms $N_{a}$ inside the 2D Sinai oscillators with $\varepsilon=E_{F}\approx\omega{N_{a}}^{1/2}$ assuming $\omega=\omega_{x}\approx\omega_{y}$. Therefore we have $\varepsilon\propto M^{1/2}\sim{N_{a}}^{1/2}$, $\Delta_{1}\sim\hbar\omega/{N_{a}}^{1/2}$ and $V_{\rm mean}\sim\alpha_{s}\hbar\omega/\ell^{3/2}\sim\alpha_{s}\hbar\omega/{N_{a}}^{b/2}$ (see (19)). Hence, the obtained exponent $b\approx 0.78$ corresponds to the number of one-particle components $\ell\sim{N_{a}}^{b/3}\sim{N_{a}}^{0.25}\sim{n_{x}}^{0.5}$. At the moment we do not have a clear explanation for this numerical dependence. This dependence corresponds to $g=\Delta_{1}/V_{\rm mean}\sim\ell^{3/2}/(\alpha_{s}{N_{a}}^{1/2})\sim 1/(\alpha_{s}{N_{a}}^{1/8})$. For such a dependence we obtain that the DTC border in 2D takes place for an excitation energy $\delta E>\delta E_{ch}\sim g^{2/3}\Delta_{1}\sim\hbar\omega/({\alpha_{s}}^{2/3}{N_{a}}^{7/12})$. Thus the thermalization can take place at rather low energy excitations above the Fermi energy with $\Delta_{1}<\delta E\ll E_{F}$.

V References