Measurement-Induced Dynamical Quantum Thermalization

1. Introduction

Thermodynamics, as a phenomenological theory, has proven to be most useful since its inception. It was constructed to describe macroscopic systems in contact with each other, where concepts like heat exchange, entropy, and temperature are defined in an operative way. In this context, thermodynamic equilibrium emerges empirically as the state that is reached after a sufficiently long time when a system is left to itself by maximizing its entropy. Such an equilibrium state can be realized by heat exchange with a thermal heat bath at a given temperature or in an isolated, classical system at constant energy by multiple scattering of the microscopic entities comprising the macroscopic system. In an isolated quantum system, however, the problem of thermalization is non-trivial, because quantum-mechanical time evolution of isolated or Hamiltonian systems is unitary: A statistically pure quantum state remains pure and, therefore, its entropy is equal to zero for all times. More generally, in an isolated quantum system, the total entropy is constant in time and cannot be maximized. In sharp contrast, experiments on trapped, ultracold atomic gases, which are nearly ideally isolated from the environment, do show exponential relaxation to thermal behavior. This constitutes the problem of quantum thermalization.

One of the most prominent scenarios for understanding thermal behavior of isolated quantum systems is the eigenstate thermalization hypothesis (ETH), put forward independently by Deutsch and Srednicki [1,2]. It asserts that in an isolated many-body system, the long-time average of the quantum-mechanical expectation value $\langle \Psi \left(t\right)|\widehat{A}|\Psi \left(t\right)\rangle$ of an observable quantity $\widehat{A}$ with respect to a statistically pure state $|\Psi \left(t\right)\rangle$ is indistinguishable from its expectation value with respect to an energy eigenstate at a typical energy and, thus, indistinguishable from the microcanonical thermal average. The ETH became a focus of quantum research after its numerical verification in hard-core boson gases in two-dimensional lattices [3,4], although it was found not to be valid in other systems [5,6]. The applicability of the ETH relies on a number of restrictive preconditions. Most notably, the energy distribution of the Hamiltonian eigenstates comprising the pure quantum state $|\Psi \left(t\right)\rangle$ must be narrow compared to its mean energy, in order to fulfill the microcanonical condition. In addition, the off-diagonal elements of $\langle \Psi \left(t\right)|\widehat{A}|\Psi \left(t\right)\rangle$ in the energy eigenbasis must effectively vanish. This, however, occurs generally only on the long-time average by destructive interference, which is slow due to the narrow frequency spectrum, in contrast to the experimentally observed, fast thermalization dynamics of ergodic systems. See Ref. [7] for a detailed discussion and the relation to ergodicity. Typically, the ETH can be validated only for local observables (correlation functions) [8].

In the present work, we explore the thermalization mechanism of measurement-induced, dynamical (heat) bath generation (DBG), put forward in Refs. [7,9]. It is based on the fact that in a sufficiently complex system1 it is not possible to completely determine the state vector, i.e., to measure all quantum numbers of the system by any realistic experiment. Thus, the measurement subdivides the Hilbert space into an observed and a non-observed sector, such that the latter acts as a canonical reservoir for the observed subsystem. As a consequence, the entanglement entropy of the reduced state vector will be maximized even though the total system remains in a pure state. This concept was introduced by Goldstein et al. by the term canonical typicality [10], building on earlier ideas of intrinsic decoherence in a quantum system by Zurek [11,12,13] and by Milburn et al. [14]. Rigorous results on thermalization due to canonical typicality and its relation to the ETH in discrete, translationally invariant systems were proven in Ref. [15]. The work of Posazhennikova et al. [7,9] extended this concept of typicality from stationary, thermal states to the time-dependent approach to equilibrium, then dubbed as DBG, albeit the calculations were done for overall grand-canonical systems, i.e., effectively open quantum systems.

We examine the validity of the predictions of the DBG by numerically exact time evolution of closed quantum many-body systems, where the total energy and particle number are conserved, and monitoring the temporal approach to thermal equilibrium behavior. As a generic, non-integrable system, we consider a Bose gas of up to $N=25$ particles, confined in a trap comprised of 5 single-particle levels and with an arbitrarily strong two-particle repulsion U. For this system, the natural observables are the (time-dependent) occupation numbers of the individual single-particle levels as well as the occupation-number correlation functions between different levels. We reach exponentially long evolution times by iteratively exponentiating the time-evolution operator. The exponentiation procedure requires computing the time-evolution operator for the entire Hilbert space, not only the components comprising the state vector at hand, as in the adaptive, time-dependent density-matrix renormalization group (t-DMRG) [16]. Although this does limit the system size, it turns out that even our few-particle systems of $N\le 25$ are complex enough to show clear thermal behavior after sufficiently long time. This indicates, in particular, that the thermodynamic limit of infinite system size is controlled by the Hilbert space dimension rather than the number of particles. We further find that the entanglement entropy of the observed subsystem with the non-observed one approaches a maximum with a bi-exponential time-dependence, independent of the energy distribution of the initial state. We show that this entropy maximum represents a thermal state by fitting the fluctuation-dissipation relation to the single-particle spectra as well as the level-occupation correlation functions. As an important result, we find that not only the local correlation functions thermalize in this sense, but also the non-local correlation functions of occupation numbers of different levels.

The paper is organized as follows. In the following Section 2 we give a pedagogical description of the partitioning of an arbitrary quantum system by the measurement, characterized by the reduced density matrix, whose entanglement entropy may be maximized although the entire system remains in a pure quantum state. In Section 3 we introduce, as a generic model for the concrete calculations, a trapped, interacting Bose gas and outline our computational method of long-time evolution. Section 4 contains results for the entanglement entropy of the observed subsystem and analyzes its approach to a maximum as a function of time for initial pure states with a narrow (microcanonical) as well as a wide (canonical) energy distribution. In Section 5 and Section 6, we present the thermalization dynamics of the boson occupation numbers in specific single-particle levels and the correlation functions between the occupation numbers of different levels, non-local in level space, respectively. We conclude with a discussion and outlook of the general applicability of our results in Section 7.

2. Measurement and Entropy in Isolated Quantum Systems

The density matrix of an isolated quantum system with the Hamiltonian $\widehat{H}$ reads

$$\begin{array}{c}\hfill \widehat{\rho }\left(t\right)=\sum _{\alpha ,\beta }{P}_{\alpha \beta }|{\psi }_{\alpha }\left(t\right)\rangle \langle {\psi }_{\beta }\left(t\right)|,\end{array}$$

(1)

where the sum runs over the orthonormal basis states $|{\psi }_{\alpha }\left(t\right)\rangle$ of the complete many-body system. In the Schrödinger picture, their time dependence, $|{\psi }_{\alpha }\left(t\right)\rangle =\widehat{U}\left(t\right)|{\psi }_{\alpha }\left(0\right)\rangle$, is given by the time-evolution operator $\widehat{U}={\mathrm{e}}^{-i\widehat{H}t}$ ($\hslash =1$ here and in the following). As a consequence, the von Neumann entropy of the entire system,

$$\begin{array}{c}\hfill S\left(t\right)=-\mathrm{tr}\{\widehat{\rho }\left(t\right)ln\widehat{\rho }\left(t\right)\}=-\mathrm{tr}\{\widehat{U}\left(t\right)\widehat{\rho }ln\widehat{\rho }{\widehat{U}}^{-1}\left(t\right)\}=S(t=0),\end{array}$$

(2)

does not depend on time due to cyclic permutation under the trace. In particular, the entropy of a pure state, whose density matrix has only one non-zero eigenvalue, is zero for all times. The ETH uses the fact that in the energy eigenbasis of the many-body system, $|{\psi }_{\alpha }\left(t\right)\rangle ={\mathrm{e}}^{-i{E}_{\alpha }t}|{\psi }_{\alpha }\left(0\right)\rangle$, the off-diagonal elements of $\widehat{\rho }\left(t\right)$ in Eq. (1) oscillate with the frequency difference of the respective eigenstates. As a consequence, the off-diagonal elements of $\widehat{\rho }\left(t\right)$ average to zero in the long-time integral, while the diagonal elements are constant. In the long-time average, the expectation value of an observable $\widehat{A}$ can then be written in the energy eigenbasis as the sum over diagonal elements,

$$\begin{array}{c}\hfill \langle A\rangle =\mathrm{tr}\left\{\widehat{\rho }\widehat{A}\right\}=\sum _{\alpha }\langle {\psi }_{\alpha }|\widehat{\rho }|{\psi }_{\alpha }\rangle \langle {\psi }_{\alpha }|\widehat{A}|{\psi }_{\alpha }\rangle .\end{array}$$

(3)

With the assumption of ergodicity (equivalence of temporal and ensemble average), this can be seen as a thermal average of the observable $\widehat{A}$, see [1,2,7] for details. However, the off-diagonal elements decay as $1/t$ only after long-time t integration, and a faster diagonalization of the density matrix is possible only on average over all off-diagonal elements. The dynamical approach to equilibrium from an arbitrary initial state is, therefore, unclear in this scenario.

Here we explore the scenario of measurement-induced partitioning. If a quantum system is sufficiently complex, more precisely, the Hilbert space dimension is sufficiently large, it is not possible to completely determine the system’s state vector by any given experiment. The measurement of an observable $\widehat{A}$ then partitions the complete Hilbert space into the subspace spanned by the quantum numbers determined by the measurement, termed system$\mathcal{S}$, and the subspace of not observed quantum numbers which we call bath or reservoir$\mathcal{R}$, as illustrated in Figure 1. The observable $\widehat{A}$ can, thus, be decomposed as $\widehat{A}={\widehat{A}}_{\mathcal{S}}\otimes {\mathbb{1}}_{\mathcal{R}}$ where ${\widehat{A}}_{\mathcal{S}}$ acts in the observed subspace only. Denoting arbitrary orthonormal basis sets of the system and of the reservoir by ${\mathbf{B}}_{\mathcal{S}}=\left\{|{\alpha }_{\mathcal{S}}\rangle \right\}$ and ${\mathbf{B}}_{\mathcal{R}}=\left\{|{\beta }_{\mathcal{R}}\rangle \right\}$, respectively, the complete density matrix reads,

$$\begin{array}{cc}\hfill \widehat{\rho }\left(t\right)& =\sum _{{\alpha }_{\mathcal{S}},{\beta }_{\mathcal{R}},{\alpha }_{\mathcal{S}}^{\prime },{\beta }_{\mathcal{R}}^{\prime }}{c}_{{\alpha }_{\mathcal{S}}{\beta }_{\mathcal{R}}}{c}_{{\alpha }_{\mathcal{S}}^{\prime }{\beta }_{\mathcal{R}}^{\prime }}^{*}|{\alpha }_{\mathcal{S}}{\beta }_{\mathcal{R}}\rangle \langle {\alpha }_{\mathcal{S}}^{\prime }{\beta }_{\mathcal{R}}^{\prime }|,\hfill \end{array}$$

(4)

with the direct product states $|{\alpha }_{\mathcal{S}}{\beta }_{\mathcal{R}}\rangle =|{\alpha }_{\mathcal{S}}\rangle \otimes |{\beta }_{\mathcal{R}}\rangle$, and the sums run over the complete system or reservoir basis, respectively. Since the observable $\widehat{A}$ acts only in the observed subspace, $\widehat{A}={\widehat{A}}_{\mathcal{S}}\otimes {\mathbb{1}}_{\mathcal{R}}$, its expectation value, valid for all times t, can be written as

$$\begin{array}{c}\hfill \langle A\rangle \left(t\right)={\mathrm{tr}}_{\mathcal{S}}\left\{{\widehat{\rho }}_{\mathcal{S}}\left(t\right)\widehat{A}\right\},\end{array}$$

(5)

where the reduced system density matrix and its entanglement entropy read, respectively,

$$\begin{array}{cc}\hfill {\rho }_{\mathcal{S}}\left(t\right)& ={\mathrm{tr}}_{\mathcal{R}}\left\{\widehat{\rho }\left(t\right)\right\}=\sum _{{\tilde{\beta }}_{\mathcal{R}}}\sum _{\begin{array}{c}{\alpha }_{\mathcal{S}},{\beta }_{\mathcal{R}}\\ {\alpha }_{\mathcal{S}}^{\prime },{\beta }_{\mathcal{R}}^{\prime }\end{array}}{c}_{{\alpha }_{\mathcal{S}}{\beta }_{\mathcal{R}}}{c}_{{\alpha }_{\mathcal{S}}^{\prime }{\beta }_{\mathcal{R}}^{\prime }}^{*}\langle {\tilde{\beta }}_{\mathcal{R}}|\widehat{U}\left(t\right)|{\alpha }_{\mathcal{S}}{\beta }_{\mathcal{R}}\rangle \langle {\alpha }_{\mathcal{S}}^{\prime }{\beta }_{\mathcal{R}}^{\prime }|{\widehat{U}}^{-1}\left(t\right)|{\tilde{\beta }}_{\mathcal{R}}\rangle ,\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill S_{\mathcal{S}}\left(t\right)& =-{\mathrm{tr}}_{\mathcal{S}}\{{\widehat{\rho }}_{\mathcal{S}}\left(t\right)ln{\widehat{\rho }}_{\mathcal{S}}\left(t\right)\}.\hfill \end{array}$$

(7)

Importantly, if the Hamiltonian $H=H_{\mathcal{S}}+H_{\mathcal{R}}+H_{\mathcal{SR}}$ couples system and bath, e.g., by particle exchange induced by $H_{\mathcal{SR}}$, cyclic permutation of the time-evolution operator $\widehat{U}\left(t\right)$ in Eq. (6) is not allowed. This is due to the noncommutativity of the action of $\widehat{U}\left(t\right)$ and the incoherent summation over bath states, ${\sum }_{{\tilde{\beta }}_{\mathcal{R}}}$ in Eq. (6). Consequently, for an arbitrary initial state, the entanglement entropy will generally be time-dependent. It is the purpose of the present paper to analyze by numerical evaluations, if and how it may approach a maximum and lead to thermal behavior of physical quantities in the subsystem $\mathcal{S}$. In the following, we will drop the subscript $\mathcal{S}$ on the entanglement entropy, since the time-independent, total entropy will not be further used.

3. Model and Time-Evolution Method

As a concrete model for the numerical time evolution, we analyze a generic, isolated, non-integrable many-body system of N bosons in a trap made up of $M=5$ single-particle states,

$$\begin{array}{c}\hfill \widehat{H}=\Delta \sum _{i=1}^M(i-1)b_i^{†}{b}_i^{\phantom{†}}+J\sum _{i=1}^M\sum _{j\ne i}{b}_i^{†}{b}_j^{\phantom{†}}+U\sum _{i=1}^M{b}_i^{†}{b}_i^{†}{b}_i^{\phantom{†}}{b}_i^{\phantom{†}}+U^{\prime }\sum _{\begin{array}{c}i,j,l,m\\ \mathrm{not}\mathrm{i}=\mathrm{j}=\mathrm{l}=\mathrm{m}\end{array}}^M{b}_i^{†}{b}_j^{†}{b}_l^{\phantom{†}}{b}_m^{\phantom{†}},\end{array}$$

(8)

where U, $U^{\prime }$ denote the intra-level and inter-level interactions, repectively, and J is the transition amplitude between any pair of levels. $\Delta$ is the single-particle level spacing, assumed to be non-degenerate and equidistant, corresponding to a one-dimensional, harmonic trap. In the present paper, we will show results for the parameter values $\Delta /J=10$, $U/J=1$, $U^{\prime }/J=0.1$, and a conserved, total particle number $N=25$. For this parameter set, the average ratio of consecutive level spacings of the many-body eigenergies is $\langle r\rangle =0.53$, close to the value for the Gaussian orthogonal ensemble (GOE) of random matrices of $\langle r\rangle =0.5307$[17,18], indicating indeed chaoticity of our many-body system. We also performed calculations for different parameter values, including $N<25$, and the particular choice is not important for the physics described here. Natural observable quantities for our model system will be the time-dependent, local occupation numbers of the single-particle levels, $\langle n_i\rangle \left(t\right)$, as well as their non-local correlation functions, $\langle n_i\left(t\right)n_j(t+\tau )\rangle$.

To propagate an initial state to sufficiently long times t, necessary to observe the approach of the entanglement entropy to a stationary, maximum value, we use the method of repeated exponentiation of the short-time evolution operator ${\widehat{U}}_0:=\widehat{U}\left(\delta t\right)=exp\left(-i\widehat{H}\delta t\right)$. For a short time interval $\delta t$, ${\widehat{U}}_0$ is expanded up to order $k_{\max }$,

$$\begin{array}{c}\hfill \widehat{U}\left(\delta t\right)=e^{-i\widehat{H}\delta t}=\sum _{k=0}^{k_{\max }}{\displaystyle \frac{{(-i\delta t)}^k}{k!}}{\widehat{H}}^k+\mathcal{O}\left[{\left(\delta t\right)}^{k_{\max }+1}\right],\end{array}$$

(9)

where we choose $\delta t$ such that $\delta t·\max \left(\widehat{H}\right)\lesssim {10}^{-1}$, with $\max \left(\widehat{H}\right)$ the maximum matrix element of $\widehat{H}$ in the entire, large, but finite-dimensional Hilbert space. For $\widehat{U}\left(\delta t\right)$, this yields a precision of $\mathcal{O}\left({10}^{-k_{\max }}\right)$. We usually choose $k_{\max }=4$. To reach long evolution times, ${\widehat{U}}_0$ is recursively exponentiated according to

$$\begin{array}{cc}\hfill {\widehat{U}}_1& =\widehat{U}\left(n\delta t\right)={\left({\widehat{U}}_0\right)}^n,\hfill \\ \hfill {\widehat{U}}_2& =\widehat{U}\left(n^2\delta t\right)={\left({\widehat{U}}_1\right)}^n,\hfill \\ \hfill {\widehat{U}}_3& =\widehat{U}\left(n^3\delta t\right)={\left({\widehat{U}}_2\right)}^n,\mathrm{etc}.\hfill \end{array}$$

(10)

In this way, exponentially long evolution times $\Delta t=n^r\delta t$ are realized at the numerical cost of $\left(r·n\right)$ matrix multiplications, growing only linearly with the number of recursions r. The precision does not significantly deteriorate for the relevant number of multiplications. A desired, high temporal resolution at long times $\Delta t$ can be achieved by applying lower-order time-evolution operators, after this long time has been reached via high-order evolution operators. We specifically use this approach in the calculation of two-time expectation values. In contrast, single-time expectation values are evaluated on a dense time series utilizing smaller steps with $r\approx 10$ recursions of $n=2$. Note that this procedure is possible because the time evolution operator does not depend on the system’s initial state and can be multiplicatively applied, in contrast to adaptive methods like the adaptive t-DMRG. Conversely, for the recursive multiplication, the evolution operators $U_r$ must be calculated for the entire Hilbert space. For our system of up to $N=25$ bosons in $M=5$ single-particle levels, the Hilbert space dimension is $\approx 24000$, and the matrix multiplications are well feasible.

An important choice is that of the initial state. In the present study, we consider only pure states of the entire system $\mathcal{S}\otimes \mathcal{R}$, as those states pose the thermalization problem in the most extreme form: the evolution from a zero-entropy state to thermal behavior at long times. We note in passing that, in this case, the entanglement entropy of the system $S_{\mathcal{S}}$ is identical to the entropy of the rest $\mathcal{R}$ [19]. Among pure states, we investigated two classes: states whose energy spectrum is narrow, $\Delta E\ll \langle E\rangle$, and states with broad energy spectrum, $\Delta E\lesssim \langle E\rangle$. The spectrum is defined as the distribution of expansion coefficients of the initial state $|{\psi }_0\rangle$ into the (discrete) energy eigenbasis via $|{\psi }_0\rangle ={\sum }_E|E\rangle \langle E|{\psi }_0\rangle ={\sum }_E{c}_E|E\rangle$. Narrow-spectrum states are approximate energy eigenstates, i.e., coherent superpositions of many-body eigenstates with nearly equal energy, which may be compatible with the ETH [7]. For the broad-spectrum states, we chose occupation-number eigenstates of single-particle levels, which are not Hamiltonian eigenstates due to the transition amplitude J and the inter-level interaction $U^{\prime }$. Such states do not fulfill the ETH conditions. The spectra of representative narrow-band and broad-band states are shown in Figure 2.

4. Entanglement Entropy

Here we present the entanglement entropy $S\left(t\right)$ for the subsystem $(\mathcal{S}:i\ge 3)$ comprized of the levels $i\ge 3$ of our system (i.e., the occupation numbers of the levels $i=3,4,5$ are the measured observables), as computed using the method described in Section 3.

4.1. Broad energy spectrum

Figure 3 shows the time evolution starting from the occupation-number state with $n_1=N=25$, $n_i=0$ for $i>1$, as the broad-spectrum, non-equilibrium initial state (cf. Figure 2, bottom panel). Clearly, the entanglement entropy of this initial state vanishes, $S(t=0)=0$. Entanglement is then generated on a fast time scale for times below $tJ\approx 1$, followed by a slow, asymptotic approach to a stationary, maximum value attained at $tJ\approx 10$ (see inset of Figure 3). Plotting the deviation from the long-time limit on a logarithmic scale in the inset of Figure 3 reveals two logarithmic slopes, that is, bi-exponential behavior with fitted time constants ${\tau }_{1}J\approx 1/4$ and ${\tau }_{2}J\approx 3/2$, see also Ref. [9].

4.2. Narrow Energy Spectrum

For comparison, we show in Figure 4 the time evolution for a microcanonical state as the initial state, i.e., a state with the narrow energy distribution of Figure 2 (upper panel). The initial entanglement entropy for this initial state is already significantly greater than zero, because each energy eigenstate comprising the microcanonical ensemble extends over the entire system, $i=1,2,3,4,5$. The time evolution of the entanglement entropy in the subsystem $(\mathcal{S}:i\ge 3)$ is essentially static (upper panel). The same is true for the occupation numbers of the individual single-particle levels $n_i$ (lower panel). This static behavior is expected since the microcanonical state contains essentially only a single frequency $\langle E\rangle$. The associated time-dependent phase factors cancel in the system density-matrix ${\widehat{\rho }}_{\mathcal{S}}$ as well as in the observable expectation values $n_i$. Note that, in this microcanonical state, the distribution of level occupation numbers $P\left(n_i\right)$ is not predefined by a temperature but instead fixed by the total energy $\langle E\rangle$ of the system. A state of this nature is well described by ETH. At the same time, such a state is essentially static and cannot describe the dynamic approach to a maximum entropy state, since it effectively contains one single frequency only. By contrast, a state with broad energy distribution is not within the scope of the ETH, while the measurement-induced partitioning into subsystem $\mathcal{S}$ and effective reservoir $\mathcal{R}$ covers both regimes, including the dynamical maximization of the entanglement entropy.

5. Thermalization Dynamics of Local Occupation Numbers and Spectra

5.1. Thermalization Dynamics of Local Occupation Numbers

The approach to a steady state and maximum of the entanglement entropy is a necessary but not sufficient condition for thermalization. In this section, we study how physical observables indeed attain a thermal distribution for sufficiently long times. In our trapped Bose gas, the natural, local observables (i.e., on an individual single-particle level) are the occupation numbers $n_i$ and their energy-resolved spectra. Note that measuring or calculating individual occupation numbers defines a different, single-level subsystem for each measurement, while the other levels are traced out.

Thus, a direct comparison with the entropy evolution in the subsystem $(\mathcal{S}:i\ge 3)$ can be made only for the subsystem occupation number $n_{\mathrm{sys}}=n_3+n_4+n_5$. Its time evolution is shown in Figure 5 (dashed line) along with the occupation numbers of the individual levels, $n_i$. It is seen that for short times, $Jt\simeq 1$, the occupation numbers vary strongly, but settle to stationary values in the long-time limit. The approach of $n_{\mathrm{sys}}$ to the long-time limit, shown in Figure 6 on a logarithmic scale, is exponential with a time scale of $\approx 10Jt$, similar to the entanglement entropy (Figure 3). From Figure 5 we also see that the long-time, individual occupation numbers systematically decrease with increasing level number i. We will demonstrate in the following that this indeed corresponds to a thermal distribution with a single temperature by analyzing the energy-resolved single-particle spectra.

5.2. Long-Time Thermalization of Local Spectra

In order to assess thermalization of physical quantities, it is necessary to compute their energy-dependent spectra and to examine whether they obey a thermal distribution. Here, we are concerned with local, single-particle spectra, while energy-dependent, non-local density correlations will be analyzed in Section 6.

Single-particle spectra are defined using the two-time lesser and greater Green functions (in the notation of Ref. [20]),

$$\begin{array}{cc}\hfill G_{ij}^{<}(t_1,t_2)& =-i\langle b_j^{†}\left(t_2\right)b_i^{\phantom{†}}\left(t_1\right)\rangle ,\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill G_{ij}^{>}(t_1,t_2)& =-i\langle b_i^{\phantom{†}}\left(t_1\right)b_j^{†}\left(t_2\right)\rangle ,\hfill \end{array}$$

(12)

from which the Keldysh and spectral functions, respectively, are obtained in the time domain as,

$$\begin{array}{cc}\hfill G_{ij}^K(t_1,t_2)& =G_{ij}^{>}(t_1,t_2)+G_{ij}^{<}(t_1,t_2),\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill A_{ij}(t_1,t_2)& =i[G_{ij}^{>}(t_1,t_2)-G_{ij}^{<}(t_1,t_2)]\hfill \end{array}$$

(14)

The energy-dependent Keldysh and spectral functions, $G_{ij}^K(t,E)$, $A_{ij}(t,E)$, are then obtained by transforming to center-of-motion (CoM) and relative time coordinates, $t=(t_1+t_2)/2$, $\tau =t_1-t_2$, respectively, and Fourier transforming with respect to the time difference $\tau$,

$$\begin{array}{c}\hfill A_{ij}(t,E)=\int \mathrm{d}\tau A_{ij}(t_1,t_2){\mathrm{e}}^{iE\tau },\end{array}$$

(15)

and similar for $G_{ij}^K(t,E)$. The total Keldysh and spectral functions of the system are then obtained by tracing over the level indices,

$$\begin{array}{c}\hfill G^K(t,E)=\sum _i{G}_{ii}(t,E),A(t,E)=\sum _i{A}_{ii}(t,E),\end{array}$$

(16)

which is independent of the choice of basis taken in level space, i.e., diagonalization is not required. In thermal equilibrium (and thus dropping the CoM time), $A\left(E\right)$ and $i{G}^K\left(E\right)$ are strictly real, positive semidefinite for $E\ge 0$, and obey the fluctuation-dissipation relation formulated for the Keldysh function [20],

$$\begin{array}{c}\hfill G^K\left(E\right)=-iA\left(E\right)coth\left({\displaystyle \frac{E}{2k_{B}T}}\right)=-iA\left(E\right)\left(2n_B(E,T)+1\right).\end{array}$$

(17)

This follows from the fluctuation-dissipation theorem for the equilibrium two-point expectation values in the frequency domain [20] ($\beta =1/k_{\mathrm{B}}T$),

$$\begin{array}{cc}\hfill {\langle \widehat{A}\left(\tau \right)\widehat{B}\left(0\right)\rangle }_E& =e^{-\beta E}{\langle \widehat{B}\left(0\right)\widehat{A}\left(\tau \right)\rangle }_E,\hfill \end{array}$$

(18)

for arbitrary operators $\widehat{A}$ and $\widehat{B}$ in the Heisenberg picture. The single-particle distribution function $n_B(E/T)$ can be extracted from Eq. (17) in terms of the ratio $G^K\left(E\right)/A\left(E\right)$.

To evaluate the two-time expectation values of Green functions in Eqs. (11), (), we write the Heisenberg operators in terms of the time evolution operators, $b\left(t_1\right)={\widehat{U}}^{-1}\left(t_1\right)b\left(0\right)\widehat{U}\left(t_1\right)$, etc., perform an appropriate, cyclic permutation under the trace to cast the time dependence on the states in the density matrix, and then evolve all the many-body state vectors appearing in the expression for $G_{ij}^{\lessgtr }(t_1,t_2)$ (before and after the action of destruction or creation operators) to the desired, long time using the algorithm described in Section 3. A high time resolution with respect to the difference time $\tau =t_1-t_2$ is obtained by subsequently applying time-evolution operators of lower order, ${\widehat{U}}_k$, $k<n$, cf. Eqs. (10). Finally, the trace over bath, then system states is taken for each pair $t_1$, $t_2$ with fixed CoM time t and Fourier transformed with respect to $\tau$ to obtain the spectra at time t. The results for $G^K\left(E\right)$ and A(E) for a generic system of $N=25$ bosons with the broad-spectrum initial state shown in Figure 2, bottom panel, are shown in Figure 7 at a CoM time of $Jt=100$ along with the sum of Lorentzian fits to each of the 5 distinct spectral peaks in $i{G}^k\left(E\right)$ and $A\left(E\right)$. The slightly negative values of $A\left(E\right)$, not allowed in equilibrium, result from numerical unprecisions of the Fourier transform, Eq. (15), and from the fact that the ∞-time limit cannot truly be reached numerically, so that there may be slight deviations from equilibrium. Figure 8 shows the spectral distribution $n_B(E,T)$ function where each data point at energy $E_i$ is extracted from Eq. (17) as $n_B(E_i,T)=(i{\overline{G}}_i^K/{\overline{A}}_i)-1/2$, where ${\overline{G}}_i^K$, ${\overline{A}}_i$, $i=1,\dots ,5$, are the peak weights of the Lorentzian fit for each of the peaks in Figure 7, and $E_i$ are the corresponding peak positions. The solid line is a fit of the thermal Bose-Einstein distribution function $n_{\mathrm{B}}(E/T)=1/[exp(E/T)-1]$. Note that the chemical potential $\mu =0$, since the particle number in the measured subsystem ${\mathcal{S}}_i$ is not fixed but determined by the system dynamics. Thus, the temperature T is the only fit parameter. The numerical results agree quantitatively remarkably well with a thermal distribution for all data points, with only a large error bar for level number $i=2$, resulting from the non-Lorentzian shape of this peak in the numerical result for $A\left(E\right)$. Note that the measurement (or calculation) of each level-occupation number creates a different partitioning of the total system, so that each level may, in principle, have its own temperature. However, for sufficiently long times, all these temperatures attain a global, equal value, as our single-parameter fit shows. See also CoM-time dependent results in Section 6. The value of the fitted temperature is determined by the expectation value of the total energy of the initial state. We conclude that the partitioning of our generic, isolated Bose gas into observed subsystems induces a dynamical approach to thermalization of the occupation number distribution at long times.

6. Thermal Behavior of Non-Local Density Correlations

Our analysis can be extended to observable quantities non-local in level space, that is, occupation-number correlation functions between different single-particle levels, ${\chi }_{ij}(t,E)={\langle {\widehat{n}}_i\left(t_1\right){\widehat{n}}_j\left(t_2\right)\rangle }_E$, where ${\widehat{n}}_i\left(t\right)=b^{†}\left(t\right)b\left(t\right)$ is the occupation number operator on level i. Thermal behavior of non-local correlation functions is usually not obtained in the ETH scenario. Here, we calculate these quantities using the procedure described in the previous section. We probe their thermal behavior by testing them against the fluctuation-dissipation theorem in Eq. (17). Figure 9 shows the comparison of the left- and right-hand side of Eq. (18) for spectra of typical occupation-density correlation functions, ${\chi }_{01}(t,E)$, ${\chi }_{10}(t,E)$, for the same parameter set as in Section 5, using $\beta =1/k_{\mathrm{B}}T$ as the only fit parameter. While for short CoM evolution times t, the deviations are rather large at all energies E (note the logarithmic scale), the agreement at longer times is very good up to $E/J\approx 60$. For even larger energies $E/J$, the density correlations have decayed to such small values that the deviations between both curves become comparable to the omnipresent and necessary thermal fluctuations. Hence, the non-local density correlation functions thermalize as a function of time, induced by measurement-induced partitioning. The extracted temperature of $T/J\approx 200$ of the non-local density correlation agrees well with the temperature of the local level-occupation distribution of Figure 8.

To exemplify the dynamical approach to thermal behavior, Figure 10 shows the temperature as extracted from the fluctuation-dissipation fit (Eq. (18)) for three different correlation functions for a large range of CoM times (solid lines, left axis). For short CoM times, $Jt\lesssim 3$, a single-parameter fit of temperature is not possible, indicating that the measured subsystem is far from equilibrium. For larger CoM times, the uncertainty of the single-temperature fit becomes of the order of $3\%$ (dashed lines, right axis), indicating that a momentary temperature at time t is well defined. However, at long and intermediate times, the solid lines still deviate from each other significantly more than the standard deviation of the fitting procedure (dashed lines) for different correlation functions. While this could indicate that different correlation functions may be at different, individual temperatures, the fitted mean temperature values (solid lines) fluctuate in time more than their instantaneous standard deviations (dashed lines). Within a time average over a characteristic fluctuation time, all three temperatures agree with each other, so that global thermal behavior is reached, but for exceedingly long times. The appearance of different transient temperatures will be a subject of further research.

7. Summary and Outlook

We presented measurement-induced, dynamical bath generation (DBG) [7], as a general thermalization mechanism for isolated quantum systems, which is based on ideas of canonical typicality [10] and intrinsic decoherence [11,12,13,14]. In particular, we analyzed the dynamical approach to thermal behavior by numerical long-time evolution of an interacting, trapped Bose gas as a generic, non-integrable quantum system. The measurement (or calculation) of an observable quantity partitions a quantum system into observed and non-observed Hilbert subspaces in such a way that the states of the non-observed Hilbert subspace, the bath, are traced over like in a grand-canonical reservoir. This tracing leads to a non-vanishing entanglement entropy of the measured subsystem with the bath despite the fact that the state vector of the entire system obeys unitary time evolution and remains a unique state for all times.

Taking a generic, non-energy eigenstate with a broad energy distribution as the initial state of the time evolution, we showed numerically that, in the long-time limit, the entanglement entropy approaches a global maximum in a bi-exponential way: The entanglement entropy reaches a near-maximum value at a rate of about one interaction constant J, and then evolves to the global maximum at a much slower rate. In the long-time limit, the deviations from the global maximum are dominated by thermal, statistical fluctuations. Correspondingly, all quantities considered, local observables as well as non-local correlation functions, reach thermal equilibrium distributions. Their temperatures, as extracted from the dissipation-fluctuation theorem, agree with each other within numerical accuracy. Here, locality and non-locality refer, in a more general sense, to the diagonal and off-diagonal elements of an observable not only in position, but, e.g., in level space.

When the entire system’s state vector is comprised of an energetically narrow spectrum of many-body eigenstates, the entanglement entropy and local occupation numbers have essentially no dynamics as expected for the microcanonical ensemble. In addition to the DBG scenario, the eigenstate thermalization hypothesis (ETH) may apply to this case in that the long-time average of the level occupation numbers is equal to the microcanonical ensemble average for an ergodic system [1,2].

These results allow for a number of general conclusions. The DBG scenario appears to have no conceptual restrictions on its applicability: In contrast to the ETH, it describes the time-dependent approach to thermal behavior. This is true even for an initial state with the lowest possible entropy, a pure state, and works for broad as well as for narrow energy distributions in the initial state. Furthermore, we showed that the DBG is not limited to the thermalization of local (in level space) but also applies to non-local correlation functions. The only condition on the applicability appears to be that the system is sufficiently complex and not integrable, that is, the Hilbert space has sufficiently large dimension and does not factorize into disconnected sectors by the Hamiltonian dynamics. In particular, we conjecture that it is, in fact, the large Hilbert space dimension that ensures a well-defined thermodynamic limit, rather than the particle number. This is because our calculated expectation values are well defined in that their temporal variations are much smaller than the mean value in the long-time, stationary limit, even though the particle numbers $N\le 25$ of our system are not macroscopic. The bi-exponential approach of the entanglement entropy to a global maximum found in our calculations suggests an intermediate, near-equilibrium state. It will be a subject of further research in how far this may be described by time-dependent temperatures of subsystems at intermediate times.