Hamiltonian Dynamics of a Sum of Interacting Random Matrices

TL;DR

This paper analytically investigates the dynamics of observables in a quantum system modeled by two interacting random matrices, revealing how non-relaxing components are determined by the Hamiltonian's spectrum and affecting late-time correlators.

Contribution

It provides an exact analytical framework for understanding non-relaxing components in quantum observables using random matrix models and diagrammatic resummation techniques.

Findings

Exact late-time correlator values for Hamiltonian H=A+B

Power-law relaxation governed by H's spectrum, independent of A

Modifications in OTOCs due to non-relaxing parts of A

Abstract

In ergodic quantum systems, physical observables have a non-relaxing component if they "overlap" with a conserved quantity. In interacting microscopic models, how to isolate the non-relaxing component is unclear. We compute exact dynamical correlators governed by a Hamiltonian composed of two large interacting random matrices, $H=A+B$. We analytically obtain the late-time value of $\langle A(t) A(0) \rangle$; this quantifies the non-relaxing part of the observable $A$. The relaxation to this value is governed by a power-law determined by the spectrum of the Hamiltonian $H$, independent of the observable $A$. For Gaussian matrices, we further compute out-of-time-ordered-correlators (OTOCs) and find that the existence of a non-relaxing part of $A$ leads to modifications of the late time values and exponents. Our results follow from exact resummation of a diagrammatic expansion and hyperoperator techniques.

Figures

Tables

Table 1. Table 1: Summary of functions and result for a few examples we checked numerically.

$A$	$B$	$R_A\left(z\right)$	$f_H\left(z\right)$
GUE	GUE	$\lambda z$	$\frac{1}{2}\left(z-\sqrt{z^2-4}\right)$	$\lambda$
MP	MP	$\frac{z}{1-z}$	$\frac{z+1-\sqrt{z^2-2z-7}}{2\left(z+2\right)}$	$\frac{1}{2}$
BIN	BIN	$\frac{z}{z^2-{\lambda }^2}$	$\frac{z}{\sqrt{\left(z^2-{\left(\lambda +\mu \right)}^2\right)\left(z^2-{\left(\lambda -\mu \right)}^2\right)}}$	$\max (\frac{1}{2},1-\frac{{\mu }^2}{2{\lambda }^2})$

Equations

$$H=A+B$$

$${⟨A(t)A(0)⟩}_c\underset{t\to \infty }{\sim }{A}_{\infty }^2+C_{A,H}{\left|⟨e^{iHt}⟩\right|}^2$$

$$\frac{1}{2}⟨[A(t),A(0){]}^2⟩$$

$$⟨\circ ⟩:=\frac{1}{N}{\mathbb{E}}_{A,B}\left[\text{Tr}\left(\circ \right)\right]$$

$$A=U^{†}\Lambda U$$

$$G(t)=⟨A(t)A(0)⟩=⟨e^{iHt}A{e}^{-iHt}A⟩$$

$$e^{iHt}=\oint \frac{dz}{2\pi i}\frac{e^{izt}}{z-H}$$

$$G(t)$$

$$=\oint \frac{dz}{2\pi i}\frac{dw}{2\pi i}{e}^{i(z-w)t}G(z,w)$$

$$G(z,w)$$

$$=\sum _{m=0}^{\infty }\sum _{n=0}^{\infty }\frac{1}{z^{n+1}}\frac{1}{w^{m+1}}⟨H^{n}A{H}^{m}A⟩$$

$$⟨H^{n}A{H}^{m}A⟩=⟨{\left(A+B\right)}^{n}A{\left(A+B\right)}^{m}A⟩$$

$$⟨X_1\dots X_n⟩=:\sum _{\pi \in \text{NC}(n)}\prod _{b\in \pi }\left\{X_{i_1}\dots X_{i_{|b|}}\right\}$$

$$f_H(z):=⟨\frac{1}{z-H}⟩=\sum _{n=0}^{\infty }\frac{⟨H^n⟩}{z^{n+1}}.$$

$$\underset{\delta \to 0^{+}}{\lim }\frac{1}{x\pm i\delta }=P\frac{1}{x}\mp i\pi \delta (x)$$

$$f(x\pm i\delta )$$

$$=\varphi (x)\mp i\pi \rho (x)$$

$$\varphi (x):=⟨P\frac{1}{x-H}⟩=P\int d\lambda \frac{\rho (\lambda )}{x-\lambda }$$

$$f_H(z)$$

$${\Sigma }_H(z)$$

$${\Sigma }_A(z)\equiv R_A(f(z))\equiv \sum _{p=0}^{\infty }\left\{A^{n+1}\right\}{f}_H(z{)}^n$$

$${\Sigma }_A(x\pm i\delta )$$

$$\equiv {\Sigma }_A^R(x)\mp i{\Sigma }_A^I(x)$$

$${\pi }_A(z,w)$$

$$=\frac{{\Sigma }_A(z)-{\Sigma }_A(w)}{f(z)-f(w)}$$

$${\pi }_H(z,w)={\pi }_A(z,w)+{\pi }_B(z,w)$$

$$h(z,w)=\frac{f(z)f(w)}{1-{\pi }_H(z,w)f(z)f(w)}$$

$$h(z,w)=-\frac{f(z)-f(w)}{z-w}$$

$$T(z,w)$$

$$=\frac{f(w){\Sigma }_A(z)-f(z){\Sigma }_A(w)}{f(z)-f(w)}$$

$$\begin{array}{rl}G(z,w)& =h(z,w){\left(T(z,w)+{\Sigma }_A(z)+{\Sigma }_A(w)\right)}^2\\ & +\sum _{n=0}^{\infty }\left\{A^{n+2}\right\}f(z)f(w)\sum _{t=0}^{n}f(z{)}^{t}f(w{)}^{n-t}\end{array}$$

$$\boxed{\begin{array}{rl}G(w,z)& =f(z)f(w)\frac{{\Sigma }_A(w)-{\Sigma }_A(z)}{f(w)-f(z)}\\ & -\frac{f(z)-f(w)}{z-w}{\left(\frac{f(z){\Sigma }_A(z)-f(w){\Sigma }_A(w)}{f(z)-f(w)}\right)}^2\end{array}}$$

$$\{A\}=\{B\}=0,$$

Full text

1. I Introduction

2. II Exact Analysis

3. II.1 Frequency space representation of correlators

4. II.2 Diagrammatic Solution

5. II.2.1 Propagator

6. II.2.2 1PI Diagrams

7. II.2.3 2PI Diagrams

8. II.2.4 Box Resummation

9. II.2.5 Triangle Resummation

10. II.3 Two–point correlation function

11. II.4 Example: $A,B$ Gaussian

12. III Late Time Asymptotics

13. III.1 General Case

14. III.2 Example: Asymptotics for a $\sqrt{x}$ Edge

15. IV Finite temperature modifications

16. V The Gaussian Rotation Approach

17. VI Numerical Confirmation

18. VII Discussion

19. A Short time expansion

I Introduction

Consider an isolated quantum system whose Hamiltonian can be decomposed into two terms

[TABLE]

Under what conditions does the observable $A$ have a non-relaxing component? There are many more or less exotic mechanisms for $\langle A(t)A(0)\rangle_{c}$ to approach a non-zero constant at late time: $A$ could commute with $H$, or the whole system could be localizedAnderson (1958); Gornyi et al. (2005); Basko et al. (2006); Serbyn et al. (2013); Huse et al. (2014); Nandkishore and Huse (2015), scarredBernien et al. (2017); Turner et al. (2018); Lin and Motrunich (2019); Choi et al. (2019), shatteredKhemani and Nandkishore (2019), or otherwise mistreated by theorists. However, even in the simple case in which $H$ is ergodic and $A$ does not commute with it, the operator $A$ appears to be “part” of the conserved energy and we expect it to have a non-relaxing component.

This intuition is misleading. For any operator $O$ and any Hamiltonian $H$, one could write $H=(H-O)+O$ and formally reproduce the decomposition in Eq. (1) without learning anything. For Eq. (1) to be nontrivial, there must be some physical sense in which the decomposition is defined. In this paper, we consider the case in which $A$ and $B$ are large independent random matrices of dimension $N$ with arbitrary spectra. They need not obey Wigner’s celebrated semicircle law, but the eigenvectors of $A$ must be in generic position with respect to those of $B$. While this choice appears technical, it is a natural model if $H$, $A$ and $B$ each satisfy the eigenstate thermalization hypothesis (ETH) and $A$ and $B$ are built out of physically distinct local operators.

Our primary result is an exact integral representation of the dynamical correlator $\langle A(t)A(0)\rangle_{c}$ at large $N$ (Eq. (42)). From this representation we derive the asymptotic structure at late time,

[TABLE]

This form has two important features: first, the late-time constant $A_{\infty}^{2}$ quantifies the non-relaxing component of $A$. It has a definite integral representation in terms of the densities of states of $A$ and $H$ (see Eq. (59)). In general, $A_{\infty}^{2}$ is not zero –this accords with the intuition that $A$ is ‘part’ of the energy– however it is not given by a simple trace overlap.

Second, the $\langle A(t)A(0)\rangle_{c}$ approaches its late-time value with a form governed by the characteristic function of $H$ alone. The particular operator $A$ only enters through the constant rescaling $C_{A,H}$. This implies that the late time power laws come from the singularities in the density of states of $H$ and do not depend on any detailed structure of $A$. For the many random matrix ensembles where the density of states of $H$ has square-root singularities –the semicircle law is one example– Eq. (2) predicts $1/t^{3}$ decay. Still, other possibilities exist: if $A$ and $B$ are projectors, for example, the resulting Hamiltonian has different singularities which produce slower $1/t$ decay.

We set up the calculation at infinite temperature, but once we have those results the full temperature dependence can be derived introducing Boltzmann weights in appropriate places (Sec.IV). In the zero temperature limit the power law approach of Eq. (2) is modified and typically the decay exponents are reduced by a factor of $2$. This behavior is reminiscent of the SYK modelKitaev (2015); Kitaev and Suh (2018); Maldacena and Stanford (2016); Cotler et al. ; Bagrets et al. (2016) at times $t>N$, the number of fermions, but $t<e^{Ns}$, the level spacing.

In the special case where $A$ and $B$ are both Gaussian random matrices, and thus satisfy the semi-circle law, these results simplify significantly: the late time value $A_{\infty}^{2}$ is given by a trace overlap between $A$ and $H$ (Eq. (83)), as one might have guessed, and the asymptotic form Eq. (2) holds exactly for all times $t$. We derive these simplified forms both from our general diagrammatic formalism, and from a more direct hyperoperator approach.

There has been a lot of recent interest in computing out of time ordered correlators (OTOCs) as a rough characterization of quantum chaosKitaev (2015); Kitaev and Suh (2018); Nahum et al. (2018); Maldacena et al. (2016); Khemani et al. (2018); von Keyserlingk et al. (2018); Vijay and Vishwanath ; Chan et al. (2018). For Gaussian $A$ and $B$ we use the hyperoperator technique to compute the OTOC,

[TABLE]

and find that it decays asymptotically more slowly than if it did not contain a conserved piece.

It is worthwhile noting that much recent analytical progress on dynamics in many-body systems has been made in the setting of random unitary circuits.Chandran and Laumann (2015); Nahum et al. (2017, 2018); von Keyserlingk et al. (2018) These models do not naturally have conservation laws, and introducing the extra structure to create themKhemani et al. (2018); Rakovszky et al. (2018) makes calculations significantly harder. On the contrary, Hamiltonian models automatically come with a conserved energy, but deriving exact results there is difficult in the absence of further structure like conformal symmetry, integrability, or a large-$N$ limit. Our work falls in this last class: we treat a Hamiltonian system exactly, at the price of introducing large-$N$ random matrices.

All of our results are at infinite temperature and averaged over the random matrix ensemble; for the rest of the paper, the “expectation value” symbol means

[TABLE]

The normalization is such that the $N\times N$ identity matrix has expectation equal to $1$.

The ensemble for $A$ and $B$ is fairly general: we just need their densities of states $\rho_{A}$ and $\rho_{B}$ to be well defined in the limit $N\to\infty$, and their eigenspaces to be in generic position. One way to realize such matrices is

[TABLE]

where $\Lambda$ is a diagonal matrix and $U$ is a Haar unitary.

The paper is organized as follows: in Sec.II we introduce a set of diagrammatic tools and use them to derive the exact frequency space two–point function. Using this result we study (Sec.III) the long time behavior of $\langle A(t)A(0)\rangle$ and identify the structure mentioned in Eq.(2). We then discuss (Sec.IV) how these results are modified at finite temperature. In Sec.V we introduce the direct hyperoperator approach for Gaussian matrices. We use it to reproduce the results of the previous sections and compute the OTOC. Finally, in Sec.VI we compare the predictions of our analysis with numerics, obtaining a satisfactory agreement.

II Exact Analysis

II.1 Frequency space representation of correlators

The correlation function we focus on in this section is

[TABLE]

It is convenient to reformulate this problem in frequency space, using the Cauchy representation of the time evolution operator:

[TABLE]

where the integral is over any contour that encloses the full spectrum of $H$ (see fig. 1). At any finite Hilbert space size $N$ this contour is closed, and in the thermodynamic limit it can be closed at infinity.

Using the integral representation just introduced, we can rewrite eq. 6 as

[TABLE]

which defines the frequency space correlator:

[TABLE]

The last formal manipulation involving the geometric series shows that we can reduce the problem of determining $G(z,w)$ to the calculation of mixed moments of $A$ and $H$:

[TABLE]

For practical purposes it is better to write the calculation in terms of the matrices $A$ and $B$, as they are the independent objects in the problem.

If $A$ and $B$ are Gaussian random matrices, the evaluation of the moments is readily accomplished using Wick’s theorem. For the more general (non-semicircle) matrices we would like to treat here, the free cumulant expansion plays the role of a generalized Wick’s theorem for non-commuting matrices. Given a set of noncommuting random variables $X_{i}$, we can recursively define their free cumulants (denoted with curly braces) through the formula

[TABLE]

where NC(k) denotes the set non–crossing partitions of $k$ objects, and $b$ is a block in the partition $\pi$. See Ref. Morampudi and Laumann, for more pedagogical details of this formalism. For aficionados of planar perturbation theory, we note that the free cumulants of a random matrix $X$ are closely related to the fully renormalized vertices of $X$.

Restating the problem in terms of free cumulants makes things simpler: in the large $N$ limit all the mixed free cumulants of independent random matrices vanishNica and Speicher , which simplifies enormously the calculation of moments like the one of eq.12. When two random matrices have this property, we say that the two variables are freely independent.

II.2 Diagrammatic Solution

The calculation of moments of large independent random matrices lends itself to a convenient diagrammatic representation through the free cumulant expansion (eq. 13). The diagrammatics are similar to those of Refs. Brézin et al., ; Brézin and Zee, although our interpretation in terms of free cumulants is somewhat more recent Nica and Speicher ; Morampudi and Laumann .

We are computing trace moments, so it is natural to represent them with a circular diagram. Factors of $A$ will be represented by the insertion of a full dot $(\bullet)$ along the circle and $B$ factors by insertions of an empty one ($\circ$).

Notice that two of the $A$ insertions are not like the others: they do not come from the expansion of $H^{n}$. Each diagram is then naturally split in two halves by those special $A$ insertions, one half associated to the $H^{n}/z^{n+1}$ term and the other one to the $H^{m}/w^{m+1}$ term.

Associate single lines (—) to factors of $1/z$ and $1/w$: lines in the upper semicircle represent factors of $1/z$ and lines in the lower semicircle factors of $1/w$.

A vertex with double line legs (=) that connects a set of insertions represents the free cumulant of those operators. Notice that all the insertions participating in a given vertex must be of the same kind, as mixed free cumulants of $A$ and $B$ vanish.

Let us summarize the diagrammatic rules we have just introduced:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Finally, in the free cumulant expansion we only need noncrossing partitions, so these diagrams are planar. The planar nature of these diagrams makes them easy to classify, which will be very useful in sections II.2.1 and II.2.4 where we compute a perturbative resummation of infinitely many diagrams.

Our aim is the following: to compute $G(z,w)$ we must sum over all planar diagrams with at least two $A$ insertions and any number of $A$ and $B$ insertions in the upper and lower halves. There are no extra combinatorial factors.

As an example of the application of these rules, let us calculate a diagram that will be included in the expansion of $G(z,w)$:

[TABLE]

In this simple case there are only vertices with one or two legs, but in general a vertex can have any number of legs.

Without loss of generality, for the rest of the paper we will assume that $\left\{A\right\}=\left\{B\right\}=0$. This does not change the dynamics, as it corresponds to a constant energy shift, but makes the diagrammatics much simpler.

II.2.1 Propagator

The full propagator of $H$, $f_{H}(z)$, is

[TABLE]

Below, we suppress the subscript $H$ when there is no risk of confusion. As usual, the complex analytic features of the propagator encode the spectrum of $H$. For example, the mean spectral density $\rho(x)$ can be obtained by inserting the Sokhotski–Plemelj formula

[TABLE]

into the definition of the propagator,

[TABLE]

Here, $\phi(x)$ is the Hilbert transform of $\rho(x)$:

[TABLE]

II.2.2 1PI Diagrams

We write a diagrammatic representation of eq. 15 by introducing a thick line to represent $f(z)$:

[TABLE]

This series can be organized into a Schwinger-Dyson equation,

[TABLE]

where we have introduced two kinds of one–particle irreducible diagrams (1PI) – those with an $A$ or a $B$ as their outermost insertion. We remind the reader that 1PI diagrams are the amputated diagrams which cannot be disconnected by cutting a single $1/z$ line. We solve Eq. (21) by introducing the ‘self-energy’ $\Sigma_{H}(z)$,

[TABLE]

Comparing to Eq. (21), we have

[TABLE]

where

[TABLE]

Let us pause and make several connections between planar perturbation theory and free probability theory explicit. The self-energies in Eq. (25) are given algebraically by

[TABLE]

and similarly for $\Sigma_{B}$. That is, the self energy $\Sigma_{A}$ is given by the free cumulant generating function $R_{A}$ evaluated at $f_{H}(z)$. Thus, additivity of the self-energy in planar perturbation theory is equivalent to the additivity of the free cumulants of freely independent random variables Zee .

The Cartesian decomposition of $f(x)$ (eq.18) induces a similar structure in $\Sigma_{A}$ (eq.26):

[TABLE]

where $\Sigma_{A}^{R}(x)$ and $\Sigma_{A}^{I}(x)$ are real functions.

II.2.3 2PI Diagrams

It is useful at this point to introduce the analog of the two–particle irreducible diagrams familiar in the context of field theory: they are (amputated) diagrams that cannot be separated into two disconnected pieces by cutting at most one $z$ and one $w$ line. Notice that all the vertices at the boundary of a 2PI diagram must be of the same type, either $A$ or $B$.

We will denote $\pi_{A}(z,w)$ the sum of all 2PI diagrams with the most external interaction of type $A$:

[TABLE]

and define analogously $\pi_{B}(z,w)$. The series corresponding to eq. 29 is

[TABLE]

The linearity of the self–energy contributions guarantees the linearity of this quantity as well:

[TABLE]

II.2.4 Box Resummation

In the diagrammatic expansion of $G(z,w)$ we need a set of disconnected diagrams that is closely related to the 2PI we just introduced: they are composed of a $1/z$ line, a $1/w$ line and any number of vertices possibly connecting those two lines. We call the sum of this class of diagrams the box, and represent it with a gray shaded area. 111The actual shape does not matter, as long as it has four sides: one side must contain only $1/z$ lines and the opposite one only contains $1/w$ lines The resummation of these diagrams proceeds in a way analogous to what we did for the propagator.

Before writing the self–consistent equation we show a few of the diagrams that contribute to the box:

[TABLE]

where we can imagine that the first line contains all the disconnected diagrams and the second line all the connected ones.

The disconnected diagrams can be summed using the results of sec. II.2.1, while the connected diagrams can be written in terms of $\pi_{H}$ and the box itself:

[TABLE]

which can be solved for $h(z,w)$ to give

[TABLE]

We can further simplify this expression using eq. 31 and eq. 22, so that after some algebra we obtain

[TABLE]

II.2.5 Triangle Resummation

There is another pattern that is important for the computation of the two–point function: (amputated) diagrams in which one of the two special insertions of $A$ participates in a vertex with at least three legs, and it connects to both a $z$ and a $w$ line. The sum of all diagrams of this kind is

[TABLE]

which can be written in a simple form using eq. 26:

[TABLE]

As a simple consistency check, notice that if $A$ is from one of the Gaussian ensembles we have $\Sigma_{A}(z)=\{A^{2}\}f(z)$ and thus $T(z,w)$ vanishes identically. This is consistent with the diagrammatic statement that Gaussian matrices do not have vertices with more than two legs.

II.3 Two–point correlation function

We are now ready to compute the two–points function $G(z,w)$: as stated in sec. II.2, we must sum over all the circular planar diagrams with at least two $A$ insertions.

Using the objects $\Sigma_{A}$, $h$, and $T$ we can sum all diagrams in which the two special $A$ insertions do not connect to each other.

We organize the remaining diagrams by the kind of interaction the two $A$ participate in:

[TABLE]

which is the diagrammatic representation of

[TABLE]

and the series appearing in the last term can be summed as usual using eq. 26:

[TABLE]

This is our main result about the two–point function: once we fix the probability distributions of $A$ and $B$, we can compute Novak and LaCroix ; Zee $f(z)$ and $\Sigma_{A}(z)$, so they can be considered inputs to the problem.

II.4 Example: $A,B$ Gaussian

If $A$ and $B$ are sampled from the GUE ensemble, all their free cumulants vanish after the second one. Without loss of generality we can set

[TABLE]

while a convenient choice for the second free cumulant is

[TABLE]

with $0\leq\lambda\leq 1$. For these random variables the self–energy series terminates after just one term:

[TABLE]

In the infinite $N$ limit, independent GUE variables are also freely independent Nica and Speicher , so we can use eq. 23 and eq. 22 to compute the propagator for $H$:

[TABLE]

and plugging $f$ and $\Sigma_{A}$ in eq. 42 we have

[TABLE]

The real time correlator is (eq. 9):

[TABLE]

which has a closed form 222The factor $1/\lambda$ on the left–hand side is included for convenience: with the choices of this section $\left<A(0)A(0)\right>/\lambda=1$. in terms of the Bessel function $J_{1}$:

[TABLE]

III Late Time Asymptotics

III.1 General Case

The integral in eq. 9 can not be expressed in terms of elementary functions except in a few special cases. In this section we prove that for a large class of Haar–invariant ensembles we have the asymptotic result

[TABLE]

where $C_{A,H}$ is a real constant and

[TABLE]

This expression lets us compute the long–time value of the correlator, and contains the non–trivial result that the decay is controlled by the Hamiltonian only.

To derive eq. 50, we deform the integration contours in eq. 9 to run infinitesimally close to the real axis:

[TABLE]

and for convenience we define the symbol

[TABLE]

so that after parametrization the integral representation of the two–point correlator becomes

[TABLE]

Comparing with Eq. (42), it is clear that contributions to the long–time value of this integral only come from Dirac $\delta(x)$ terms in $f(z)$, and the divergence in $h(x+i\delta,y-i\delta)$ as $x$ approaches $y$ (see Eq.(36)) .

The first kind of constant is trivial, so we will assume that

[TABLE]

and focus on the second kind.

The asymptotic behavior of the inverse Fourier transform $\mathcal{F}^{-1}\{\circ\}(t)$ is determined by singularities in frequency spaceLighthill , with smaller positive powers $\alpha_{i}>0$ corresponding to slower real time decay:

[TABLE]

For many common choices of $A$ and $B$, including the Gaussian, the Orthogonal PolynomialBrézin et al. , and the Wishart ensemble, all the singularities in $\rho_{H}(x)$ are of the form

[TABLE]

so we obtain the asymptotics of eq. 54 expanding the integrand in powers of $\rho$ and integrating term by term. The expansion coefficients depend on $x$ and $y$, so they could in principle modify the singular behavior, but in hindsight we realize that this not the case.

Keeping terms up to $O(\rho(x)\rho(y))$ gives

[TABLE]

The first term is time–independent after integration, and gives an approximate value of $\left<A(\infty)A(0)\right>$, but we can do better: using eq. 16 before expanding in powers of $\rho$ and collecting all terms proportional to $\delta(x-y)$, we obtain the exact expression

[TABLE]

where we have reintroduced the subscript in $\rho_{H}(x)$ to make clear that it is the spectral density of $H$, not $A$.

III.2 Example: Asymptotics for a $\sqrt{x}$ Edge

The square–root singularity at the edges of the GUE spectral density is very common, so we believe it is useful to analyze in detail this case. Using the same conventions as in section II.4, eq. 9 reads

[TABLE]

The integrand has four branching points in $z=\pm 2$ and $w=\pm 2$ (see eq. 46), so we cut the complex plane and deform the integration contour as presented in fig. 2.

The first term in eq. 60 is factorized and can be estimated using standard 333We remind the reader that $\sqrt{-1}=\exp{(-i\pi/2)}$ with our choice of branch cut. methods:

[TABLE]

which gives

[TABLE]

Using eq. 46, we can rewrite the remaining terms in eq. 60 as

[TABLE]

During the calculation we can drop all analytic terms from the sum, as they integrate to zero on a closed contour. The factorized term gives the same result as eq. 63, and one can prove that

[TABLE]

Even without going into the details of integration, we can recognize that it must be constant, by taking the time derivative under the integral sign:

[TABLE]

In conclusion, we have proved that

[TABLE]

which is consistent with the exact result eq.49, since

[TABLE]

IV Finite temperature modifications

In order to shift from infinite to finite temperature $\beta^{-1}$, observe that the exact analysis of $G(w,z)$ is unmodified, so the results we have derived about $G(t)$ immediately transfer mutatis mutandis: it is sufficient to replace $t\to t+i\beta$ in appropriate locations, which induces some thermal reweighting of the integrals. Here, we point out the key modifications that need to be made in Secs. II and III for finite $\beta$.

At finite temperature, the dynamical autocorrelator is

[TABLE]

At large $N$, the concentration of measure allows one to split the disorder average between numerator and denominator. Defining the partition function

[TABLE]

the correlator becomes

[TABLE]

where the average on the right is the usual disordered averaged trace from Eq. (4).

Retracing the steps of Sec.II, we find that all the information about temperature disappears from $G(z,w)$, and is only contained in the integration measure:

[TABLE]

The exponential factor $e^{-\beta z}$ does not cause convergence problems if the density of states has bounded support, as in most random matrix ensembles.

In particular, the late time asymptotics of $\langle A(t)A(0)\rangle_{\beta}$ are still governed by the decomposition Eq. (LABEL:eq:Gxy_approx) of the exact $G(w,z)$. Plugging Eq. (LABEL:eq:Gxy_approx) into Eq. (73), we find the analog of Eq.(2):

[TABLE]

The late time constant $A_{\infty,\beta}^{2}$ is the same as Eq.(59) reweighted by the Boltzmann weight,

[TABLE]

Somewhat more interesting is the approach to the constant: for any finite $\beta$ the power law is the same as in infinite temperature case, while in the limit $\beta\to\infty$ the decay is slower since

[TABLE]

where $E_{0}$ is the ground state energy. Thus, if the finite temperature relaxation follows a power law $1/t^{\alpha}$, the zero temperature system relaxes as $1/t^{\alpha/2}$.

V The Gaussian Rotation Approach

If $A$ and $B$ are both Gaussian, there is an alternative approach that lets us compute any real–time correlator based on an orthogonal transformation.

When $A$ and $H$ are freely independent, we can easily compute the two–point function

[TABLE]

using the non–crossing rules between $A$ and $H$, but the problem is of course that we are interested in the case when they are not.

In section II we kept $A$ fixed and expanded $H$ in terms of $A$ and $B$, which do have a non–crossing rule. Here we do the opposite: we hold $H$ fixed, and seek a change of variables that turns $A$ into something freely independent with $H$. This is easily done in the Gaussian case, but it is not clear how to construct such a transformation for general ensembles.

As in sec.II.4, we use traceless matrices with second moments

[TABLE]

and it is convenient to extract the $\lambda$ dependence defining the unit variance variables

[TABLE]

to make the algebra in the rest of the section a little cleaner.

We define the variable $C$ through the orthogonal transformation

[TABLE]

which makes it Gaussian and independent with $H$. In the $N\to\infty$ limit, independent Gaussian variables become also freely independentNica and Speicher , so we have found a variable with the requested non–crossing rule with $H$.

The calculation of the correlator at this point is straightforward: we express $A$ as

[TABLE]

and using the free cumulant expansion we get

[TABLE]

which correctly reproduces the result we obtained through the diagrammatic formalism (eq. 49).

Using this method we can actually compute any correlator of Gaussian operators easily. For example, the out-of-time-order correlator of $A$ with itself:

[TABLE]

The free cumulant expansion gives in this case,

[TABLE]

This expression is in good agreement with numerics, as we can see from fig. 3. The short time behavior is compatible with what is known in literatureVijay and Vishwanath

[TABLE]

but we see a modification of the exponent in the long–time power law: if $A$ is independent of $H$, previous workVijay and Vishwanath finds a $t^{-4}$ decay to the infinite time value. The partial conservation of $A$ leads to a slower $1/t^{3}$ decay:

[TABLE]

VI Numerical Confirmation

In this section we check the predictions of eq. 42 in a few interesting cases. Table 1 contains a summary of the relevant functions and the resulting long–time constant, while fig. 4 compares the analytical prediction with the numerical results.

In all three examples, $A$ and $B$ are sampled from the same ensemble simply because the resulting expressions are cleaner, but this is not necessary: our results works just as well in the mixed case.

GUE

We used the conventions of section II.4.

Marchenko–Pastur (MP)

To sample an MP matrix $A$ we first sample a standard $GUE$ matrix $M$ with

[TABLE]

and then we compute

[TABLE]

The same procedure is repeated for the matrix $B$.

Binary Matrices (BIN)

These are matrices of the form

[TABLE]

where $U_{i}$ are Haar–random $N\times N$ unitaries and $D$ is a diagonal matrix filled with half $-1$ and half $+1$ values.

While the result eq. 42 is still valid for these matrices, the approximation that leads to eq. LABEL:eq:Gxy_approx breaks down: $\rho_{H}(x)$ has $x^{-1/2}$ edges (see the second column in Table 1). This means that the value of the long–time constant eq.59 is correct, as confirmed by numerics, but determining the approach requires more work.

We instead numerically compute the characteristic function $\left<e^{iHt}\right>$ and fit

[TABLE]

to the numerics. Figure 4 shows that this approximation is not as clean as in the other two cases, but after a short time the error settles to $1/N$, which is the best we can hope for.

This suggests that even though the approximation of eq. LABEL:eq:Gxy_approx is not valid in this case, the approach to the constant is still determined by the characteristic function of $H$.

VII Discussion

There are two directions along which it would be interesting to extend the current work to local, finite dimensional, ergodic quantum systems.

First, if $A$ is a local operator and $B$ is a sum of local operators which itself satisfies ETH, then we expect $H=A+B$ to satisfy ETH and the observable $A$ to be partially conserved. This implies that $\langle A(t)A(0)\rangle\to c/L$ where $L$ is the size of the extended system and $c$ is a constant quantifying how conserved $A$ is. This can be computed explicitly using Eq. (59) and the results tested against ETH systems.

More technically challenging is to extend the analysis here to chains of locally interacting random matrices where one might hope to compute the energy diffusion constant explicitly from the dynamical correlators. Here the exact resummations available in the random matrix case are complicated by the locality structure of the chain. Some technical steps along this axis have been developed in Ref. Morampudi and Laumann, .

The order of limits is important: $N\to\infty$ must be taken before $t\to\infty$. At finite $N$, we expect corrections of order $1/N$ to the late time value of $\langle A(t)A(0)\rangle_{c}$, though we have not computed them. They can be calculated perturbatively by resumming diagrams that tessellate a torus with a hole. For $A$ independent of $H$, the $1/N$ corrections can be computed non-perturbatively from dephasing the spectral representation:

[TABLE]

where $\alpha$ runs over the energy eigenbasis. However, for $A$ part of $H$, this ‘diagonal ensemble’ calculation is not straightforward, as the $\ket{\alpha}$ are correlated with $A$. Indeed, these correlations must produce both the $O(1/N^{0})$ late time value which we have computed and any $O(1/N)$ corrections.

Acknowledgements.

The authors would like to thank A. Chandran and A. Polkovnikov for stimulating discussions. C.R.L. acknowledges support from the NSF through grant PHY-1752727. This work was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611, and at the Galileo Galilei Institute in Florence. Any opinion, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the NSF.

Appendix A Short time expansion

In the GUE case, the two–point correlation function can also be computed summing its short time expansion: we write $G(t)$ as a power series in the Liouvillian hyperoperator $\mathcal{L}\coloneqq[H,\cdot]$.

[TABLE]

and evaluate explicitly $\left<\left(\mathcal{L}^{n}\tilde{A}\right)\tilde{A}\right>$.

These moments exhibit a clear pattern OEIS Foundation Inc. : for odd powers of $\mathcal{L}$ the expression vanishes, while for even powers we have

[TABLE]

where $C_{n}$ is the n–th Catalan number. The short time series is then

[TABLE]

and it can be summed:

[TABLE]

which agrees with the result presented in the main text.