Probing Thermalization through Spectral Analysis with Matrix Product Operators

Yilun Yang , Sofyan Iblisdir, J. Ignacio Cirac, and Mari Carmen Bañuls 1,4 Max-Planck-Institut für Quantenoptik, Hans-Kopfermann-Str. 1, D-85748 Garching, Germany Departament de Física Quàntica i Astronomia & Institut de Ciències del Cosmos, Universitat de Barcelona, 08028 Barcelona, Spain Departamento de Análisis y Matemática Aplicada, Universidad Complutense de Madrid & Instituto de Ciencias Matemáticas, 28040 Madrid, Spain Munich Center for Quantum Science and Technology (MCQST), Schellingstr. 4, D-80799 München, Germany

The study of one-dimensional quantum many-body systems has motivated the emergence of a number of techniques, based on tensor network states (TNS). More concretely, they use matrix product states (MPS) and matrix product density operators (MPDO) [1][2][3][4][5] to approximate the ground states, low-lying excitations, thermal states, as well as time evolution. These methods have enabled the indepth study of a multitude of models and the analysis of relevant physical phenomena.
The success of such techniques is rooted in the ability of MPS and related ansatzes to accurately describe states that fulfill an area law of entanglement [6,7], satisfied (or only slightly violated) by many of the problems mentioned above [8]. There are, however, important open questions that such techniques cannot easily solve. In particular, excited states at finite energy density are difficult to approximate, except in very particular cases [9,10], as they generically display volume law entanglement and, additionally, are embedded in highly dense spectral regions, which severely hinders the convergence of the algorithms. Out-of-equilibrium dynamics is also problematic: under time evolution a volume law often emerges, that makes an MPS approximation inadequate, except for short times. As a consequence, it is virtually impossible for standard MPS techniques to address the fundamental questions of equilibration and thermalization of relatively large closed quantum systems.
A few alternative tensor network algorithms have tried to overcome these problems by avoiding the explicit representation of the states [11][12][13][14]. Although they extend the applicability of the toolbox and allow access to additional dynamical quantities in some scenarios, the fundamental goal of accessing the longtime behavior in a general case, and thus deciding the appearance of equilibration or thermalization, has not been achieved.
Here we introduce a new powerful tool to fill in these gaps. Our method is based on the use of matrix product operators (MPO) to approximate a family of generalized densities of states, and provides a means to directly address thermalization. More concretely, we combine TNS and the kernel polynomial method (KPM) [15] in a general scheme that provides access not only to the full density of states (DOS) of a given many-body Hamiltonian [16], but also to energy functions that are intimately related to the out-ofequilibrium dynamics, including the local density of states (LDOS). With these functions it is possible to probe the eigenstate thermalization hypothesis (ETH) [17,18] across the spectrum, and to verify the thermalization of initial states without explicitly simulating the time evolution.
Generalized DOS.-Let us consider a quantum manybody Hamiltonian with spectral decomposition H ¼ P k E k jkihkj. We are interested in energy functions of the form where O is any operator and δðxÞ is the Dirac delta function. We will aim an approximation to where δ M is a smooth function such that lim M→∞ δ M ¼ δ.  [17,18,[21][22][23], which postulates that, regarding physical observables [24], energy eigenstates look thermal, i.e., they have expectation values close to those of an equilibrium ensemble with a temperature set to get the same mean energy [25]. If ETH holds we thus expect the estimate O M ðEÞ to be a smooth function of energy, and to be equal to the thermal value at the same mean energy, where βðEÞ is the corresponding temperature. Probing this relation constitutes a weak test of ETH. Second, we can use the estimates (ii) and (iii) to approximate the long-time averaged expectation valueŌ ¼ lim T→∞ ð1=TÞ R T 0 dthΨðtÞjÔjΨðtÞi which, if the spectrum is not degenerate, is given by the expectation value in the diagonal ensemble,Ō ¼ O Diag ðΨÞ ≡ P k hkjOjkijhkjΨij 2 . Under the nondegeneracy condition, hkjOjki ¼ OðE k Þ, and we can write If the system thermalizes, the longtime value will be thermal, so we expect for hEi ¼ hΨjHjΨi. Hence it is possible to probe the thermalization of individual initial states without the need to explicitly simulate time evolution. Instead, we can estimate, as we detail in the following, the expectation value of any local observable in the diagonal ensemble for initial states that can be written as an MPS, and compare this result to the expectation value in the Gibbs ensemble for which the mean energy is hEi (which for local Hamiltonians can be efficiently approximated using MPS tools). Notice that if the energy spectrum has degeneracies, it is still possible to estimate the longtime averagedŌ with our method, and perform this comparison, although with a higher computational cost [26].
Finally, the LDOS encodes information about the evolution of a state under H at arbitrarily long times. Indeed, its Fourier transform (assuming H is constant in time) is the survival probability, FðtÞ ≡ jhΨð0ÞjΨðtÞij 2 , which is sensitive to all time regimes of the evolution [27], The decay of the survival probability after a quench presents different regimes, and shows sensitivity towards ergodicity and thermalization [27][28][29][30].
Chebyshev expansions.-The basis of our numerical strategy is the expansion of the Dirac delta function in terms of Chebyshev polynomials T n , defined by the recurrence relation [15] T 0 ðxÞ ¼ 1; Any piecewise continuous function fðxÞ with x ∈ ½−1; þ1 admits such expansion [15,31,32], and can be approximated by a truncated sum The moments μ n ¼ R 1 −1 fðxÞT n ðxÞdx are the coefficients of the full expansion, while the γ n are introduced by the KPM to improve the quality of the truncated approximation, and depend on the order of the truncation M, but not on f [26].
Using the expansion for the delta function [15] for each term in Eq. (1), g M ðE; OÞ can be written in the form (9), with moments whereH ¼ H=ν þ ΔE is the rescaled and potentially shifted Hamiltonian, such that the spectrum ε k ¼ E k =ν þ ΔE is strictly contained in ½−1; 1 [26].
PHYSICAL REVIEW LETTERS 124, 100602 (2020) We can construct fixed bond dimension MPO approximations to the polynomials T ðDÞ n ðHÞ ≈ T n ðHÞ for any Hamiltonian H that is itself expressed as an MPO. Starting from T 0 ðHÞ ¼ 1 and T 1 ðHÞ ¼H (both exact MPO), we apply the recurrence relation between Chebyshev polynomials (8). This increases the bond dimension, so at each step we approximate the result with the maximum D allowed using standard TNS techniques [1], T The case of the LDOS allows for a more efficient implementation, since in that case the traces to be evaluated reduce to the single expectation value hΨjT ðDÞ n ðHÞjΨi. Thus, instead of each full polynomial, it is enough to find an MPS approximation of the vectors resulting from applying them, jt n i ≡ T n ðHÞjΨi, which satisfy the same recurrence relation. This reduction of Chebyshev expansions to states was used in Refs. [33,34] to estimate spectral function.
The sources of error in our approach are mainly the truncation of the Chebyshev expansion to a finite M, which limits the energy resolution, and the truncation error from the MPO approximation with finite D (see Ref. [26] for a detailed analysis). Note that because of the largely varying DOS across energy regions (in particular for local models as considered here, the DOS is Gaussian in the thermodynamic limit [35][36][37]), the precision of g M ðE; OÞ estimated with the above procedure worsens near the edges of the spectrum as discussed below. We can alleviate this problem by applying separate expansions to the Hamiltonian projected onto the different energy intervals, H → θðH − E cut ÞH (for high) or θðE cut − HÞH (for low energies) [38]. Since the step function θðxÞ can also be approximated using the KPM, this construction can be realized within our numerical method [26].
Models.-We have applied the method to two quantum spin chains with open boundary conditions (the scheme can be also used for periodic chains). The first is the Ising model, in general nonintegrable, except in the limits g ¼ 0 (classical) or h ¼ 0 (transverse field Ising model). This Hamiltonian has been profusely studied in the context of quantum quenches. Nontrivial dynamics has been observed and investigated in the nonintegrable regime [39][40][41][42][43], in particular, for the parameters that we consider, ðJ; g; hÞ ¼ ð1; −1.05; 0.5Þ. For comparison, we analyze also the integrable point (1,0.8,0). Second, we consider the PXP model, where P i ¼ ð1 − σ ½i z Þ=2. This kinetically constrained model was recently realized in a Rydberg atom chain experiment [44], and the observation of persistent revivals for particular initial configurations has triggered intense theoretical investigation about quantum scars as a possible mechanism to prevent thermalization [45][46][47][48][49].
Thermalization probes.-To probe thermalization in the Ising model, we consider three initial states that we call jXþi, jYþi, jZþi, defined as translationally invariant products of totally polarized spins in the corresponding directions. Their LDOS for the integrable and nonintegrable Ising models, and the DOS of both Hamiltonians are shown in Fig. 1 for a chain of N ¼ 80 particles. The results for the DOS are very precise already for moderate bond dimensions D and truncation parameter M. We have also found that the accuracy converges faster with D for our scheme than for previous results, see Sec. III of the Supplemental Material [26].
In the nonintegrable case we analyze the thermalization of O ¼ σ hold. Notice that if we had not used different (θ projected) Chebyshev expansions for different energy sectors, the results do not converge in the outer parts of the spectrum (pink lines in the figures).
We also probe thermalization for the initial states mentioned above by checking relation (6). For each of the jXþi, jYþi, and jZþi states, the figures show the result of evaluating the right-hand side of Eq. (5) for the observable O analyzed in the corresponding model versus E ¼ hΨjHjΨi. If the state thermalizes, and the approximation (5) is good enough, we expect that the result agrees with the thermal expectation value (black curve) at the same mean energy.
In the nonintegrable case (left panel of Fig. 2), the error bars are compatible with thermalization for the three states. This is particularly interesting for the jXþi and jZþi states, for which numerical simulations are not able to reach thermalization times [39,40], but there are arguments for eventual thermalization [42]. The significantly larger error bar for the jZþi state is related to the closeness of this state to the edge of the spectrum, which makes it sensitive to the discrete character of the latter, as evidenced also in the corresponding LDOS (lower left panel in Fig. 1). In the integrable case, the value of the most energetic of the states, jZþi, is not compatible with the assumption of thermal equilibrium, even with error bars. In this case, if the system equilibrates, we expect it to be to a generalized Gibbs ensemble [50]. Nevertheless, our observation cannot be taken as a test of such effect, because the estimate (5), in a case with degeneracies in the spectrum, does not necessarily correspond to the expectation value in the long time limit [26].
A similar analysis for the PXP model is shown in Fig. 3 (left panel) for a system of 40 spins. We compare the microcanonical estimate (3) for the operator O ¼ σ ½N=2 z (yellow line and symbols) to the thermal value (black line), and observe that the agreement is best close to the center of the spectrum, but values increasingly differ (error bars considered) towards the edges. This observation is compatible with exact diagonalization results [46] (for much smaller systems) that predict the existence of ETH-violating eigenstates (scar states) in all regions of the spectrum. Closer to the edges of the spectrum, the ratio of scar states with respect to ETH ones becomes non-negligible, which explains the more evident breaking of ETH in these regions.
We next consider two initial states jZ 2 i ≡ j↑↓↑↓…i and jZ 3 i ≡ j↑↓↓↑↓↓…i, for which unexpectedly long-lived oscillations have been experimentally observed in Rydberg atoms [44]. It has been recently postulated [45][46][47][48] that the slow dynamics of these states is due to their large overlap with a few scar states. These states lie nevertheless in the middle of the spectrum, E ¼ 0, where there is an exponentially large degeneracy, so that Eq. (5) (shown in the figure and compatible with the thermal value within error bars) is not necessarily a good estimate of the long-time limit.
The survival probability of these states, in contrast, does show the peculiarities of the real time dynamics of these two states. As shown in the right panel of Fig. 3, the fidelities of both states show periodic revivals. We observe that for both jZ 2 i and jZ 3 i, the height of the peaks seems to decrease exponentially with the system size, which we choose to be multiples of 6 [46]. The locations of the peaks are more robust (see inset), at times t Z 2 ≈ 3πn=2 and t Z 3 ≈ 9πn=8, for n ∈ Z, in agreement with the prediction in Ref. [46].
Discussion.-We have presented a technique, based on Chebyshev expansions and MPS algorithms, to compute generalized densities of states, and have shown how it can be used to directly probe thermalization in one-dimensional quantum many-body models.
Crucially, our method does not encounter the volume law of entanglement typically exhibited by high energy eigenstates, because it does not target such states individually,  Fig. 2); for initial states jZ 2 i, jZ 3 i (red and blue symbols), at E ¼ 0, (5) agrees with the thermal value. Right: Survival probability of both states as a function of time for different sizes. Revival times are almost independent of system size, and agree with predictions in Ref. [46] (dashed lines in the inset). but mixtures of them as operators, for which not the entanglement entropy but the mutual information affects the approximability. In particular, thermal equilibrium states can be efficiently described by MPOs [51,52]. We aim at averages of intensive quantities in microcanonical ensembles, which can also be described by such mixtures (see Ref. [26] for a quantitative analysis of the bond dimension needed for a given level of accuracy as a function of the system size).
We consider a broad range of potential extensions and applications to be analyzed in the future. First of all, our calculations for spin models can be easily extended to disordered, quasiperiodic, long-range interacting, bosonic, fermionic, and even two-dimensional systems. Beyond Hamiltonians, our scheme carries over to any sort of MPO, and could be useful to explore the spectral properties of Lindbladians, random MPO, or others. Regarding the survival probability, we are looking at extensions of our scheme that would allow us to monitor the evolution of (local) observables, even at finite temperature, or to detect so-called correlation holes at intermediate timescales [53].