Canonical Density Matrices from Eigenstates of Mixed Systems

One key issue of the foundation of statistical mechanics is the emergence of equilibrium ensembles in isolated and closed quantum systems. Recently, it was predicted that in the thermodynamic (N→∞) limit of large quantum many-body systems, canonical density matrices emerge for small subsystems from almost all pure states. This notion of canonical typicality is assumed to originate from the entanglement between subsystem and environment and the resulting intrinsic quantum complexity of the many-body state. For individual eigenstates, it has been shown that local observables show thermal properties provided the eigenstate thermalization hypothesis holds, which requires the system to be quantum-chaotic. In the present paper, we study the emergence of thermal states in the regime of a quantum analog of a mixed phase space. Specifically, we study the emergence of the canonical density matrix of an impurity upon reduction from isolated energy eigenstates of a large but finite quantum system the impurity is embedded in. Our system can be tuned by means of a single parameter from quantum integrability to quantum chaos and corresponds in between to a system with mixed quantum phase space. We show that the probability for finding a canonical density matrix when reducing the ensemble of energy eigenstates of the finite many-body system can be quantitatively controlled and tuned by the degree of quantum chaos present. For the transition from quantum integrability to quantum chaos, we find a continuous and universal (i.e., size-independent) relation between the fraction of canonical eigenstates and the degree of chaoticity as measured by the Brody parameter or the Shannon entropy.


Introduction
As first recognized by Ludwig Boltzmann [1,2] "molecular" chaos lies at the core of the foundation of classical statistical mechanics. Only when the phase space of an isolated mechanical system is structureless can the motion be safely assumed to be ergodic and the equal a priori probability for phase space points on the energy hypersurface, the basic tenet of the microcanonical ensemble, is realized. Moreover, chaotic dynamics is "mixing", thereby enforcing the approach to the thermal equilibrium state from "almost all" out-of-equilibrium initial conditions. While any large isolated system is expected to be described by a microcanonical ensemble, any well-defined small subsystem thereof that is only allowed to exchange energy with the remainder of the large system (referred to as bath or environment in the following) is described by the canonical ensemble. The phasespace density of the subsystem is weighted by the Boltzmann factor e −βH s , where H s is the Hamilton function of the subsystem, β = 1/k B T with T the temperature imprinted by the environment and k B the Boltzmann constant. However, when the phase space of the system is not chaotic but rather dominated by regular motion on KAM tori [3,4], neither ergodicity nor mixing is a priori assured, and thermalization of an initial non-equilibrium state may be elusive. The implicit assumption of classical equilibrium statistical mechanics is that in the limit of a large number of degrees of freedom, chaos is generic for any interacting many-particle system.
How those concepts translate into quantum physics has remained a topic of great conceptual interest and lively debate [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]. Renewed interest is stimulated by the experimental accessibility of ultracold quantum gases [22][23][24][25][26][27][28], trapped ions [29], and nanosystems [30,31], where many of the underlying concepts became quantitatively accessible in large but finite quantum systems in unprecedented detail. The foundation of thermalization of quantum systems has been pioneered by von Neumann in terms of the quantum ergodic theorem [32][33][34][35][36][37]. Accordingly, the entropy is an increasing function of time and expectation values of generic macroscopic observables for pure states formed by coherent superposition of states within microscopic energy shells converge to that of the microcanonical ensemble provided that the energy spectrum of the system is strictly non-degenerate. Recently, this description of thermal equilibrium states was extended to the notion of canonical typicality [37][38][39][40]. Accordingly, starting from almost any pure state formed by a coherent superposition of energy eigenstates of a large isolated many-body system with eigenenergies within a given energy shell [E, E + ∆E] of macroscopically small thickness ∆E, the reduction to a small subsystem by tracing out the degrees of freedom of the bath will yield the same reduced density matrix one would obtain from the reduction of the microcanonical density matrix for the entire system. If the bath is sufficiently large and the interactions between the bath and the subsystem sufficiently weak, the reduced density matrix corresponds to the standard canonical density matrixρ s = e −βĤ s /Tr[e −βĤ s ] withĤ s the Hamilton operator of the subsystem. The proof of this canonical typicality invokes the intrinsic randomness of the expansion coefficients of the pure state in terms of entangled subsystem-bath states. The latter assumption goes back to the notion of intrinsic quantum complexity of entangled states in large systems put forward already by Schrödinger [41].
An alternative approach to thermalization is tied to the eigenstate thermalization hypothesis (ETH) [5][6][7][8] first put forward by Landau and Lifshitz [42], stating that basic properties of statistical mechanics can emerge not only from ensemble averages but from typical single wavefunctions. However, the condition under which such an equivalence may emerge has remained open. The more recent formulation of the ETH [5][6][7][8] invokes the notion of quantum chaos and Berry's conjecture. Characteristics of quantum chaos were originally identified in few-degrees of freedom systems whose classical limit exhibits chaos [43][44][45][46][47][48][49]. Nowadays, the notion of quantum chaos is invoked more generally for systems that display the same signatures such as energy level distributions predicted by random matrix theory (RMT) [43,48,50,51] or randomness of wavefunction amplitudes [5,52] even when a well-defined classically chaotic counterpart is not known. The ETH conjectures that for chaotic systems, the diagonal matrix elements of any generic local observable taken in the energy eigenstate basis are smooth functions of the total energy while the off-diagonal elements are exponentially decreasing randomly fluctuating variables with zero mean [6][7][8]. If the ETH is valid for a specific system, individual eigenstates show thermal properties upon reduction to a small subsystem. The ETH has been shown to hold for a large variety of systems without a classical analogue [24,[53][54][55][56][57][58][59][60][61][62][63][64]. Deviations from the ETH have been observed for local observables in finite systems of hard-core bosons and spin-less fermions [57,58,65] when the energy level distribution deviates from the Wigner-Dyson level statistics of RMT characteristic for chaotic systems.
In the present paper, we explore the quantitative relationship between thermal properties of reduced density matrices (RDMs) emerging from single isolated eigenstates of the entire system and quantum chaos. More specifically, we want to address the question: Is for large but finite systems quantum chaos a conditio sine qua non for the emergence of the Gibbs ensemble, i.e., the canonical ensemble of the subsystem, from eigenstates of the entire system? Or is quantum entanglement and complexity in these systems itself sufficient to render the reduced density matrix of a small subsystem canonical? To this end, we determine the fraction of canonical density matrices emerging upon reduction from the entire set of eigenstates. We explore the existence of a quantitative relationship between the fraction of eigenstates that upon reduction lead to canonical eigenstates, also termed fraction of canonical eigenstates, and the degree of quantum chaos of the entire system. We unravel the connection between this eigenstate canonicity and quantum chaos by exact diagonalization of a large yet finite mesoscopic quantum system. We emphasize that this measure addresses isolated energy eigenstates of the many-body system, in contrast to coherent superpositions of energy eigenstates from a given energy shell of finite width with random expansion coefficients as invoked in the well-established notion of canonical typicality [37][38][39][40].
As a prototypical case in point, we consider an itinerant impurity embedded in a spinpolarized Fermi-Hubbard system. Unlike impurity models for disordered systems [66], our model is fully deterministic. All key ingredients for the realization of the present system, i.e., discrete lattice, tunable interactions, and impurity can be experimentally realized with ultracold fermionic atoms (see, e.g., [67][68][69][70][71][72]). In the present scenario, the impurity serves as a probe or "thermometer" in the isolated many-body quantum system providing an unambiguous subsystem-bath decomposition with tunable coupling strength between subsystem and bath. Moreover, our system features a tunable transition from quantum chaos to quantum integrability without invoking any extrinsic stochasticity or disorder [66]. The fact that the subsystem consists of a distinguishable particle has a number of distinct advantages: The reduced density matrix of the probe is uniquely defined and its properties are basis-independent. No choice of a specific basis for the probe such as the independentparticle basis is involved. Moreover, its canonical RDM approaches a Maxwell-Boltzmann rather than a Fermi-Dirac distribution for indistinguishable fermions. Its thermal state is thus characterized by a single equilibrium parameter, the temperature T, without the need for introducing a chemical potential µ, thereby improving the numerical accuracy of the test of canonicity. We measure the proximity of the reduced density matrix of the impurity to the canonical density matrix and identify a direct and size-independent correlation between the fraction of canonical eigenstates and quantum chaos.
The paper is structured as follows. In Section 2, we introduce our impurity-Fermi-Hubbard model which serves as a prototypical (sub)system-environment model system. Quantitative measures for quantum chaos are introduced in Section 3. The mapping of spectral properties of this isolated many-body system onto thermal states of the impurity within the framework of the microcanonical and canonical ensembles are discussed in Section 4. The distance in Liouville space between the reduced density matrix of the impurity and a generic canonical density matrix is analyzed and the relation between the fraction of canonical eigenstates and quantum chaos is established in Section 5. Concluding remarks are given in Section 6.

The Fermi-Hubbard Model with Impurity
We investigate a variant of the single-band one-dimensional Fermi-Hubbard model which is particularly well suited to study entanglement and quantum correlations between subsystem and its environment or bath. The bath is represented by spin-polarized fermions enforcing single occupancy of sites by bath particles while the distinguishable impurity can occupy any site. Accordingly, the Hamiltonian of the total system is given bŷ where the Hamiltonian of the subsystem, i.e., the impurity (I), iŝ while the Hamiltonian of the bath iŝ The interaction between the subsystem and the bath is given bŷ The operatorsâ j andâ † j (b j andb † j ) are the creation and annihilation operators of the impurity (bath particles) on site j with the anticommutation The operatorsn j =â † jâ j andN j =b † jb j correspond to the number operators of impurity and bathn j |j = n j |j andN j |j = N j |j with occupation numbers n j and N j of site j, respectively. J I (J B ) describes the hopping matrix elements of the impurity (bath particles). The bath particles interact with each other by a nearestneighbor interaction with strength W BB while the impurity interacts with the bath particles via an on-site interaction with strength W IB . The Hubbard chain has M s sites with Dirichlet boundary conditions imposed at the edges. An additional very weak external background potential (V J I , J B ) with on-site matrix element V(j) (j = 1, . . . , M s ) is applied, for which we use a linear (n = 1) or quadratic (n = 2) function in order to remove residual geometric symmetries such that the irreducible state space coincides with the entire state space and symmetry related degeneracies are lifted. Alternative impurity models were recently suggested for the investigation of the ETH [64]. We solve the system via exact diagonalization to determine all eigenstates and eigenenergies of the entire system. The dimension of the Hilbert space of the system is d H = M s ( M s N B ), where N B is the number of bath particles. We consider typical half-filling configurations with N B ≈ M s /2. The largest M s considered is M s = 15 resulting in a Hilbert space dimension of d H = 96,525 for N B = 7. We set J I = J B = J which also defines the unit of energy (J = 1) in the following. The key advantage of the present model is that it allows to control and tune the properties of the bath separately by varying W BB while keeping fixed the properties of the subsystem whose reduced density matrix we probe. This clear-cut subsystem-bath decomposition allows for the unambiguous probing of the emergence of canonical density matrices, thereby avoiding any ad hoc separation by "cutting out" of the subsystem which then requires the grand canonical density matrix for an open quantum system since both energy and particles can be exchanged [65]. Moreover, its thermal state is unambiguously characterized by T rather than by T and µ as for indistinguishable fermions, thereby improving the numerical reliability of the performed tests.
The present system should be realizable for ultracold fermionic atoms trapped in optical lattices [28,[69][70][71][72][73][74]. All key ingredients required for its realization including tunable interactions and impurity-bath mixtures are available in the toolbox of ultracold atomic physics. We note that tuning the nearest-neighbor interaction W BB between the atoms in optical lattices to large values in the regime of strong correlations, W BB /J B 1, still poses an experimental challenge which might be overcome in the near future.

Measures of Quantum Chaos
The present single-band Fermi-Hubbard model does not possess an obvious classical counterpart whose phase space consists of regions of regular and/or chaotic motion. Lacking such direct quantum-classical correspondence, quantum integrability and quantum chaos in the present system is identified by signatures of the quantum system that have been shown to probe chaotic and regular motion in systems where quantum-classical correspondence does prevail. Several measures of quantum chaos have been proposed that are based on either properties of eigenstates or of the spectrum [14,49,58,[75][76][77][78][79]. As will be shown below, by tuning W BB , we can continuously tune the entire system from the limit of quantum integrability to the limit of quantum chaos across the transition region of a mixed quantum system in which integrable and chaotic motion coexist and explore its impact on the fraction of eigenstates which upon reduction lead to canonical density matrices. The influence of the continuous transition from quantum integrability to quantum chaos on the thermal state of the subsystem will be explored with the help of the present prototypical system.

Spectral Measures
Starting point for analyzing and quantifying quantum chaos by means of spectral statistics is the cumulative spectral function also called the staircase function where E α are the energy eigenvalues of the entire system, and Θ is the Heaviside step function. Its spectral derivative is the density of states (DOS) Examples for N(E) and Ω(E) of the present system are shown in Figure 1.
The smoothed "average" spectral staircase functionN(E) fitted to a polynomial of order 10, also shown in Figure 1, provides the reference for spectral unfolding required for certain measures of quantum fluctuations about the (classical) mean. Accordingly, the unfolded energy spectrum is given by e α =N(E α ). For systems for which quantumclassical correspondence holds,N(E) corresponds to the classical phase space volume in units of Planck's constant h and Ω(E) to the microcanonical energy shell. We note that the saturation of N(E) observed with increasing E (Figure 1a) or, likewise, the bell-shaped curve for the DOS (Figure 1b) decreasing at large E is in the present case a consequence of the single-band approximation of the Fermi-Hubbard model (Equation (1)) and, more generally, appears for systems with a spectrum bounded from above. For realistic macroscopic systems, N(E) and Ω(E) should generically increase monotonically with E. As discussed in more detail below, this non-generic decrease of the density of states observed for the present as well as for other finite and mesoscopic systems has implications for the ensuing thermal properties.
The probability density P(s) of the nearest-neighbor level spacings (NNLS), s = e α+1 − e α , features distinctively different shapes for quantum integrable and quantum chaotic systems. While for integrable systems, the NNLS have been predicted by Berry and Tabor [50] to feature an exponential (or Poissonian) distribution P P (s) = exp (−s), for chaotic systems it closely follows random matrix theory [43]. In our case of a time-reversal symmetric system, the corresponding random-matrix ensemble is the Gaussian orthogonal ensemble (GOE) which has been shown (see e.g., [46]) to closely follow the Wigner-Dyson distribution (or Wigner surmise) given by A complementary spectral measure first proposed by Gurevich and Pevzner [80] and applied to quantum chaos [81,82] has the advantage that it does not require spectral unfolding but can be applied to the spectral raw data, i.e., the restricted gap ratios r α where . The distribution of restricted gap ratios has been shown to obey for 3 × 3 GOE matrices the analytical prediction For chaotic systems, this prediction remains very accurate even for large systems [82]. In the limit of quantum integrable systems, the distribution of restricted gap ratios is given by (see [81]) The search for generic spectral measures for the transition regime between the quantum integrable and quantum chaotic limit has remained an open problem. For systems possessing a classical counterpart with a mixed phase space in which integrable and chaotic motion coexist, several models for the NNLS have been proposed [83][84][85][86][87]. Empirically, one of the best fits to spectral data for mixed systems has been provided by a heuristic ansatz suggested by Brody [88] which allows for a one-parameter smooth interpolation of the NNLS distribution in the transition region between the quantum integrable and quantum chaotic limit, where the Brody parameter γ characterizes the transition from the integrable (γ = 0) to the chaotic limit (γ = 1) and b follows from the normalization as The Brody parameter can be viewed as a measure of the strength of level repulsion between neighboring levels of the quantum system. For mixed few-degrees of freedom systems with a classical analogue, γ could be identified as a measure for the chaotic fraction of classical phase space [86,89]. Moreover, γ has also been found to be directly proportional to the degree of phase-space (de)localization of eigenstates as measured by their Husimi distribution [77]. The parameterization of the transition from quantum integrability to quantum chaos in terms of a variable exponent γ has the salient feature that even for very small but finite γ, 0 < γ 1, P B (0) = 0, reflecting the fact that any perturbation of quantum integrability immediately causes level repulsion and suppresses the probability density for any exact degeneracy. We recall that non-degeneracy is one of the key prerequisites of von Neumann's quantum ergodic theorem [32]. We further note that the Hasegawa distribution [84] sometimes provides an even more accurate fit to the NNLS distribution (see, e.g., [90]), however, at the price of a second adjustable parameter.
To determine γ, we fit Equation (12) to the data for P(s) (Figure 2). The quality of the fit is evaluated through the χ 2 -function which measures the deviation of the distribution of nearest-neighbor spacings P(s) from the Brody distribution P B (s) using a bin size of ∆s. As an additional measure for the uncertainty of γ, we use the fact that the Brody parameter can be alternatively determined from a fit to the integral ds P(s ) rather than to P(s) itself. The small differences found between the two fits can be used as a measure for the numerical error. For W BB = 0 and for the linear tilt of the external potential V(j) (Equation (5)), we observe an excess of (near-) degenerate states as compared to the prediction of the exponential (Poisson) distribution in the first bin at s = 0 with ∆s = 0.01. This hints at the presence of an only weakly broken symmetry which disappears when using a quadratic tilt. For reasons of consistency, we employ for all W BB a linear tilt in the following. Neglecting the first bin in the fitting procedure for W BB = 0, we obtain γ = 0.005 and, overall, a very good agreement with the Poisson distribution ( Figure 2a). As the intra-bath interaction is varied from W BB = 0 to W BB = 1, we observe a continuous transition from a near Poissonian to an approximate Wigner-Dyson NNLS distribution (Figure 2a-c). The Brody parameter monotonically increases from γ 0.005 at W BB = 0 to γ 0.9 at W BB = 1. We note that after reaching a plateau at γ 0.93 near W BB = 3, the Brody parameter decreases again for W BB > 5 and vanishes in the strongly correlated limit of W BB 1. The decrease of the Brody parameter for large W BB results from clustering of the energy spectrum in the strongly interacting regime. The bath fragments into clusters of particles with the interactions between separate clusters suppressed. Thus, a partially ordered system emerges reducing the degree of quantum chaoticity. We will focus in the following on the parameter range W BB ≤ 1 within which the transition from a nearly quantum integrable to a nearly fully quantum chaotic system occurs.
For the two limiting cases of quantum integrability (W BB → 0) and quantum chaos (W BB → 1) of the present Fermi-Hubbard system, we can also apply the predictions for the restricted gap ratio distribution (Equations (10) and (11)) . We find for these two limiting cases very good agreement between the prediction and the data (Figure 3), confirming that the identification of quantum integrability and quantum chaos is independent of the particular choice of the spectral measure.   (10)) and for integrable spectra (Equation (11)).
For the first moment of the restricted gap ratio distribution, we find r = 0.5284 for W BB = 1 agreeing to within 0.5% with the GOE expectation value for asymptotically large matrices r GOE = 0.5307 [82]. Conversely, for W BB = 0, we find r = 0.3811 in very good agreement with the prediction for a Poisson distribution r P = 0.3863. As there is presently no interpolation function W(r) available for the transition between the quantum integrable limit (Equation (11)) and the quantum chaotic limit (Equation (10)), we will focus in the following on the Brody distribution for the NNLS as spectral measure for the transition regime.

Measures for Wavefunctions
As an alternative to spectral measures, on can also explore and quantify chaos through the complexity of the eigenstates. According to Berry's conjecture, the eigenstates of a chaotic system feature randomly distributed amplitudes over an appropriate basis, e.g., in quantum billiards, they correspond to randomly distributed plane waves [75]. Following this conjecture, a large number of such measures have been proposed. They include the statistical distribution of eigenvectors [9,12,46,91], the configuration-space probability distribution [92], the configuration-space self-avoiding path correlation function [52], the Wigner function-based wavefunction autocorrelation function [75], the inverse participation ratio [93], the Shannon entropy [94], and the phase space localization measured in terms of the information entropy encoded in the Husimi distribution [77]. One limitation for the quantitative significance of most of these measures (with the possible exception of [77]) is their dependence on the chosen basis of representation. For systems that can be continuously tuned from integrable to chaotic, the eigenstates of the integrable limit suggest themselves as a convenient basis to monitor the transition to chaos [9,12,58]. For manybody systems, the eigenstates of the mean-field Hamiltonian often provide the reference basis for measuring quantum chaoticity [14]. In the following, we use the eigenstates |ψ 0 α of the integrable system with W BB = 0 as a basis for determining the statistical distribution of eigenvectors. From the amplitudes c α α = ψ 0 α |ψ α and probabilities |c α α | 2 , we calculate the Shannon entropy [94] for each eigenstate |ψ α We observe that for W BB = 1, the Shannon entropy as a function of E α forms an inverted parabola-like function with remarkably small eigenstate-to-eigenstate fluctuations (Figure 4). At the apex near the center of the spectrum, S α reaches a maximum S max close to the GOE limit S GOE ≈ ln 0.48d H [58] with d H the dimension of the Hilbert space (Figure 4a). States in the tails of the spectrum show strong deviations from this limit as the eigenstates in this region are less complex and do not fulfill the ETH [58]. Best agreement with GOE predictions can therefore be expected near the center of the spectrum at α ≈ d H /2 with the highest density of states. For smaller W BB (Figure 4b-d), the Shanon entropy reveals a significantly diminished complexity of the eigenstates indicated by a reduced S max and, at the same time, drastically increased state-to-state fluctuations. Probing the generic features of the wavefunctions, we will use the dependence of the scaled Shannon entropȳ as an alternative wavefunction-based measure of quantum chaoticity complementing the Brody parameter γ as spectral measure. Numerically, we determine S max by averaging over small intervals of energy and calculating the maximum of the resulting smooth curve. Empirically, we find that the dependence of the Brody parameter γ, i.e., the degree of quantum chaoticity on the interaction parameter of the bath particles, γ(W BB ) ( Figure 5) can be accurately approximated by with γ 0 = 0.88 and W 0 BB = 0.15. While a monotonic increase is intuitively expected, the origin of this particularly simple functional form remains to be understood. Remarkably, the evolutions of γ andS as a function of W BB closely mirror each other, thereby representing two independent measures of the degree of quantum chaoticity during the transition from integrability to chaos. Overall, the agreement between γ andS is very good. Residual differences can be viewed as a measure for the residual uncertainty in the quantitative determination of the degree of the eigenstate quantum chaoticity.  Figure 4 and correspond to the scaled standard deviation around S max . Error of the fit to the Brody distribution as measured by the square root of the χ 2 function (Equation (14)) (gray line and right y-axis).

The Reduced Density Matrix of the Impurity
The impurity embedded in the Fermi-Hubbard system serves as a "thermometer", i.e., as a sensitive probe of the thermal state of the interacting many-body system. We aim at exploring the emergence of thermal properties of the impurity when the entire (subsystem and bath) system is in a given pure and stationary eigenstate ofĤ with energy E α and vanishing state entropy (or von Neumann entropy S vN = 0). Such an isolated large quantum system can be viewed as the limiting case of the quantum microcanonical ensemble where the width of the energy shell ∆E vanishes, i.e., ∆E → 0. Unlike other approaches, it does not invoke any coarse-graining over a macroscopically small but finite width of the energy shell nor any random interactions. For such a quantum system without any a priori built-in statistical randomness, we pose the following question: Starting from a given isolated eigenstate of the entire system, under which conditions will the reduced density matrix of the impurity correspond to a canonical density matrix, i.e., the thermometer will be accurately represented by a Gibbs ensemble or, for short, be in a Gibbs state? If such a thermal state emerges, what will be its temperature T, or its inverse temperature β = 1/k B T? We refer to this process as emergence of a thermal equilibrium state rather than the frequently used term "thermalization" as the latter (implicitly) implies a time-dependent approach to an equilibrium state starting from an out-of-equilibrium (statistical or pure) initial state that represents a coherent superposition of different energy eigenstates. We neither invoke any ensemble average over states from the microcanonical energy shell of finite thickness ∆E nor do we invoke wave packet dynamics of a nonstationary state of the entire system.
For finite isolated systems, in particular, systems with a bounded spectrum such as the present Fermi-Hubbard model, the extraction of proper thermodynamic (or thermostatic) variables from the microcanonical ensemble requires special care. As has been recently demonstrated [95,96], the alternative definitions of the entropy used as the fundamental thermodynamic potential for the microcanonical ensemble yield, in general, inequivalent results. The standard definition [97] attributed to Boltzmann with Ω(E) the DOS of the entire closed system, implies an inverse temperature that may violate certain thermodynamic relations for mesoscopic systems with a bounded spectrum [95,96]. As shown more than 100 years ago [97,98], the Gibbs entropy defined by results in an inverse temperature that is free of such inconsistencies. From Equations (19) and (21), it follows that the two inverse temperature definitions are interrelated through the specific heat C [95] β Boltzmann = (1 − k B /C)β Gibbs (22) with C = (∂T Gibbs /∂E) −1 and T Gibbs = β −1 Gibbs /k B . Only for systems with a small specific heat of the order of k B or smaller, differences between β Boltzmann and β Gibbs become noticeable. This is in particular the case for systems with a bounded spectrum. While β Boltzmann (E) features negative values as soon as the density of states Ω(E) = N (E) decreases (Equation (19)), β Gibbs (E) remains always positive semi-definite (Equation (21)). Figure 6 presents a comparison between β Gibbs and β Boltzmann for the present Fermi-Hubbard model with an impurity where we have applied the microcanonical thermodynamic relations for β Boltzmann and β Gibbs (Equations (19) and (21)) to the numerically determined spectral data (Figure 1) of the entire system over a wide range of energies E. The two inverse temperatures closely follow each other in parallel with β Gibbs shifted upwards relative to β Boltzmann as long as Ω (E) > 0. For larger E when β Boltzmann turns negative, the discrepancies increase as β Gibbs remains positive for all E.
Alternatively, the entire system can be assigned an inverse temperature β c by treating the system as a canonical ensemble. Accordingly, the energy E can be expressed in terms of the canonical expectation value where Z c = Tr exp (−β cĤ ) is the canonical partition function andĤ is the Hamiltonian of the entire system (see Equation (1)). For a given E, Equation (23) yields an implicit relation for β c also shown in Figure 6. Obviously, for this finite system, β c is close to β Boltzmann . In the thermodynamic limit, we would expect β c = β Boltzmann . In spite of the fact that the size of our system is still far from the thermodynamic limit (N → ∞), the agreement between different thermodynamic ensembles is already remarkably close. Deviations appear primarily near the tails of the density of states and are larger in the region of negative β Boltzmann where the DOS decreases rather than increases with E.  (20), blue) as well as the canonical expectation value (Equation (23), dashed black). The energy is restricted to the interval [E min , E peak + E FWHM /2, ] with E min the lower bound where the DOS of the entire system is ≥15% of its peak value at E peak , and E FWHM the full-width-at-half-maximum of the DOS. Bath-bath interaction strength W BB = 1 and impurity-bath interaction W IB = 1 (see Figure 1b).
The conceptually interesting question now arises which of these temperatures, if any, will be imprinted on the impurity upon an exact calculation of its reduced density matrix by tracing out all bath degrees of freedom from a given single exact eigenstate of a the isolated many-body system, and without invoking any a priori assumption of the microcanonical ensemble.
To address this question, we start from the density operator for any pure energy eigenstate |ψ α of the entire system given by the projector |ψ α ψ α |. Consequently, the reduced density matrix (RDM) of the impurity follows from tracing out all bath degrees of freedom, which will, in general, depend on the parent state |ψ α it is derived from. We explore now the generic properties of D (I) α independent of the particular parent state. Specifically, we investigate whether a given D yielding natural orbitals |η m,α with natural occupation numbers n m,α [99]. We emphasize that within the present approach, the RDMs D (I) α and their eigenvalues, the occupation numbers n m,α , which characterize the thermal state, are a priori uniquely determined and not influenced by the choice of any (approximate) basis. Compared to previous investigations, this is one distinguishing feature of the present study of the thermal state emerging from an isolated deterministic many-body system. RDMs have been previously employed in studies of disordered fermionic systems [100][101][102].
m,α = η m,α |Ĥ I |η m,α , which, in turn, should be close to the eigenstates ofĤ I . Moreover, the resulting value for β extracted from the fit to the exponential distribution allows the identification of the inverse temperature uniquely characterizing the thermal distribution. For a finite-size system with an impurity and a bath with an order of magnitude of 10 particles and finite impurity-bath coupling, the residual interaction of the impurity with the bath is not negligible and should therefore be included to improve the numerical accuracy. We account for the residual impurity-bath interaction on the level of the meanfield (MF) or Hartree approximation [14]. Accordingly, the energies (I) of the impurity appearing in the Boltzmann factor include a correction term where the MF interaction operator in site-representation reads the reduced one-body density of residual bath particles at the site j when the entire system is in state |ψ α . In Equation (28), the partial trace over all but one (N B − 1) bath particles and the impurity (I) is denoted by Tr N B −1,I . The energy fluctuations provide a measure for the proximity of the natural orbitals of the RDM to the eigenstates of the (perturbed) single-particle Hamilton operator of the subsystem,Ĥ I,eff =Ĥ I +Ŵ (IB) MF,α . The energy fluctuations (Equation (29)) vanish only when the natural orbitals |η m,α with which the matrix elements in Equation (29) are evaluated do coincide with the eigenstates ofĤ I,eff . Therefore, the variance ∆¯ (I) m,α can serve as a distance measure of the natural orbitals from eigenstates of the impurity Hamiltonian operator. The MF correction in Equation (27) follows from the Liouville-von Neumann equation for the reduced system where the interaction with the bath consists of the MF term and a collision operator. The collision operator describes the correlations between the impurity and the bath particles and contains the so-called two-particle (subsystem-bath) cumulant ∆ 12 . We numerically monitor the validity of the MF approximation through the magnitude of the two-particle correlation energy determined by ∆ 12 . Consistently, we find that for all many-particle states |ψ α which reduce to a near-canonical RDM for the impurity, the correlation energy is negligible compared to the MF energy thereby justifying Equation (26). Of course, in the limit of weak impurity-bath coupling, the MF correction (Equation (27)) becomes negligible as well.
A representative example for the spectrum of the impurity RDM, i.e., the occupation number distribution of natural orbitals of the impurity RDM emerging from a single energy eigenstate of the entire system with state index α = 4364 (with α sorted by energy) and energy eigenvalue E α = −2.396 lying on the tail of the DOS with positive β for W BB = 1, is shown in Figure 7.
Indeed, a Boltzmann distribution ∝ e −β¯ (I) m,α characterizing the canonical density matrix is observed. Moreover, the fit to an exponential yields β ≈ 0.58 in close agreement with β Boltzmann = 0.58 predicted by Equation (19) for the inverse temperature within the microcanonical ensemble (see also Figure 6) and reproduces the distribution of occupation numbers very well. It also agrees with β c predicted by Equation (23) where the entire system is treated as a canonical ensemble. We note that the Boltzmann-like decay of the diagonal elements would remain qualitatively unchanged when neglecting the MF correction in Equation (26) but the fit to β would deteriorate. Thus, from the reduction of state α = 4364, we have verified that a canonical density matrix emerges.  (29)). The blue solid line corresponds to the best exponential fit yielding the exponent β ≈ 0.58 in agreement with β Boltzmann deduced for this state from Equation (19). The inset shows the same plot on a logarithmic scale.
On a conceptual level, the present results confirm the analysis by Dunkel and Hilbert [95] who showed that the recently observed experimental single-particle population distribution in an isolated finite cold-atom system [22] is governed by β Boltzmann . Thus, the canonical density matrix of a small system emerging from tracing out bath variables is characterized by the inverse Boltzmann temperature β Boltzmann rather than by β Gibbs . Consequently, level inversion in a small system in thermal contact with a bath, in particular, spin systems [103,104], can be properly characterized by negative β Boltzmann . The point to be noted is that while β Boltzmann describes the canonical density matrix, the use of β Gibbs is required for consistency in thermodynamic relations such as the Carnot efficiency [95,96]. In the following, we present the numerical results for the canonical density matrix of the impurity in terms of β Boltzmann which we denote, from now on, for notational simplicity by β. We point out that β can be straightforwardly transformed into β Gibbs using Equation (22) and that none of the conclusions to be drawn in the following are altered by this transformation.

Eigenstate Canonicity and Degree of Quantum Chaoticity
The demonstration of the emergence of a canonical density matrix from a particular eigenstate |ψ α (α = 4364) of the entire system invites now the following questions: Is the reduction to a canonical density matrix generic, i.e., will it emerge for almost all |ψ α ? Is this appearance related to the quantum chaos present in the underlying many-body system? On a more quantitative footing: For how many of the eigenstates will a canonical density matrix emerge and does this number depend on the degree of quantum chaos of the system?
We explore these questions by determining the fraction of many-body eigenstates reducing to a canonical density matrix of the impurity, referred to in the following as eigenstate canonicity, as a function of the exact total energy E α for the complete set of eigenstates α of the entire system and for varying bath-bath interaction W BB . The corresponding degree of quantum chaoticity of the entire system is measured by either the Brody parameter (Equation (12)) or the Shannon entropy (Equation (16)). Striking differences in the approach to the thermal state with inverse temperature β appear which are controlled by the Brody parameter γ (or Shannon entropyS): At W BB = 1, when the system is chaotic as indicated by a Brody parameter γ ≈ 0.9 (or scaled Shannon entropȳ S = 0.9), a thermal distribution with a well-defined inverse temperature β, consistent with the (micro)canonical ensemble prediction (Equations (19) and (23)), emerges for an overwhelming fraction of states with the exception of states in the tails of the spectrum where the DOS is strongly suppressed (Figure 8a). The large deviations in the tails are consistent with the corresponding deviations ofS in the same spectral region (Figure 4a). With decreasing W BB and, correspondingly, decreasing γ orS, an increasing fraction of states yields values of β that are far from the thermal ensemble prediction. Moreover, the quality of the fit to a canonical density matrix measured by the variance of ∆β and indicated by the color coding of Figure 8 drastically deteriorates. In other words, for a significant fraction of states, the emerging RDMs do not conform with the constraints of a canonical density matrix.
In order to quantify the decomposition of the Hilbert space into the subspace of states |ψ α whose reduction to the subsystem yields a canonical density matrix and into the complement whose reduction fails to yield such a thermal state, we introduce a threshold for the variance of the inverse temperature ∆β th above which we consider the eigenstate canonicity to be failing. We then calculate for all states |ψ α the fraction of emerging canonical density matrices satisfying ∆β ≤ ∆β th . Of course, the resulting fraction of states will depend on the precise value of ∆β th chosen. We have determined these fractions for thresholds ranging from ∆β th = 5 × 10 −3 to 1.5 × 10 −2 . Changes of the fractions due to variation of ∆β th are indicated by the vertical error bars in Figure 9. An unambiguous trend of a monotonic increase of the fraction of canonical density matrices with chaoticity is emerging, obviously unaffected by the choice of ∆β th . This fraction representing Gibbs states, denoted in the following by G, monotonically increases with quantum chaoticity as parameterized by either the Brody parameter γ, G(γ), or alternatively by the scaled Shannon entropy, G(S) (Figure 9). Since γ andS both increase monotonically with the bath interaction W BB (see Figure 5), this implies also a monotonic relationship G(W BB ). The conceptually important observation emerging from Figure 9 is that the degree of canonicity of the RDM, G(γ), undergoes a continuous transition from the quantum-integrable (γ → 0) to the quantum-chaotic limit (γ → 1). The strength of level repulsion in the NNLS parameterized by γ directly determines the probability of finding the RDM of the impurity represented by a Gibbs ensemble.
The approach of the RDM of the impurity to the Gibbs ensemble with MF,α ] can be also directly observed in the spatial site representation (j 1 , j 2 ) of the RDM of the impurity (Figure 10b).
We illustrate the RDM in the site representation for two energetically nearest-neighbor states (α = 13,637 and α = 13,638) when the system is in the transition regime between integrable and non-integrable (in the present case, W BB = 0.1). We quantify the approach to D Gibbs α through the density matrix site correlation function where j|D (I) α |j is the RDM of the impurity (Equation (24)) in the site basis. While the state α = 13,638 results in a nearly diagonal RDM in the site basis (Figure 10a) with rapidly decaying site correlations closely following the prediction for a Gibbs ensemble (Equation (30)), the adjacent state α = 13,637 yields a RDM with significant off-diagonal entries, extended site correlations, and strong deviations from Equation (30). Thus, the emergence of a thermal density matrix in the transition regime between quantum integra-bility and quantum chaos displays strong state-to-state fluctuations and is not a smooth function of the energy E α .  As quantitative measure for the distance of a given RDM from the Gibbs ensemble, we use the trace-class norm with ||M|| 1 = Tr √ M † M the largest of the Schatten p-norms (p = 1). For hermitian positivesemidefinite matrices of unit trace, the Schatten 1-norm is bounded by 0 ≤ ||M 1 − M 2 || 1 ≤ 2. We observe for RDMs derived from all eigenstates of the entire system ( Figure 11) an overall reduction of distances ∆D α from a canonical density matrix with increasing W BB . For W BB = 1 (Figure 11a) in the (near) quantum chaotic limit, the vast majority of impurity RDMs have a distance of 0.15 from an ideal Gibbs ensemble (apart from those reduced from many-body states in the tail regions of the spectrum with low DOS). The distribution of ∆D α mirrors the distribution of Shannon entropies (Figure 4). We note that for the present finite quantum system, we find that the distance measured by the Schatten 1-norm has a lower bound of ∆D α 0.05. As the Schatten 1-norm is sensitive to small deviations in both diagonal and off-diagonal elements, these deviations are due to residual fluctuations (Equation (29), Figure 7) of the natural orbitals of the impurity which are expected to vanish in the thermodynamic limit N → ∞. Indeed, plotting the value of the smallest distance (∆D α ) min as a function of the dimension of the Hilbert space of the system d H for three numerically feasible system sizes indicates that the minimal distance vanishes in the thermodynamic limit as d −1/3 H (inset Figure 11a). With decreasing W BB , e.g., W BB = 0.1 in Figure 11b, the mean distance of RDMs from a Gibbs state significantly increases and, moreover, the spread becomes much larger reflecting, again, the behavior of the Shannon entropy (Figure 4c). The emergence of canonical density matrices, i.e., of Gibbs states for almost all |ψ α in the quantum chaotic limit (γ 1 orS 1) can be viewed as a rather specific manifestation and extension of the ETH [5][6][7][8]. The local observable in this case is the RDM of the impurity, D (I) α , itself. Its diagonal elements are, indeed, a smooth function of the total energy E α as predicted by ETH but now, more specifically, Boltzmann-distributed ∝ e −β α¯ (I) m,α over impurity states with the inverse temperature imprinted by E α . The present analysis covers, in addition, also the transition regime between quantum integrability and quantum chaos (0 < γ < 1) where, in general, the ETH does not apply. A canonical density matrix may still emerge but now only for a decreasing fraction of eigenstates of the finite large system. The size of this fraction G is predicted by the degree of quantum chaoticity as measured by the Brody parameter γ or Shannon entropyS (Figure 9).
The direct relation between the emergence of the canonical density matrix for a small subsystem from eigenstate reduction and the quantum chaoticity of the large system it is embedded in, established here for a finite quantum system, raises the conceptual question as to the extension of this connection to the thermodynamic (N → ∞) limit. Clearly, this question cannot be conclusively addressed by the present method of exact diagonalization. Nevertheless, we can provide evidence to this effect by exploring the scaling with system size still within computational reach. We first establish that the degree of quantum chaoticity as measured by the Brody parameter γ (or the Shannon entropȳ S) indeed increases with system size at fixed strength of the interaction W BB that breaks quantum integrability (Figure 12). . The observed increase of quantum chaoticity with system size is qualitatively in line with properties of classical chaos: In a mixed phase space with surviving local regular structures such as tori, their influence on phase space dynamics is rapidly diminishing with increasing phase space dimension a prominent example of which is Arnold diffusion [3,4]. This increase of quantum chaoticity with system size at fixed interaction strength turns out to be key for the emergence of a universal, i.e., (nearly) system-size-independent, interrelation between the fraction of canonical eigenstates and the degree of quantum chaoticity. Both the Brody parameter γ as well as the fraction of density matrices complying with the Gibbs ensemble increase with system size at fixed bath interaction strength. As a consequence, a near universal, i.e., size-independent, relation G(γ) between the fraction of (approximate) Gibbs states and the degree of quantum chaos as measured by γ emerges (Figure 13).
The data for different combinations of values of W BB and M s fall on the same curve. A very similar relation would emerge for G(S) as a function of the scaled Shannon entropy. We have thus established the remarkable feature that the fraction of canonical eigenstates, i.e., the likelihood that a subsystem is in a Gibbs state when the large but finite system is in a pure energy eigenstate with zero von Neumann entropy is controlled and can be tuned by γ (orS) and, in turn, by the degree of level repulsion in the quantum many-body system which is controlled by γ.

Conclusions and Outlook
In this work, we have explored the emergence of a thermal state (or Gibbs ensemble) of a small (sub)system in contact with a bath when the combined large but finite deterministic quantum system is isolated and in a well-defined energy eigenstate. As prototypical case, we have considered an impurity embedded in an interacting spin-polarized Fermi-Hubbard many-body bath which facilitates a clear-cut subsystem-bath decomposition and a tunable transition of the entire system from quantum integrability to quantum chaos. By tracing out the bath degrees of freedom, we have investigated how many of the resulting reduced density matrices of the subsystem represent a canonical density matrix. We have shown that the probability for finding a canonical density matrix monotonically increases with the degree of quantum chaos. The degree of quantum chaos is identified here by both the energy-level statistics as well as by the randomness of the eigenstates as measured by the Shannon entropy. The likelihood for the emergence of thermal states is thus found to be controlled by the degree of quantum chaoticity as parameterized by the Brody parameter or the Shannon entropy. Even though our simulations are limited to finite-size systems, the present results for varying system sizes suggest that the relation between the fraction of eigenstates of the isolated many-body system whose reduction to a small subsystem yields a reduced canonical density matrix and the degree of quantum chaoticity is universal, i.e., size-independent. Each many-body eigenstate represents the fine-grained version of the energy shell of the microcanonical ensemble of the entire impurity-bath system. This connection between the fraction of canonical eigenstates and quantum chaoticity thus offers a direct quantum analogue to the role of classical chaos which Boltzmann invoked in deducing the classical (micro-)canonical ensemble. One can view this as an example of classicalquantum correspondence to this cornerstone of the foundation of statistical mechanics. The statistical ensemble properties can already emerge for isolated energy eigenstates without invoking any randomness, e.g., coarse-graining over a macroscopically thin energy shell or superposition of many eigenstates of the isolated large system as frequently employed. The emergence of statistical ensemble properties from the reduction of pure states was already early anticipated by Landau and Lifshitz [42] and later related to quantum chaos [14]. The present study establishes a direct quantitative relationship between the degree of canonicity and the degree of quantum chaos, in particular, also covering the transition regime from quantum integrability to quantum chaos.
The present results are also expected to have implications for the topical issue of thermalization in finite quantum systems [24,26,27,30]. In this paper, we intentionally avoided this notion and, instead, focused on thermal equilibrium states as we deduce the canonical density matrix from stationary energy eigenstates bypassing any explicit time dependence of the dynamics. Thermalization of an initial non-equilibrium state is, by contrast, a fundamental probe of the time evolution of quantum many-body systems. Up to now, one primary focus has been on quantum quenches, the relaxation of out-off equilibrium initial states. Their time evolution has typically shown a transition from an exponential decay for weakly perturbed many-body systems to a Gaussian decay in the strongly coupled limit, however, without an unambiguous correlation to quantum chaos [105,106]. The Shannon entropy was found to linearly increase with time before reaching saturation [107]. For disordered systems, an initial rapid decay followed by a slow power-law relaxation of occupation numbers has been observed [66,102]. The extension of the present study to a non-equilibrium initial state of a deterministic many-body system would yield the time evolution of the entire one-body RDM, and its eigenvalues and eigenvectors, the time dependence of which remains to be explored. Moreover, the dependence of the relaxation dynamics of the RDM on the choice of the initial state for systems in the transition regime between quantum integrability and quantum chaos (i.e., for intermediate values of the Brody parameter γ) is of particular interest. Most importantly, will quantum chaos play an analogous role for the process of mixing as classical chaos does for classical non-equilibrium dynamics and the relaxation to equilibrium? The origin and properties of such "quantum mixing" remain a widely open question.

Data Availability Statement:
The presented data can be obtained from the corresponding author upon a reasonable request.