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 *tE* ∼ ln(1/¯*h*)/*h* [17,19]. Here, *h*¯ 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 *tE*, 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* > *tE*, the decay of quantum correlations stops and they remain on the level of quantum fluctuations being proportional to *h*¯ [26–28]. Since the level of quantum fluctuations is proportional to *h*¯, the classical diffusive spreading over the momentum is affected by quantum corrections only on a significantly larger diffusive time scale *tD* ∝ 1/¯*h*2 *tE* ∝ ln(1/¯*h*) [17,19,26–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 *δ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*(*s*) 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 Δ1. In metallic quantum dots, this case with *U*/Δ1 ≈ 1/*g* corresponds to a large conductance of a dot *g* = *ETh*/Δ1 1, where *ETh* = *h*¯ /*tD* is the Thouless energy with *tD* being a diffusion spread time over the dot [35–37]. In this case, the main question is about critical interaction strength *U* or excitation energy *δ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 *Uc*, called the Åberg criterion [40], compares the typical two-body matrix elements with the energy level spacing Δ*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* < *Uc* ∼ Δ*c* and the RMT statistics sets in for *U* > *Uc* ∼ Δ*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 *δE* above the critical border *δEch* determined from the above criterion [41]:

$$
\delta E \succ \delta E\_{\text{ch}} \approx \text{g}^{2/3} \Delta\_1 \quad , \text{ g} = \Delta\_1 / \mathcal{U}. \tag{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 Δ2 ∼ Δ*c* being parametrically larger than the three-particle level spacing Δ3 Δ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 *δE* < *δEch*) to dynamically thermalized individual eigenstates at relatively strong interactions (presumably for *δE* > *δEch*). Thus, for the TBRIM with fermions, the validity of DTC for individual eigenstates at *U* > *Uc* ∼ Δ*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–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–59]. In fact, the SYK model, in its fermionic formulation [56], corresponds to the TBRIM considered with a conductance close to zero *g* → 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–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 \mathbf{x}^2 + \omega\_y^2 \mathbf{y}^2) + V\_d(\mathbf{x}, \mathbf{y}).\tag{2}$$

Here, the first two terms describe the 2D oscillator with frequencies *ωx*, *ωy* and the last term gives the potential wall of elastic disk of radius *rd*. We choose the dimensionless units with mass *m* = 1, frequencies *ωx* = 1, *ωy* = √2 and disk radius *rd* = 1. The disk center is located at (*xd*, *yd*)=(−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*, *px*]=[*y*, *py*] = *ih*¯ [48]. For the quantum problem, we use the value of the dimensionless Planck constant *h*¯ = 1 so that the ground state energy is *Eg* = 1.685. In the following, the energies are expressed in atomic like units of energy *Eu* = *h*¯ *ωx* (for our choice of Sinai oscillator parameters, we also have *Eu* = *h*¯ *ωx* = *h*¯ 2/(*mrd* 2)) [48] with the typical size of oscillator ground state being equal to the disk radius: *a*0 = Δ*xosc* = (*h*¯ /*mωx*)1/2 = *rd*.

**Figure 1.** Color plot of one-particle eigenstates *ϕk*(*x*, *y*) of the Sinai Hamiltonian in coordinate plane (*x*, *y*) with −7.6 ≤ *x* ≤ 7.6 and −5.4 ≤ *y* ≤ 5.4 for orbital numbers *k* = 1 (ground state) (**a**), *k* = 6 (**b**), *k* = 11 (**c**) and *k* = 16 (**d**). The numerical values of the color bar apply to the signed and nonlinearly rescaled wave function amplitude: sgn[*ϕk*(*x*, *y*)] |*ϕk*(*x*, *y*)/*ϕ*max| 1/2, where *ϕ*max is the maximum of |*ϕk*(*x*, *y*)| and the exponent 1/2 provides amplification of regions of small amplitude.

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*(*ε*) = *ε*2/(2 √2) − *ε*/2 [48]. Thus, the one-particle density of states *ρ*1(*ε*) and corresponding level spacing Δ1 are:

$$\rho\_1(\varepsilon) = \frac{dk}{d\varepsilon} = \frac{\varepsilon}{\sqrt{2}} - \frac{1}{2} \quad , \quad \Delta\_1 = \frac{1}{\rho} \approx \frac{\sqrt{2}}{\varepsilon} \approx \frac{0.84}{\sqrt{k}}.\tag{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 *rc* is significantly smaller than the disk radius *rd* and the typical size of the wave function, i.e., *rc rd*. 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 *rd* ∼ 1 μm = 10−4 cm, while, for Li and other alkali atoms, we have *rc* ∼ 3 × 10−7 cm [64,65]. Of course, in the limit of small *rc rd*, 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* ≤ *rc* and being zero for *r* > *rc*, where we simply choose *rc* = 0.2*rd*, which corresponds well to the short range interaction regime. The precise value of *rc* is not very important since a slight modification *rc* → *r*¯*c* can be absorbed in a modified amplitude according to *U* → *U*¯ = *U*(*r*¯*c*/*rc*)2, a relation we verified numerically for various values of *r*¯*c* < *rd*. We numerically checked that, at a used *rc* = 0.2*rd* value, we are in the regime when the interaction matrix elements are proportional to *rc* 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 *εk* and corresponding eigenstates *ϕk*(**r**) of the Sinai oscillator (2). As repulsive interaction potential *v*(**r**), we choose the short ranged box function *v*(**r**) = *U* if |**r**| ≤ *rc* = 0.2 (since *rd* = 1) and *v*(**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 ≤ *L* ≤ *M* spinless fermions takes the form:

$$H = \sum\_{k=1}^{M} \varepsilon\_k \, c\_k^\dagger c\_k + \sum\_{i$$

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

$$\mathcal{V}\_{\text{ij,kl}} = \mathcal{V}\_{\text{ij,kl}} - \mathcal{V}\_{\text{ij,kl}} \quad , \quad \mathcal{V}\_{\text{ij,kl}} = \int d\mathbf{r\_1} \int d\mathbf{r\_2} \,\phi\_i^\*(\mathbf{r\_1}) \,\phi\_j^\*(\mathbf{r\_2}) \,\upsilon(\mathbf{r\_1} - \mathbf{r\_2}) \,\eta\_k(\mathbf{r\_1}) \,\eta\_l(\mathbf{r\_2}) \tag{5}$$

and *c*† *k* , *ck* 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 *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: *V*¯ *ij*,*kl* = *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 *Vij*,*kl* according to (5) allows for also reducing the *k*, *l* sums to *k* < *l*. Furthermore, the ordering of the two fermion operators *cl ck* in (4) is also important and necessary to obtain positive expectation values if the interaction is repulsive. The anti-symmetrized matrix elements *Vij*,*kl* correspond to a *M*2 × *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 *Vij*,*ij* → *Vij*,*ij* − (1/*M*2) ∑*k*<*l Vkl*,*kl* to ensure that this matrix has a vanishing trace ( One can easily show that the trace of the *M*2 × *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(\{n\_k\}) = \sum\_{k=1}^{M} \varepsilon\_k n\_{k\prime} \tag{6}$$

where {*nk*} represents a configuration such that *nk* ∈ {0, 1} and ∑*k nk* = *L*. The associated eigenstates are the basis states where each orbital is either occupied (if *nk* = 1) or unoccupied (if *nk* = 0) and, in this work, we will denote these states in the usual occupation number representation: |*nM* ··· *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 *nk* act as quasi-random numbers) with mean and variance (see also Equation (A.4) of Ref. [46]):

$$E\_{\text{mean}} = L\overline{\varepsilon} \quad , \quad \sigma\_0^2 = \frac{L(M-L)}{M-1} \left(\overline{\varepsilon^2} - \overline{\varepsilon}^2\right) \quad , \quad \overline{\varepsilon^n} = \frac{1}{M} \sum\_{k=1}^M \varepsilon\_k^n \quad , \quad n = 1, 2. \tag{7}$$

Therefore, the many-body level spacing ΔMB or inverse Heisenberg time at the band center *E* = *E*mean is given by ΔMB = 1/*t*H = √2*π*(*σ*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 *ε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, *εn* ≈ 2 *εn M*/(*n* + 2) and *ε*2 − *ε* 2 ≈ *ε*2 *M*/18 ≈ √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 Δ*c* = √2*π* [*σ*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{\varepsilon^2} - \overline{\varepsilon}^2\right) \quad \Rightarrow \quad \frac{\sigma\_0^2(L=2)}{\sigma\_0^2} = \frac{2(M-2)}{L(M-L)}\tag{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/Δ*c* = *K*/[ √2*π σ*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 *εk* values).

According to the Åberg criterion [38,39,41], the onset of chaotic mixing happens for typical interaction matrix elements *U* comparable to Δ*c*. Therefore, we compute the quantity *V*mean = 2 |*Vij*,*kl*|2 (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*mean ≈ 2 exp(ln |*Vij*,*kl*|), which excludes this scenario. Using this quantity, we introduce the dimensionless Åberg parameter and the critical interaction amplitude *Uc* by *A* = *V*mean/Δ*c* = *U*/*Uc* such that *A* = 1 if *U* = *Uc*. We expect [38,39,41] to be the onset of strong/chaotic mixing at *A* 1 and a perturbative regime for *A* 1, while, at *A* = 1, we have the critical interaction strength *U* = *Uc*. The value of *Uc* depends on the parameters *L*, *M*, *σ*0 and the overlap of the one-particle eigenstates according to (5). To obtain some useful analytical expression of *Uc*, we note that the quantity *V*mean, numerically computed for 4 ≤ *M* ≤ 30, can be quite accurately fitted by *V*mean ≈ 3 × 10−4 *U*/*εM*. Furthermore, we remind readers of the expression Δ*c* = (1/*K*) 2 4*π*(*M* − 2)(*ε*2 − *ε* 2)/(*M* − 1), which can be simplified in the limit *M* 1 and *L* 1, such that *K* ≈ (*M* − *L*)2*L*2/4, resulting in: Δ*c* = 4/3 √2*π εM*/[(*M* − *L*)2 *L*2]. Here, we also used the found expression above *ε*2 − *ε* 2 ≈ *ε*2 *M*/18. From this, we find that *Uc* = Δ*cU*/*V*mean ≈ *C M*/[(*M* − *L*)2 *L*2] with a numerical constant *C* ≈ 16 × 104√*π*/9 ≈ 3.15 × 104, where we also used *ε*2 *M* ≈ 2 √2 *M* according to (3). Below, we will give more accurate numerical values of *V*mean, Δ*c* and *Uc* for the parameter choice of *M* and *L* numerically relevant in this work.

We note that this estimate for *A* = *U*/*Uc* applies to energies close to the many body band center of *H* and that, for energies away from the band center, the value of Δ*c* is enhanced, thus reducing the effective value of *A*. Furthermore, we computed *V*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 ≤ *A* ≤ 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*mean = 3.865 × 10−5 *U* and, from (8), that Δ*c* = √2*π* [*σ*0(*L* = 2)/*K*] = 6.1706 × 10−3, where the quantities *εn* were exactly computed from the numerical orbitals energies *εk*. From this, we find that *Uc* = Δ*c U*/*V*mean ≈ 159.65. This expression is more accurate than the more general analytical estimate for arbitrary *M* 1 and *L* 1 given in the last section (which would provide *Uc* ≈ 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*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 *σ*2(*A*) of the Gaussian density of states by a fit of the integrated density of states *P*(*E*) using

$$P(E) = (1 + \text{erf}[q(E)])/2 \quad , \quad q(E) = (E - E\_{\text{mean}})/[\sqrt{2}\,\sigma(A)].\tag{9}$$

such that *P* (*E*) is a Gaussian of width *σ*(*A*) and center *E*mean (see Appendix A of Ref. [46] for more details). From this, we find the behavior:

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

where *α* is a constant depending on *M* and *L*; for *M* = 16, *L* = 7 the fit values of *σ*0 and *α* are *σ*0 = 3.013 ± 0.009 and *α* = 0.00877 ± 0.00010. It is also possible to determine *σ*(*A*) using the expression *σ*2(*A*) = TrFock % (*H* − *E*mean**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 *Em*. This provides the same behavior as (10) with the very similar numerical values *σ*0 = 3.013 ± 0.007 and *α* = 0.00858 ± 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* ≤ 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*(*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*2 mean. In this case, one can show theoretically [46] that *σ*2(*A*) = *σ*2 0 + *K*2*V*2 mean, 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*mean = *A*Δ*c* = *A* √2*π σ*0(*L* = 2)/*K*, we indeed obtain (10) with *α* = *α*th = 4*π*(*M* − 2)*K*2/[*K*2(*L*(*M* − *L*)]. For *M* = 16, *L* = 7, we find *σ*0 = 3.0279 (see (7)) and *α*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* ≈ 2*L* 1, we have *K*2 ≈ *K* ≈ *L*4/4 and *α*th ≈ 32*π*/*L*5, which is parametrically small for very large *L*. Therefore, there is a considerable range of values 1 < *A* < 1/ √*α* 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* ≤ 10 such that the effect of the increase of the spectral width *σ*(*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 [47]. As in [45,46], we determine the occupation numbers *nk* = *c*† *k ck* for such a state, as well as the corresponding fermion entropy *S* [47] and the effective total one-particle energy *E*1*p* by :

$$S = -\sum\_{k=1}^{M} \left( n\_k \ln n\_k + (1 - n\_k) \ln(1 - n\_k) \right) \quad , \quad E\_{1p} = \sum\_{k=1}^{M} \varepsilon\_k n\_k \tag{11}$$

based on the assumption of weakly interacting fermions. In the regime of modest interaction *A* 5 (for *M* = 16, *L* = 7), corresponding to a constant spectral width *σ*(*A*) ≈ *σ*0, we have typically *E*1*p* ≈ *E*ex (for exact eigenstates of *H*) or *E*1*p* ≈ *H* (for other states). If the given state is thermalized, its occupation numbers *nk* should be close to the theoretical Fermi–Dirac filling factor *n*(*εk*) with *n*(*ε*) = 1/(1 + exp[*β*(*ε* − *μ*)]), where inverse temperature *β* = 1/*T* and chemical potential *μ* are determined by the conditions:

$$L = \sum\_{k=1}^{M} n(\varepsilon\_k) \quad , \quad E = \sum\_{k=1}^{M} \varepsilon\_k n(\varepsilon\_k) . \tag{12}$$

Here, *E* is normally given by *E*1*p*, but one may also consider the value *E*ex (or *H*) 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*th(*E*) using (11) with *nk* being replaced by *n*(*εk*) (where *β*, *μ* are determined from (12) for the energy *E*).

The many-body states with energies above *E*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 ≤ *E* ≤ 39 ≈ *E*mean (for *M* = 16, *L* = 7). In Figure 2, we compare the thermalized Fermi–Dirac occupation number *n*(*ε*) with the occupation numbers *nk* for two eigenstates at level numbers *m* = 123 (1354, with *m* = 1 corresponding to the ground state) with approximate energy eigenvalue *E* ≈ 32 (*E* ≈ 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.

**Figure 2.** Orbital occupation number *nk* versus orbital energies *εk* (black stars) of individual eigenstates at level numbers *m* = 123 (**a**,**c**,**e**), 1354 (**b**,**d**,**f**), and Åberg parameter *A* = 0.35 (**a**,**b**), *A* = 3.5 (**c**,**d**), *A* = 10 (**e**,**f**), (with *m* = 1 corresponding to the ground state). The thin blue (thick red) curves show the theoretical Fermi–Dirac occupation number *n*(*ε*) = 1/(1 + exp[*β*(*ε* − *μ*)]), where inverse temperature *β* and chemical potential *μ* are determined from (12) with *E* = *E*1p (*E* = *E*ex). The horizontal green lines correspond to the constant value 0.5 whose intersections with the red or blue curves provide the positions of the chemical potential. In this and all subsequent figures, the orbital number is *M* = 16, the number of particles is *L* = 7 and the corresponding dimension of the many body Hilbert space is *d* = 11440. Table 1 gives for each of these levels the values of *E*ex, *E*1p, *S*, *S*th, *β*, *μ* and for both energies for the latter three parameters.

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 *nk* values are rounded to 1 (or 0) if *nk* > 0.5 (*nk* < 0.5). For *m* = 1354, the values of *nk* 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*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*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*1p as reference energy to compute temperature and chemical potential. The temperature obtained from *E*ex is too small because here the increase of *σ*(*A*) is already quite strong and *E*ex rather strongly deviates from *E*1p. In addition, the value of *S*th using *E*ex does not match *S*. Obviously, at stronger interaction values, it is necessary to use *E*1p to test the thermalization hypothesis of a given state.

Figure 3 shows the mutual dependence between the three parameters *β*, *μ* on *E* when solving the conditions (12). The chemical potential as a function of *β* = 1/*T* is rather constant except for smallest values of *β* where *μ* ∼ 1/*β* with a negative prefactor. One can actually easily show from (12) that in the limit *β* → 0 the chemical potential does not depend on *εk* and is given by *μ* = − ln[1 + (*M* − 2*L*)/*L*]/*β* providing a singularity if *L* = *M*/2 with negative (positive) prefactor for *L* < *M*/2 (*L* > *M*/2) and *μ* = 0 for *L* = *M*/2. The temperature (*β*−1) vanishes for *E* close to the lower energy border and diverges for *E* close to the band center *E*mean.

**Figure 3.** Dependence of chemical potential *μ* on inverse temperature *β* = 1/*T* (**a**) and of *β* = 1/*T* on energy *E* (**b**) where *β* and *μ* are determined from (12) for a given energy *E*.

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* ≤ 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].

**Figure 4.** Histogram of unfolded level spacing statistics (blue line) for the exact energy eigenvalues *Em* of *H* (using the lower half of the spectrum with 1 ≤ *m* ≤ *d*/2). The different panels correspond to the Åberg parameter values *A* = 0.2 (**a**), *A* = 0.5 (**b**), *A* = 3.5 (**c**), *A* = 10 (**d**). The unfolding is done using the integrated density of states (9), where *q*(*E*) is replaced by a fit polynomial of degree 5. The Poisson distribution *p*Pois(*s*) = exp(−*s*) (black line) and the Wigner surmise *p*Wig(*s*) = *π* 2 *s* exp(−*π* 4 *s*2) (green line) are also shown for comparison.

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*1p or *E*ex with the theoretical thermalized entropy *S*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*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*th(*E*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*1p) are very close to the theoretical line with the average ratio *S*/*S*th(*E*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*ex and *E*1p are still rather close. For *A* = 10, all states are well thermalized (when using the energy *E*1p) while the data points for *E*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*1p (but more concentrated on energy values closer to the center with larger entropy values and larger temperatures), while, for *E*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).

**Figure 5.** Dependence of the fermion entropy *S* on the effective one-particle total energy *E*1*p* (blue cross symbols) and the exact many-body energy *E*ex (red plus symbols). The green curve shows the theoretical entropy *S*th(*E*) obtained from the Fermi–Dirac occupation numbers as explained in the text. The different panels correspond to the Åberg parameter values *A* = 0.2 (**a**), *A* = 0.5 (**b**), *A* = 3.5 (**c**), *A* = 10 (**d**).

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 *nk* (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*(*ε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.

**Figure 6.** Color density plot of the orbital occupation number *nk* in the plane of energy *E* and orbital index *k*. (**a**) *nk* values of exact eigenstates of *H* with Åberg parameter *A* = 3.5; (**b**) thermalized Fermi–Dirac occupation number *n*(*εk*), where *β* and *μ* are determined from (12) as a function of total energy *E*. The occupation number *nk* is averaged over all eigenstates (**a**) or several representative values of *E* (**b**) inside a given energy cell. The energy interval 29 ≤ *E* ≤ 39 corresponds roughly to the lower half of the spectrum (at *M* = 16, *L* = 7) for states with positive temperature and is similar to the energy interval used in Figures 4 and 5. The color bar provides the translation between *nk* values and colors (red for maximum *nk* = 1, green for *nk* = 0.5 and blue for minimum *nk* = 0).

### *3.6. Thermalization of Quantum Time Evolution*

The question arises how or if a time dependent state |*ψ*(*t*)>= exp(−*iHt*)|*ψ*(0)>, obeying the quantum time evolution with the Hamiltonian *H* and an initial state |*ψ*(0)> being a non-interacting eigenstate |*nM* ··· *n*2*n*1)> (with all *nk* ∈ {0, 1} and ∑*k nk* = *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) |*φ*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) |*φ*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) |*φ*3>= |0000011000110111> and (iv) |*φ*4>= |1000100011001011>. The states |*φ*3> and |*φ*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 *nk* to 1 (or 0) if *nk* > 0.5 (*nk* < 0.5) (states of top panels in Figure 2). The approximate energies (6) of these four states are *E* ≈ 30 (|*φ*1>), *E* ≈ 32 (|*φ*2> and |*φ*3>) and *E* ≈ 35 (|*φ*4>).

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

$$
\Delta t = \frac{t\_{\rm H}}{d} = \frac{1}{\sqrt{2\pi \,\sigma^2(A)}},
\tag{13}
$$

where *t*H is the Heisenberg time (at the given value of *A*), *d* the dimension of the Hilbert space and *σ*(*A*) the width of the Gaussian density of states given in (10). The quantity Δ*t* is the shortest physical time scale of the system (inverse of the largest energy scale) and obviously for *t* Δ*t* the unitary evolution operator is close to the unit matrix multiplied by a uniform phase factor: exp(−*iHt*) ≈ exp(−*iE*mean*t*) **1** since the eigenvalues *Em* of *H* satisfy |*Em* − *E*mean| *σ*(*A*). We expect that any signification deviation of |*ψ*(*t*)> with respect to the initial condition |*ψ*(0)> happens at *t* ≥ Δ*t* (or later in case of very weak interaction). Furthermore, by analyzing the time evolution in terms of the ratio *t*/Δ*t* the results do not depend on the global energy scale of the spectral width. The longest time scale is the Heisenberg time *t*H ≈ 104Δ*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 *nk*(*t*) =<*ψ*(*t*)|*c*† *k ck* |*ψ*(*t*)> and present them in a color density plot in the plane (*k*, *t*/Δ*t*). In addition, at the time value used last, we compute the effective total one-particle energy *E*1p using the relation (11) (note that *E*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 *β*, chemical potential *μ* and the thermalized Fermi–Dirac filling factor *n*(*εk*) 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 *nk* values). The latter is shown in Figure 6 and applies also to all subsequent figures with color density plots for *nk* 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 |*φ*1> and the two Åberg parameter values *A* = 1 and *A* = 3.5 using a linear time scale with integer multiples of Δ*t* and for *t* ≤ 2000 Δ*t* ≈ *t*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 Δ*t* for *t* < 1000 Δ*t*, and a modest decay for *t* > 1000 Δ*t*. At the same time, the first orbitals above the Fermi sea are slightly excited. At final *t* = 2000 Δ*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 *nk* with 8 ≤ *k* ≤ 11, there are later and more modest excitations with a periodic time structure. Here, the final state at *t* = 2000 Δ*t* is also not thermalized, but it is closer to thermalization as for the case *A* = 1.

**Figure 7.** Color density plot of the orbital occupation number *nk* in the plane of orbital index *k* and time *t* for the time dependent state |*ψ*(*t*)>= exp(−*iHt*)|*ψ*(0)> with initial condition |*ψ*(0)>= |*φ*1>= |0000100000111111>. The time values are integer multiples of the elementary quantum time step Δ*t* = *t*H/*d* = 1/[ √2*π σ*(*A*)] where *t*H is the Heisenberg time (at the given value of *A*). The bar behind the vertical white line with the label "th" shows the theoretical thermalized Fermi–Dirac occupation numbers *n*(*εk*) where *β* and *μ* are determined from (12) using the energy *E* = *E*1p of the state |*ψ*(*t*)> at the last time value *t* = 2000 Δ*t*. The two panels correspond to the Åberg parameter *A* = 1 (**a**), *A* = 3.5 (**b**). For the translation of colors to *nk* values, the color bar of Figure 6 applies.

The linear time scale used in Figure 7 is not very convenient since it cannot capture a rapid decay/increase of *nk* well 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 Δ*t* ≤ *t* ≤ 106 Δ*t* ≈ 102 *t*H. Note that, in these figures, the different *nk* 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 *nk*, there will be, for larger time values, a quasi random selection of different time positions with respect to the period.

**Figure 8.** Color density plot of the orbital occupation number *nk* in the plane of orbital index *k* and time *t* for the time dependent state |*ψ*(*t*) >= exp(−*iHt*)|*ψ*(0) >. The time axis is shown in logarithmic scale with time values *tn* = 10 (*n*/100)−1 Δ*t* and integer *n* ∈ {0, 1, ... , 700} corresponding to 0.1 ≤ *tn*/Δ*t* ≤ 106. The elementary quantum time step Δ*t* is the same as in Figure 7. The bar behind the vertical white line with the label "th" shows the theoretical thermalized Fermi–Dirac occupation numbers *n*(*εk*) where *β* and *μ* are determined from (12) using the energy *E* = *E*1p of the state |*ψ*(*t*)> at the last time value *t* = 106 Δ*t*. The additional longer tick below the *t*-axis right next to the tick for 103 gives the position of the maximal time value *t*/Δ*t* = 2000 of Figure 7. The different panels correspond to the initial state |*ψ*(0)>= |*φ*1>= |0000100000111111> (**a**,**c**,**e**) or |*ψ*(0)>= |*φ*2>= |0010100000011111> (**b**,**d**,**f**) and Åberg parameter values *A* = 1 (**a**,**b**), *A* = 3.5 (**c**,**d**), *A* = 10 (**e**,**f**). For the translation of colors to *nk* values, the color bar of Figure 6 applies.

In Figure 8, the time evolution for the initial states |*φ*1> and |*φ*2> is shown for the Åberg values *A* = 1, 3.5, 10. For |*φ*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 Δ*t* up to 106 Δ*t* and for *A* = 3.5 the initial decay of *n*6 and *n*12 happens at *t* ≈ 10 Δ*t*. For |*φ*1> at *A* = 10, the decay starts at *t* ≈ 3 Δ*t* and an approximate thermalization happens at *t* > 40 Δ*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 |*φ*2> at *A* = 1, the decay of excited orbitals 12 and 14 starts at *t* ≈ 100 Δ*t* and saturates at *t* ≈ 1000 Δ*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 Δ*t*. The final state at *t* = 106 Δ*t* is not thermalized even though some orbitals have *nk* values close to thermalization. For |*φ*2> at *A* = 3.5, the decay of excited orbitals 12 and 14 starts at *t* ≈ 10 Δ*t* and, for *t* > 300 Δ*t*, there is thermalization (but requiring some time average as for |*φ*1> at *A* = 10). Interestingly, at intermediate times 10 Δ*t* < *t* < 100 Δ*t*, the high orbitals 13 and 16 are temporarily slightly excited and decay afterwards rather quickly to their thermalized values. For |*φ*2> at *A* = 10, the decay of excited orbitals 12 and 14 starts even at *t* ≈ 3 Δ*t* and thermalization seems to set in at *t* > 30 Δ*t*.

**Figure 9.** As in Figure 8 but for the initial states |*ψ*(0)>= |*φ*3>= |0000011000110111> (**a**,**c**,**e**) and |*ψ*(0)>= |*φ*4>= |1000100011001011> (**b**,**d**,**f**) (with the same *A* values as in Figure 8 for each row). These initial states can be obtained from the eigenstates of *H* for *A* = 0.35 at level numbers *m* = 123 or 1354, respectively, by rounding the occupation numbers to 1 (or 0) if *nk* > 0.5 (*nk* < 0.5) (see also top panels of Figure 2).

Figure 9 is similar to Figure 8 except for the initial states |*φ*3> and |*φ*4> which have occupation numbers *nk* ∈ {0, 1} obtained by rounding the *nk* 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* ≈ 300 Δ*t* (*t* ≈ (10 − 20) Δ*t* or *t* ≈ (2 − 3) Δ*t*) for *A* = 1 (*A* = 3.5 or *A* = 10, respectively). There is no thermalization for both states at *A* = 1 (but some *nk* values are close to thermalized values), approximate thermalization for *A* = 3.5 and |*φ*3> and good thermalization for *A* = 3.5 and |*φ*4> as well as *A* = 10 (both states).

Using the time dependent values *nk*(*t*), one can immediately determine the corresponding entropy *S*(*t*) using (11). At *t* = 0, we have obviously *S*(0) = 0, since, for all four initially considered states, we have perfect occupation number values of either *nk* = 0 or *nk* = 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*th(*E*1p) (with *E*1p determined for the state |*ψ*(*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 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*(*t*) for the two initial states |*φ*1>, |*φ*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*th(*E*1p). We note that *E*1p is not conserved at strong interactions and that its initial value *E*1p ≈ 30 (*E*1p ≈ 35) at *t* = 0 evolves to *E*1p ≈ 33.5 (*E*1p ≈ 37) at large times for |*φ*1> (|*φ*4>) corresponding roughly to *S* ≈ *S*th(*E*1p) ≈ 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*1p values. For *A* = 3.5 and |*φ*4>, there is also saturation of *S* to its thermalized value. For *A* = 3.5 and |*φ*1>, there seems to be an approximate saturation at a quite low value *S* ≈ 6 but with periodic fluctuations in the range 6 ± 0.3. For *A* = 1 and |*φ*4>, there is a quite late and approximate saturation with some fluctuations that are visible for *t* > 104 Δ*t* and with *S* ≈ 10 ± 0.2. For *A* = 1 and |*φ*1>, there is a late periodic regime for *t* > 103 Δ with a quite large amplitude *S* ≈ 3 ± 1 and with *S*max ≈ 4 significantly below the thermalized entropy *S*th(*E*1p) ≈ 5.5. The panels using a normal (instead of logarithmic) time scale with *t* ≤ 200 Δ*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*.

**Figure 10.** Time dependence of the entropy *S*, computed by (11), of the state |*ψ*(*t*) >= exp(−*iHt*)|*ψ*(0) > for the Åberg parameter values *A* = 1 (red lines), *A* = 3.5 (green lines), *A* = 10 (blue lines) and initial states |*ψ*(0)>= |*φ*1>= |0000100000111111> (**a**,**c**); |*ψ*(0)>= |*φ*4>= |1000100011001011 > (**b**,**d**); thick colored lines show numerical data of *S*(*t*) and thin horizontal colored lines show the thermalized entropy *S*th(*E*1p) with *E*1p being determined from |*ψ*(*t*) > at *t* = 106 Δ*t*; panels (**a**,**b**) use a linear time axis: 0 ≤ *t* ≤ 200 Δ*t*; panels (**c**,**d**) use a logarithmic time axis: 0.1 Δ*t* ≤ *t* ≤ 106 Δ*t*; Δ*t* is the elementary quantum time step (see also Figure 6).

The periodic (or quasi-periodic) time dependence of *nk*(*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, ...) of exact eigenstates of *H* contribute mostly in the expansion of |*ψ*(*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* Δ*t* but which immediately becomes linear for *t* - Δ*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 Δ*t* in the definition of Γ*c*, i.e., Γ*c* is defined as the ratio of the initial slope *S* (*t*) over the global spectral bandwidth ∼ *σ*(*A*) ∼ 1/Δ*t*) the quantity Γ*c* = *dS*(*t*)/*d*(*t*/Δ*t*) = Δ*tS* (*t*) for *t* = *ξ*Δ*t*, where *ξ* - 1 is a numerical constant of order one. To determine Γ*c* practically, we perform first the fit *S*(*t*) = *S*¯ ∞ (1 − exp[−*γ*¯1(*t*/Δ*t*)]) for 0 ≤ *t*/Δ*t* ≤ 100 and use the exponential decay rate *γ*¯1 to perform a refined fit *S*(*t*) = *S*∞ (1 − exp[−*γ*1(*t*/Δ*t*) − *γ*2(*t*/Δ*t*)2]) for the interval 0 ≤ *t*/Δ*t* ≤ 5/*γ*¯1. From this, we determine Γ*c* = *S*∞ *γ*1, which is rather close to *S*¯ ∞ *γ*¯1 for *A* ≤ 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 *ξ* between 1 and roughly 1/*γ*1,

which is indeed useful to smear out some oscillations in the initial increase of *S*(*t*) for smaller values *A* ≤ 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 Γ*c* on the parameter *A* for our four initial states. At first sight, one observes a behavior Γ*c* ∝ *A*2 for *A* 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 Γ*c* ∝ *Ap* for *A* ≤ 2 provide exponents close to *p* ≈ 1.75 − 1.85) and it turns out that these deviations correspond to a logarithmic correction: Γ*c* = *f*(*A*)=(*C*1 − *C*2 ln[*g*(*A*)]) *g*(*A*) with *g*(*A*) = *A*2 (for fits with *A* ≤ 1) or with *g*(*A*) = *A*2/(1 + *C*3 *A*2) (for fits with all *A* values).

**Figure 11.** Dependence of the initial slope Γ*c* = Δ*t S* (*t*) (at small time values *t* ∼ Δ*t*) of the time dependent entropy on the Åberg parameter *A* in double logarithmic scale. Δ*t* is the elementary quantum time step (see also Figure 6). The practical determination of Γ*c* is done using the fit *S*(*t*) = *S*∞ (1 − exp[−*γ*1(*t*/Δ*t*) − *γ*2(*t*/Δ*t*)2]) which provides Γ*c* = *S*∞ *γ*1. The different data points correspond to the four different initial states used in Figures 8 and 9. The dashed line corresponds to the power law behavior ∝ *A*2 and the light blue line corresponds to the fit Γ*c* = *f*(*A*)=(*C*1 − *C*2 ln[*g*(*A*)]) *g*(*A*) with *g*(*A*) = *A*2/(1 + *C*3*A*2) and fit values *C*1 = 0.107 ± 0.009, *C*2 = 0.0081 ± 0.0023, *C*3 = 0.092 ± 0.017 for the initial state |*ψ*(0)>= |*φ*1>= |0000100000111111> corresponding to the red plus symbols. Fit values for the other initial states can be found in Table 2. The full black line corresponds to *f*0(*A*) = *f*(*A*)*C*3=0 = [*C*1 − *C*2 ln(*A*2)] *A*2. The simpler fit Γ*c* = *f*0(*A*) in the range 0.025 ≤ *A* ≤ 1 provides the values *C*1 = 0.107 ± 0.002 and *C*2 = 0.0078 ± 0.0005 which are identical (within error bars) to the values found by the more general fit Γ*c* = *f*(*A*) for the full range of *A* values.

**Table 2.** Values of the fit parameters *C*1, ... , *C*5 for the initial states |*φ*1 >, ... , |*φ*4 > used for the analytical fits of Γ*c* (Γ*F*) in Figure 11 (see also Figure 12 and Figure 13).

To understand this behavior, we write for sufficiently small times *nk*(*t*) ≈ 1 − *δnk*(*t*) (if *nk*(0) = 1) or *nk*(*t*) ≈ *δnk*(*t*) (if *nk*(0) = 0), where *δnk*(*t*) is the small modification of *nk*(*t*). Time dependent perturbation theory suggests that *δnk*(*t*) ∼ (*t*/Δ*t*)2 for *t* Δ*t* and *δnk*(*t*) ≈ *ak A*2 *t*/Δ*t* for *t* - Δ*t* such that still *δnk*(*t*) 1 with coefficients *ak* 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

*δn*2 *k*, we obtain: *S* ≈ − ∑*k* (*δnk* ln *δnk* − *δnk*) and Γ*c* = Δ*t S* (*t*) ≈ −Δ*t* ∑*k δn k*(*t*)ln *δnk*(*t*) with *t* = *ξ* Δ*t*. Since *δn k*(*t*) ≈ *akA*2/Δ*t*, we find indeed the behavior :

$$\Gamma\_{\mathfrak{c}} = \left[\mathbb{C}\_1 - \mathbb{C}\_2 \ln(A^2)\right] A^2 \quad , \quad \mathbb{C}\_1 = -\sum\_k a\_k \langle \ln(a\_k \xi) \rangle \quad , \quad \mathbb{C}\_2 = \sum\_k a\_{k'} \tag{14}$$

where · · · indicates an average over some modest values of *ξ* - 1. The precise values of *ak* 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 → *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 Γ*c* = *f*(*a*) for the full range of *A* values.

### *3.7. Survival Probability and Fermi's Golden Rule*

The knowledge of the time dependent states |*ψ*(*t*)> allows us also to compute the decay function *p*dec(*t*) = | <*ψ*(0)|*ψ*(*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* Δ*t*, we expect a quadratic decay: 1 − *p*dec(*t*) ≈ (*H* − *E*mean)2*t* 2 ≈ const.(*t*/Δ*t*)2 with · · · being the quantum expectation value with respect to |*ψ*(0)> and a numerical constant 1 since 1/Δ*t* represents roughly the spectral width of *H*. For *t* - Δ*t*, but such that 1 − *p*dec(*t*) 1, we have, according to Fermi's golden rule: 1 − *p*dec(*t*) = Γ*F* (*t*/Δ*t*), where Γ*F* is the decay rate (Again, for practical reasons and similarly to Γ*c*, we incorporate in the definition of Γ*F* the time scale Δ*t*, i.e., Γ*F* = Δ*t* × usual decay rate found in the literature and meaning that Γ*F* is defined as the ratio of the usual decay rate over the global spectral bandwidth) of the state.

To determine Γ*F* numerically, we apply the fit: *p*dec(*t*) = *C* exp(−ΓF *t*/Δ*t*) in two steps. First, we use the interval 1 ≤ *t*/Δ*t* ≤ 50 and, if 5/Γ*F* < 50, corresponding to a rapid decay (which happens for larger values of *A*), we repeat the fit for the reduced interval 1 ≤ *t*/Δ*t* ≤ 5/Γ*F*. The choice of the Amplitude *C* = 1 and the condition *t* ≥ Δ*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 |*φ*1> and the Åberg values *A* = 3.5 and *A* = 10. In both cases, the shown maximal time value *t*max = 50 Δ*t* (if *A* = 3.5) or *t*max ≈ 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 ≤ *t*/Δ*t* ≤ 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 ≤ *t*/Δ*t* ≤ 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*dec(*t*) does not decay completely at large times and which typically happens in the perturbative regime *A* 1.

**Figure 12.** Decay function *p*dec(*t*) = | <*ψ*(0)|*ψ*(*t*)> | 2 obtained numerically from |*ψ*(*t*)> with the initial state |*ψ*(0)>= |0000100000111111> (thin red line) and the fit *p*dec(*t*) = *C* exp(−ΓF *t*/Δ*t*) (thick green line) for the two Åberg values *A* = 3.5 (**a**) and *A* = 10 (**b**). The fit values are *C* = 0.959 ± 0.011, ΓF = 0.0483 ± 0.0015; (**a**)*C* = 1.339 ± 0.015, ΓF = 0.369 ± 0.005; (**b**) corresponding to the decay times Γ−1 F = 20.7 (**a**), 2.71; (**b**). Δ*t* is the elementary quantum time step (see also Figure 6).

Figure 13 shows the dependence of Γ*F* on *A* for the usual four initial states together with the fit Γ*F* = *f*(*A*) = *C*4*A*2/(1 + *C*5*A*2) for the data with initial state |*φ*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 Γ*F* ∝ *A*2 is highly accurate (with no logarithmic correction). Similarly to Γ*c*, there is only a slight dependence of the values of Γ*F* and the fit values on the choice of initial state.

**Figure 13.** Dependence of the decay rate Γ*F* corresponding to Fermi's gold rule on the Åberg parameter *A* in double logarithmic scale. The practical determination of Γ*F* is done using the exponential fit function of Figure 12 for the numerically computed decay function *p*dec(*t*). The different data points correspond to the four different initial states used in Figures 8 and 9. The dashed line corresponds to the power law behavior ∝ *A*2 and the full black line corresponds to the fit Γ*F* = *f*(*A*) = *C*4*A*2/(1 + *C*5*A*2) and fit values *C*4 = 0.0048 ± 0.0001, *C*5 = 0.0054 ± 0.0003 for the initial state |*ψ*(0) >= |*φ*1 >= |0000100000111111> corresponding to the red plus symbols. Fit values for the other initial states can be found in Table 2. Δ*t* is the elementary quantum time step (see also Figure 6).

Theoretically, we expect according to Fermi's Golden rule that: Γ*F* ≈ (Δ*t*) 2*πV*2 Fock *ρc*(*E*), where *V*2 Fock = TrFock(*V*2)/(*Kd*) = *σ*2 0 *αA*2/*K* according to the discussion below (10) and *ρc*(*E*) is the effective two-body density of states for states directly coupled by the interaction such that *ρc*(*E*mean) = 1/Δ*c* (see discussion below (7)). We note that *V*Fock is the typical interaction matrix element in Fock space which is slightly larger than *V*mean (see the theoretical discussion above for the computation of the coefficient *α* used in (10) and Appendix A of [46]). The factor Δ*t* is due to our particular definition (Again, for practical reasons and similarly to Γ*c*, we incorporate in the definition of Γ*F* the time scale Δ*t*, i.e., Γ*F* = Δ*t* × usual decay rate found in the literature and meaning that Γ*F* is defined as the ratio of the usual decay rate over the global spectral bandwidth) of decay rates. The expression of Γ*F* is actually also valid for larger values of *A* provided we use the density of states *ρc*(*E*) in the presence of interactions that provides an additional factor 1/ √ 1 + *αA*2 according to (10). Therefore, at the band center, we have: 2*π* Δ*t ρc* = *K*/[*σ*0(*L*) *σ*0(*L* = 2)(1 + *αA*2)], which gives, together with (8):

$$
\Gamma\_F = \frac{\sigma\_0(L)}{\sigma\_0(L=2)} \left( \frac{aA^2}{1 + aA^2} \right) = \sqrt{\frac{L(M-L)}{2(M-2)}} \left( \frac{aA^2}{1 + aA^2} \right). \tag{15}
$$

For *M* = 16 and *L* = 7, the square root factor is 1.5 and we have to compare 1.5*α* ≈ 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 *α* ≈ 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 Δ*t* being 1/Γ*F* ≈ 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 Figures 8 and 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(\mathbf{x}, y, t) = \varsigma \psi(t) |\Psi^{\dagger}(\mathbf{x}, y)\,\Psi(\mathbf{x}, y)|\psi(t) > \quad , \quad \Psi^{(\dagger)}(\mathbf{x}, y) = \sum\_{k} \varrho\_{k}^{(\ast)}(\mathbf{x}, y) \,\mathfrak{c}\_{k}^{(\dagger)} \,\tag{16}$$

where *ϕk*(*x*, *y*) is the one-particle eigenstate of orbital *k*, with some examples shown in Figure 1. Here, Ψ(†)(*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 *ρ*(*x*, *y*, *t*), only the values *k* ≤ *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(\mathbf{x}, \mathbf{y}; t) = \sum\_{k,l=1}^{M} \varphi\_k(\mathbf{x}, \mathbf{y}) \, \varphi\_l(\mathbf{x}, \mathbf{y}) \, n\_{kl}(t) \quad , \quad n\_{kl}(t) = \varsigma \psi(t) |c\_k^\dagger c\_l| \, \psi(t) > \tag{17}$$

where *nkl*(*t*) is the density matrix in orbital representation generalizing the occupation numbers *nk*(*t*) which are its diagonal elements. Due to the complex phases of |*ψ*(*t*)> (when expanded in the usual basis of non-interacting many-body states), the density matrix is complex valued but hermitian: *n*∗ *kl*(*t*) = *nlk*(*t*). Therefore, its anti-symmetric imaginary part does not contribute in *ρ*(*x*, *y*; *t*). We have numerically evaluated (17) and we present in Figure 14 color plots of the density matrix and the spatial density *ρ*(*x*, *y*; *t*) for *A* = 3.5, the initial state |*ψ*(0)>= |*φ*2> and four time values *t*/Δ*t* = 0.1, 30, 100, 1000. Since the density *ρ*(*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 Δ*ρ*(*x*, *y*; *t*) = *ρ*(*x*, *y*; *t*) − *ρ*(*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*/Δ*t* = 0.1, density matrix and spatial density are essentially identical to the initial condition at *t* = 0. For Δ*ρ*, 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*/Δ*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*/Δ*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 Δ*ρ* is roughly close to the case at *t*/Δ*t* = 30 but with some significant differences.

**Figure 14.** Time dependent density matrix |*nkl*(*t*)| (**a**), spatial density *ρ*(*x*, *y*, *t*) (**b**), spatial density difference with respect to the initial condition Δ*ρ*(*x*, *y*, *t*) = *ρ*(*x*, *y*, *t*) − *ρ*(*x*, *y*, 0) (**c**) all computed from |*ψ*(*t*)> for the initial state |*ψ*(0)>= |*φ*2 >= |0010100000011111> and the Åberg parameter *A* = 3.5. Panels in column (**a**) correspond to (*k*, *l*) plane with *k*, *l* ∈ {1, ... , 16} being orbital index numbers. Panels in columns (**b**,**c**) correspond to the same rectangular domain in (*x*, *y*) plane as in Figure 1. The five rows of panels correspond to the time values *t*/Δ*t* = 0.1, 30, 100, 1000 and the thermalized case (label "th") where the density matrix is diagonal with entries being the thermalized occupations numbers *nkk* = *n*(*εk*) at energy *E* = 32.9 (typical total one-particle energy of |*ψ*(*t*)> for *t*/Δ*t* ≥ 1000). The numerical values of the color bar represent values of |*nkl*| (**a**), (*ρ*/*ρ*max)1/2 (**b**), sgn(Δ*ρ*)(|Δ*ρ*|/Δ*ρ*max)1/2 (**c**) where *ρ*max or Δ*ρ*max are maximal values of *ρ* or |Δ*ρ*|, respectively.

Finally, at *t*/Δ*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 Δ*ρ* has a similar but different structure as for *t*/Δ*t* = 100 or *t*/Δ*t* = 30.

Apparently, at intermediate times 20 ≤ *t*/Δ*t* ≤ 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 × 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 *ρ*(*x*, *y*; *t*) ≈ ∑*k ϕ*2 *k* (*x*, *y*) *nk*(*t*) in terms of the time dependent occupation numbers. Only for intermediate time scales with more visible quantum coherence (more off-diagonal elements *nkl*(*t*) = 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 Δ*ρ* 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 *ρ* when *t* is slightly changed by a few multiples of Δ*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 *nkk* = *n*(*εk*) at energy *E* = 32.9, which is the typical total one-particle average energy of |*ψ*(*t*)> for the long time limit *t*/Δ*t* ≥ 1000 showing that, at *t*/Δ*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\_{\text{corr}}(\mathbf{x}, y; \mathbf{x}\_0, y\_0; t) = \epsilon \varphi(t) |\Psi^\dagger(\mathbf{x}, y)\Psi(\mathbf{x}\_0, y\_0)| \\ \psi(t) > = \sum\_{k,l=1}^{M} \varphi\_k(\mathbf{x}, y) \, \varrho\_l(\mathbf{x}\_0, y\_0) \, n\_{kl}(t) \tag{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 *ϕ*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 *nkl*(*t*) = 0 for *k* = *l*. In Figure 15, we present density plots of absolute value, real and imaginary part of *ρ*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 *ρ*corr,th(*x*, *y*; *x*0, *y*0) vanishes (giving a blue panel due to zero values).

**Figure 15.** Time dependent spatial density correlator *ρ*corr(*x*, *y*; *x*0, *y*0; *t*) shown in the same rectangular domain in (*x*, *y*) plane as in Figure 1. The columns correspond to absolute value (**a**), real part (**b**), imaginary part (**c**). The initial point is given by (*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 *ϕ*1(*x*, *y*) visible in panel (**a**) of Figure 1. Initial state, Åberg parameter and meaning of row labels are as in Figure 14. The numerical values of the color bar represent values of sgn(*u*)|*u*/*u*max| 1/2 where *u* is absolute value (**a**), real part (**b**), imaginary part (**c**), of *ρ*corr and *u*max is the maximal value of |*u*|. The data for thermalized case and imaginary part is completely zero (blue panel in bottom right corner) since the spatial density correlator is for the thermalized case purely real.

There are significant time dependent fluctuations of *ρ*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*/Δ*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 *ω* ≈ *ωx*,*y*,*z* ≈ 2*π*10*Hz* with *a*0 = 7*h*¯ /(*mω*) = 6.5*μm* and oscillator level spacing *Eu* = *h*¯ *ω* ≈ 0.5*nK* (nanoKelvin). The typical scattering length is *as* ≈ 3*nm* being small compared to *a*0. The atomic density is *ρ*0 = 1/(*a*0)3 ≈ 4 × 109*cm*−3. Since *as a*0, the two-body interaction is of *δ*-function type with *v*(**r1** − **r2**)=(4*πh*¯ 2*as*/*m*)*δ*(**r1** − **r2**) [66,72]. In our numerical simulations, *δ*(**r**) is replaced by a function *H*(*rc* − |**r**|)/(*C* 2*drd c* ) with a small *rc*, volume *C* of the unit sphere in *d* dimensions and with the Heaviside function *H*(*x*) = 1 (or 0) for *x* ≥ 0 and *H* = 0 for *x* < 0. Hence, our parameter *U* introduced in Section 3 corresponds to *U* = (*C* 2*drd c* )(4*πh*¯ 2*as*/*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 *Us* = 4*πh*¯ 2*ρ*0*as*/*m* = 4*π*(*as*/*a*0)*h*¯ *ω* [66,72] so that *Us*/¯*hω* ∼ 6 × 10−3. Compared to sodium, the mass of Li atoms is approximately three times smaller so that, for the same *ω*, we have *a*0 ≈ 10 μm and *Us*/¯*hω* ≈ 4 × 10−3. We think that the scattering length can be significantly increased via the Feshbach resonance, allowing for reaching effective interaction values *Us*/¯*hω* ∼ 1 being similar to the value *A* ∼ 3 used in our numerical studies with the onset of dynamical thermalization.

Usually, a 3D trap with fermionic atoms can capture about *Na* ∼ 105 atoms with *ω* ≈ *ωx* ∼ *ωy* ∼ *ωz* ∼ 2*π*10*Hz*. Following the result (1), it is interesting to determine the DTC border dependence on *Na* 1 for a Sinai oscillator with *rd* ∼ 1*μm* ∼ *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 - ≤ *Na* components (e.g., in the basis of oscillator eigenfunctions). The Fermi energy of the trap is then *EF* = *h*¯(*Naωxωyωz*)1/3 ≈ *h*¯ *ωNa* 1/3 [52,53]. Assuming that all these components have random amplitudes of a typical size 1/ √-, we then obtain an estimate for a typical matrix element of two-body interaction between one-particle eigenstates

$$
\hbar \mathcal{U}\_2 \approx a\_s \hbar \omega / \ell^{3/2}, \quad a\_s = 4\pi (a\_s / a\_0), \quad a\_0 = \sqrt{\hbar / m \omega}.\tag{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 *EF*, we have the one-particle level spacing Δ1 = *dEF*/*dNa* ≈ *h*¯ *ω*/(3*Na* 2/3). Hence, the effective conductance appears in in (1) is *g* = Δ1/*U*2 ≈ -3/2/(3*αsNa* 2/3). Thus, from (1), we obtain the dynamical thermalization border for excitation energy *δE* in a 3D Sinai oscillator trap with *Na* fermionic atoms:

$$
\delta E > \delta E\_{\rm ch} \approx \Lambda\_1 \mathbf{g}^{2/3} \approx 2\ell \Lambda\_1 / (\mathbf{u}\_s^{2/3} \mathbf{N}\_d^{4/9}) \sim N\_d^{-5/9} \Lambda\_1 / \mathbf{u}\_s^{2/3} \sim \hbar \omega / (\mathbf{u}\_s^{2/3} \mathbf{N}\_d^{1/9}) \sim E\_F / (\mathbf{u}\_s^{2/3} \mathbf{N}\_d^{4/9}).\tag{20}
$$

It is assumed that *δE EF*. Here, the last three relations are written in an assumption that - ∼ *Na*. Thus, at large *Na* values and not too small *αs*, the critical energy border *δEch* for dynamical thermalization is rather small compared to *EF*. However, still *δEch* Δ1. Here, we used the maximal value for the number of components - ∼ *Na*. It is possible that, in reality, can be significantly smaller than *Na*. However, the determination of the dependence -(*Na*) 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 *EF* (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* → −*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.

**Author Contributions:** All authors equally contributed to all stages of this work.

**Funding:** 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.

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

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Appendix A. Two-Body Matrix Elements near the Fermi Energy**

In this Appendix, we present numerical results for the dependence of the quantity *V*mean(*ε*) = 2 *V*2 *ij*,*kl* as a function of *ε* where the average is done *only for orbitals with energies εn close to ε* (for *n* ∈ {*i*, *j*, *k*, *l*}), i.e.: |*εn* − *ε*| ≤ Δ*ε* with Δ*ε* = 2. We note that this is different from the quantity *V*mean used in Section 3 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 *ε* (which will be identified with the Fermi energy *EF*) is that these transitions are dominant in the presence of the Pauli blockade near the Fermi level.

We remind readers that, according to the discussion of Sections 2 and 3, the matrix elements *Vij*,*kl* were computed for an interaction potential of amplitude *U* for |**r1** − **r2**| < *rc* (with the radius *rc* = 0.2*rd* = 0.2) and being zero for |**r1** − **r2**| ≥ *rc*. Furthermore, they have been anti-symmetrized and a diagonal shift *Vij*,*ij* → *Vij*,*ij* − (1/*M*2) ∑*k*<*l Vkl*,*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 *ε* close to the maximal orbital energy *εM*). Due to this, we considered two values of *M* = 30 and *M* = 60.

The numerically obtained dependence is shown in Figure A1 and is well described by the fit *V*mean/*U* = *a*/*εb* with *a* = 1.56 × 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 *rc* = 0.2*rd*, which leads to a decrease of the effective interaction strength being proportional to *rc* 2.

**Figure A1.** Dependence of the average two-body matrix element *V*mean rescaled by the amplitude of interaction strength *U* on one-particle energy *ε* for two-body interaction transitions in a vicinity of Fermi energy *EF* = *ε*; green symbols are for number of one-particle orbitals *M* = 30, red symbols are for *M* = 60; the blue line shows the fit *V*mean/*U* = *a*/*εb* with *a* = 0.000156 ± 2 × 10−6; *b* = 0.781 ± 0.005.

We note that the Fermi energy is determined by the number of fermionic atoms *Na* inside the 2D Sinai oscillators with *ε* = *EF* ≈ *ωNa* 1/2 assuming *ω* = *ωx* ≈ *ωy*. Therefore, we have *ε* ∝ *M*1/2 ∼ *Na* 1/2, Δ1 ∼ *h*¯ *ω*/*Na* 1/2 and *V*mean ∼ *αsh*¯ *ω*/-3/2 ∼ *αsh*¯ *ω*/*Na b*/2 (see (19)). Hence, the obtained exponent *b* ≈ 0.78 corresponds to the number of one-particle components - ∼ *Na b*/3 ∼ *Na* 0.25 ∼ *nx* 0.5. At the moment, we do not have a clear explanation for this numerical dependence. This dependence corresponds to *g* = Δ1/*V*mean ∼ -3/2/(*αsNa* 1/2) ∼ 1/(*αsNa* 1/8). For such a dependence, we obtain that the DTC border in 2D takes place for an excitation energy *δE* > *δEch* ∼ *g*2/3Δ1 ∼ *h*¯ *ω*/(*αs* 2/3*Na* 7/12). Thus, the thermalization can take place at rather low energy excitations above the Fermi energy with Δ1 < *δE EF*.

### **References**

c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Quantum-Heat Fluctuation Relations in Three-Level Systems Under Projective Measurements**

**Guido Giachetti 1,2,\*, Stefano Gherardini 3, Andrea Trombettoni 4,5 and Stefano Ruffo 1,2,6**

Received: 20 January 2020; Accepted: 27 February 2020; Published: 13 March 2020

**Abstract:** We study the statistics of energy fluctuations in a three-level quantum system subject to a sequence of projective quantum measurements. We check that, as expected, the quantum Jarzynski equality holds provided that the initial state is thermal. The latter condition is trivially satisfied for two-level systems, while this is generally no longer true for *N*-level systems, with *N* > 2. Focusing on three-level systems, we discuss the occurrence of a unique energy scale factor *β*eff that formally plays the role of an effective inverse temperature in the Jarzynski equality. To this aim, we introduce a suitable parametrization of the initial state in terms of a thermal and a non-thermal component. We determine the value of *β*eff for a large number of measurements and study its dependence on the initial state. Our predictions could be checked experimentally in quantum optics.

**Keywords:** fluctuation theorems; nonequilibrium statistical mechanics; quantum thermodynamics

### **1. Introduction**

Fluctuation theorems relate fluctuations of thermodynamic quantities of a given system to equilibrium properties evaluated at the steady state [1–4]. This statement finds fulfillment in the Jarzynski equality, whose validity has been extensively theoretically and experimentally discussed in the last two decades for both classical and quantum systems [5–17].

The evaluation of the relevant work originated by using a coherent modulation of the system Hamiltonian has been the subject of intense investigation [18–24]. A special focus was devoted to the study of heat and entropy production, obeying of the second law of thermodynamics, the interaction with one or more external bodies, and/or the inclusion of an observer [25–34].

In this respect, the energy variation and the emission/absorption of heat induced—and sometimes enhanced—by the application of a sequence of quantum measurements was studied [17,35–41]. In such a case, the term quantum heat has been used [41,42]; we will employ it as well in the following to refer to the fact that the fluctuations of energy exchange are induced by quantum projective measurements performed during the time evolution of the system. As recently discussed in Ref. [17], the information about the fluctuations of energy exchanges between a quantum system and an external environment may

be enclosed in an energy scaling parameter that only depends on the initial and the asymptotic (for long times) quantum states resulting from the system dynamics.

In Ref. [41], the effect of stochastic fluctuations on the distribution of the energy exchanged by a quantum two-level system with an external environment under sequences of quantum measurements was characterized and the corresponding quantum-heat probability density function was derived. It has been shown that, when a stochastic protocol of measurements is applied, the quantum Jarzynski equality is obeyed. In this way, the quantum-heat transfer was characterized for two-level systems subject to projective measurements. Two-level systems have the property that a density matrix in the energy basis (as the one obtained after a measurement of the Hamiltonian operator [18]) can be always written as a thermal state and, therefore, the Jarzynski equality has a 1 on its right-hand side. Therefore, a natural issue to be investigated is the study of quantum-heat fluctuation relations for *N*-level systems, e.g., *N* = 3, where this property of the initial state of a two-points measurement scheme of being thermal is no longer valid. So, it would be desirable, particularly for the case of a large number of quantum measurements, to study the properties of the characteristic function of the quantum heat when initial states cannot be written as thermal states.

With the goal of characterizing the effects of having arbitrary initial conditions, in this paper, we study quantum systems described by finite dimensional Hilbert spaces, focusing on the case of three-level systems. We observe that finite-level quantum systems may present peculiar features with respect to continuum systems. As shown in Ref. [23], even when the quantum Jarzynski equality holds and the average of the exponential of the work equals the free energy difference, the variance of the energy difference may diverge for continuum systems, an exception being provided by finite-level quantum systems.

In this paper we analyze, using numerical simulations, (i) the distribution of the quantum heat originated by a three-level system under a sequence of *M* projective measurements in the limit of a large *M*, and (ii) the behavior of an energy parameter *β*eff, such that the Jarzynski equality has 1 on its right-hand side, always in the limit of a large *M*. We also discuss the dependence of *β*eff on the initial state, before the application of the sequence of measurements.

### **2. The Protocol**

Let us consider a quantum system described by a finite dimensional Hilbert space. We denote with *H* the time-independent Hamiltonian of the system that admits the spectral decomposition

$$H = \sum\_{k=1}^{N} E\_k |E\_k\rangle\langle E\_k| \,\,\,\,\,\tag{1}$$

where *N* is the dimension of the Hilbert space. We assume that no degeneration occurs in the eigenstates of *H*.

At time *t* = 0−, just before the first measurement of *H* is performed, the system is supposed to be in an arbitrary quantum state described by the density matrix *ρ*0 s.t. [*H*, *ρ*0] = 0. This allows us to write

$$\rho\_0 = \sum\_{k=1}^{N} c\_k |E\_k\rangle\langle E\_k| \,\,\,\,\,\tag{2}$$

where 1 ≥ *ck* ≥ 0 ∀*k* = 1, . . . , *N* and ∑*N k ck* = 1.

Then, we assume that the fluctuations of the energy variations, induced by a given transformation of the state of the system, are evaluated by means of the so-called two-point measurement (TPM) scheme [18]. According to this scheme, a quantum projective measurement of the system hamiltonian is performed

both at the initial and the final times of the transformation. This hypothesis justifies the initialization of the system in a mixed state, as given in Equation (2). By performing a first projective energy measurement, at time *t* = 0+, the system is in one of the states *ρn* = |*EnEn*| with probability *pn* = *En*|*ρ*0|*En*, while the system energy is *En*. Afterwards, we suppose that the system *S* is subject to a number *M* of consecutive projective measurements of the generic observable

$$O = \sum\_{k=1}^{N} \Omega\_k |\Omega\_k\rangle\langle\Omega\_k| \; , \tag{3}$$

where Ω*k* and |Ω*k* denote, respectively, the outcomes and the eigenstates of *O*. According to the postulates of quantum mechanics, the state of the system after one of these projective measurements is given by one of the projectors |Ω*n*Ω*n*|. Between two consecutive measurements, the system evolves with the unitary dynamics generated by *H*, i.e., *U*(*τi*) = *e*−*iHτi* , where *h*¯ has been set to unity and the waiting time *τi* is the time difference between the (*i* − 1)th and the *i* th measurement of *O*.

In general, the waiting times *τi* can be random variables, and the sequence (*τ*1, ... , *τM*) is distributed according to the joint probability density function *p*(*τ*1, ... , *τM*). The last (i.e, the *M*th) measurement of *O* is immediately followed by a second projective measurement of the energy, as prescribed by the TPM scheme. By denoting with *Em* the outcome resulting from the second energy measurement of the scheme, the final state of the system is *ρm* = |*EmEm*| and the quantum heat *Q* exchanged during the transformation is thus given by

$$Q = E\_{\mathfrak{M}} - E\_{\mathfrak{N}} \,. \tag{4}$$

As *Q* is a random variable, one can define the characteristic function

$$G(\varepsilon) \equiv \left\langle e^{-\varepsilon Q} \right\rangle = \sum\_{m,n} p\_{m|n} p\_n e^{-\varepsilon (E\_m - E\_n)} \, , \tag{5}$$

where *pm*|*n* denotes the probability of obtaining *Em* at the end of the protocol conditioned to have measured *En* at the first energy measurement of the TPM scheme. If the initial state is thermal, i.e., *ρ*0 = *e*−*βH*/*Z*, then one recovers the Jarzynski equality stating that

$$G(\beta) = \left\langle e^{-\beta Q} \right\rangle = 1 \,\,. \tag{6}$$

Let us notice that *G*() is a convex function such that *G*(0) = 1 and *G*(±∞) → +∞, as discussed in Ref. [17]. Hence, as long as *∂G ∂* (0) = 0, one can unambiguously introduce the parameter *β*eff = 0 defined by the relation

$$G(\beta\_{\rm eff}) = 1 \,, \tag{7}$$

which formally plays the role of an effective inverse temperature. Focusing on three-level systems, in the following, we will study the characteristic function *G*() and the properties of such a parameter *β*eff. For comparison, we first pause in the next subsection to discuss what happens for two-level systems.

#### *Intermezzo on Two-Level Quantum Systems*

We pause here to remind the reader of the results for two-level systems. The state of any two-level system, diagonal on the Hamiltonian basis, is a thermal state for some value of *β*. Of course, if the state is thermal, the value of *β*eff trivially coincides with *β*. In particular, in Ref. [41], the energy exchanged between a two-level quantum system and a measurement apparatus was analyzed, with the assumption that the repeated interaction with the measurement device can be reliably modeled by a sequence of projective measurements occurring instantly and at random times. Numerically, it has been observed that,

as compared with the case of measurements occurring at fixed times, the two-level system exchanges more energy in the presence of randomness when the average time between consecutive measurements is sufficiently small in comparison with the inverse resonance frequency. However, the quantum-heat Jarzynski equality, related to the equilibrium properties of the transformation applied to the system, is still obeyed, as well as when the waiting times between consecutive measurements are randomly distributed and for each random realization. These results are theoretically supported by the fact that, in the analyzed case, the dynamical evolution of the quantum system is unital [11,12]. A discussion on the values of the parameter *β*eff, extracted from experimental data for nitrogen-vacancy (NV) centers in diamonds subject to projective measurements in a regime where an effective two-level approximation is valid was recently presented in Ref. [17].

In Figure 1, we plot the quantum-heat characteristic function *e*−*βQ* as a function of the parameter *c*1 that appears in the decomposition of the initial state *ρ*0 with respect to the energy eigenstates |*E*1 and |*E*2

$$
\rho\_0 = c\_1 |E\_1\rangle\langle E\_1| + c\_2 |E\_2\rangle\langle E\_2|\,,\tag{8}
$$

where *c*2 = 1 − *c*1. The function is plotted for three values of the parameter |*a*| 2, used to parametrize the eigenstates {|Ω1, |Ω2} of *O* as a function of the energy eigenstates of the system, i.e.,

$$|\Omega\_1\rangle = a|E\_1\rangle - b|E\_2\rangle \quad \text{and} \quad |\Omega\_2\rangle = b|E\_1\rangle + a|E\_2\rangle \,. \tag{9}$$

with |*a*| 2 + |*b*| 2 = 1 and *a*∗*b* = *ab*∗. As a result, one can observe that *e*−*βQ* = 1 for the value of *c*1 ensuring that *ρ*0 = *e*−*βH*/*Z*. Further details can be found in Ref. [41].

**Figure 1.** Quantum-heat characteristic function *e*−*βQ* for a two-level quantum system as a function of *c*1 in Equation (8) for three values of |*a*| 2, which characterizes the initial state. The function is obtained from numerical simulations performed for a system with a Hamiltonian *H* = *J*(|01| + |10|) subject to a sequence of *M* = 5 projective measurements. The latter are separated by a fixed waiting time *τ* = 0.5, averaged over 2000 realizations, with *E*1,2 = ±1 and *β* = 3/2. Units are used with ¯*h* = 1 and *J* = 1.

The *N* > 2 case is trickier: Since, in general, the initial state is no longer thermal, it is not trivial to determine the value of *β*eff, and its dependence on the initial condition is interesting to investigate. In this regard, in the following, we will numerically address the *N* = 3 case by providing results on the

asymptotic behavior of the system in the limit *M* 1. For the sake of simplicity, from here on, we assume that the values of the waiting times *τi* are fixed for any *i* = 1, ... , *M*. However, our findings turn out to be the same for any choice of the marginal probability distribution functions *p*(*τi*), up to "pathological" cases such as *p*(*τ*) = *δ*(*τ*). So, having random waiting times does not significantly alter the scenario emerging from the presented results.

### **3. Parametrization of the Initial State**

In this paragraph, we introduce a parametrization of the initial state *ρ*0 for the *N* = 3 case, which will be useful for a thermodynamic analysis of the system.

As previously recounted, for a two-level quantum system in a mixed state, as given by Equation (2), it is always formally possible to define a temperature. In particular, one can implicitly find an effective inverse temperature, making *ρ*0 a thermal state, by solving the following equation for *β* ∈ R

$$e^{-\beta(E\_2 - E\_1)} = \frac{c\_2}{c\_1} \,. \tag{10}$$

In the *N* = 3 case, however, we have two independent parameters in (2) (the third parameter indeed enforces the condition Tr[*ρ*0] = 1), and thus, in general, it is no longer possible to formally define a single temperature for the state. Here, we propose the following parametrization of *c*1, *c*2, *c*3 that generalizes the one of Equation (10).

We denote as *partial* effective temperatures the three parameters *b*1, *b*2, *b*3, defined through the ratios of *c*1, *c*2, and *c*3 *c*2

$$\frac{c\_2}{c\_1} = e^{-b\_1(E\_2 - E\_1)}, \qquad \frac{c\_3}{c\_2} = e^{-b\_2(E\_3 - E\_2)}, \qquad \frac{c\_1}{c\_3} = e^{-b\_3(E\_1 - E\_3)}, \tag{11}$$

such that, for a thermal state, *bk* = *β*, ∀*k* = 1, 2, 3. The three parameters are not independent, as they are constrained by the relation *c*2

*c*1 *c*3 *c*2 *c*1 *c*3 = 1 , (12)

which gives in turn the following equality

$$b\_1(E\_2 - E\_1) + b\_2(E\_3 - E\_2) + b\_3(E\_1 - E\_3) = 0\,\text{.}\tag{13}$$

By introducing Δ1 = *E*2 − *E*1, Δ2 = *E*3 − *E*2, and Δ3 = *E*1 − *E*3, Equation (13) can be written as

$$\sum\_{k=1}^{3} b\_k \Delta\_k = 0 \,, \tag{14}$$

where by definition

3 ∑ *k*=1 Δ*k* = 0 . (15)

Thus, as expected, the thermal state is a solution of the condition (14) for any choice of *E*1, *E*2, *E*3. This has also a geometric interpretation. In the space of the Δ*k*, *k* = 1, 2, 3, Equation (14) becomes an orthogonality condition between the Δ*k* and the *bk* vectors, while (15) defines a plane that is orthogonal to the vector (1, 1, 1). When *bk* is proportional to (1, 1, 1), the orthogonality condition is automatically satisfied and one finds a thermal state. This suggests that, in general, one can conveniently parametrize *bk* in terms of both the components that are orthogonal and parallel to the plane ∑3 *k*=1 Δ*k* = 0. Such terms

have the physical meaning of the thermal and non-thermal components of the initial state. Formally, this means that we can parametrize each *bk* through the fictitious inverse temperature *β* and a deviation *α*, i.e.,

$$\beta(b\_1, b\_2, b\_3) = \beta(1, 1, 1) + \frac{a}{v} (\Delta\_3 - \Delta\_2, \Delta\_1 - \Delta\_3, \Delta\_2 - \Delta\_1) \,. \tag{16}$$

where *v* acts as a normalization constant

$$
\sigma^2 = 3 \left( \Lambda\_1^2 + \Lambda\_2^2 + \Lambda\_3^2 \right) \,. \tag{17}
$$

Hence, taking into account the normalization constraint, the coefficients *ck* are given by

$$\varepsilon\_1 = \frac{1}{1 + e^{-b\_1 \Delta\_1} + e^{b\_3 \Delta\_3}}, \qquad \varepsilon\_2 = \frac{1}{1 + e^{-b\_2 \Delta\_2} + e^{b\_1 \Delta\_1}}, \qquad \varepsilon\_3 = \frac{1}{1 + e^{-b\_3 \Delta\_3} + e^{b\_2 \Delta\_2}}, \tag{18}$$

or, in terms of the parameters *α* and *β*,

$$\mathcal{L}\_1 = \frac{1}{\overline{Z}} \exp\left[-\beta E\_1 + \frac{a}{\upsilon}(E\_2 - E\_3)^2\right],\\ \mathcal{c}\_2 = \frac{1}{\overline{Z}} \exp\left[-\beta E\_2 + \frac{a}{\upsilon}(E\_3 - E\_1)^2\right],\\ \mathcal{c}\_3 = \frac{1}{\overline{Z}} \exp\left[-\beta E\_3 + \frac{a}{\upsilon}(E\_1 - E\_2)^2\right],\tag{19}$$

where *Z*˜ is a pseudo-partition function ensuring the normalization of the *ck*'s

$$Z = Z(a, \beta) \equiv e^{-\beta E\_1 + \frac{\mathbf{g}}{v}(E\_2 - E\_3)^2} + e^{-\beta E\_2 + \frac{\mathbf{g}}{v}(E\_3 - E\_1)^2} + e^{-\beta E\_3 + \frac{\mathbf{g}}{v}(E\_1 - E\_2)^2}.\tag{20}$$

Let us provide some physical intuition about the parameters *α* and *β*: For *α* = 0, we recover a thermal state, whereby *c*1 > *c*2 > *c*3 if *β* > 0, or vice versa if *β* < 0. On the other hand, the non-thermal component *α* can be used to obtain a non-monotonic behavior of the coefficients *ck*. For example, for *β* = 0, since (*E*3 − *E*1)2 is greater than both (*E*3 − *E*2)2 and (*E*1 − *E*2)2, one finds that *c*2 > (<) *c*1, *c*3 if *α* > (<) 0.

As a final remark, it is worth noting that we can reduce the dimension of the space of the parameters. In particular, without loss of generality, one can choose the zero of the energy by taking *E*2 = 0 (and then *E*3 > 0, *E*1 < 0), or we can reduce our analysis to the cases with *β* > 0. As a matter of fact, the parametrization is left unchanged by the transformation {*β* → −*β*, *Ek* → −*Ek*}, with the result that the case of *β* < 0 can be explored by simply considering *β* > 0 in the fictitious system with *E k* = −*Ek* (here, the choice of *E*2 = 0 guarantees that this second case can be simply obtained by substituting *E*1 with *E*3).

### **4. Large** *M* **Limit**

Here, we numerically investigate the behavior of a three-level system subject to a sequence of *M* projective quantum measurements with a large *M* (asymptotic limit) and where *τ* is not infinitesimal. From here on, we adopt the language of spin-1 systems, and we thus identify *O* with *Sz*.

In the asymptotic limit, the behavior of the system is expected not to depend on the choice of the evolution Hamiltonian, with the exception that at least one of the eigenstates of *Sz* is also an energy eigenstate. In such a case, indeed, if the energy outcome corresponding to the common eigenstate is obtained by the first measurement of the TPM scheme, then the evolution is trivially deterministic, as the system is locked in the measured eigenstate.

Choosing a generic observable (with no eigenstates in common with *H*), numerical simulations (cf. Figure 2) suggest that our protocol leads the system to the completely uniform state. The latter can be interpreted as a canonical state with *β* = 0 (notice that this result holds in the situation we are analyzing, with a finite dimensional Hilbert space). The system evolves with Hamiltonian *H* = *ω*1*Sz* + *ω*2*Sx*, where the energy units are chosen such that *ω*1 = 1 and *ω*2 = 1 2 . It is initialized in the state *ρ*0 with

{*c*1 = 0.8, *c*2 = 0.01, *c*3 = 0.19}, corresponding to *α* ≈ −2, 32 and *β* ≈ 1, 96 , and we performed *M* = 20 projective measurements of the observable *O* = *Sz* separated by the time *τ* = 1.

**Figure 2.** Histogram of the initial (dashed) and final (solid) energy outcomes for the TPM scheme described in the text performed over 3 × 105 realizations. While the initial state is non-uniform, the final state is practically uniform over the three energy levels.

This numerical finding allows us to derive an analytic expression of the quantum-heat characteristic function. In this regard, as the final state is independent of the initial one for a large *M*, the joint probability of obtaining *En* and *Em*, after the first and the second energy measurement respectively, is equal to

$$p\_{\rm mu} = \frac{1}{3}c\_{\rm u} \,. \tag{21}$$

Hence,

$$G(\varepsilon) = \left\langle e^{-\varepsilon Q} \right\rangle = \frac{1}{\mathfrak{Z}} \sum\_{m,n=1}^{3} c\_m e^{-\varepsilon (E\_m - E\_n)} = \frac{1}{\mathfrak{Z}} \sum\_{n=1}^{3} e^{-\varepsilon E\_n} \sum\_{m=1}^{3} c\_m e^{\varepsilon E\_m}.\tag{22}$$

As a consequence, *G* can be expressed in terms of the partition function *Z*(*β*) and of the pseudo-partition function introduced in Equation (20), i.e.,

$$G(\epsilon; \mathfrak{a}, \boldsymbol{\beta}) = \frac{Z(\epsilon)}{Z(0)} \frac{\mathbb{Z}(\mathfrak{a}, \boldsymbol{\beta} - \epsilon)}{\mathbb{Z}(\mathfrak{a}, \boldsymbol{\beta})} \,. \tag{23}$$

Regardless of the choice of the system parameters, the already known results are straightforwardly recovered, i.e., *G*(0) = 1 and *G*(*β*) = 1 for *α* = 0 (initial thermal state). In Figure 3, our analytical expression for *G*(*ε*) is compared with its numerical estimate for two different Hamiltonians, showing a very good agreement.

We remark that the distribution of *ρ* after the second energy measurement could also be obtained by simply imposing the maximization of the von Neumann entropy. This is reasonable, since the measurement device is macroscopic and can provide any amount of energy. As a final remark, notice that, in the numerical findings of Figure 3, the statistics of the quantum-heat fluctuations originated by the system respect the same ergodic hypothesis that is satisfied whenever a sequence of quantum measurements is performed on a quantum system [43–45]. In particular, in Figure 3, one can observe that the analytical

expression of *G*() for a large *M* (i.e., in the asymptotic regime obtained by indefinitely increasing the time duration of the implemented protocol) practically coincides with the numerical results obtained by simulating a sequence with a finite number of measurements (*M* = 20) but over a quite large number (3 × 105) of realizations. This evidence is quite important, because it means that the quantum-heat statistics is homogeneous and fully take into account even phenomena occurring with very small probability in a single realization of the protocol.

**Figure 3.** Comparison of the analytic expression (22) of the asymptotic (large-*M*) quantum-heat characteristic function *G*() (blue solid lines) with the numerical results averaged over 3 × 105 realizations (red dots). The initial state is the same as in Figure 2 and, again, *O* = *Sz*. In panel (**a**), the Hamiltonian is the same as in Figure 2, while in panel (**b**) the Hamiltonian is *H* = *ω*1*S*2 *z* + *ω*2*Sx*, with *ω*1 = 2*ω*2 = 1.

### **5. Estimates of** *β***eff**

In this section, we study the behavior of *β*eff, i.e., the nontrivial solution of *G*(*β*eff) = 1, as a function of the initial state (parametrized by *α* and *β*) and of the energy levels of the system. Let us first notice that, by starting from Equation (22), obtaining an analytical expression for *β*eff in the general case appears to be a very non-trivial task. Thus, in Figure 4, we numerically compute *β*eff as a function of *α* (the non-thermal component of *ρ*0) for different values of *β*. Three representative cases for the energy levels are taken, i.e., {*E*1 = −1, *E*2 = 0, *E*3 = 3}, {*E*1 = −1, *E*2 = 0, *E*3 = 1}, and {*E*1 = −3, *E*2 = 0, *E*3 = 1}, respectively. This choice allows us to deal both with the cases *E*3 − *E*2 > *E*2 − *E*1 and *E*3 − *E*2 < *E*2 − *E*1. The choice of the energy unit is such that the smallest energy gap between *E*3 − *E*2 and *E*2 − *E*1 is set to one. As stated above, we consider *β* > 0; the corresponding negative values of the inverse temperature are obtained by taking *E k* = −*Ek* with *β* = −*β*. As expected, for *α* = 0, we have *β*eff = *β*, regardless of the values of the *Ek*'s.

In the next two subsections, we continue discussing in detail the findings of Figure 4, presenting the asymptotic behaviors for large positive and negative values of *α*.

**Figure 4.** Behavior of *β*eff as a function of *α* for different values of *β* ∈ [0, 2.5]. We have chosen: (**a**) {*E*1 = −1, *E*2 = 0, *E*3 = 3}, (**b**) {*E*1 = −1, *E*2 = 0, *E*3 = 1}, and (**c**) {*E*1 = −3, *E*2 = 0, *E*3 = 1}, respectively.

### *5.1. Asymptotic Behavior for a Large Positive α*

From Figure 4, one can deduce that, for large positive values of *α* (corresponding to having as the initial density operator the pure state *ρ*0 = |*E*2*E*2|), *β*eff → *β*¯ eff, which only depends on *E*1 and *E*3. This asymptotic value *β*¯ eff is positive if *E*3 − *E*2 > *E*2 − *E*1, negative if *E*3 − *E*2 < *E*2 − *E*1, and zero when *E*3 − *E*2 = *E*2 − *E*1. To better explain the plots in Figure 4, let us consider the analytic expression of *G*(). In this regard, for a large *α* and finite *β*, we can write

$$Z(u, \beta) \approx e^{-\beta E\_2 + \frac{4}{v} (E\_3 - E\_1)^2} \,, \tag{24}$$

so that, by using Equation (23), the condition *G*(*β*eff) = 1 reads as *e*−*β*¯ eff(*E*1−*E*2) + *e*−*β*¯ eff(*E*3−*E*2) = 2, or, setting *E*2 = 0,

$$e^{-\vec{\beta}\_{\text{off}}E\_1} + e^{-\vec{\beta}\_{\text{off}}E\_3} = \mathbf{2} \,. \tag{25}$$

Notice that, if *E*3 = −*E*1 the only solution of Equation (25) is *β*¯ eff = 0, while a positive solution appears for *E*3 > −*E*1 and a negative one if *E*3 < −*E*1, thus confirming what was observed in the numerical simulations. Moreover, by replacing *E*1 → −*E*3 and *E*3 → −*E*3, the value of *β*¯ eff changes its sign.

Now, without loss of generality, let us fix the energy unit so that *E*1 = −1. The behavior of *β*¯ eff as a function of *E*3 is shown in Figure 5. We observe a monotonically increasing behavior of *β*¯ eff up to a constant value for *E*3 |*E*1| = 1. Once again, this value can be analytically computed from Equation (25), which, for a large value of *E*3, gives *β*¯ eff = ln 2. Putting together all of the above considerations and restoring the energy scales, the limits of *β*¯ eff are the following

$$-\frac{\ln 2}{E\_3 - E\_2} < \vec{\beta}\_{\rm eff} < \frac{\ln 2}{E\_2 - E\_1} \,. \tag{26}$$

The lower and the upper bounds of *β*¯ eff are also shown in Figure 5, in which *E*1 = −1 and *E*2 = 0.

### *5.2. Asymptotic Behavior for a Large Negative α*

From Figure 4, one can also conclude that, for large negative values of *α*, the behavior of *β*eff is linear with *α*

$$
\beta\_{\rm eff} \approx r \alpha \,, \tag{27}
$$

with *r* > 0 if *E*3 − *E*2 > *E*2 − *E*1, *r* = 0 when *E*3 − *E*2 = *E*2 − *E*1, and *r* is negative otherwise. This divergence is easily understood: In fact, the limit *α* → −∞ (for finite *β*) corresponds to the initial state *ρ*0 = |*E*1*E*1| when *E*3 − *E*2 < *E*2 − *E*1 or *ρ*0 = |*E*3*E*3| if *E*3 − *E*2 > *E*2 − *E*1. On the other hand, those states (thermal states with *β*eff = *β* = ±∞) are also reached in the limits *β* → ±∞ with a finite *α*. This simple argument does not imply the linear divergence of *β*eff as in Equation (27), nor does it provide insights about the value of *r*, which, however, can be derived from Equation (23). Although the calculation makes a distinction on the sign of *r* depending on whether *E*3 − *E*2 is greater or smaller than *E*2 − *E*1, the result is independent of this detail. In particular, by considering the case *E*3 − *E*2 > *E*2 − *E*1 (*r* > 0) and taking into account the divergence of *β*eff = *rα*, we find in the *α* → −∞ regime that the characteristic function *G*(*β*eff) has the following form:

$$G(\beta\_{\rm eff}) = \frac{1}{3} + \text{const} \times e^{-a|\Delta\_3| \left[r - \frac{(\Delta\_1 - \Delta\_2)}{v}\right]} \,. \tag{28}$$

Hence, in order to ensure that *G*(*β*eff) = 1 3 in the limit *α* → −∞, the following relation has to be satisfied

$$r = \frac{E\_1 + E\_3 - 2E\_2}{v} \,. \tag{29}$$

The numerical estimate of *rv* as a function of *E*3 is shown in Figure 5. The numerical results confirm the linear dependence of *β*eff as a function of *α*.

**Figure 5.** (**a**) Behavior of the asymptotic value *β*¯ eff for a large positive *α* (*α* = 20) as a function of *E*3 (solid blue line) with *E*1 = −1 and *E*2 = 0. We compare the curve with its limiting (lower and upper) values, defined in Equation (26) (dash-dotted red lines). (**b**) Behavior of the asymptotic slope *r*, rescaled for *v*, for a large negative *α* (*α* = −20) as a function of *E*3. In both cases, *E*2 = 0 and *E*1 = −1.

### *5.3. Limits of the Adopted Parametrization*

The parametrization introduced in Section 3 is singular in correspondence of the initial states *ρ*0 with one or more coefficients *ck* equal to zero. In this regard, as remarked above, initial pure states can be easily obtained in the limits *β* → ±∞, *α* finite (corresponding to *ρ*0 = |*E*1*E*1| and *ρ*0 = |*E*3*E*3|, respectively) and *α* → +∞, *β* finite (that provides *ρ*0 = |*E*2*E*2|). Instead, initial states with only a coefficient *ck* equal to zero, namely

$$\rho\_0 = q|E\_1\rangle\langle E\_1| + (1-q)|E\_2\rangle\langle E\_2| \ , \quad \rho\_0 = q|E\_1\rangle\langle E\_1| + (1-q)|E\_3\rangle\langle E\_3| \ , \quad \rho\_0 = q|E\_2\rangle\langle E\_2| + (1-q)|E\_3\rangle\langle E\_3| \ , \tag{30}$$