^a^ainstitutetext: Department of Physics, Stanford University, Stanford, California 94305, USA

Operator Krylov complexity in random matrix theory

Abstract

Krylov complexity, as a novel measure of operator complexity under Heisenberg evolution, exhibits many interesting universal behaviors and also bounds many other complexity measures. In this work, we study Krylov complexity $\mathcal{K}(t)$ in Random Matrix Theory (RMT). In large $N$ limit: (1) For infinite temperature, we analytically show that the Lanczos coefficient $\{b_{n}\}$ saturate to constant plateau $\lim\limits_{n\rightarrow\infty}b_{n}=b$ , rendering a linear growing complexity $\mathcal{K}(t)\sim t$ , in contrast to the exponential-in-time growth in chaotic local systems in thermodynamic limit. After numerically comparing this plateau value $b$ to a large class of chaotic local quantum systems, we find that up to small fluctuations, it actually bounds the $\{b_{n}\}$ in chaotic local quantum systems. Therefore we conjecture that in chaotic local quantum systems after scrambling time, the speed of linear growth of Krylov complexity cannot be larger than that in RMT. (2) For low temperature, we analytically show that $b_{n}$ will first exhibit linear growth with $n$ , whose slope saturates the famous chaos bound. After hitting the same plateau $b$ , $b_{n}$ will then remain constant. This indicates $\mathcal{K}(t)\sim e^{2\pi t/\beta}$ before scrambling time $t_{*}\sim O(\beta\log\beta)$ , and after that it will grow linearly in time, with the same speed as in infinite temperature. We finally remark on the effect of finite $N$ corrections.

Keywords:

Krylov complexity, Random Matrix Thoery, Scrambling

1 Introduction

Operator complexity describes the phenomenon that a simple operator $O$ becomes complex under Heisenberg evolution $O(t)$ in chaotic local quantum systems. Many complexity measures Roberts2017 ; Cotler2017chaos ; Jefferson2017 ; Roberts2018 ; Yang2018 ; Khan2018 ; Qi2019 ; Lucas2019 ; Balasubramanian2020 ; Balasubramanian2021 make this manifest. An intuitive one is the operator size and OTOC Roberts2018 ; Qi2019 ; Lucas2019 , which counts the average size of spatial support for an operator written in the basis made from tensor product of simple operators. Since operator size is bounded by the size of the whole system, then one would expect that size would cease to grow and saturate soon after it spreads onto the whole system at the time scale of scrambling time $t_{*}$ Maldacena2016bound . However, since the Heisenberg evolution proceeds, one would expect that the operator is still growing more and more complex, simply not in the sense of the size of its spatial support.

Behaviour of Krylov complexity in generic chaotic local systems.

Krylov complexity Parker2019 , denoted as $\mathcal{K}(t)$ , proposed as a measure of operator complexity, achieves this goal. Namely, it successfully captures the entire time scale and different stages of complexity dynamics. In generic chaotic local quantum systems and starting from a simple operator, the dynamics of $\mathcal{K}(t)$ often exhibit three stages:

1.

Stage one: exponential-in-time growth. $\mathcal{K}(t)\sim e^{2\alpha t}$ . This exponential growth will extend up to scrambling time $t_{*}\sim\log(S)$ and $\mathcal{K}(t)$ reaches the value of order $O(S)$ , where $S$ is the number of degree of freedom. The complexity growth in this stage is mainly due to spreading a simple operator to the whole system, therefore describing the same physics as in operator size. The temperature-dependent exponential coefficient $\alpha$ is related and actually bounds the Lyapunov exponential of OTOC Maldecena2016 ; Maldacena2016bound .
2.

Stage two: linear-in-time growth. $\mathcal{K}(t)\sim vt$ . Up to scrambling time in stage one, the growth of ‘complexity in real space’ has finished. In this linear growth region, $\mathcal{K}(t)$ probes the growth of ‘complexity in Hilbert space’. Roughly speaking, $\mathcal{K}(t)$ in this stage characterizes the fact that $O(t)$ explores more and more linearly independent basis that it has not been occupied before in the operator vector space as large as $O(e^{2S})$ . This stage will last up to a time scale $O(e^{2S})$ where $\mathcal{K}(t)$ reaches the value of order $O(e^{2S})$ Rabinovici2021 , the dimension of operator vector space, as claimed.
3.

Stage three: remain constant. $\mathcal{K}(t)\sim\text{Const}$ . At the end of linear-in-time stage, $O(t)$ has explored the edge of Krylov subspace and occupied as much linearly independent operator basis as it can, therefore it will remain constant afterward. This constant value is precisely $\frac{1}{2}K$ Rabinovici2021 , which is the half of the Krylov subspace dimension $K\sim O(e^{2S})$ . The prefactor $\frac{1}{2}$ means that $O(t)$ is evenly distributed on every basis of Krylov subspace, indicating that it is fully scrambled. There might be possible Poincaré recurrence at the time of order $O(e^{e^{2S}})$ corresponding to some doubly-non-perturbative physics that is beyond our scope Rabinovici2021 .

Refer to caption — Figure 1: Schematic plot of the main result of this work. In (a) and (b), we show the behavior of Krylov complexity $\mathcal{K}(t)$ and Lanczos coefficient $\{b_{n}\}$ in Random Matrix Theory in large $N$ limit with finite or infinite temperature. The plateau value equating unity comes from normalization of radius of Wigner semicircle to be unity. In (c) and (d), we compare the random matrix results with generic chaotic local quantum systems.

Motivation of our work.

In this paper, we mainly focus on the linear-in-time stage, and sometime on exponential-in-time stage. Our motivation is that, at the conjunction of stage one and two, $O(t)$ has just finished its ‘scrambling in real space’, and begun to start its ‘scrambling in Hilbert space’. Therefore in stage two, $O(t)$ completely abandoned any spatial locality structure, or equivalently, the local Hamiltonian $H$ no longer looks local from the perspective of $O(t)$ . Therefore, in this stage, $H$ would act like a chaotic Hamiltonian without spatial locality, resembling a typical random matrix.

This is our main reasoning and motivation for studying Krylov complexity in Random Matrix Theory (RMT), which intuitively should dominate the behavior of complexity growth after scrambling time. This is consistent with the fact that in stage one, the exponential growth in time crucially comes from the argument of the locality of $O(t)$ and $H$ Parker2019 . Without locality, we would expect the growth rate to decrease to linear.

Another important motivation is that the linear growth of complexity characterizes the expanding volume of the Einstein-Rosen bridge in the dual holographic gravity susskind2014computational . We expect our RMT analysis would shed light on bulk gravitational systems.

Main results and structure of the paper.

In this paper, we obtain several analytical results of Krylov complexity in a specific random matrix ensemble, called GUE (Gaussian Unitary Ensemble) in large $N$ limit, thanks to the large $N$ factorization.

In section 2, we review the construction of Krylov basis which originates from an orthogonalization procedure, the definition of Krylov complexity, and associated Lanczos coefficient $\{b_{n}\}$ .

In section 3.1, we study the Krylov complexity of arbitrary traceless initial operator $O$ , and we show that if we properly scale the radius of Wigner semicircle to unity, for Lanczos coefficient $\{b_{n}\}$ , we have:

\lim\limits_{n\rightarrow\infty}\lim\limits_{N\rightarrow\infty}b_{n}=1

(1)

This indicates that in large $N$ limit, $\mathcal{K}(t)$ in RMT grows linearly in time, in contrast to the chaotic local quantum systems in thermodynamic limit, whose K-complexity exhibits an exponential growth (in thermodynamic limit, scrambling time $t_{*}$ goes to infinity, therefore chaotic local system only have stage one, without stage two and three). We remark that this does not mean K-complexity in RMT is smaller than that in chaotic local systems. This is because, in RMT theory, we usually scale the Hamiltonian such that its eigenenergy is of order $O(1)$ . But in chaotic local systems, the energy spectrum is extensive which is of order $O(S)$ . Since the energy scale is conjugate to the time scale, we should properly rescale our RMT result.

Indeed, in section 3.2, we numerically compared the Lanczos coefficient $\{b_{n}\}$ in generic chaotic local systems in stage two with that from properly rescaled RMT, and find that the former is generally of the same order and smaller than the latter up to fluctuation. This indicates that the stage two growth in chaotic local systems originates the onset of random matrix behavior, as claimed in the motivation. Due to various empirical findings from numerics and theoretic reasoning, we therefore conjecture that stage two growth in chaotic local systems is bounded by scaled RMT growth. This is schematically summarized in figure 1(c), (d).

Nevertheless, the complexity dynamics of RMT is by itself interesting. In section 3.3, still working in large $N$ limit, we deviate from infinite temperature and study K-complexity in RMT at finite and especially low temperature. We analytically show that $\{b_{n}\}$ will first grow linear in $n$ with slope $\alpha$ proportional to temperature $T$ , until it reaches the plateau value $b_{n}=1$ , the same as in the case of infinite temperature. Through numerical simulations, we fix $\alpha=\pi\beta^{-1}$ ( $\beta=T^{-1}$ is inverse temperature). Translating the behavior of $\{b_{n}\}$ to the dynamics of K-complexity, we find that $\mathcal{K}(t)$ will first grow exponentially with the Lyapunov exponent saturating the chaos bound Maldecena2016 ; Maldacena2016bound ; Parker2019 , and after a scrambling time of order $O(\beta\log\beta)$ , $\mathcal{K}(t)$ switch to linear growth in time, with the velocity the same as if it’s in infinite temperature. This is summarized schematically in figure 1(a), (b).

In section 4 we remark on the finite $N$ corrections and in section 5, we summarize our result and discuss some possible future directions.

Related works.

In reference Bhattacharyya2023 ; Erdmenger2023 , the authors studied the ‘state’ Krylov complexity generated by random matrix Hamiltonian. Namely, they expand $|\psi(t)\rangle=e^{-iHt}|\psi_{0}\rangle$ on the Krylov basis which is constructed from the Lanczos orthogonalization of $\{H^{n}|\psi_{0}\rangle\}$ . Despite having the same name, in our work we actually study ‘operator’ Krylov complexity generated by Heisenberg evolution of RMT, which is a different quantity. In reference Kar2022 , the authors considered an ‘energy-refined’ version of Krylov complexity and its relation with RMT universality class. Though bearing similar motivations to us, they studied different quantities from ours.

2 Review of Krylov complexity

In this section, we briefly review the definition of Krylov complexity and construction of Krylov basis Parker2019 . Experts on this topic may skip this section and directly dive to the next one with new results.

Given an Hermtian initial operator $O\in\mathcal{B}(\mathcal{H})$ which is an bounded operator on Hilbert space $\mathcal{H}$ with dimension $N$ , its Heisenberg evolution is given by $O(t)=e^{iHt}Oe^{-iHt}$ , where $H$ is the Hamiltonian governing the time evolution of the quantum systems. Since we are interested in the dynamics of $O(t)$ , it is convenient to work in the vector space $\mathcal{B}(\mathcal{H})$ of operator, with dimension $N^{2}$ . This is achieved through operator-to-state mapping: $O\rightarrow|O)$ .

In order to have a notion of normalized and orthogonal vectors, we need to equip this operator linear space with an inner product structure $(A|B)$ , defined as

(A|B)\equiv\begin{cases}\operatorname{tr}[A^{\dagger}B],\ &\beta=0\\ \operatorname{tr}[e^{-\beta H/2}A^{\dagger}e^{-\beta H/2}B],\ &\beta\neq 0\end{cases}

(2)

for infinite temperature and finite temperature respectively. The motivation for inserting thermal density matrix into finite temperature inner product is because it resembles the thermal Wightman two-point function Parker2019 .

The Heisenberg evolution is also mapped to a superoperator, called Liouvillian $\mathcal{L}$ , defined as the commutators of Hamiltonian: $\mathcal{L}O\equiv[H,O]$ . One can easily show that such an inner product would render Livioullian a Hermitian super-operator, namely $(A|\mathcal{L}B)=(\mathcal{L}A|B)$ . By Baker-Hausdoff-Campell formula, the time evolution is thus expanded in the basis of $\{|\mathcal{L}^{n}O)\}$ :

|O(t))=e^{i\mathcal{L}t}|O)=\sum\limits_{n=0}^{\infty}\frac{(it)^{n}}{n!}|\mathcal{L}^{n}O)

(3)

As we claimed in section 1, the idea of Krylov complexity is that can quantify the complexity ‘beyond real space’: the complexity in Hilbert space. Namely, we expect the growth of $\mathcal{K}(t)$ to characterize the process in which $O(t)$ explores more and more linearly independent basis in $\mathcal{B}(\mathcal{H})$ that it has not occupied before. Certainly $\{|\mathcal{L}^{n}O)\}$ is linearly independent in general, however, they are neither normalized nor orthogonal. Therefore, to make this idea explicit, we need to perform an orthogonalization procedure on $\{|\mathcal{L}^{n}O)\}$ and generate an orthonormal basis $\{|O_{n})\}$ .

This procedure is the Lanczos algorithm. Starting from $|O)$ , define a normalized vector $|O_{0})=(O|O)^{-1/2}|O)$ . The second step is to define $|O_{1})=b_{1}^{-1}\mathcal{L}|O_{0})$ where $b_{1}=(\mathcal{L}O_{0}|\mathcal{L}O_{0})^{1/2}$ . Then inductively define:

$\displaystyle\|A_{n})$	$\displaystyle=\mathcal{L}\|O_{n-1})-b_{n-1}\|O_{n-2})$	(4)
$\displaystyle b_{n}$	$\displaystyle=(A_{n}\|A_{n})^{1/2}$
$\displaystyle\|O_{n})$	$\displaystyle=b_{n}^{-1}\|A_{n})$

The output of this algorithm is a set of orthonormal basis $\{|O_{n})\}$ satisfying $(O_{n}|O_{m})=\delta_{nm}$ and a set of real and positive coefficient $\{b_{n}\}$ called Lanczos coefficient. One also notice that since $O$ is Hermitian, then $|O_{n})$ is also a Hermitian operator which can be expressed as a degree- $n$ -polynomial of $\mathcal{L}$ acting on $|O)$ .

One also notices that the matrix of super-operator $\mathcal{L}$ written in the basis $\{|O_{n})\}$ is a tridiagonal real and Hermitian matrix with zero diagonal entries and its off-diagonal entries equating $\{b_{n}\}$ , namely:

L_{mn}=(O_{m}|\mathcal{L}|O_{n})=\begin{bmatrix}0&b_{1}&0&0&\cdots\\ b_{1}&0&b_{2}&0&\cdots\\ 0&b_{2}&0&b_{3}&\cdots\\ 0&0&b_{3}&0&\ddots\\ \vdots&\vdots&\vdots&\ddots&\ddots\end{bmatrix}

(5)

We also briefly remark that the dimension of Krylov space (spanning of $\{|O_{n})\}$ ) denoted as $K$ , is of the same order but slightly smaller than the dimension of $\mathcal{B}(\mathcal{H})$ . This is because Hamiltonian commutes with any function of itself. In other words, $\{|H^{n})\}$ spans the kernel space of $\mathcal{L}$ . Therefore even if initial operator $|O)$ may have non zero occupation on $\text{ker}(\mathcal{L})$ , all other $\{\mathcal{L}^{n}|O),n\geq 1\}$ are perpendicular to $\text{ker}(\mathcal{L})$ . Therefore at most we have $K\leq N^{2}-\text{dim}(\text{ker}(\mathcal{L}))+1$ . A careful examination shows that Rabinovici2021 $K\leq N^{2}-N+1$ . Nevertheless, in this paper, we will not probe the physics near the edge of Krylov space.

Now we are ready to study the dynamics of $O(t)$ by expanding $|O(t))$ in the basis $\{O_{n}\}$ with coefficient $\varphi_{n}(t)$ :

|O(t))=\sum\limits_{n=0}^{K-1}\varphi_{n}(t)|O_{n})

(6)

Then $\varphi_{n}(t)$ satisfies the following Schrodinger-like equation:

i\partial_{t}\varphi_{n}=-b_{n+1}\varphi_{n+1}-b_{n}\varphi_{n-1}

(7)

which describes a single particle initially at the left-end and then starts hopping on the one dimensional (semi-infinite, if we imagine $K$ is large) Krylov chain with non-uniform position-dependent hopping coefficients $\{b_{n}\}$ .

The Krylov complexity is then defined as the average position of this wave packet:

\mathcal{K}(t)\equiv\sum\limits_{n=0}^{K-1}n|\varphi_{n}(t)|^{2}

(8)

This definition is very intuitive compared with our daily life experience of defining complexity: To handle a complex task (say, build a house), we try to decompose it into a sequence of simple tasks (say, adding bricks), which we manage to deal with step by step. Then we use the number of steps to denote the overall complexity (say, the number of bricks). Many complexity measures in quantum systems follow this spirit. For instance, in quantum computing we often use circuit complexity to quantify a complex quantum state $|\psi\rangle$ , which is defined as the number of elementary gates that are needed to assemble into a large unitary rotating simple state $|\psi_{0}\rangle$ (often a product state with zero entanglement entropy) to $|\psi\rangle$ . For our case of Krylov complexity, the simple operation is just the Heisenberg commutator $\mathcal{L}=[H,\cdot]$ and K-complexity counts the average number of operations we needed to make a simple $O$ to be as complex as $O(t)$ . This analogy is summarized in table 1.

	Simple initial object	Simple operation	Complex final object	Complexity measure
Circuit complexity	$\|\psi_{0}\rangle$ product state	$\{U_{i}\}$ sets of universal gates	$\|\psi\rangle$ highly entangled	number of operations
Krylov complexity	$O$ local operator	$\mathcal{L}=[H,\cdot]$ Heisenber commutator	$O(t)$ highly scrambled	average number of operation

Table 1: Analogy between circuit complexity (as a state complexity measure shown above, but can also be readily generalized to measure operator complexity) and Krylov complexity (as an operator complexity measure). The universal set of gates

\{U_{i}\}

are often one-qubit gates and two-qubit gates. One choice is all single qubit rotation together with CNOT gate.

We also remark on the difference between circuit complexity and Krylov complexity. For the former, the choice of ‘simple operation’, namely a set of universal one-qubit/two-qubit gates, has ambiguity. But for the latter, the choice of commutator is very natural, since the dynamics is generated by it in the first place.

Going back to $\mathcal{K}(t)$ , let’s find out the relation between $\{b_{n}\}$ and $\mathcal{K}(t)$ through some intuitive examples. The first example is that $b_{n}\equiv b$ is constant. In this case, the hopping coefficient is uniform so we expect the center of the wave packet travel in a constant velocity proportional to $b$ . Therefore $K(t)$ growth linear in time.

The next example of interest is to consider $b_{n}=\alpha n+\gamma$ . Before that, we first want to take the continuous limit of the wave equation. Define $\phi_{n}(t)=i^{-n}\varphi_{n}(t)$ to make $\phi_{n}(t)$ satisfy a real wave equation $\partial_{t}\phi_{n}(t)=-b_{n+1}\phi_{n+1}+b_{n}\phi_{n-1}$ with initial condition $\phi_{n}(0)=\delta_{n0}$ . We define a continuous position variable $x=\epsilon n$ , where $\epsilon$ plays the role of lattice constant. Then $\epsilon\rightarrow 0$ limit results in the following continuous equation Jian2021 :

(\partial_{t}+2\alpha x\partial_{x}+\alpha)\phi(x,t)=0

(9)

Given initial distribution $\phi_{0}(x)$ , we obtain the unique solution:

\phi(x,t)=e^{-\alpha t}\phi_{0}(xe^{-2\alpha t})

(10)

This shows that up to an overall damping factor $e^{-\alpha t}$ which keeps the normalization of wavefunction the same, the initial wave packet is stretched with factor $e^{2\alpha t}$ . Therefore the average position $\langle x\rangle_{(t)}=e^{2\alpha t}\langle x\rangle_{(0)}$ will also grow exponentially, indicating a $\mathcal{K}(t)\sim e^{2\alpha t}$ .

We also notice that $\mathcal{K}(t)$ is fully determined by two-point auto-correlation function $C(t)\equiv(O|O(t))$ . This is because $\mathcal{K}(t)$ is determined by $\{b_{n}\}$ , which is in turn determined by a set of moments $\{\mu_{2n}\equiv(O|\mathcal{L}^{2n}|O)\}$ , the derivatives of $C(t)$ at $t=0$ . It is also worth remarking the relation between $\{\mu_{2n}\}$ and $\{b_{n}\}$ is highly non-linear Parker2019 :

\mu_{2n}=(L^{2n})_{00}=\sum\limits_{\text{path}}b_{i_{1}}b_{i_{2}}\cdots b_{i_{2n}}

(11)

where $L$ is the tridiagonal matrix of $\mathcal{L}$ written in Krylov basis, defined in equation (5). Here, a Dyck path is defined by starting from $n=0$ and return $n=0$ with exactly $2n$ steps. The number of such paths with $2n$ steps is given by the Catalan number $C_{n}=\frac{(2n)!}{(n+1)!n!}$ . The relation from $\{\mu_{2n}\}$ to $\{b_{n}\}$ is also complicated Parker2019 :

b_{1}^{2}\cdots b_{n}^{2}=\operatorname{det}(\mu_{i+j})_{0\leq i,j\leq n}

(12)

Given $\{b_{n}\}$ , in order to gain some intuition about the scaling behavior of $\{\mu_{2n}\}$ wrt $n$ , which is important for our discussion in later sections, there is an explicit inequality bound for $\mu_{2n}$ if we assume $b_{n}$ is non-decreasing with $n$ :

C_{n}(b_{1}^{2}b_{2}^{2}\cdots b_{n}^{2})\geq\mu_{2n}>(b_{1}^{2}b_{2}^{2}\cdots b_{n}^{2})

(13)

This is because $(b_{1}^{2}b_{2}^{2}\cdots b_{n}^{2})$ is the maximal weight among all Dyck paths. For physically relevant cases, $b_{n}$ schematically has the common form $b_{n}\sim\alpha n^{\delta}$ . $\delta=1,\frac{1}{2},0$ respectively correspond to Parker2019 : chaotic local systems at stage one; integrable systems at stage one; free systems at stage one or chaotic local systems at stage two. Therefore, considering $n$ to be large but finite, we may apply Stirling formula to equation (13):

2\delta n\log n+\left(2\log 2+2\log\alpha-2\delta\right)n+O(\log n)\geq\log\mu_{2n}>2\delta n\log n+\left(2\log\alpha-2\delta\right)n+O(\log n)

(14)

If we further assume that $\log\mu_{2n}=A_{1}n\log n+A_{2}n+O(\log n)$ with the same scaling, then the leading-order coefficient $A_{1}$ completely fix $\delta$ and next-leading-order coefficient $A_{2}$ provide bound on $\alpha$ .

Ever since the concept of Krylov complexity has been proposed Parker2019 , it has attracted considerable attention from various communities since it exhibits many universal behaviors in chaotic local systems described in detail in section 1. People studied this quantity in many contexts ranging from condensed matter models focusing on scrambling behaviors Barbón2019 ; Dymarsky2020 ; Tianrui2020 ; Chen2021 ; Noh2021 ; Fabian2022 ; Kim2022 ; Hörnedal2022 ; Bhattacharjee2022 ; Sinong2022 ; Zhongying2022 ; Fan2022 ; Rabinovici2021 ; Rabinovici20221 ; Rabinovici20222 ; Erdmenger2023 ; Kar2022 , open quantum systems Liu2023 ; Bhattacharya2022 , field theory Dymarsky2021_CFT ; Caputa2021 ; caputa2021geometry ; Adhikari2022 ; ADHIKARI2023116263 , and holographic gravity Avdoshkin2020 ; Barbón2020 ; Magán2020 ; MUCK2022115948 ; Banerjee2022 ; Kar2022 ; Jian2021 . Now we are ready to add one more analytically tractable example, in RMT, to this exciting field with flourishing development.

3 Large $N$ limit of K-complexity in RMT

3.1 Infinite temperature

In this section, we analytically calculate Krylov complexity of RMT in large $N$ limit at infinite temperature. Our strategy is to first calculate $\{\mu_{2n}\}$ using RMT techniques. Ideally one would apply equation (12) to get an expression of $\{b_{n}\}$ , but it seems hopeless. Our strategy is to focus on the scaling behavior of $\mu_{2n}$ for large but finite $n$ , then appeal to equation (11) and inequality (14) to estimate the scaling behavior of $\{b_{n}\}$ with respect to $n$ .

For readers’ convenience, we first recall some definitions about Krylov complexity:

\mu_{2n}=(O|\mathcal{L}^{2n}|O),\ \mathcal{L}=[H,\cdot]

(15)

where the inner product in infinite temperature is defined as $(A|B)=\operatorname{tr}\left[A^{\dagger}B\right]$ . We also notice that $\mu_{2n+1}=0$ since $\mathcal{L}^{2n+1}O$ is anti-Hermitian therefore its trace with Hermitian matrix $O$ is zero. Using identity:

\mathcal{L}^{2n}O=\left([H,\right)^{2n}O]=\sum\limits_{k=0}^{2n}C_{2n}^{k}(-1)^{k}H^{k}OH^{2n-k}

(16)

we derive the formal equation for $\mu_{2n}$ :

\mu_{2n}=\sum\limits_{k=0}^{2n}C_{2n}^{k}(-1)^{k}\operatorname{tr}[OH^{k}OH^{2n-k}]

(17)

For simplicity, we consider the simplest type of RMT, which is an $N\times N$ random matrix from GUE (Gaussian Unitary Ensemble). In GUE, the probabilistic distribution of $H$ is symmetric under any unitary transformation: $P(H)=P(UHU^{\dagger})$ . Therefore, the moment $\mu_{2n}(O)=\mu_{2n}(UOU^{\dagger}),\ \forall U$ , which means that it finally depends only on the spectrum of $O$ .

We can compare the case in Bhattacharyya2023 ; Erdmenger2023 , where the author there studied the state complexity rather than operator complexity generated by GUE Hamiltonian. In their case, the initial state independence exist after average: $|\psi\rangle\rightarrow U|\psi\rangle$ means all initial states give the same Lanczos coefficient.

To proceed, we use the above mentioned fact that $P(H)=P(UHU^{\dagger})$ , therefore we can do the following replacement:

\mathbb{E}\operatorname{tr}[OH^{k}OH^{2n-k}]=\mathbb{E}\operatorname{tr}[OUH^{k}U^{\dagger}OUH^{2n-k}U^{\dagger}],\ \forall U

(18)

where notation $\mathbb{E}$ means taking average over GUE Hamiltonian. Since the above equation is valid for arbitrary $U$ , we can first keep $H$ fixed and average over $U$ under Haar measure, using Weingarten function Roberts2017 ; Cotler2017chaos :

	$\displaystyle\int dU\ U_{i_{1}j_{1}}U_{i_{2}j_{2}}U^{}_{i_{3}j_{3}}U^{}_{i_{4}j_{4}}=$	$\displaystyle\frac{\delta_{i_{1}i_{3}}\delta_{i_{2}i_{4}}\delta_{j_{1}j_{3}}\delta_{j_{2}j_{4}}+\delta_{i_{1}i_{4}}\delta_{i_{2}i_{3}}\delta_{j_{1}j_{4}}\delta_{j_{2}j_{3}}}{N^{2}-1}$		(19)
		$\displaystyle-\frac{\delta_{i_{1}i_{3}}\delta_{i_{2}i_{4}}\delta_{j_{1}j_{4}}\delta_{j_{2}j_{3}}+\delta_{i_{1}i_{4}}\delta_{i_{2}i_{3}}\delta_{j_{1}j_{3}}\delta_{j_{2}j_{4}}}{N^{3}-N}$		(19)

Then we arrive at the following useful formula:

	$\displaystyle\int dU\cdot\operatorname{tr}[OUAU^{\dagger}OUBU^{\dagger}]=\frac{1}{N^{2}-1}\bigg{(}$	$\displaystyle\operatorname{tr}[O^{2}]\operatorname{tr}[A]\operatorname{tr}[B]+\operatorname{tr}[O]^{2}\operatorname{tr}[AB]$		(20)
		$\displaystyle-N^{-1}\operatorname{tr}[O]^{2}\operatorname{tr}[A]\operatorname{tr}[B]-N^{-1}\operatorname{tr}[O^{2}]\operatorname{tr}[AB]\bigg{)}$		(20)

For simplicity we can consider a traceless normalized operator with $\operatorname{tr}[O]=0,\ \operatorname{tr}[O^{2}]=(O|O)=1$ , then the moment $\mu_{2n}$ is given by:

$\displaystyle\mathbb{E}\ \mu_{2n}$	$\displaystyle=\mathbb{E}\ \frac{1}{N^{2}-1}\sum\limits_{k}C_{2n}^{k}(-1)^{k}\bigg{(}\operatorname{tr}[H^{k}]\operatorname{tr}[H^{2n-k}]-N^{-1}\operatorname{tr}[H^{2n}]\bigg{)}$	(21)
	$\displaystyle\approx\mathbb{E}\ N^{-2}\sum\limits_{k}C_{2n}^{k}(-1)^{k}\operatorname{tr}[H^{k}]\operatorname{tr}[H^{2n-k}]$
	$\displaystyle=\mathbb{E}N^{-2}\sum_{i,j}(\lambda_{1}-\lambda_{j})^{2n}$

where $\lambda_{i}$ is the eigenvalue of $H$ . For later convenience, in the second line, we perform approximation by only retaining the leading order terms in large $N$ limit.

As a side remark, the second line of equation (21) can be derived using another approach. Instead of insisting that $O$ is a fixed operator, we can otherwise consider $O$ to be also a random Hermitian matrix under independent GUE distribution from $H$ :

\overline{O_{ij}O_{mn}}=N^{-2}\delta_{in}\delta_{jm}

(22)

where the variance $N^{-2}$ is determined by normalization condition:

\overline{(O|O)}=\overline{\operatorname{tr}[O^{2}]}=\sum\limits_{i,j=1}^{N}\overline{O_{ij}O_{ji}}=1

(23)

With the GUE distribution of $O$ we have identity $\overline{\operatorname{tr}[OAOB]}=N^{-2}\operatorname{tr}[A]\operatorname{tr}[B]$ . Applying this to equation (17), we have:

\overline{\mu_{2n}}=N^{-2}\sum\limits_{k}C_{2n}^{k}(-1)^{k}\operatorname{tr}[H^{k}]\operatorname{tr}[H^{2n-k}]

(24)

which is similar to equation (21).

Nevertheless, we will proceed by adopting that $O$ is a fixed operator. Starting from equation (21), we see that $\mathbb{E}\mu_{2n}$ is identically the same for any traceless initial operator $O$ , and only depend on the spectrum of $H$ :

$\displaystyle\mathbb{E}\ \mu_{2n}$	$\displaystyle=\mathbb{E}\ N^{-2}\sum\limits_{k}C_{2n}^{k}(-1)^{k}\operatorname{tr}[H^{k}]\operatorname{tr}[H^{2n-k}]$	(25)
	$\displaystyle=\mathbb{E}\ N^{-2}\sum\limits_{k,i,j}C_{2n}^{k}(-1)^{k}\lambda_{i}^{k}\lambda_{j}^{2n-k}$
	$\displaystyle=\mathbb{E}\ N^{-2}\sum\limits_{i,j}(\lambda_{i}-\lambda_{j})^{2n}$
	$\displaystyle=N^{-2}\int d\lambda d\lambda^{\prime}\cdot(\lambda-\lambda^{\prime})^{2n}\rho(\lambda,\lambda^{\prime})$

where $\rho(\lambda,\lambda^{\prime})$ is the averaged two-point spectral density of RMT.

Here we shall pause for a moment and make two relevant remarks.

1.

Order of taking random average: The order of average we take here is the same as that when people study Krylov complexity in SYK model Maldecena2016 ; Sachdev2015 . There, people first calculated random averaged two-point auto-correlation function $\mathbb{E}\ C(t)=\mathbb{E}(O|O(t))$ . Then by taking derivatives we arrive at disorder averaged $\mathbb{E}\mu_{2n}$ . The Krylov coefficient $b_{n}=f(\mathbb{E}\mu)$ is obtained subsequently using equation (12), which we schematically denote as $f(\cdot)$ . Of course, we can directly perform average over $\mathbb{E}b_{n}$ , and since the fluctuation in RMT is controllably small in large $N$ limit, we expect the order of average to not affect the mean value,i.e., $\mathbb{E}f(\mu)=f(\mathbb{E}\mu)$ , even if $f(\cdot)$ is a highly nonlinear map. Nevertheless, in the right panel of figure 2, we evaluate the statistical variance of $\{b_{n}\}$ from numerical simulation of random matrix, and indeed verified that the variance is suppressed by some positive power of $N^{-1}$ . Therefore, for notational symplicity, we will sometimes drop the average symbol $\mathbb{E}$ .
2.

Average over initial operator $O$ : In some recent papersJafferis20231 ; Jafferis20232 , people find that if the ensemble of $O$ and $H$ are sampled independently, the averaged $n$ -point correlator of $O$ cannot reproduce the known results of RMT for $n\geq 4$ . Here in our story, every data comes from the two-point auto-correlation function of $O$ , therefore this possible subtlety will not show up in our story. Nevertheless, to avoid any possible ambiguity, we will proceed by adopting that $O$ is a fixed operator.

Going back to the main theme, we shall proceed by first introducing some notation and knowledge in random matrix theory. The spectral one point function $\rho(\lambda)$ and two point function $\rho(\lambda,\lambda^{\prime})$ in GUE is defined by:

	$\displaystyle\rho(\lambda)=\mathbb{E}\sum_{i}\delta(\lambda-\lambda_{i})=R_{1}(\lambda)$		(26)
	$\displaystyle\rho(\lambda,\lambda^{\prime})=\mathbb{E}\sum_{i,j}\delta(\lambda-\lambda_{i})\delta(\lambda^{\prime}-\lambda_{j})=\delta(\lambda-\lambda^{\prime})R_{1}(\lambda)+R_{2}(\lambda,\lambda^{\prime})$		(27)

where $R_{1},R_{2}$ are one-point and two-point spectral correlators, following the standard notation in RMT mehta2004random . For spectrum one point function $R_{1}(\lambda)$ , it is given by the famous Wigner semicircle in large $N$ limit:

R_{1}(\lambda)=\begin{cases}2\pi^{-1}N\sqrt{1-\lambda^{2}}&\text{for}\ |\lambda|<1\\ 0&\text{for}\ |\lambda|>1\end{cases}

(28)

where we scaled the probability distribution of GUE to be $P(H)\propto\exp(-2N\operatorname{tr}H^{2})$ in order to set the radius of Wigner semicircle to be unity.

For spectrum two-point function $R_{2}(\lambda,\lambda^{\prime})$ , it has disconnected part $R_{1}(\lambda)R_{2}(\lambda^{\prime})$ and connected part $T_{2}(\lambda,\lambda^{\prime})$ related by:

R_{2}(\lambda,\lambda^{\prime})=R_{1}(\lambda)R_{1}(\lambda^{\prime})-T_{2}(\lambda,\lambda^{\prime})

(29)

The generic result for $T_{2}(\lambda,\lambda^{\prime})$ is complicated, and its simplification appears only when we zoom into the center of the Wigner circle with a few windows with width determined by mean level spacing $O(N^{-1})$ , i.e., we take the limit where $N\rightarrow+\infty,N\lambda\rightarrow\text{Const},N\lambda^{\prime}\rightarrow\text{Const}$ :

T_{2}(\lambda,\lambda^{\prime})\approx(2\pi^{-1}N)^{2}\left[\frac{\sin 2N(\lambda-\lambda^{\prime})}{2N(\lambda-\lambda^{\prime})}\right]^{2},\ \text{when}\ N\rightarrow+\infty,N\lambda,N\lambda^{\prime}\rightarrow\text{Const}

(30)

Plugging everything into the expression we have:

\mu_{2n}\approx 4\pi^{-2}\int d\lambda d\lambda^{\prime}\cdot(\lambda-\lambda^{\prime})^{2n}\left\{\sqrt{1-\lambda^{2}}\sqrt{1-\lambda^{\prime 2}}-\left[\frac{\sin 2N(\lambda-\lambda^{\prime})}{2N(\lambda-\lambda^{\prime})}\right]^{2}\right\}

(31)

We first omit the connected part in large $N$ limit. So the leading order in large $N$ would be:

$\displaystyle\mu_{2n}$	$\displaystyle\approx 4\pi^{-2}\int_{-1}^{1}d\lambda d\lambda^{\prime}\sqrt{1-\lambda^{2}}\sqrt{1-\lambda^{\prime 2}}(\lambda-\lambda^{\prime})^{2n}$	(32)
	$\displaystyle=4\pi^{-2}\int_{0}^{\pi}d\theta d\theta^{\prime}\sin^{2}\theta\sin^{2}\theta^{\prime}(\cos\theta-\cos\theta^{\prime})^{2n}$
	$\displaystyle=4\pi^{-2}\sum\limits_{k=0}^{2n}C_{2n}^{k}(-1)^{k}\left[\int_{0}^{\pi}d\theta\sin^{2}\theta\cos^{k}\theta\right]\left[\int_{0}^{\pi}d\theta\sin^{2}\theta\cos^{2n-k}\theta\right]$
	$\displaystyle=4\pi^{-2}\sum\limits_{k=0}^{2n}C_{2n}^{k}(-1)^{k}\left[F_{k}-F_{k+2}\right]\left[F_{2n-k}-F_{2n-k+2}\right]$

where the integral $F_{k}$ is elementary:

F_{k}=\int_{0}^{\pi}d\theta\cos^{k}\theta=\begin{cases}\frac{\Gamma\left(\frac{1}{2}\right)\Gamma\left(\frac{k}{2}+\frac{1}{2}\right)}{\Gamma\left(\frac{k}{2}+1\right)}&\text{for}\ k\in\text{even}\\ 0&\text{for}\ k\in\text{odd}\end{cases}

(33)

Plugging the expression of $F_{k}$ , then we arrived at the analytical form of $\mu_{2n}$ in leading order of large $N$ :

\mu_{2n}=2\pi^{-1}\sum\limits_{m=0}^{n}\frac{(2n)!}{(2m)!(2n-2m)!}\frac{\Gamma\left(m+\frac{1}{2}\right)\Gamma\left(n-m+\frac{1}{2}\right)}{\Gamma\left(m+2\right)\Gamma\left(n-m+2\right)}

(34)

It seems hopeless to obtain a closed expression for this summation, therefore we try to seek the asymptotic behavior of $\mu_{2n}$ when $n\gg 1$ . To do this, we appeal to the Stirling formula $\log\Gamma(z)\approx z\log z-z+O(\log z)$ :

$\displaystyle\mu_{2n}$	$\displaystyle\approx 2\pi^{-1}\sum\limits_{m=0}^{n}\exp\left\{-n\left[2x\log x+2(1-x)\log(1-x)\right]-3\log n-\frac{3}{2}\log x-\frac{3}{2}\log(1-x)\right\}$	(35)
	$\displaystyle\approx 2\pi^{-1}n^{-2}\int_{-1}^{1}dx\exp\left(-nf(x)+O(\log n)\right),\ \text{where}\ x\equiv\frac{m}{n},\ f(x)\equiv 2x\log x+2(1-x)\log(1-x)$
	$\displaystyle\approx 2\pi^{-1}n^{-2}\int_{-\infty}^{+\infty}dx\exp\left[-nf\left(1/2\right)-\frac{1}{2}nf^{\prime\prime}\left(1/2\right)(x-1/2)^{2}\right]$
	$\displaystyle\approx\exp\left[(2\log 2)n+O(\log n)\right]$

In the first line, we use Stirling formula; in the second we approximate summation by integral; in the third line we perform integral using saddle point approximation, namely expanding around minimal of $f(x)$ and perform Gaussian integral, and in the last line we keep contribution only up to order $O(n)$ on the exponential. This is because, in the first line of approximation, the precision is already only up to order $n$ , since if we want to keep to $O(\log n)$ , we should use a finer Stirling formula $\log\Gamma(z)\approx z\log z-z-\frac{1}{2}\log z+O(1)$

To summarize, we have an exponential scaling behavior of $\mu_{2n}$ :

\log\mu_{2n}\approx n\cdot 2\log 2+O(\log n)

(36)

Next, we notice that this exponential behavior of $\mu_{2n}$ would translate to the Constant behaviour of $b_{n}\equiv b$ for large $n$ . From equation (11), we see that if $b_{n}$ equals constant, we only need to count the number of paths:

\mu_{2n}=C_{n}b^{2n}=\frac{(2n)!}{(n+1)!(n)!}b^{2n}

(37)

where $C_{n}$ is the Catalan number which counts the number of paths Parker2019 . In order to compare the scaling behaviour of $\mu_{2n}$ for large but finite $n$ , we again appeals to Stirling formula:

\log\mu_{2n}\approx{n\cdot 2\log 2+2n\log b+O(\log n)}

(38)

\Longrightarrow b=1

(39)

This suggests that $b_{n}\approx 1$ will approach a constant value when $n$ is large!

We confirm this analytical prediction by direct numerical simulation of GUE random matrices, as presented in the left panel of figure 4. Numerically we see that $b_{n}$ is already close enough to 1 even when $n=10$ . Since there is no scale at all for RMT at large $N$ and infinite temperature, we may well represent $b_{n}=1$ for $n$ larger than some order one value. This is the reason why we draw the behavior of $b_{n}$ as a horizontal straight line in the schematic plot in figure 1(a).

As we have already pointed out in section 2, $b_{n}$ equating constant indicates that in large $N$ limit, $\mathcal{K}(t)$ in RMT grows linearly in time, in sharp contrast to the chaotic quantum system in thermodynamic limit, whose K-complexity grow exponentially in time, as mentioned in section 1. We remark that this does not mean K-complexity in RMT is smaller than that of chaotic local systems. This is because, in RMT theory, we usually scale the Hamiltonian such that its eigenenergy is of order $O(1)$ . But in chaotic local systems, the energy is an extensive value which is of order $O(S)$ , where $S$ is the number of degree of freedom. Since the energy scale is conjugate to the time scale, we should properly rescale our RMT result in order to compare with chaotic local systems.

We also remark that only the fact that it approaches constant is important, not the specific value of that constant. One can show directly from the Lanczos algorithm with infinite temperature inner product, that $b_{n}$ is proportional to the overall scaling of Hamiltonian. So, here deriving $b_{n}$ to approach 1 is simply a coincidence that we scale the radius semicircle to unity.

As a side remark, to make sure the linear growth of complexity on constant- $b$ -plateau, we need to ensure that the fluctuation of $b_{n}$ is small enough in order to free from Anderson localization on Krylov chain Rabinovici20221 ; Rabinovici20222 . Luckily, this is also satisfied in large $N$ limit, where in the right panel of figure 2 we show that the statistical variance of $\{b_{n}\}$ is suppressed by some positive power of $N^{-1}$ .

3.2 Numerically comparing with generic chaotic local systems.

As explained in section 1, one motivation for studying Krylov complexity in RMT is that for a generic chaotic Hamiltonian with locality, the complexity growth after scrambling time $t_{*}$ is governed by random matrix behavior. This is because, at $t_{*}$ , $O(t)$ has just finished its ‘scrambling in real space’, and begun to start its ‘scrambling in Hilbert space’. Therefore in stage two, $O(t)$ completely abandoned any spatial locality structure, or equivalently, the local Hamiltonian $H$ no longer looks local from the perspective of $O(t)$ . Therefore, in this stage, $H$ would act like a chaotic Hamiltonian without spatial locality, resembling a typical random matrix.

In this section, we numerically study the Lanczos coefficient $\{b_{n}\}$ of various initial operator $O$ in several chaotic local systems in various dimensions, including Ising model in one or two dimensions with both transverse field and longitudinal field, SYK model Maldecena2016 ; Sachdev2015 , one-dimensional spinless fermion model with nearest-neighbor and next-nearest-neighbor hopping and density-density interaction, and PXP model as an effective model describing blockated Rydberg atom array Bernien2017 . Their Hamiltonian is explicitly given by:

H_{\text{Ising}}=-\sum\limits_{\langle ij\rangle}\sigma_{i}^{x}\sigma_{j}^{x}+\sum\limits_{i}\left(g\sigma^{z}_{i}+h\sigma_{i}^{x}\right)

(40)

H_{\text{SYK}}=\sum\limits_{i<j<k<l}J_{ijkl}\chi_{i}\chi_{j}\chi_{k}\chi_{l},\ \overline{J^{2}_{ijkl}}=\frac{(4-1)!}{N_{m}^{4-1}}J^{2}

(41)

H_{\text{fermion}}=-\sum\limits_{i}\left(t_{1}c_{i+1}^{\dagger}c_{i}+t_{2}c_{i+2}^{\dagger}c_{i}+\text{h.c.}\right)+\sum\limits_{i}\left(V_{1}n_{i+1}n_{i}+V_{2}n_{i+2}n_{i}\right)

(42)

H_{\text{PXP}}=P\left[\sum\limits_{i}\left(\Omega\sigma_{i}^{x}+U\sigma_{i}^{z}\right)\right]P

(43)

In SYK model, we use $N_{m}$ to denote the number of Majoranas. In PXP model, $P$ is a projection operator on to the subspace that any nearest neighbor sites cannot be both spin up.

To compare with RMT result, we should properly rescale $H_{\text{RMT}}$ to $\eta S\cdot H_{\text{RMT}}$ such that $\eta S\sim O(S)$ . Therefore the plateau value of $b_{n}$ would rescale to $\eta S$ , which is the same order as $b_{n}$ from $H_{\text{chaotic}}$ . This rescaling also shifts the radius of Wigner semicircle from unity to $\eta S$ , the same order as the extensive energy from $H_{\text{chaotic}}$ . However, the order one coefficient $\eta$ is model-dependent and exhibits ambiguity. In order to obtain a model-independent rescaling ansatz, we define the rescaling factor to be such that the diameter of Wigner semicircle equals the width of the spectrum of $H_{\text{chaotic}}$ . Therefore, the plateau value of RMT is scaled to be $b=(E_{\max}-E_{\min})/2$ .

Numerical results in figure 3 shows that up to some erratic fluctuations, $b=(E_{\max}-E_{\min})/2$ is indeed slightly larger than $\{b_{n}\}$ from $H_{\text{chaotic}}$ , with SYK model almost saturate RMT result. Therefore, we may carefully conjecture that in any chaotic local quantum systems, the infinite temperature Lanczos coefficient $\{b_{n}\}$ is bounded by $b_{n}\leq\eta_{0}(E_{\max}-E_{\min})/2$ , with $\eta_{0}\sim O(1)$ a model-independent constant.

3.3 Finite temperature

For finite temperature, recall the inner product is defined by Wightman function:

(A|B)=\operatorname{tr}\left[e^{-\beta H/2}A^{\dagger}e^{-\beta H/2}B\right]

(44)

Similar calculation gives the expression of $\mu_{2n}$ :

\mu_{2n}=N^{-2}\sum\limits_{i,j}(\lambda_{i}-\lambda_{j})^{2n}e^{-\beta(\lambda_{i}+\lambda_{j})/2}

(45)

Using the same RMT technique in large $N$ limit:

	$\displaystyle\mu_{2n}$	$\displaystyle=4\pi^{-2}\int_{-1}^{1}d\lambda d\lambda^{\prime}\sqrt{1-\lambda^{2}}\sqrt{1-\lambda^{\prime 2}}(\lambda-\lambda^{\prime})^{2n}e^{-\beta(\lambda+\lambda^{\prime})/2}$		(46)
		$\displaystyle=4\pi^{-2}\sum\limits_{k=0}^{2n}C_{2n}^{k}(-1)^{k}\left[F_{k}-F_{k+2}\right]\left[F_{2n-k}-F_{2n-k+2}\right]$		(46)

where the integral $F_{k}$ is given by the derivative of the modified Bessel function of the first kind:

F_{k}=\int_{0}^{\pi}d\theta\cos^{k}\theta e^{-\beta\cos\theta/2}=\pi(-1)^{k}I_{0}^{(k)}(\beta/2)

(47)

with $I_{0}(z)=\pi^{-1}\int d\theta e^{-z\cos\theta}$ is the zeroth order of modified Bessel function of first kind, and $I_{0}^{(k)}(z)$ denotes its $k$ -th derivative. Using the property that $I_{0}(z)=I_{0}(-z)$ and the identity $\partial^{n}(fg)=\sum C_{n}^{k}f^{(k)}g^{(n-k)}$ , we can perform summation over $k$ and arrive at a simplified result:

\mu_{2n}=4\partial_{x}^{2n}\left[\left(I_{0}(x)-I_{0}^{(2)}(x)\right)^{2}\right]\bigg{|}_{x=\beta/2}

(48)

Of course one can use the expansion $I_{0}(z)=\sum\frac{1}{(k!)^{2}}(z/2)^{2k}$ to study the behaviour of $\nu_{2n}$ in general temperature, but to get a feeling, we can study low temperature region when $\beta$ is large. We use the asymptotic behaviour of Bessel function:

I_{0}(x)\sim\frac{e^{x}}{\sqrt{2\pi x}},\ x\gg 1

(49)

So we have $(I_{0}(x)-I_{0}^{(2)}(x))^{2}\propto e^{2x}x^{-3}\sim e^{2x}$ , where we ignored the $\log\beta$ correction on the exponential. Therefore:

\mu_{2n}\sim\partial_{x}^{2n}\left[e^{2x}\right]\bigg{|}_{x=\beta/2}=2^{2n}e^{\beta}\Longrightarrow\mu_{2n}\sim\exp(n\cdot 2\log 2)

(50)

Interestingly, under this simple approximation, one may observe that finite temperature seemingly doesn’t affect the scaling behaviour of $n$ , compared to infinite temperature. This also means $b_{n}=1=\text{Const}$ , following the same argument starting from equation (36) to (39).

Although the above calculation is in low temperature limit, we check numerically the exponential scaling behaviour of $\mu_{2n}$ actually appears for finite $\beta$ and large $n$ , namely we first take $n$ goes to infinity while keeping $\beta$ finite, see figure 5.

We also numerically calculate $b_{n}$ by direct simulation of RMT. The result is shown in figure 5. We also confirm that the $b_{n}=1$ region appears at large $n$ and finite $\beta$ . For finite $n$ but large $\beta$ , there is a region where $b_{n}$ linearly increases with $n$ .

In order to explain this linear growing region of $b_{n}$ , we directly consider the integral of $\mu_{2n}$ :

\mu_{2n}=4\pi^{-2}\int_{0}^{\pi}d\theta d\theta^{\prime}\cdot\sin^{2}\theta\sin^{2}\theta^{\prime}(\cos\theta-\cos\theta^{\prime})^{2n}e^{-\beta(\cos\theta+\cos\theta^{\prime})/2}

(51)

We want to evaluate this integral by taking low temperature ( $\beta\rightarrow\infty$ ). In this way we change dummy variable $\varphi=\pi-\theta,\varphi^{\prime}=\pi-\theta^{\prime}$ and expand the integrand for small $\varphi,\varphi^{\prime}$ :

\mu_{2n}\approx 2\pi^{-2}e^{\beta}2^{-2n}\int_{-\pi}^{\pi}d\varphi d\varphi^{\prime}(\varphi^{2}-\varphi^{\prime 2})^{2n}e^{-\beta(\varphi^{2}+\varphi^{\prime 2})/4}

(52)

where we also extend the integral range from $[0,\pi]$ to $[-\pi,\pi]$ using inversion symmetry of integrand. For low temperature, this Gaussian distribution is located at zero with width $O(\beta^{-1/2})$ , parametrically small compared to integral range $[-\pi,\pi]$ . Therefore we can safely extend the integral range to infinity. By changing variable to $\overline{\varphi}\equiv\varphi+\varphi^{\prime},\delta\varphi=\varphi-\varphi^{\prime}$ , this Gaussian integral is evaluated as:

	$\displaystyle\mu_{2n}$	$\displaystyle\approx\pi^{-1}e^{\beta}2^{-2n+2}\beta^{-1}\left(\left\langle\overline{\varphi}^{2n+4}\right\rangle_{\sigma}\left\langle\delta\varphi^{2n}\right\rangle_{\sigma}+\left\langle\overline{\varphi}^{2n}\right\rangle_{\sigma}\left\langle\delta\varphi^{2n+4}\right\rangle_{\sigma}-2\left\langle\overline{\varphi}^{2n+2}\right\rangle_{\sigma}\left\langle\delta\varphi^{2n+2}\right\rangle_{\sigma}\right)$		(53)
		$\displaystyle=\pi^{-1}\beta^{-2n-3}2^{3}e^{\beta}\left(\frac{(2n+2)!}{(n+1)!}\right)^{2}\frac{1}{2n+1}$		(53)

where $\langle\cdot\rangle_{\sigma}$ is the Gaussian average with variance $\sigma=4\beta^{-1}$ . Going from the first line to the second line, we use $\langle z^{2n}\rangle_{\sigma}=\sigma^{n}\frac{(2n)!}{2^{n}n!}$ .

Using Stirling formula for large but finite $n$ , we find that the scaling behaviour of $\log\mu_{2n}$ with respect to $n$ :

\log\mu_{2n}\approx 2n\log n+(4\log 2-2+2\log\beta^{-1})n+O(\log n)

(54)

Assuming ansatz $b_{n}=\alpha n^{\delta}$ and using inequality (14), we can fix $\delta=1$ and provide a bound for $\alpha$ :

4\beta^{-1}>\alpha>2\beta^{-1}

(55)

This is indeed consistent with numerical simulation in the right panel of figure 5, where we found $b_{n}=\pi\beta^{-1}n$ , with $\alpha=\pi\beta^{-1}$ , which is indeed within the bound.

We notice that from the continuous limit of the operator wave function on Krylov chain, as analyzed in section 2, we immediately conclude that the Krylov complexity $\mathcal{K}(t)$ at low temperature will exhibit the following two-stage growth:

1.

Exponential growth before scrambling time. $\mathcal{K}(t)\sim e^{\frac{2\pi}{\beta}t}$ . We also notice that $\mathcal{K}(t)$ is interpreted as the average position of the wave function on the semi-infinite Krylov chain, therefore this exponential growth will proceed until the average position reaches the plateau, where the linear growth of $b_{n}$ ceases, namely when $\mathcal{K}(t_{*})\sim O(\beta)$ . This shows the scrambling time scale $t_{*}$ is of order $O(\beta\log\beta)$ .
2.

Linear growth after scrambling time. $\mathcal{K}(t)\sim t$ . When $t\gtrsim t_{*}$ , the average position of the wave packet will travel to the region where $b_{n}$ is constant. Therefore it becomes a usual wave packet traveling ahead with constant velocity.

We also notice that at low temperature and before scrambling time, the exponential growth of Krylov complexity saturates the chaos bound of Lyaponov exponent $\lambda_{L}=\frac{2\pi}{\beta}$ Maldacena2016bound .

On the one hand, this is natural since the original literature on Krylov complexity Parker2019 has proved that K-complexity can bound many other complexity measures like OTOC and operator size at infinite temperature, and for finite temperature the authors there conjectured the same bound to be true and verified on SYK model.

On the other hand, the initial linear grow of $b_{n}$ in RMT should come as a surprise, since the linear growing behavior of $b_{n}$ , though generic in chaotic many body systems, actually originates from the locality of Hamiltonian Parker2019 . But clearly a random matrix Hamiltonian lacks locality. One interpretation is that the spectrum of RMT at low energy, which scales as $\rho(E)\propto\sqrt{E}$ agrees with the spectrum of a quantum black hole in JT gravity (see reference Mertens2023 for a recent review) or SYK model Maldecena2016 in low energy, which scales as $\rho_{JT}(E)\propto\sinh(2\pi\sqrt{2CE})\sim\sqrt{E}$ Mertens2023 . Matching of low energy density of state suggests that RMT has similar behavior with black hole in the view of fast scrambler YasuhiroSekino2008 .

To justify that this low temperature chaos bound is universal, one necessary condition is that the coefficient of exponential growth in time, equivalently the slope of $b_{n}$ wrt to $n$ , should not change when we vary the overall scaling of Hamiltonian while keeping temperature to be fixed. We notice that this is not possible in infinite temperature case. We recall the fact that $\{b_{n}\}$ is simply the matrix elements of Liouville superoperator $\mathcal{L}=[H,\cdot]$ writing in an orthonormal Krylov basis $\{|O_{n})\}$ . Therefore, naively we should expect $\{b_{n}\}$ to be proportional to the overall scaling of $H$ . This is indeed true in infinite temperature, since the inner product $(A|B)=\operatorname{tr}[A^{\dagger}B]$ used to normalize the basis $\{|O_{n})\}$ do not depend on $H$ , therefore one can show by induction that $\{|O_{n})\}$ do not depend on the overall scaling of $H$ . However, this is not the case for finite temperature, where we choose Wightman function to define the inner product $(A|B)_{\beta}=\operatorname{tr}[e^{-\beta H/2}A^{\dagger}e^{-\beta H/2}B]$ , which explicitly depend on $H$ . As a result the basis $\{|O_{n})\}$ depends on the overall scaling of $H$ in a complicated way and $\{b_{n}\}$ no longer proportional to it.

From right panel of figure 4 we find that for finite temperature, the constant plateau part of $\{b_{n}\}$ is indeed affected by overall scaling, resembling infinite temperature behavior, while the initial slope remains intact from overall scaling! This perfectly justified the fact that the saturation of RMT to chaos bound is not a coincidence.

One interesting thing is that for RMT in $N\rightarrow\infty$ , we obtain a finite scrambling time $t_{*}^{\text{RMT}}\sim O(\beta\log\beta)$ , which is in contrast with chaotic local quantum systems $t_{*}^{\text{chaotic}}\sim O(\beta\log S)$ , diverging in thermodynamic limit $S\rightarrow\infty$ . Similar to the situation discussed in section 3.2, the reason is merely the overall scaling of Hamiltonian. From the right panel of figure 5, we can of course prolong the scrambling time by scaling $H_{\text{RMT}}\rightarrow O(S)\cdot H_{\text{RMT}}$ , in order to match the scrambling time. Since $S\sim\log N$ , this rescaling is equivalent to changing the variance of RMT probability distribution: $\log P(H)\sim-2N\operatorname{tr}H^{2}\longrightarrow\log P(H)\sim-\frac{2N}{\log^{2}N}\operatorname{tr}H^{2}$ , which is strange.

In RMT with scale controlled by probability distribution $P(H)\propto\exp(-2N\operatorname{tr}H^{2})$ , the energy spectrum is of order $O(1)$ , which remains finite when $N\rightarrow\infty$ . For chaotic local quantum systems, it energy spectrum diverges since energy is extensive.

4 Remarks on finite $N$ corrections

A general remark would be that the level repulsion in connected part would have more impact for $\mu_{2n}$ with larger $n$ . This is because level repulsion would make $\lambda_{i}-\lambda_{j}$ larger. We also notice that, as $n$ gets larger, $x^{2n}$ have increasingly more weight on large- $x$ region.

To observe this, we first collect the exact result of $R_{1}(\lambda)$ and $T_{2}(\lambda,\lambda^{\prime})$ without any approximation on $N$ :

	$\displaystyle R_{1}(\lambda)=\sqrt{2N}^{1}K_{N}\left(\sqrt{2N}\lambda,\sqrt{2N}\lambda\right)$		(56)
	$\displaystyle T_{2}(\lambda,\lambda^{\prime})=\sqrt{2N}^{2}\left[K_{N}\left(\sqrt{2N}\lambda,\sqrt{2N}\lambda^{\prime}\right)\right]^{2}$		(57)

with the kernel function given in terms of normalized wave-function of harmonic oscillator $\varphi_{j}(x)=(2^{j}j!\sqrt{\pi})^{-1/2}\exp(-x^{2}/2)H_{j}(x)$ :

K_{N}(x,y)=\sum\limits_{j=0}^{N-1}\varphi_{j}(x)\varphi_{j}(y)

(58)

Then from the exact expression of $\mu_{2n}$ :

\mu_{2n}=N^{-2}\int d\lambda a\lambda^{\prime}(\lambda^{2})(\lambda-\lambda^{\prime})^{2n}\left[R_{1}(\lambda)R_{1}(\lambda^{\prime})-T_{2}(\lambda,\lambda^{\prime})\right],

(59)

we obtain some numerical results shown in figure 6. We first indeed see that when $N$ becomes large, the result approaches infinite $N$ limit. In the right panel of figure 6, we see that even for finite $N$ , the large $n$ behavior of $\mu_{2n}$ is still exponential in $n$ . From equation (11), a simple interpretation is that $\mu_{2n}$ would be exponential in the largest eigenvalue of $L$ .

We can also infer information of ${b_{n}}$ from the exponential growth. Notice that this does not contradict with the fact that $b_{n}$ in large but finite $n$ , namely $n\leq N^{2}$ would become zero, due to the termination of Krylov subspace. For example, for $N=20$ ’s case, we necessarily have $b_{400}=0$ , and from equation (12) we only need information of $\mu_{2n}$ at most up to $\mu_{400}$ , so all the long exponential growth( $\mu_{2n>400}$ ) in figure 6 actually do not provide independent information, they only provide information on $b_{\max}$ :

\max\limits_{0\leq n\leq N^{2}}[b_{n}]\sim\lim\limits_{n\rightarrow+\infty}\left[\frac{\log(\mu_{2n})}{2n}\right]

(60)

So, the decrease and termination Rabinovici2021 of $b_{N^{2}}=0$ actually hide in the subleading term on the exponential.

The argument of equation (60) is in the following. For large $n$ , $\mu_{2n}$ ’s dependence on ${b_{n}}$ is actually like a path integral:

	$\displaystyle\mu_{2n}=\sum\limits_{\text{path}}e^{S(\text{path})},\ \text{path}=(i_{1},i_{2},...,i_{2n}),\ \text{with}\ i_{1}=i_{2n}=0$		(61)
	$\displaystyle S(\text{path})=\sum\limits_{k=1}^{2n}\log b_{i_{k}}\approx 2n\oint dx\ \log(b(x))=2n\int_{0}^{1}dt\ \log(b(x(t))),\ \text{where}\ t=k/2n$		(62)
	$\displaystyle\Longrightarrow e^{\log(\mu_{2n})}\approx\int\mathcal{D}[x(t)]\ \exp\left[2n\int_{0}^{1}dt\log(b(x(t)))\right]\bigg{\|}_{x(0)=x(1)=0,\ \dot{x}(t)<\epsilon}$		(63)

For large $n$ , this path integral is dominated by its saddle point, i.e., the particular path with the largest action. For finite $N$ , $b_{n}$ would first grow and then decrease to zero Parker2019 , meaning that there is necessarily a $b_{\text{max}}$ , therefore, the saddle path would be ‘wandering around the site with $b_{\text{max}}$ ’, which gives the argument of equation (60).

In the limit $N\rightarrow\infty$ , the decreasing slope of $b_{n}$ would go to zero as $N^{-2}$ Parker2019 , then there are many paths with near the same action, and the saddle point approximation may not be valid.

In conclusion, although $b_{n}$ and $\mu_{2n}$ have an exact one-to-one mapping rule, this mapping is complicated. The large-scale structure of $b_{n}$ is hidden in the subleading terms and early- $n$ -behaviour of $\mu_{2n}$ , which is relatively hard to extract theoretically.

Averaged Krylov complexity and Spectrum Form Factor

As an aside, an interesting observation is that $\mu_{2n}$ in RMT is actually the moment of the superoperator. At infinite temperature, from equation (21) we observe:

\mu_{2n}=N^{-2}\sum\limits_{i,j}(\lambda_{i}-\lambda_{j})^{2n}=N^{-2}\operatorname{tr}(\mathcal{L}^{2n})=N^{-2}(-1)^{n}\left(\frac{\partial}{\partial t}\right)^{2n}\operatorname{tr}(e^{i\mathcal{L}t})\bigg{|}_{t=0}

(64)

which means that it is related to the higher derivative of SSF (Spectrum Form Factor) or analytically continued partition function. For finite temperature, from equation (45) we obtain:

	$\displaystyle\mu_{2n}$	$\displaystyle=(-1)^{n}\left(\frac{\partial}{\partial t}\right)^{2n}\left[\sum e^{-\beta(\lambda_{i}+\lambda_{j})/2+it(\lambda_{i}-\lambda_{j})}\right]\bigg{\|}_{t=0}$		(65)
		$\displaystyle=(-1)^{n}\left(\frac{\partial}{\partial t}\right)^{2n}\big{\|}Z(\beta/2,t)\big{\|}^{2}\bigg{\|}_{t=0}$		(65)

This observation also shows that the following quantity contains the same spectral information in random matrix theory:

\mathcal{K}(t)\longleftrightarrow\{b_{n}\}\longleftrightarrow\{\mu_{2n}\}\longleftrightarrow Z(t)

(66)

Therefore, we may expect that the decay-dip-ramp-plateau structure Cotler2017 in SSF $(t)$ may translate correspondingly to some behavior in time dependence of complexity $\mathcal{K}(t)$ . A possible obstacle is that the discontinuity of $Z(t)$ at dip-time and plateau-time, which is highly non-perturbative physics, would be hard to extract from its derivatives at $t=0$ , which is perturbative.

5 Conclusions

This study delves into the realm of Random Matrix Theory (RMT) to explore Krylov complexity $\mathcal{K}(t)$ . In the large $N$ limit at infinite temperature, we analytically demonstrate the saturation of the Lanczos coefficient ${b_{n}}$ to a constant plateau $\lim\limits_{n\rightarrow\infty}b_{n}=b$ , resulting in a linear increase in complexity $\mathcal{K}(t)\sim t$ . Comparison of this plateau value $b$ with a diverse set of chaotic local quantum systems suggests its bounding effect on ${b_{n}}$ in such systems. Hence, we conjecture that, after scrambling time, the linear growth rate of Krylov complexity in chaotic local systems cannot exceed that in RMT.

At low temperature, we analytically establish that $b_{n}$ initially undergoes linear growth with $n$ , reaching a plateau characterized by the renowned chaos bound. Upon reaching this common plateau $b$ , $b_{n}$ stabilizes, indicating $\mathcal{K}(t)\sim e^{2\pi t/\beta}$ before the scrambling time $t_{*}\sim O(\beta\log\beta)$ . Subsequently, a linear growth phase ensues, mirroring the behavior at infinite temperature.

We conclude by addressing the influence of finite $N$ corrections, and the relation between K-complexity and spectrum form factor. Our results confirmed and strengthened previous understanding of RMT as a fast scrambler in the view of complexity growth.

Acknowledgements.

I would thank Chang Liu and Hui Zhai for previous collaboration on related topics. I would also thank Yiming Chen, Yingfei Gu, Chang Liu, Xiao-Liang Qi, Hui Zhai, Pengfei Zhang, Tian-Gang Zhou, and Yi-Neng Zhou for helpful discussion. I would especially thank professor Hui Zhai, not only for leading me to the topic of Krylov complexity but also for his encouragement and support throughout this project. The majority of this work was finished when I was studying at Institute for Adanced Study, Tsinghua University, and I would like to thank for their hospitality.

References

(1) D.A. Roberts and B. Yoshida, Chaos and complexity by design, Journal of High Energy Physics 2017 (2017) 121.
(2) J. Cotler, N. Hunter-Jones, J. Liu and B. Yoshida, Chaos, complexity, and random matrices, Journal of High Energy Physics 2017 (2017) 48.
(3) R.A. Jefferson and R.C. Myers, Circuit complexity in quantum field theory, Journal of High Energy Physics 2017 (2017) 107.
(4) D.A. Roberts, D. Stanford and A. Streicher, Operator growth in the syk model, Journal of High Energy Physics 2018 (2018) 122.
(5) R.-Q. Yang, Complexity for quantum field theory states and applications to thermofield double states, Phys. Rev. D 97 (2018) 066004.
(6) R. Khan, C. Krishnan and S. Sharma, Circuit complexity in fermionic field theory, Phys. Rev. D 98 (2018) 126001.
(7) X.-L. Qi and A. Streicher, Quantum epidemiology: operator growth, thermal effects, and syk, Journal of High Energy Physics 2019 (2019) 12.
(8) A. Lucas, Operator size at finite temperature and planckian bounds on quantum dynamics, Phys. Rev. Lett. 122 (2019) 216601.
(9) V. Balasubramanian, M. DeCross, A. Kar and O. Parrikar, Quantum complexity of time evolution with chaotic hamiltonians, Journal of High Energy Physics 2020 (2020) 134.
(10) V. Balasubramanian, M. DeCross, A. Kar, Y.C. Li and O. Parrikar, Complexity growth in integrable and chaotic models, Journal of High Energy Physics 2021 (2021) 11.
(11) J. Maldacena, S.H. Shenker and D. Stanford, A bound on chaos, Journal of High Energy Physics 2016 (2016) 106.
(12) D.E. Parker, X. Cao, A. Avdoshkin, T. Scaffidi and E. Altman, A universal operator growth hypothesis, Phys. Rev. X 9 (2019) 041017.
(13) J. Maldacena and D. Stanford, Remarks on the sachdev-ye-kitaev model, Phys. Rev. D 94 (2016) 106002.
(14) E. Rabinovici, A. Sánchez-Garrido, R. Shir and J. Sonner, Operator complexity: a journey to the edge of krylov space, Journal of High Energy Physics 2021 (2021) 62.
(15) L. Susskind, Computational complexity and black hole horizons, 1402.5674.
(16) A. Bhattacharyya, S.S. Haque, G. Jafari, J. Murugan and D. Rapotu, Krylov complexity and spectral form factor for noisy random matrix models, Journal of High Energy Physics 2023 (2023) 157.
(17) J. Erdmenger, S.-K. Jian and Z.-Y. Xian, Universal chaotic dynamics from krylov space, Journal of High Energy Physics 2023 (2023) .
(18) A. Kar, L. Lamprou, M. Rozali and J. Sully, Random matrix theory for complexity growth and black hole interiors, Journal of High Energy Physics 2022 (2022) 16.
(19) S.-K. Jian, B. Swingle and Z.-Y. Xian, Complexity growth of operators in the syk model and in jt gravity, Journal of High Energy Physics 2021 (2021) 14.
(20) J.L.F. Barbón, E. Rabinovici, R. Shir and R. Sinha, On the evolution of operator complexity beyond scrambling, Journal of High Energy Physics 2019 (2019) 264.
(21) A. Dymarsky and A. Gorsky, Quantum chaos as delocalization in krylov space, Phys. Rev. B 102 (2020) 085137.
(22) T. Xu, T. Scaffidi and X. Cao, Does scrambling equal chaos?, Phys. Rev. Lett. 124 (2020) 140602.
(23) C.-F. Chen and A. Lucas, Operator growth bounds from graph theory, Communications in Mathematical Physics 385 (2021) 1273.
(24) J.D. Noh, Operator growth in the transverse-field ising spin chain with integrability-breaking longitudinal field, Phys. Rev. E 104 (2021) 034112.
(25) F.B. Trigueros and C.-J. Lin, Krylov complexity of many-body localization: Operator localization in Krylov basis, SciPost Phys. 13 (2022) 037.
(26) J. Kim, J. Murugan, J. Olle and D. Rosa, Operator delocalization in quantum networks, Phys. Rev. A 105 (2022) L010201.
(27) N. Hörnedal, N. Carabba, A.S. Matsoukas-Roubeas and A. del Campo, Ultimate speed limits to the growth of operator complexity, Communications Physics 5 (2022) 207.
(28) B. Bhattacharjee, X. Cao, P. Nandy and T. Pathak, Krylov complexity in saddle-dominated scrambling, Journal of High Energy Physics 2022 (2022) 174.
(29) P. Caputa and S. Liu, Quantum complexity and topological phases of matter, Phys. Rev. B 106 (2022) 195125.
(30) Z.-Y. Fan, Universal relation for operator complexity, Phys. Rev. A 105 (2022) 062210.
(31) Z.-Y. Fan, The growth of operator entropy in operator growth, Journal of High Energy Physics 2022 (2022) 232.
(32) E. Rabinovici, A. Sánchez-Garrido, R. Shir and J. Sonner, Krylov localization and suppression of complexity, Journal of High Energy Physics 2022 (2022) 211.
(33) E. Rabinovici, A. Sánchez-Garrido, R. Shir and J. Sonner, Krylov complexity from integrability to chaos, Journal of High Energy Physics 2022 (2022) 151.
(34) C. Liu, H. Tang and H. Zhai, Krylov complexity in open quantum systems, Phys. Rev. Res. 5 (2023) 033085.
(35) A. Bhattacharya, P. Nandy, P.P. Nath and H. Sahu, Operator growth and krylov construction in dissipative open quantum systems, Journal of High Energy Physics 2022 (2022) 81.
(36) A. Dymarsky and M. Smolkin, Krylov complexity in conformal field theory, Phys. Rev. D 104 (2021) L081702.
(37) P. Caputa and S. Datta, Operator growth in 2d cft, Journal of High Energy Physics 2021 (2021) 188.
(38) P. Caputa, J.M. Magan and D. Patramanis, Geometry of krylov complexity, 2109.03824.
(39) K. Adhikari and S. Choudhury, C ${\cal c}$ osmological k ${\cal k}$ rylov c ${\cal c}$ omplexity, Fortschritte der Physik 70 (2022) .
(40) K. Adhikari, S. Choudhury and A. Roy, Krylov complexity in quantum field theory, Nuclear Physics B 993 (2023) 116263.
(41) A. Avdoshkin and A. Dymarsky, Euclidean operator growth and quantum chaos, Physical Review Research 2 (2020) .
(42) J.L. Barbón, J. Martín-García and M. Sasieta, Momentum/complexity duality and the black hole interior, Journal of High Energy Physics 2020 (2020) 169.
(43) J.M. Magán and J. Simón, On operator growth and emergent poincaré symmetries, Journal of High Energy Physics 2020 (2020) 71.
(44) W. Mück and Y. Yang, Krylov complexity and orthogonal polynomials, Nuclear Physics B 984 (2022) 115948.
(45) A. Banerjee, A. Bhattacharyya, P. Drashni and S. Pawar, From cfts to theories with bondi-metzner-sachs symmetries: Complexity and out-of-time-ordered correlators, Phys. Rev. D 106 (2022) 126022.
(46) S. Sachdev, Bekenstein-hawking entropy and strange metals, Phys. Rev. X 5 (2015) 041025.
(47) D.L. Jafferis, D.K. Kolchmeyer, B. Mukhametzhanov and J. Sonner, Matrix models for eigenstate thermalization, Phys. Rev. X 13 (2023) 031033.
(48) D.L. Jafferis, D.K. Kolchmeyer, B. Mukhametzhanov and J. Sonner, Jackiw-teitelboim gravity with matter, generalized eigenstate thermalization hypothesis, and random matrices, Phys. Rev. D 108 (2023) 066015.
(49) M. Mehta, Random Matrices, ISSN, Elsevier Science (2004).
(50) H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler et al., Probing many-body dynamics on a 51-atom quantum simulator, Nature 551 (2017) 579.
(51) T.G. Mertens and G.J. Turiaci, Solvable models of quantum black holes: a review on jackiw–teitelboim gravity, Living Reviews in Relativity 26 (2023) 4.
(52) Y. Sekino and L. Susskind, Fast scramblers, Journal of High Energy Physics 2008 (2008) 065.
(53) J.S. Cotler, G. Gur-Ari, M. Hanada, J. Polchinski, P. Saad, S.H. Shenker et al., Black holes and random matrices, Journal of High Energy Physics 2017 (2017) 118.

$\displaystyle\|A_{n})$	$\displaystyle=\mathcal{L}\|O_{n-1})-b_{n-1}\|O_{n-2})$	(4)
$\displaystyle b_{n}$	$\displaystyle=(A_{n}\|A_{n})^{1/2}$
$\displaystyle\|O_{n})$	$\displaystyle=b_{n}^{-1}\|A_{n})$