Universality and Emergent Classicality in Krylov Space

Introduction— Universality is a central organizing principle in modern physics. It refers to the phenomenon where a class of macroscopic systems with drastically different microscopic details obey the same scaling laws for coarse-grained observables. Renormalization-group fixed points provide some of its most striking examples. Another example is the level-spacing statistics of generic nonintegrable quantum systems, described by random-matrix ensembles. In this Letter, we present a notion of universality for late-time dynamics in Krylov space that applies to a broad class of physical systems. It describes the ballistic motion of a quantum state in a distinguished orthogonal basis. This motion is determined by the energy distribution of the state through a universal relation. We call this basis the Krylov basis following standard terminology. The Krylov basis relevant for describing time-evolution $|\psi(t)\rangle:= e^{-iHt}|\psi(0)\rangle$ is constructed through the recursion

with boundary conditions $b_{-1}=a_{-1}=0$ and $|0\rangle:= |\psi(0)\rangle$, such that the basis states $\{|n\rangle\}$ are orthonormal. The Hamiltonian $H$ acts like a tridiagonal matrix in this basis, i.e., $M_{mn}:=\langle m|H|n\rangle=b_{n-1}\delta_{m,n-1}+a_n\delta_{m,n}+b_n\delta_{m,n+1}$. The Lanczos coefficients $b_n$, $a_n$ are fixed by orthogonality, in terms of $\langle\psi|H^k|\psi\rangle$ for $k$ up to $2n+2$. The Krylov basis is entirely determined by the survival amplitude

or equivalently, its Fourier transform

which we call the spectral function. In this basis the state is expanded as

where $\varphi_n(t)$ are determined in terms of $C(t)$. The moments of the distribution $\{|\varphi_n(t)|^2\}$ are given by powers of the “complexity” operator $\hat n:=\sum_n n|n\rangle\langle n|$ in the state $|\psi(t)\rangle$:

The Krylov basis has emerged as a framework for characterizing quantum dynamics [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. This activity was stimulated by [12] which argued that operator wavefunctions spread universally in the Krylov basis under Heisenberg evolution, and that $\langle\hat n(t)\rangle$ grows exponentially for local operators in generic nonintegrable systems. The mean $\langle\hat n(t)\rangle$, called Krylov complexity is a simple measure of operator growth obtained from the operator autocorrelation function. In general, Krylov complexity shows rich phenomenology across several time scales [13, 14].

Statement— For a broad class of systems where $\Phi(E)$ has bounded support, $\langle\hat n(t)\rangle$ is expected to exhibit a parametrically long regime of linear growth [15, 16, 17, 18]. Extended linear growth is also a central feature of holographic proposals relating quantum complexity to the growth of Einstein-Rosen bridges [19, 20]. In finite systems, a maximum energy directly leads to $\Phi(E)$ vanishing beyond a cutoff. In sufficiently large systems, this regime may persist for times exponential in system size, before a slow drift sets in and eventually leads to saturation. For infinite systems with bounded spectral support, the linear growth may persist indefinitely.

We shift the energies by a constant so they lie in $[-E_0,E_0]$. We assume scenarios where $\Phi(E)$ is exactly, or after coarse-graining, smooth and continuous, vanishing for $|E|>|E_0|$ with no gaps inside $[-E_0,E_0]$. The generating function of $\langle\hat n^k(t)\rangle$ was studied in [21] which found that at late times

where the slopes $v_k$ have the spectral representation

which can also be rewritten in the time-domain

if $E_0$ is known, but only $C(t)$ is directly accessible. Ref. [21] derived (6) using the propagation generated by $\hat{n}$ in the energy basis. The main technical ingredient in the derivation is the asymptotic form of $|n\rangle$ in the energy basis (22), that follows from the orthogonal polynomial asymptotics found by Szegő [22]. The $O(t^{k-1})$ bound on corrections requires further assumptions but is expected to be generic. The asymptotics (22) also offer a way to directly approximate $\varphi_n(t)$ as explored in Sec. SM A. While completing this paper, we became aware of independent derivations of (7) in recent works [23, 24], for general $k$ and for $k=1$, respectively.

A simpler physical argument is the following. In our conventions, $\lim_{n\to\infty} b_n=E_0/2$ and $\lim_{n\to\infty} a_n=0$. The Lanczos coefficients approach their asymptotic value at a scale $n_c$. If this approach is algebraic $\sim n^{-\gamma}$ for $\gamma>0$, $n_c(\epsilon)$ can be defined as $\left|b_{n_c(\epsilon)}-\frac{E_0}{2}\right|=\epsilon$. Sometimes $n_c$ can be estimated as the $n$ at which the saddle point of $\int dE\,E^{2n}\Phi(E)$ reaches $|E|=E_0$, which depends on the detailed profile of $\Phi(E)$.

The equations of motion for the Heisenberg picture $\hat{n}(t)$ are

expressed in terms of $B:=i[M,\hat{n}]$, $F:=-[M,[M,\hat{n}]]$. We denote the asymptotic form of operators at indices $n\gg n_c$ using the subscript $n_c$. $M_{n_c}$ is a constant symmetric tridiagonal matrix with off-diagonal entries $E_0/2$ and zero along the diagonal. $B_{n_c}$ is $i$ times the antisymmetrized $M_{n_c}$. This implies

where in writing (10) we neglect boundary terms and deviations of the Lanczos coefficients from their asymptotic values. Since the wavepacket moves towards larger $n$, we take the positive root

Once the wavepacket $\varphi_n(t)$ has propagated into this region $n\gg n_c$, (10), (11) can be used to approximate the solution of (9) within expectation values. Vanishing of $F_{n_c}$ implies the ballistic growth in (6). Using energy conservation of $|\psi(t)\rangle$, the asymptotic relation (11) gives rise to the spectral relation (7).

The onset time $t_{\mathrm{on}}^{(k)}$ for the validity of (6) depends on $k$, but this effect is small [21], so we restrict to the single scale $t_{\mathrm{on}}$. Subleading corrections in (6) encode details of the transient propagation and the approach of the Lanczos coefficients to their asymptotic values. The width of $\Phi(E)$ sets the scale for Lanczos coefficients; if this scale does not depend on the system size but $E_0$ does then $t_{\mathrm{on}}$ grows with the system size.

Emergent classicality— When (6) is valid for sufficiently large times $\langle\hat{n}^k(t)\rangle/(v_k t^k)$ approaches 1 and the Krylov distribution is simply and universally characterized by the spectral representation of $|\psi\rangle$ through (7). In this regime, there is a simple relation between the width of $|\psi(t)\rangle$ in the Krylov and energy bases. In many contexts, the width of $\Phi(E)$ is subextensive but $E_0$ is extensive. If $\Phi(E)$ is sharply peaked about its mean $\mu_1$, expanding $(E_0^2-E^2)^{k/2}$ about $\mu_1$ gives $v_k\sim(E_0^2-\mu_1^2)^{k/2}$ with corrections that are suppressed by the system size. $|\psi(t)\rangle$ becomes effectively classical in the Krylov basis in the following sense. $\langle\hat{n}^k(t)\rangle\approx\langle\hat{n}(t)\rangle^k$, which means that $\hat{n}$ becomes a classical macroscopic variable. Smooth functions of $\hat{n}$ become deterministic, evaluated at late times. The Krylov probabilities move along a ballistic trajectory determined by (6).

Surprisingly, the phases in $\varphi_n(t)$ also become deterministic. Despite the obstruction to defining a global phase operator conjugate to $\hat{n}$, a simple workaround is to use the shift operator $T=\sum_n|n+1\rangle\langle n|$ [25]. The identification $T=e^{ip}$ is valid when $|\psi(t)\rangle$ has negligible support close to $n=0$. $\arg\langle T(t)\rangle$ can be interpreted as the mean $\langle p\rangle$. Because $M_{n_c}+B_{n_c}\approx E_0T$, the $p$ distribution is approximately $t$-independent and

with analogous diagnostics for fluctuations of $p$. When $\Phi(E)$ is sharply peaked, so is the $p$ distribution and $\langle p\rangle$ can be approximated by $\arccos(\mu_1/E_0)$, see Figure 3. This, together with $|\langle T(t)\rangle|\approx1$ implies an approximately linear winding phase for $\varphi_n(t)$ and can be seen as one possible mechanism underlying Krylov winding studied in [26].

Thus, at late times $|\psi(t)\rangle$ describes a well-defined trajectory $(\sqrt{E_0^2-\mu_1^2}t,\arccos(\mu_1/E_0))$ for the semiclassical variables $n$ and $p$. This is similar to a coherent state with large mean occupation number and a well-defined phase. The number-phase uncertainty is respected because variance of $\hat{n}$ is not small in absolute terms. In fact, the variance of $\hat{n}$ grows as $t^2$ but it is the relative width

that converges to $(\sqrt{v_2-v_1^2})/v_1$. The approach of $\delta n(t)$ to its $t\to\infty$ asymptotic defines the scale $t_{\mathrm{cl}}$, marking the point where the subleading terms in (6) become negligible for $k=1,2$.

There are physical scenarios where the sharpness of Krylov distribution does not hold. Take $|\psi\rangle$ uniformly populated within an energy window of width $2E_0$. $E_0$ could be taken finite to mimic a microcanonical state in a large system. Then (6) holds for exponentially long times $t$ with $O(1)$ values of $v_k$ and $\delta n(t)$ asymptotes to $\approx0.284$. Similarly $\lim_{t\to\infty}\delta n(t)\approx0.202$ for the dynamics of a fully delocalized state under a Gaussian random matrix of large size assuming that coarse-grained $\Phi(E)$ converges to the Wigner semicircle distribution. The narrowness condition is also not sufficient if $\Phi(E)$ is sharply peaked very close to $E_0$.

Conversely, a narrow $\Phi(E)$ away from the edges $\pm E_0$ is sufficient but not necessary for the wavepacket to become sharply peaked in the Krylov basis. $\delta n(t)$ decays as $\frac{1}{\sigma t}$ when $\Phi(E) \propto e^{-\frac{E^2}{2\sigma^2}}$, while $\langle\hat{n}(t)\rangle=\sigma^2t^2$.

Operator growth— Here we estimate $t_{\mathrm{on}}$ and $t_{\mathrm{cl}}$ using a physically motivated model for $\Phi(\omega)$ relevant for operator growth. Following [12], we assume that $\Phi(\omega)=\Phi(-\omega)$ decays as $\sim e^{-\frac{\pi}{2\alpha}|\omega|}$ in the bulk, and that it vanishes at an extensive maximal frequency $\omega_0$. Accordingly, our analysis relies exclusively on $\omega_0$ and $\alpha$. We assume that $\Phi(\omega)$ can be treated as a smooth function and ignore additional structure such as its precise behavior near $\omega_0$. For symmetrically ordered $C(t)$ at finite inverse temperature $\beta$, it is expected that $\alpha\sim\frac{\pi}{\beta}$. To leading order in $\alpha/\omega_0\ll1$, $v_2-v_1^2\sim\alpha^4/\omega_0^2$ and $v_1\sim\omega_0$.

Matching $b_n\approx\alpha n$ for $n\lesssim n_c$ and $b_n=\omega_0/2$ for $n\ge n_c$, we get $n_c\approx\omega_0/2\alpha$. Using the known transient exponential growth $e^{2\alpha t_{\mathrm{on}}}\approx n_c$ we get

To find $t_{\mathrm{cl}}$, we need control over the subleading terms in (6). In the End Matter we describe a phenomenological model designed with input from numerical experiments that gives

Matching (15) with its asymptotic limit $\sim\alpha^2/\omega_0^2$ gives $t_{\mathrm{cl}}\sim\omega_0^2/\alpha^3$.

As a concrete physical example, we study the two-point function (44) of the transverse magnetic field in the critical spin-1/2 XY model (also studied in [9, 27]). This is a free system with $\Phi(\omega)\sim e^{-\beta|\omega|/2}$ and fixed $\omega_0=2$.

The saturation of $\delta n(t)$ becomes numerically accessible because $\omega_0$ is a fixed $O(1)$ number. It mimics features of generic systems as $\omega_0\gg\alpha=\frac{\pi}{\beta}$ up to an overall rescaling of time. Figure 1 illustrates $t_{\mathrm{on}}$ and $t_{\mathrm{cl}}$ consistent with our model in this system.

Double-Scaled SYK (DSSYK)— The SYK model [28, 29] of $N$ Majoranas with random $q$-body interactions has a solvable double-scaling limit $N\to\infty$ and $q\to\infty$ keeping $\lambda:=\frac{2q^2}{N}=: -\log q$ fixed [30, 31, 32]. At low energies, SYK is dual to JT gravity [29], an effective theory of quantum gravity in two dimensions. Sine-dilaton gravity extends this description to DSSYK beyond the low-energy regime [33, 34].

The Krylov basis constructed from the infinite-temperature thermofield double, using disorder-averaged moments, coincides with the chord-number basis [35, 36, 37]. The chord number gives a notion of a discrete length in the bulk and helps to understand the emergence of a bulk Hilbert space, with semiclassical features when $\lambda\to0$. The dictionary between bulk length and chord number is given by [35, 36, 37, 38, 39, 40]

$\langle\hat{n}(t)\rangle_\beta$ is the Krylov complexity of complexified time evolution, incorporating the Euclidean segment required to prepare the thermofield double. To characterize the corresponding normalized Krylov distribution, we consider the moments

The $\Phi(E)$ relevant for building the Krylov basis is the SYK density of states. In the spectral representation, the Euclidean evolution amounts to replacing $\Phi(E)$ by the normalized thermal measure $Z(\beta)^{-1}e^{-\beta E}\Phi(E)$, while keeping the Krylov basis fixed, see [41] for a discussion. $\langle\tilde{\ell}(t)\rangle$ can be interpreted as the expectation value of the geodesic length of a wormhole in an effective twosided AdS$_2$ black hole geometry whose inverse temperature $\beta_{\mathrm{BH}}$ differs from the DSSYK inverse temperature $\beta$, with the two related as [33]

$\mathcal{J} > 0$ is an overall energy scale present in the SYK Hamiltonian. $\theta_s$ parameterizes the semiclassical saddle. From the Krylov perspective, $\theta_s$ is the asymptotic value of $\langle p(t)\rangle$. In the semiclassical limit,

with length fluctuations suppressed because $\langle\hat{n}^{k}(t)\rangle_{\beta}=\langle\hat{n}(t)\rangle_{\beta}^{k}$ to leading semiclassical order. A bulk prescription for $\langle\hat{n}^{k}(t)\rangle_{\beta}$ was proposed in [42], while their quantum fluctuations at infinite temperature were studied through higher-loop expansions in [23].

(7) provides a direct route to the late-time growth of all chord-number moments, including their quantum corrections. The loop expansion in $v_k$ directly follows from the loop expansion in $\lambda$ of the thermal energy moments, see (37) and the following discussion. The leading answer is

agreeing with (19) for $t\gg\beta_{\mathrm{BH}}$. In this example $t_{\mathrm{on}}\sim O(\lambda^0)$ due to extensive energy variance [21, 36]. The one-loop answer $v_k^{(1)}$ is given in (38). For $k=1$, it reproduces the large-$t$ value of one-loop $\langle\tilde{l}(t)\rangle/t$ computed in [43, 44] after analytic continuation. We show the numerical agreement with (6) in Figure 2 where we also include a one-parameter deformation widely studied in recent times in the context of end-of-the-world branes [45, 46, 47].

Discussion

The heuristics presented in this letter have the merit that they plausibly apply to large, finite systems. There is a parametrically long regime in $n$ where the Lanczos coefficients are expected to be well-described by fluctuations about smooth functions $b(x)$ and $a(x)$ sampled at $x=n/D$, where $D$ is the Krylov space dimension, exponential in system size. The constant Lanczos sequence we utilized emerges from the Taylor expansion $b\left(\frac{n_c}{D}\right)+\left(x-\frac{n_c}{D}\right)\partial_xb+\ldots$. The gradient in the Lanczos coefficients remains parametrically small for $n\ll D$, suggesting an exponentially long window with growth $\langle\hat{n}^{k}(t)\rangle\sim v_kt^k$. The size and correlations of

the fluctuations play a crucial role. For random matrix ensembles, the former decay as $1/\sqrt{D}$ [48, 49], which is favorable for the scaling (6) to persist for exponentially long times. In real systems, the fluctuations can be much larger and are known to impede propagation in the Krylov basis—like for operator dynamics in integrable systems [50]. Since the average $b_n$ is extensive, there could be much larger fluctuations and yet our results could still hold over a time window polynomial in the system size. It is an interesting problem to quantify the size of fluctuations that would be compatible with an approximate form of (7).

The more intriguing point is the possible robustness of the relation between the width of $\Phi(E)$ and the sharpness of the Krylov distribution. It could very well extend far beyond cases where either the $\langle\hat n(k) \rangle\sim v k^t$ scaling or its predicted coefficients (7) are reliable. A preliminary numerical study is presented in Sec. SM C.

In the semiclassical limit of DSSYK, $p$ along with $\lambda\hat n$ become phase-space variables of the dual theory. Our work suggests that a similar classical picture can emerge dynamically at late times in more general systems. In this regime $\hat n$ and $p$ become concentrated around a single trajectory if the state has a narrow spectral spread. This raises the question of when such an emergent classical description is physically meaningful. Our perspective complements recent studies of classicality in Krylov space using different approaches [51, 52, 53, 54, 55].

In DSSYK the chord number operator is proportional to the microscopic size operator [35] within a time regime that overlaps with the regime studied in this work. In this window, $\langle\hat n^k(t) \rangle_\beta$ describes moments of the size distribution of the complexified time evolution operator $e^{-(\beta/2+it)H}$ [56]. Thus (7) determines the leading growth of its size moments, including their fluctuations, directly from the ensemble-averaged density of states. More generally, the operator Krylov complexity moments whose late-time behavior we studied are known to upper bound operator size moments [12]. Both examples we studied suggest a deeper connection between operator size and spectral data that involves the spectral edge $E_0$. Understanding the physical mechanism behind this connection and its broader applicability is an interesting problem for the future.

Acknowledgements— I thank B. Bhattacharjee, A. Dymarsky, S. Aguilar-Gutierrez, O. Parrikar and D. Rosa for discussions that helped this work. I am supported by the São Paulo Research Foundation (FAPESP) through the grant 2024/13100-8. I acknowledge the use of Chat-GPT and Claude for language feedback and LaTeX troubleshooting.

References

End Matter

Orthogonal polynomials and Szego Asymptotics

The Krylov vectors $|n\rangle= p_n(H)|0\rangle$ where $p_n(x)$ is the orthonormal polynomial of order $n$ satisfying the same recursion as (1) with $p_0(x)=1$. If the integral

exists, the $n\gg1$ asymptotics are given by

in terms of $\cos(\theta)=E_0/E$. $\gamma(\theta)$ is also calculable in terms of $\Phi(E)$ but we ignore its precise value in this work.

The propagator associated with the amplitude $\langle E_1|e^{in\hat{\varphi}}|E_2\rangle$ underlies the physical time-evolution of the global mode in the bipartite scar construction of [57]. [21] derived (6), (7) by computing an approximation of this propagator using (22). The integral resulting from inserting this approximate propagator in the generating function $\langle\psi(t)|e^{in\hat{\varphi}}|\psi(t)\rangle$ can be evaluated in the $t\to\infty$ limit using the stationary phase method. Expanding in small $\varphi$ isolates the $v_k t^k$ piece and gives the explicit form (7).

DSSYK Complexity

After ensemble averaging, the energy spectrum of DSSYK is given by

such that $E_0=\frac{2\mathcal{J}}{\sqrt{\lambda(1-q)}}$. The density of states is also conveniently expressed in terms of $\theta$ as

where $\widetilde{\Phi}(\theta)=E_0\Phi(E_0\cos(\theta))\sin(\theta)$. $(a,q)_n$ stands for the q-Pochammer symbol $\prod_{j=1}^{n}(1-aq^{j-1})$. The corresponding Lanczos coefficients are given by

There is a simple generalization at the level of an effective chord Hamiltonian description given by the measure

and Lanczos coefficients

which reduces to (24), (25) when $\kappa=0$. We restrict to $|\kappa|<1$ to avoid bound-state contributions.

Combining (17) with (7) gives

Semiclassical Limit

The semiclassical $\lambda\to0$ limit when $\kappa=0$ is given by

where following [43] we call

The saddle point is defined by $\left.\partial_\theta F\right|_{\theta=\theta_s}=0$ giving rise to the first equation in (18). It is conveniently parametrized using

(18) and (28) imply the classical result $v_k=\left(\frac{2\mathcal{J}}{\lambda}\sin(\theta_s)\right)^k$. The one-loop $Z(\beta)$ is given by [43, 44]

The one-loop answer is accurate away from the edges, i.e., $\lambda\ll\pi-\theta_s$. Using (33) we write

where $\mu_1^{(1)}$ is the $O(1)$ contribution to the energy. To compute the $O(\lambda)$ correction to $\lambda^k v_k$, it is sufficient to keep $\mu_1$ up to $\mu_1^{(1)}$ and energy variance to $O\left(\frac{1}{\lambda}\right)$.

This can be seen by expanding $z_k(E) = (E_0^2-E^2)^{k/2}$ around $\mu_1$

where $z_k^{(r)}$ is the $r$th derivative of $z_k$. It is easy to verify that $z_k^{(r)}(E)/z_k(E)$ scales as $\lambda^r$ by switching variables to $\epsilon\to\lambda E$ and by restricting to $\theta_s$ away from the edges. From the moment-cumulant formula, it can be seen that the central moment in (37) comes with a leading factor of $\lambda^{-\lfloor r/2\rfloor}$, which means that $r$th term in (37) contributes at $\lambda^{\lceil r/2\rceil}$, taking $r \ge2$. To obtain $\lambda^k v_k$ through $O(\lambda^p)$ it suffices to truncate the expansion on the rhs of (37) at $r=2p$, while evaluating the semiclassical free energy $\log Z(\beta)$ up to $O(\lambda^{p-1})$. The one-loop result computed using (33) is

Operator growth

Estimating the timescales— The heuristic derivation in the main text suggests the approximate operator equation

valid for the physical trajectory in the relevant regime, where crucially $X$ is time-independent. Matching the solution at $t=t_{\mathrm{on}}$ gives

where the intercept can be understood as $\langle X\rangle$. Numerics confirm the scaling of intercept in (40), although the two-parameter model does not exactly give the numerical value of the intercept. (39) implies that $\langle\hat{n}(t)^2\rangle$ is a quadratic polynomial in $t$. Importantly, (39) constrains the coefficients of these polynomials in terms of those of $\langle\hat{n}(t)\rangle^2$. We write a quadratic polynomial ansatz for the variance

where $A_2=v_2-v_1^2$. $A_0$ can be approximated by evaluating (41) at $t=t_{\mathrm{on}}$ using the answer from the exponential growth regime, that gives $A_0\approx e^{4\alpha t_{\mathrm{on}}}=\left(\frac{\omega_0}{2\alpha}\right)^2$. We empirically find that $A_1<0$ and $A_2\ll|A_1|\ll A_0$, which we can attribute to small negative covariance of $B_{n_c}$ with $X$. This hierarchy of coefficients predicts approximate $\sqrt{A_0}/(v_1t)$ decay for parametrically long times and a small negative drift below $\sqrt{v_2-v_1^2}/v_1$, before $\delta n(t)$ reaches saturation. These features are confirmed in Figure 1 and Figure S2.

The XY Spin Chain— We consider isotropic XY Hamiltonian in terms of Pauli matrices $\sigma_j^\alpha$ at site $j$:

with periodic boundary conditions, and the $N\to\infty$ limit of the (unnormalized) two-point function

which was evaluated in [27, 58] as

Taking the Fourier transform gives the convolution $\Phi(E)=(g_\beta*g_\beta)(E)$ in terms of

where $\Theta$ is the Heaviside function. The claim $\omega_0=2$ follows directly from the support of $g_\beta(x)$.

Numerics using truncated Jacobi matrix— For numerical computations of complexity, we truncate the matrix $M$ at $n_{\mathrm{cut}}$. Consistency of the answers obtained using the cutoff is enforced by ensuring that $\varphi_n(t)\approx0$ (up to numerical precision) for the sampled $t$, for $n$ much smaller than $n_{\mathrm{cut}}$.

For (44), it is straightforward to sample moments with enough accuracy such that $n$ is much larger than $n_{\mathrm{cut}}$ and the $b_n$ coefficients have converged to $\omega_0/2$ up to small $<10^{-3}$ fluctuations. To access late times to see the saturation of $\delta n(t)$, we padded the numerically computed $M$ with $\omega_0/2$ along the off-diagonal. We verified that this padding does not change the time-dependence of $\langle\hat{n}(t)\rangle$ to any observable effect in this regime.

Krylov amplitudes from saddle point

The asymptotic form (22) provides a good starting point for computing $\varphi_n(t)$

using the stationary phase method, where we remind that $\widetilde{\Phi}(\theta)=E_0\Phi(E_0\cos\theta)\sin(\theta)$. We expand the $\cos(n\theta+\gamma)$ and keep

where we reuse $F(\theta)$ to mean

Assuming that $\partial_\theta\gamma(\theta)$ isn’t negative with a very large modulus, the integrand is localized around $\partial_\theta F(\theta)=0$ for $t>0$. Now we will further drop $\partial_\theta\gamma(\theta)$ as a first approximation. This is not necessarily justified, but is done for the sake of simplicity. Then the stationary points are defined by the equation

which has solutions $\theta_1=\arcsin(n/E_0t)$ and $\theta_2=\pi-\theta_1$. Evaluating $\partial_\theta^2F(\theta)$ at the stationary points gives $r_1:=\partial_\theta^2F(\theta)|_{\theta=\theta_1}=-\sqrt{(E_0t)^2-(n)^2}$ and $r_2=-r_1$. The stationary point answer gives

in terms of $\chi=-n\theta_1-\sqrt{E_0^2t^2-n^2}+\frac{\pi}{4}$. The overall phase in (S5) should not be trusted because $\gamma$ was neglected. (S5) diverges as $n=E_0t$. This signals the breakdown of the Gaussian approximation as the two saddles $\theta_1$ and $\theta_2$ merge at $\frac{\pi}{2}$. Expanding $F(\theta)$ to cubic order around $\theta=\frac{\pi}{2}$ gives the expansion near $n=E_0t$

To get a simple analytical approximation for (S6), we assume that $\widetilde{\Phi}(\theta)$ is slowly varying around $\theta=\frac{\pi}{2}$, so that we can approximate the integral by only keeping $\sqrt{\widetilde{\Phi}(\frac{\pi}{2})}$ in the Taylor expansion. Then the integral over $\theta$ can be cast into the standard Airy integral by extending integration domain to $\mathbb{R}$. This gives the answer near $n\sim E_0t$

where $\xi=\left(\frac{2}{E_0t}\right)^{\frac{1}{3}}(E_0t-n)$. (S7) resembles the continuum limit answers obtained by [15]. The oscillatory asymptotics of (S7) matches on to (S5). It is easy to see this for the prefactor. To match the phase we can take and $e^{i\chi}+e^{-i\chi}\sim\cos(\chi)$ and then expand in $(E_0t-n)/E_0t$ to get $\sin\left(\frac{2}{3}\xi^{3/2}+\frac{\pi}{4}\right)$ which matches the oscillatory piece in the expansion of $\operatorname{Ai}(-\xi)$ for $\xi>0$. (S7) should be matched on to (S5) for some $n_{\mathrm{match}}$ such that $E_0t-n_{\mathrm{match}}\lesssim O(1)(E_0t)^{\frac{1}{3}}$. Together (S5), (S7) can provide a remarkably good approximation for $|\varphi_n(t)|^2$.

In Figure S1, we present an example for DSSYK. A word of caution is that numerics show that the phase $\gamma(\theta)$ can cause important corrections when taking the $\lambda\to0$ limit at fixed $t$.

Varying bandwidth

We numerically studied complexity dynamics for the measure

and $\Phi(\omega)=0$ for $\omega>\omega_0$ to substantiate the model following [37]. (S8) shows a distinct discontinuous edge-behavior. We use it to examine cases where $\alpha$ is held fixed but $\omega_0$ is allowed to vary, to complement the study of the XY model spectral function where $\alpha$ varied but $\omega_0$ was held fixed. Figure S2 confirms that $\delta_n(t)$ decays but $\frac{1}{\alpha t}$ independent of $\omega_0$, before drifting slightly below its infinite time asymptotic.

Finite size signatures

We numerically study the Hamiltonian

The three-body interaction is chosen to preserve a spectral reflection symmetry while maintaining nonintegrability. For the initial states considered, this symmetry enforces $a_n=0$ reducing overall fluctuations. For the system sizes studied, $11\leq N\leq14$, the expected linear-growth regime is not yet clearly resolved.

Nonetheless $\delta n(t)$ shows a clear trend that we illustrate with two examples

The density of states $\rho(E)$ is approximately Gaussian towards the center of the spectrum with variance growing with $N$. For $|\psi_1\rangle$, $\Phi(E)=\rho(E)$, while for $|\psi_2\rangle$, $\Phi(E)$ has a fixed width $\sigma=4$ independent of $N$. Since the spectral bandwidth grows linearly with $N$, the width of $\Phi(E)$ relative to the bandwidth scales as $N^{-1/2}$ and $N^{-1}$, respectively.

Figure S3 shows $\delta n(t)$ for a time window that could become the linear growth regime in larger systems. $\delta n(t)$ changes slowly over this window and is seen to become flatter for larger $N$. $\delta n$ is inversely correlated with the width of $\Phi(E)$ rescaled by the bandwidth. This correlation persists across different Hamiltonians and initial states, but eventually saturates when the width $\sigma$ becomes sufficiently small.