Eigenstate thermalization in the two-dimensional transverse field Ising model: II. Off-diagonal matrix elements of observables

TL;DR

This paper investigates the off-diagonal matrix elements of observables in the eigenstates of the 2D transverse field Ising model, linking quantum chaos to matrix element structure and confirming the eigenstate thermalization hypothesis.

Contribution

It provides a detailed analysis of off-diagonal matrix elements in the 2D transverse field Ising model, demonstrating the applicability of random matrix theory and ETH in this context.

Findings

Ratio of variances approaches 2 in the quantum chaotic regime

Off-diagonal matrix elements depend on eigenstate energy differences

Results support the eigenstate thermalization hypothesis

Abstract

We study the matrix elements of few-body observables, focusing on the off-diagonal ones, in the eigenstates of the two-dimensional transverse field Ising model. By resolving all symmetries, we relate the onset of quantum chaos to the structure of the matrix elements. In particular, we show that a general result of the theory of random matrices, namely, the value 2 of the ratio of variances (diagonal to off-diagonal) of the matrix elements of Hermitian operators, occurs in the quantum chaotic regime. Furthermore, we explore the behavior of the off-diagonal matrix elements of observables as a function of the eigenstate energy differences, and show that it is in accordance with the eigenstate thermalization hypothesis ansatz.

Equations

$$O_{\alpha \beta }=\mathcal{O}(\bar{E}){\delta }_{\alpha \beta }+e^{-S(\bar{E})/2}{f}_O(\bar{E},\omega )R_{\alpha \beta },$$

$$O(t)=\sum _{\alpha }|c_{\alpha }{|}^2{O}_{\alpha \alpha }+\sum _{\begin{array}{c}\alpha ,\beta \\ (\alpha \ne \beta )\end{array}}{c}_{\alpha }^{\ast }{c}_{\beta }{e}^{\mathrm{i}(E_{\alpha }-E_{\beta })t}{O}_{\alpha \beta }.$$

$$\hat{H}=-J\sum _{⟨\mathbf{i},\mathbf{j}⟩}{\hat{\sigma }}_{\mathbf{i}}^z{\hat{\sigma }}_{\mathbf{j}}^z+g\sum _{\mathbf{i}}{\hat{\sigma }}_{\mathbf{i}}^x,$$

$$A_{\alpha \beta }\equiv ⟨\alpha |\hat{A}|\beta ⟩=\sum _{i,j}⟨\alpha |i⟩⟨i|\hat{A}|j⟩⟨j|\beta ⟩=\sum _i{A}_i{c}_i^{\alpha }{c}_i^{\beta },$$

$$\text{var}{A}_{\alpha \alpha }$$

$$\text{var}{A}_{\alpha \beta }^{(\alpha \ne \beta )}$$

$$\frac{\text{var}{A}_{\alpha \alpha }}{\text{var}{A}_{\alpha \beta }^{(\alpha \ne \beta )}}=2.$$

$${\hat{S}}^{\text{F}}\equiv \frac{1}{N}\sum _{\mathbf{i},\mathbf{j}}{\hat{\sigma }}_{\mathbf{i}}^z{\hat{\sigma }}_{\mathbf{j}}^z.$$

$${\hat{S}}^{\text{NN}}\equiv \frac{1}{N}\sum _{⟨\mathbf{i},\mathbf{j}⟩}{\hat{\sigma }}_{\mathbf{i}}^z{\hat{\sigma }}_{\mathbf{j}}^z.$$

$$\overline{|O_{\alpha \beta }|}=\frac{1}{\mathcal{N}}\sum _{\begin{array}{c}\alpha ,\beta \\ (\alpha \ne \beta )\end{array}}|O_{\alpha \beta }|,$$

$${\mathcal{D}}_{{\mathbf{k}}_0}=699,600\{\begin{array}{l}350,064\{\begin{array}{l}181,012\{\begin{array}{l}93,202\\ 87,810\end{array}\\ \\ 169,052\{\begin{array}{l}84,662\\ 84,390\end{array}\text{par}\end{array}\\ 349,536\{\begin{array}{l}180,232\{\begin{array}{l}91,652\\ 88,580\end{array}\\ \\ 169,304\{\begin{array}{l}85,164\\ 84,140\end{array}\text{par}\end{array}\end{array}$$

$${\mathcal{D}}_{{\mathbf{k}}_0}=1,342,208\{\begin{array}{l}671,104\{\begin{array}{l}337,192\{\begin{array}{l}170,440\{\begin{array}{l}86,056\\ 84,384\end{array}\\ 166,752\end{array}\\ \\ 333,912\{\begin{array}{l}166,752\\ 167,160\{\begin{array}{l}84,384\\ 82,776\end{array}\end{array}\text{par}\end{array}\\ 671,104\{\begin{array}{l}337,192\{\begin{array}{l}170,440\{\begin{array}{l}86,056\\ 84,384\end{array}\\ 166,752\end{array}\\ \\ 333,912\{\begin{array}{l}166,752\\ 167,160\{\begin{array}{l}84,384\\ 82,776\end{array}\end{array}\text{par}\end{array}\end{array}$$

Full text

I Introduction

Whether quantum statistical behavior can emerge in isolated systems, which evolve unitarily as dictated by quantum mechanics, is a question that has fascinated physicists since the early days of quantum mechanics von Neumann (1929). It has taken almost a century for experimental setups with the degree of isolation and control required to address such a question to become widely available Bloch et al. (2008); Cazalilla et al. (2011). Those experimental setups, involving mostly ultracold gases trapped in ultrahigh vacuum, have begun to be used to study quantum thermalization (or the lack thereof) in a variety of settings Kinoshita et al. (2006); Gring et al. (2012); Trotzky et al. (2012); Langen et al. (2015); Schreiber et al. (2015); Kaufman et al. (2016); Clos et al. (2016); Choi et al. (2016); Smith et al. (2016); Neill et al. (2016).

In parallel with the experimental activity, computational studies have shown that equilibration of observables can occur in quantum systems even if their dynamics are unitary Rigol et al. (2007, 2006); Kollath et al. (2007); Manmana et al. (2007); Rigol et al. (2008); Cramer et al. (2008); Flesch et al. (2008); Rigol (2009a, b); Eckstein et al. (2009); Rigol and Santos (2010); Bañuls et al. (2011); Calabrese et al. (2011); Gramsch and Rigol (2012); Khatami et al. (2013); Zangara et al. (2013); Wright et al. (2014); Sorg et al. (2014); Fagotti et al. (2014); Bertini et al. (2015); Brandino et al. (2015); Balz and Reimann (2017), but only those that are nonintegrable (quantum chaotic) are generally described by traditional statistical mechanics after equilibration Rigol et al. (2007); Rigol (2016). This can be understood in the context of the eigenstate thermalization hypothesis (ETH) Deutsch (1991); Srednicki (1994); Rigol et al. (2007) (see Ref. D’Alessio et al. (2016) for a recent review), which states that the diagonal matrix elements of observables in the eigenstates of the Hamiltonian are smooth functions of the energy 111The differences between the diagonal matrix elements of observables in neighbor eigenstates are exponentially small in the system size, while the off-diagonal ones are exponentially small in the system size. For an observable $\hat{O}$, the ETH ansatz can be written as Srednicki (1999); D’Alessio et al. (2016)

[TABLE]

where $\bar{E}\equiv(E_{\alpha}+E_{\beta})/2$, $S(E)$ is the thermodynamic entropy at energy $E$, and $\omega\equiv E_{\alpha}-E_{\beta}$. $\cal O$ and $f_{O}$ are smooth functions of their arguments, while $R_{\alpha\beta}$ are random numbers with zero mean and unit variance. The connection with statistical mechanics is immediate through ${\cal O}(E)$, which is the statistical mechanics prediction for $\hat{O}$ at the mean energy $E$. Eigenstate thermalization has been observed in a variety of nonintegrable lattice models Rigol et al. (2008); Rigol (2009a, b); Rigol and Santos (2010); Biroli et al. (2010); Roux (2010); Neuenhahn and Marquardt (2012); Khatami et al. (2013); Steinigeweg et al. (2013); Sorg et al. (2014); Beugeling et al. (2014); Kim et al. (2014); Steinigeweg et al. (2014); Khodja et al. (2015); Beugeling et al. (2015); Fratus and Srednicki (2015); Mondaini et al. (2016); Luitz (2016); Luitz and Bar Lev (2016).

To understand how the ETH ansatz (1) explains thermalization in isolated quantum systems, let us consider an initial state $|\Psi_{I}\rangle=\sum_{\alpha}c_{\alpha}|\alpha\rangle$, where $|\alpha\rangle$ are the eigenstates of the Hamiltonian $\hat{H}$ generating the dynamics ($\hat{H}|\alpha\rangle=E_{\alpha}|\alpha\rangle$). Using that $|\Psi(t)\rangle\equiv\exp[-{\rm i}\hat{H}t]|\Psi_{I}\rangle=\sum_{\alpha}c_{\alpha}\exp[-{\rm i}E_{\alpha}t]|\alpha\rangle$, we set $\hbar=1$, the time evolution of $\hat{O}$, $O(t)\equiv\langle\Psi(t)|\hat{O}|\Psi(t)\rangle$, follows from

[TABLE]

The fact that off-diagonal matrix elements of observables are exponentially small in the system size ensures that, after dephasing, the second term in Eq. (2) is exponentially smaller than the first term. The smooth dependence of $O_{\alpha\alpha}$ on the energy, combined with the fact that the width of the energy distribution of physical initial states (such as those in quantum quenches involving local Hamiltonians Rigol et al. (2008)) is generally subextensive as in statistical mechanics ensembles, then ensures that the expectation values of observables after relaxation are the ones predicted by statistical mechanics D’Alessio et al. (2016).

While the expectation values of observables after equilibration are determined by diagonal matrix elements, the first term in Eq. (2), the dynamics that results in equilibration and the fluctuations of observables about their equilibrated values are determined by the exponentially small off-diagonal matrix elements (and the initial state), the second term in Eq. (2). Also, the fact that fluctuation-dissipation relations hold in isolated quantum-chaotic systems without the need of assuming thermal equilibrium is encoded in the structure of the off-diagonal matrix elements of observables Khatami et al. (2013); D’Alessio et al. (2016). Early studies of eigenstate thermalization explored the qualitative behavior of off-diagonal matrix elements of observables in near-integrable and quantum-chaotic regimes Rigol et al. (2008); Rigol (2009b); Santos and Rigol (2010a), while more recent studies have looked into their distribution and scaling of their magnitude with system size, as well as the functional form of $f_{O}(\bar{E},\omega)$, in one-dimensional chains Khatami et al. (2013); Steinigeweg et al. (2013); Beugeling et al. (2015); D’Alessio et al. (2016); Luitz and Bar Lev (2016). Still, the off-diagonal matrix elements of observables remain much less studied than the diagonal ones.

Here, we study the behavior of the off-diagonal matrix elements of few-body observables, such as the structure factor and nearest neighbor spin correlations, in the eigenstates of the two-dimensional transverse field Ising model (2D-TFIM) on the square lattice. This work extends the analysis of the same model in Ref. Mondaini et al. (2016), in which we studied quantum chaos indicators and the diagonal matrix elements of observables. The 2D-TFIM exhibits quantum-chaotic behavior, and the diagonal matrix elements of few-body observables comply with the ETH ansatz (1), for nonvanishing but finite values of the field. We note that, despite its simplicity, the 2D-TFIM hosts both a zero-temperature quantum phase transition with increasing the strength of the transverse field, and a finite-temperature one, separating an ordered and a paramagnetic phase. Remarkably, the one-dimensional TFIM was recently realized experimentally using ultracold bosonic atoms in tilted optical lattices Simon et al. (2011).

One of our main goals in this work is to show that, in the quantum-chaotic regime, the matrix elements of observables in the 2D-TFIM satisfy a striking prediction from the Gaussian orthogonal ensemble (GOE) of random matrix theory. Namely, that the ratio between the variance of diagonal and off-diagonal matrix elements of observables in very small energy windows is universal and equal to 2 D’Alessio et al. (2016). We also study the scaling of the magnitude of off-diagonal matrix elements of observables, and explore the existence and behavior of the function $f_{O}(\bar{E},\omega)$ introduced in the ETH ansatz (1).

The presentation is organized as follows. In Sec. II, we introduce the model and the numerical methods used. Section III explores the connection between the onset of quantum chaos and the GOE result for the ratio between variances of the diagonal and off-diagonal matrix elements of observables. The scaling and behavior of the off-diagonal matrix elements, and their relation to the ETH ansatz, are studied in Sec. IV. A summary of our results is presented in Sec. V.

II Model and Numerical Method

The Hamiltonian of the 2D-TFIM, assuming periodic boundary conditions, reads

[TABLE]

where $\hat{\sigma}_{\bf i}^{z}$ ($\hat{\sigma}_{\bf i}^{x}$) is the $z$ ($x$) Pauli matrix at site ${\bf i}$ of the lattice. The strength of the nearest neighbor (denoted by $\langle{\bf i},{\bf j}\rangle$ in the constrained summation) Ising exchange interactions is given by $J$. In this work, we focus on the ferromagnetic case $(J>0)$, setting $J=1$ as our energy scale. (In Ref. Mondaini et al. (2016), we considered both the ferromagnetic and antiferromagnetic cases.) Lastly, $g$ denotes the magnitude of the transverse field. The majority of the results presented refer to regular square lattices with number of sites $N=\ell_{x}\times\ell_{y}$, where $\ell_{x}$ and $\ell_{y}$ denote the linear dimensions of the system in the $x$ and $y$ directions, respectively. We also solve for a tilted lattice with 20 sites, one of the two largest lattices (lattice 20A) studied in Ref. Mondaini et al. (2016). Here, the largest lattice considered is a regular lattice with 25 sites, while the smallest one, also regular, has 15 sites. We study various system sizes to carry out finite-size scaling analyses of observables of interest.

The 2D-TFIM [Eq. (3)] on the 2D square lattice possesses a variety of symmetries, which once identified and taken into account allow one to block diagonalize the Hamiltonian. Subsectors of the Hamiltonian with no symmetries are needed to establish the existence of quantum chaotic behavior using the energy spectrum Santos and Rigol (2010b, a); Mondaini et al. (2016). We make use of the following symmetries: translation ($\hat{T}$), spin-flipping ($\hat{Z}_{2}$), mirror in $x$ ($\hat{S}_{x}$), mirror in $y$ ($\hat{S}_{y}$) and, when applicable, mirror along the $x=y$ line ($\hat{S}_{xy}$). The latter symmetry is present when $\ell_{x}=\ell_{y}$ and the parity is the same under $\hat{S}_{x}$ and $\hat{S}_{y}$. After applying these symmetries, the subsectors of the Hamiltonian are diagonalized using full exact diagonalization.

In the best case scenario, when all those symmetries are present, the original Hilbert space (of dimension ${\cal D}=2^{N}$) can be split into smaller subsectors, some of which have dimension ${\cal D}^{\prime}\approx\frac{2^{N-4}}{N}$. As discussed in Appendix A, the largest subsector we fully diagonalize has $166,752$ states. It corresponds to the zero-momentum subsector of the $5\times 5$ lattice, after applying the $\hat{Z}_{2},\,\hat{S}_{x},$ and $\hat{S}_{y}$ operations, and for parities under $\hat{S}_{x}$ and $\hat{S}_{y}$ obeying $\lambda_{\hat{S}_{x}}\cdot\lambda_{\hat{S}_{y}}=-1$. For regular lattices we focus on the zero momentum subsector [$\boldsymbol{k}=(0,0)$], while for the 20-sites tilted lattice we focus on the $\boldsymbol{k}=\frac{\pi}{5}(2,1)$ momentum subsector, in which the only remaining symmetry to be resolved is $\hat{Z}_{2}$.

III Quantum chaos indicators

III.1 Ratio of adjacent gaps

We begin studying the so-called rigidity of the spectrum, i.e., by checking whether level repulsion takes place Bohigas et al. (1984); Brody et al. (1981); Haake (1991). This is a fundamental insight from random matrix theory, which has been shown to apply to a variety of nonintegrable lattice models D’Alessio et al. (2016).

Given that the 2D-TFIM is integrable in both the classical $(g\to 0)$ and paramagnetic $g\to\infty$ (one-site) limits, we test level repulsion when $g$ and $J$ have the same magnitude, $g=J=1$. Specifically, we compute the ratio of adjacent energy gaps Oganesyan and Huse (2007); Atas et al. (2013), $r_{n}\equiv\min\left(\delta_{n+1},\delta_{n}\right)/\max\left(\delta_{n+1},\delta_{n}\right)$, where $\delta_{n}=E_{n}-E_{n-1}$ is the difference between consecutive energy levels, and $\{E_{n}\}$ is the sorted list of eigenenergies (from the lowest to the highest) in each subsector of the Hamiltonian. Here is where the necessity of resolving all symmetries becomes apparent: subsectors of the Hamiltonian in which there are unresolved symmetries exhibit uncorrelated energy levels (and possibly extensive degeneracies), i.e., level repulsion is absent even if the Hamiltonian is quantum chaotic Santos and Rigol (2010b, a); Mondaini et al. (2016).

Figure 1(a) shows the distribution of the ratio of adjacent gaps obtained in the central one half of the spectrum 222Since the 2D-TFIM [Eq. (3)] only has short-ranged interactions and no randomness, the random matrix theory predictions for the level spacing statistics only apply away from the edges of the spectrum Brody et al. (1981); Kaplan and Papenbrock (2000); Flores et al. (2001); Santos and Rigol (2010b, a)., averaged between the equivalent symmetry irreducible subsectors with either $\lambda_{\hat{S}_{x}}\cdot\lambda_{\hat{S}_{y}}=-1$ or $\lambda_{\hat{S}_{x}}\cdot\lambda_{\hat{S}_{y}}=1$ in the $5\times 5$ lattice. The former contain ${\cal D}^{\prime}\approx 167000$ states, while the extra $\hat{S}_{xy}$ mirror symmetry in the latter results in subsectors with ${\cal D}^{\prime}\approx 85000$ states (see Appendix A). For Hamiltonians that are time-reversal symmetric, the appropriate random matrix ensemble is the GOE, whose distribution of the ratio of adjacent gaps is: $P_{\text{GOE}}(r)=(27/4)[(r+r^{2})\Theta(1-r)]/(1+r+r^{2})^{5/2}$ Atas et al. (2013). This prediction is depicted by the (green) dashed line in Fig. 1, and it is in almost perfect agreement with the distribution of the adjacent energy gaps obtained numerically. From this, one can conclude that quantum chaotic behavior takes place in the 2D-TFIM. The small differences seen between the analytic predictions and the numerical results are due to the fact that the former are exact only for $3\times 3$ matrices Atas et al. (2013); Mondaini et al. (2016).

In Fig. 2(a), we show numerical results for the average of adjacent energy gaps in the central one half of the spectrum, $\langle r\rangle_{\text{mid.}}$, as a function of the transverse field $g$ for four lattice sizes. The prediction from the GOE is $\langle r\rangle_{\text{GOE}}\approx 0.5359$ Atas et al. (2013), and is depicted in Fig. 2(a) as a horizontal dashed line. Even for the smallest lattice sizes studied, the average ratio of adjacent gaps is quite close to $\langle r\rangle_{\text{GOE}}$ for $g\simeq J$. With increasing lattice size one can see that the agreement improves for $g\simeq J$ and extends toward $g\ll 1$ and $g\gg 1$, which suggests that any nonzero but finite value of $g$ results in quantum chaotic behavior in the thermodynamic limit Santos and Rigol (2010b); Rigol and Santos (2010).

III.2 Ratio of diagonal and off-diagonal variances of matrix elements of observables

Given the quantum chaotic behavior observed in the energy spectrum, one might wonder whether other properties of random matrices are present in the 2D-TFIM. Of particular relevance to eigenstate thermalization, one can show that the variance of the diagonal and off-diagonal matrix elements of Hermitian operators in the eigenstates of random matrices are proportional to each other. For the GOE, the proportionality constant is exactly 2 D’Alessio et al. (2016).

This follows from the fact that the eigenstates of real symmetric random matrices are essentially orthonormal random vectors in arbitrary bases. Let us imagine we have a Hermitian operator $\hat{A}$, with $\hat{A}|i\rangle=A_{i}|i\rangle$. The matrix elements of $\hat{A}$ in the eigenkets of a real symmetric random matrix $\{|\alpha\rangle\}$ read

[TABLE]

where we defined $c_{i}^{\alpha}\equiv\langle i|\alpha\rangle=\langle\alpha|i\rangle$. The $c_{i}^{\alpha}$’s are Gaussian distributed with zero mean and variance equal to $1/{\cal D}$ (to leading order), where ${\cal D}$ is the dimension of the random matrix. Two results follow immediately from this: (1) $\overline{c_{i}^{\alpha}c_{j}^{\beta}}=(1/{\cal D})\delta_{\alpha\beta}\delta_{ij}$, and (2) different moments of the distribution of $c_{i}^{\alpha}$’s are related, e.g., $\overline{(c_{i}^{\alpha})^{4}}=3\overline{(c_{i}^{\alpha})^{2}}$.

We are interested in the fluctuations of the diagonal and off-diagonal matrix elements of operator $\hat{A}$ in the eigenstates of a random matrix: $\text{var\,}A_{\alpha\alpha}=\overline{A_{\alpha\alpha}^{2}}-\overline{A_{\alpha\alpha}}^{2}$ and $\text{var\,}A_{\alpha\beta}^{(\alpha\neq\beta)}=\overline{A_{\alpha\beta}^{2}}-\overline{A_{\alpha\beta}}^{2}$, respectively. Using the two results mentioned above, one gets

[TABLE]

and

[TABLE]

which means that the ratio between the variance of the diagonal and off-diagonal matrix elements of $\hat{A}$ is

[TABLE]

This is consistent with the ETH ansatz, as in the latter the eigenstate to eigenstate fluctuations of the diagonal matrix elements of observables are exponentially small in the system size, as the off-diagonal matrix elements are. The question that remains is whether the ratio of the variances of diagonal and off-diagonal matrix elements of observables in physical systems with short-range interactions and no randomness is constant away from the edges of the spectrum and equal to 2. Since in such systems the diagonal and off-diagonal matrix elements of observables are expected to have structure in $\bar{E}$ and $\omega$ [see the ETH ansatz (1)] the calculation of the variances has to be carried out within sufficiently small energy windows so that ${\cal O}(\bar{E})$ and $e^{-S(\bar{E})/2}f_{O}(\bar{E},\omega)$ are essentially constant (as the corresponding terms are in random matrices D’Alessio et al. (2016)).

To address this question, we first study the matrix elements of the ferromagnetic structure factor,

[TABLE]

This non-local few-body observable is an order parameter for the phase transitions that occur in the 2D-TFIM. Its expectation value is extensive (order 1) in the ordered (paramagnetic) phase. Since the focus of our study are energy eigenstates in the central one half of the spectrum (“high temperature” eigenstates; see Appendix B), the eigenstate expectation values of $\hat{S}^{\text{F}}$ in our calculations are $O(1)$ Mondaini et al. (2016).

In Fig. 1(b) and 1(c), we show the ratio of the variances of the diagonal and off-diagonal matrix elements of the structure factor in the energy eigenstates of the 2D-TFIM, as a function of the eigenstate index. The results were obtained on the $5\times 5$ lattice, within two symmetry irreducible subsectors with $\lambda_{\hat{Z}_{2}}=\lambda_{\hat{S}_{x}}=-\lambda_{\hat{S}_{y}}=1$ [Fig. 1(b)] and $\lambda_{\hat{Z}_{2}}=\lambda_{\hat{S}_{x}}=\lambda_{\hat{S}_{y}}=\lambda_{\hat{S}_{xy}}=1$ [Fig. 1(c)], for $g=1$. It is remarkable that, away from the edges of the spectrum, the ratios of variances fluctuate about 2. Actually, the average of the ratios of variances within the central one half of the spectrum, shown as horizontal dashed-dotted lines, is very close to 2 (closer for the largest symmetry irreducible subsector). The insets in Figs. 1(b) and 1(c) show the distribution of the ratios of variances when considering, once again, only the central one half of the spectrum (vertical dashed-dotted lines depict the average).

In Figs. 2(b), we plot the average ratio of variances of the diagonal and off-diagonal matrix elements of the structure factor as a function of the transverse field strength, for the same lattice sizes as in Fig. 2(a). The correlation between the values of $g$ for which the average ratio of variances is closest to 2 [Fig. 2(b)] and for which the average ratio of adjacent energy gaps is closest to the GOE prediction [Fig. 2(a)] is apparent. With increasing system size, one can see that the range of values of $g$ over which the numerical results are closest to random matrix theory prediction increases. We note that, when departing from the GOE predictions, the average ratio of variances increases. This is the result of an increase in the eigenstate-to-eigenstate fluctuations of the diagonal matrix elements Rigol (2009a, b), and an increase of the number of off-diagonal matrix elements that become very small Rigol (2009b); Khatami et al. (2013). The latter effect can make the ratio of variances become arbitrarily large.

Our results contrast the ones obtained by Steinigeweg et al. Steinigeweg et al. (2013) for the ratios of variances of the diagonal and off-diagonal matrix elements of current operators in nonintegrable fermionic chains. They were found to exhibit a significant dependence on the energy of the eigenstates, to differ from 2 in the center of the energy spectrum, and to differ between the spin and energy currents. In our calculations, we have found that lack of agreement with the random matrix theory prediction can be a result of finite-size effects and/or the width of the windows used to compute the variances. Finite-size effects affect different observables in different ways, and the size of the windows used to compute the variances needs to be selected with care to avoid the influence of the nontrivial structure of matrix elements of observables in physical Hamiltonians (not present in random matrix theory).

In Figs. 3(a) and 3(b), we show results of calculations identical to those reported in Figs. 1(b) and 1(c) but for the nearest-neighbor spin-spin correlation function,

[TABLE]

$\hat{S}^{\text{NN}}$ exhibits larger finite-size effects than the structure factor, as apparent from the fact that for $\hat{S}^{\text{NN}}$: (1) the fluctuations of the ratio of variances are slightly larger throughout the spectrum, and (2) the average exhibits a larger departure from 2. In Fig. 3(c), we show that if one increases the size of the windows used to compute the variances, in order to reduce the fluctuations of their ratio, then the nontrivial structure of the diagonal and off-diagonal matrix elements of observables kicks in and the average ratio of variances departs from the GOE prediction. The departure depends strongly on the lattice geometries used. See Appendix D for the equivalent of Fig. 3(c) for the ferromagnetic structure factor.

In Appendix C, we show results for the eigenstate expectation values of $\hat{S}^{\text{F}}$ and $\hat{S}^{\text{NN}}$ for different values of $g$ and system sizes. They further help gaining an understanding of finite-size effects in different observables, as well as how ETH sets in when $g$ is increased from zero and breaks down as $g$ becomes much larger than $J$.

IV Off-diagonal matrix elements and the ETH ansatz

Having unveiled random matrix theory behavior in the matrix elements of observables in microscopic energy windows, we now study the off-diagonal matrix elements of the structure factor and the nearest-neighbor spin-spin correlation function as a function of $\omega=E_{\alpha}-E_{\beta}$.

Figure 4 shows the absolute value of the off-diagonal matrix elements of the structure factor [Figs. 4(a) and 4(b)] and of the nearest-neighbor spin-spin correlation function [Figs. 4(c) and 4(d)] versus $\omega$, for $\omega>0$ ($S^{\rm{F}}_{\alpha\beta}$ and $S^{\rm NN}_{\alpha\beta}$ are symmetric) and $2|\bar{E}|/N=|E_{\alpha}+E_{\beta}|/N\leq 0.1$. We report results for the two largest lattices with square aspect ratio, namely, the $\sqrt{20}\times\sqrt{20}$ lattice [Figs. 4(a) and 4(c)] and the $5\times 5$ lattice [Figs. 4(b) and 4(d)]. The off-diagonal matrix elements of each observable are qualitatively similar in the two lattices.

In order to study the behavior of the smooth $e^{-S(\bar{E})/2}f_{O}(\bar{E},\omega)$ function [see Eq. (1)], we compute the running (or coarse-grained) average over small windows of width $\delta\omega$ for the two observables of interest. The widths of the windows are different for different lattices. They are selected such that the result of the averaging produces a smooth curve that is not sensitive to the exact value of $\delta\omega$ chosen. The results of such a coarse-graining procedure are reported in Fig. 4 as continuous black lines.

Since we have chosen $\bar{E}/N\approx 0$, the energy at infinite temperature, $e^{-S(\bar{E})/2}$ is nothing but $({\cal D}^{\prime})^{-1/2}$ where ${\cal D}^{\prime}$ is the size of the subsector studied. That the off-diagonal matrix elements in our calculations are indeed proportional to $({\cal D}^{\prime})^{-1/2}$ (and, hence, exponentially small in the system size) can be verified by computing the average

[TABLE]

where ${\cal N}$ is the number of terms contributing to the sum, and plotting it vs ${\cal D}^{\prime}$ for different lattice sizes and, within a given lattice size, for different subsectors. In Fig. 5, we present such a plot for the average between the off-diagonal matrix elements of eigenstates with $2|\bar{E}|/N\leq 0.1$, in all subsectors of the zero momentum sector of the regular lattices studied, and for $g=1$. The numerical results exhibit an excellent agreement with the expected $({\cal D}^{\prime})^{-1/2}$ behavior. This allows us to extract $f_{O}(\bar{E},\omega)$, up to a constant, from the running averages.

Figure 6 shows $f_{O}(\bar{E},\omega)$ determined this way, for the structure factor [Fig. 6(a)] and for the nearest-neighbor spin-spin correlation function [Fig. 6(b)], for one subsector of the 20-sites tilted lattice and two subsectors of the $5\times 5$ lattice. The results for each observable in the two lattices are very close to each other. For both observables, $f_{O}(\bar{E}\approx 0,\omega)$ is nearly constant for very small values of $\omega$ (this is the regime in which the random matrix predictions were tested in Sec. III.2) and decays rapidly for large $\omega$. In the latter regime, the decay of $f_{S^{\rm{F}}}(\bar{E}\approx 0,\omega)$ appears to be a composition of several exponentials, while the one of $f_{S^{\rm{NN}}}(\bar{E}\approx 0,\omega)$ is close to that of a single exponential. Exponential decays of $f_{O}(\bar{E},\omega)$ for local observables $\hat{O}$ at large $\omega$ have been observed in previous works Khatami et al. (2013); Beugeling et al. (2015); D’Alessio et al. (2016). They can be understood in terms of perturbation theory for systems with a bounded spectrum. The reason is that many-particle processes, which are suppressed exponentially, are required to connect eigenstates with large energy differences D’Alessio et al. (2016). The inset in Fig. 6(b) shows that the exponential behavior found in the $5\times 5$ lattice is robust to changes in the size of the window used in the calculations.

V Summary

Using large (within full exact-diagonalization calculations) Hilbert space sizes, we studied properties of the matrix elements of two few-body observables in the eigenstates of the 2D-TFIM. We showed that the onset of quantum chaos, identified using a level spacing analysis, results in the applicability of another random matrix theory prediction. The ratio of the variances of the diagonal and off-diagonal matrix elements of observables, calculated within very small energy windows, is constant across the spectrum (excluding the edges) and equal to 2. We also studied the behavior of the off-diagonal matrix elements as a function of $\omega=E_{\alpha}-E_{\beta}$. We showed that, in the quantum chaotic regime, their smooth part can be well described by a function $e^{-S(\bar{E})/2}f_{O}(\bar{E},\omega)$ as prescribed by the ETH ansatz. $f_{O}(\bar{E},\omega)$ was shown to be nearly constant for small values of $\omega$ and to exhibit a rapid (exponential) decay for large values of $\omega$. An interesting problem left for future studies is correlating the temporal evolution of observables at short, intermediate, and long times with the large, intermediate, and small $\omega$ behavior of $f_{O}(\bar{E},\omega)$, and unveiling the effect of the initial states in the dynamics at different time scales.

Acknowledgements.

RM is financially supported by the National Natural Science Foundation of China (NSFC) (Grant Nos. U1530401, 11674021 and 11650110441) and MR by the U.S. Office of Naval Research, Grant No. N00014-14-1-0540. The computations were performed in the Tianhe-2JK at the Beijing Computational Science Research Center (CSRC), the Institute for CyberScience at Penn State, and the Center for High-Performance Computing at the University of Southern California.

Appendix A Hilbert space subsectors

We rewrite the original Fock basis, written in terms of the eigenstates of $\hat{\sigma}^{z}_{\bf i}$, using the following:

Translations $\hat{T}$: $\hat{T}_{\boldsymbol{R}}|\psi_{\hat{T}_{\boldsymbol{R}}}\rangle=e^{-{\rm i}\boldsymbol{k}\cdot\boldsymbol{R}}|\psi_{\hat{T}_{\boldsymbol{R}}}\rangle$, where $\boldsymbol{k}$ defines the $N$ possible momentum sectors. For the regular lattices, we diagonalize only the zero momentum sector, i.e., we deal with real matrices. For the 20-sites tilted lattice, we diagonalize the $\boldsymbol{k}=\frac{\pi}{5}(2,1)$ momentum sector (see Ref. Mondaini et al. (2016)).

Spin-flip $\hat{Z}_{2}$: $\hat{Z}_{2}|\psi_{\hat{T}_{\boldsymbol{R}},\hat{Z}_{2}}\rangle=\lambda_{\hat{Z}_{2}}|\psi_{\hat{T}_{\boldsymbol{R}},\hat{Z}_{2}}\rangle$, where $\lambda_{\hat{Z}_{2}}=\pm 1$. The remaining symmetries apply only to the zero momentum sector of the regular lattices.

Mirror in $x$, $\hat{S}_{x}$: $\hat{S}_{x}|\psi_{\hat{T}_{\boldsymbol{R}},\hat{Z}_{2},\hat{S}_{x}}\rangle=\lambda_{\hat{S}_{x}}|\psi_{\hat{T}_{\boldsymbol{R}},\hat{Z}_{2},\hat{S}_{x}}\rangle$, where $\lambda_{\hat{S}_{x}}=\pm 1$.

Mirror in $y$, $\hat{S}_{y}$: $\hat{S}_{y}|\psi_{\hat{T}_{\boldsymbol{R}},\hat{Z}_{2},\hat{S}_{x},\hat{S}_{y}}\rangle=\lambda_{\hat{S}_{y}}|\psi_{\hat{T}_{\boldsymbol{R}},\hat{Z}_{2},\hat{S}_{x},\hat{S}_{y}}\rangle$, where $\lambda_{\hat{S}_{y}}=\pm 1$.

If $\ell_{x}=\ell_{y}$, then for $\lambda_{\hat{S}_{x}}=\lambda_{\hat{S}_{y}}$ one also has a mirror symmetry along the $x=y$ line, $\hat{S}_{xy}$: $\hat{S}_{xy}|\psi_{\hat{T}_{\boldsymbol{R}},\hat{Z}_{2},\hat{S}_{x},\hat{S}_{y},\hat{S}_{xy}}\rangle=\lambda_{\hat{S}_{xy}}|\psi_{\hat{T}_{\boldsymbol{R}},\hat{Z}_{2},\hat{S}_{x},\hat{S}_{y},\hat{S}_{xy}}\rangle$, where $\lambda_{\hat{S}_{xy}}=\pm 1$.

The splitting of sectors after each symmetry is given below for the two largest regular lattices that we study, namely, the $6\times 4$ and $5\times 5$ lattices. We start with the largest sector from translational symmetry, $\boldsymbol{k}=(0,0)$, and then apply $\hat{Z}_{2}$, $\hat{S}_{x}$, $\hat{S}_{y}$, and, if applicable $\hat{S}_{xy}$. The top (bottom) number in each curly brace refers to the size of the subsector with positive (negative) parity.

For the $6\times 4$ lattice the splitting of the sectors reads (ordered to show the splitting under $\hat{Z}_{2}$, $\hat{S}_{x}$, and $\hat{S}_{y}$):

[TABLE]

For the $5\times 5$ lattice the splitting of the sectors reads (ordered to show the splitting under $\hat{Z}_{2}$, $\hat{S}_{x}$, $\hat{S}_{y}$, and, if applicable $\hat{S}_{xy}$):

[TABLE]

Appendix B Many-body density of states

The many-body density of states of systems with few-body interactions is, in general, Gaussian Brody et al. (1981). This is the case in the 2D-TFIM. The many-body density of states for subsectors that are even with respect to all parity operations (when applicable) for the lattices studied in this work are presented in Fig. 7 for $g=1$, along with Gaussian fits. The inset shows that the ratio between the width of the Gaussian fits and $\ln{\cal D}^{\prime}$ ($\ln{\cal D}^{\prime}$ is proportional to $N$) vanishes as $(\ln{\cal D}^{\prime})^{-1/2}$. This means that, with increasing system size, the overwhelming majority of the eigenstates of the Hamiltonian have $E_{\alpha}/N$ increasingly close to 0, i.e., they are “infinite temperature” eigenstates. Those are the 2D-TFIM eigenstates studied in this work, as we have focused on eigenstates that are located in the central one half of the spectrum.

Appendix C Diagonal matrix elements

In Fig. 8, we plot the expectation values of the structure factor [Fig. 8(a)–8(c)] and of the nearest-neighbor spin-spin correlation function [Fig. 8(d)–8(f)] in the eigenstates of the 2D-TFIM in four regular lattices. We show results from all parity subsectors of the zero momentum sector, for three values of $g$. As one can see in Figs. 8(a), 8(b), 8(d), and  8(e), as $g$ departs from $g=0$ the expectation values of the observables become smooth functions of the energy (away from the edges of the spectrum), i.e., eigenstate thermalization occurs. If one further increases $g$, when $g\gg J$, the system approaches an integrable regime and eigenstate thermalization breaks down [Figs. 8(c) and 8(f)].

Appendix D Ratio of variances for different energy windows

In Fig. 3(c), we show how the ratio of variances of diagonal and off-diagonal matrix elements of the nearest-neighbor spin-spin correlation function depends on the size of the window used in the calculation. Figure 9 shows equivalent results for the ferromagnetic structure factor. The deviations of the results for the latter observable from the random matrix prediction exhibit a behavior that is qualitatively similar to that seen in Fig. 3(c). However, for identical window sizes, the deviations of the results for the ferromagnetic structure factor are much smaller than the ones for the nearest-neighbor spin-spin correlation function. This is expected because, as mentioned before, the latter observable exhibits stronger finite-size effects in our calculations.