# Exactly solvable models for universal operator growth

Oleksandr Gamayun, Murtaza Ali Mir, Oleg Lychkovskiy, Zoran Ristivojevic

► To cite this version:

Oleksandr Gamayun, Murtaza Ali Mir, Oleg Lychkovskiy, Zoran Ristivojevic. Exactly solv-able models for universal operator growth. Journal of High Energy Physics, 2025, 2025 (256), 10.1007/JHEP07(2025)256. [hal-05378783](https://hal.science/hal-05378783v1)

HAL Id: hal-05378783

https://hal.science/hal-05378783v1

Submitted on 23 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Exactly solvable models for universal operator growth

Oleksandr Gamayun$^{a,1}$ Murtaza Ali Mir$^{b,2}$ Oleg Lychkovskiy$^{c,3}$ Zoran Ristivojevic$^{d,4}$

$^a$London Institute for Mathematical Sciences
Royal Institution, 21 Albemarle St, London W1S 4BS, UK

$^b$Centre of Excellence, Chinar Quantum AI,
F3 Uptown Tower, Hyderpora, Srinagar, India

$^c$Castiglioncello, Italy

$^d$Laboratoire de Physique Theorique, Universit de Toulouse, CNRS
118 Route de Narbonne, 31062 Toulouse Cedex 9, France

E-mail: og@lims.ac.uk, murtaza.mir@skoltech.ru, lychkovskiy@gmail.com,
zoran.ristivojevic@utoulouse.fr

**ABSTRACT:** Quantum observables of generic many-body systems exhibit a universal pattern of growth in the Krylov space of operators. This pattern becomes particularly manifest in the Lanczos basis, where the evolution superoperator assumes the tridiagonal form. According to the universal operator growth hypothesis, the nonzero elements of the superoperator, known as Lanczos coefficients, grow asymptotically linearly. We introduce and explore broad families of Lanczos coefficients that are consistent with the universal operator growth and lead to the exactly solvable dynamics. Within these families, the subleading terms of asymptotic expansion of the Lanczos sequence can be controlled and fine-tuned to produce diverse dynamical patterns. For one of the families, the Krylov complexity is computed exactly.

ARXIV EPRINT: 2504.03435

<footnote>$^1$Corresponding author.</footnote>
# Contents

1 Introduction: recursion method and universal operator growth 1

2 Recursion method and orthogonal polynomials 4

3 A family of models based on the continuous Hahn polynomials 6

3.1 General model 6

3.2 Particular cases 7

3.3 Asymptotic behavior 8

3.4 Derivation based on the properties of the hypergeometric functions 10

3.5 Krylov complexity 10

4 A family of models with odd-even alterations in the Lanczos sequence 12

5 Correlation functions relaxing to nonzero stationary values 14

5.1 Lanczos coefficients from stationary value 14

5.2 The stationary value from the Lanczos coefficients 17

6 Discussion 18

A Continuous Hahn polynomials 19

A.1 Fourier transform 20

A.2 Transformation to the form of eq. (1.1) 20

A.3 The case $\tilde{a}_n = 0$ 21

A.4 Evaluation of $\phi_n(t)$ 22

# 1 Introduction: recursion method and universal operator growth

Describing the dynamics of quantum many-body systems is among the main objectives of condensed matter and quantum field theories. This task has a reputation of being extremely complex in general. The most well-understood tractable scenario emerges if a system under study is close in some sense to a collection of noninteracting (quasi)particles. In this case a diverse and sophisticated toolbox of perturbative techniques is available and often sufficient for an exhaustive quantitative description. The opposite case of a generic quantum many-body system far from any free-particle point has been until recently widely believed to be intractable (apart from some exceptional types of systems and special techniques, including AdS/CFT correspondence [1], low-dimensional systems solvable by matrix product state methods [2, 3], etc.). This assessment is being revised nowadays, in large part thanks to the universal operator growth hypothesis proposed in ref. [4]. The latter behavior emerges in the framework of the recursion method developed long ago [5-7].
In essence, the recursion method amounts to solving coupled Heisenberg equations of motion in the tridiagonal Lanczos basis in operator space. For our purposes, it is enough to briefly outline the outcome of the recursion method without going into derivations and technical details [4, 7, 8]. We consider a many-body system described by a local Hamiltonian $H$ and focus on the most basic object within the method – an infinite-temperature autocorrelation function. It is defined as $C(t) = \text{Tr}\{\mathcal{O}(t)\mathcal{O}(0)\}/\text{Tr}\{\mathcal{O}^2\}$, where $\mathcal{O}$ is some local observable that does not depend on time in the Schrödinger picture. In the Heisenberg picture, it depends on time and evolves according to the Heisenberg equation of motion, $\partial_t\mathcal{O}(t) = i[H, \mathcal{O}(t)]$. Here and in the following we set the reduced Planck constant to one, $\hbar = 1$. In order to compute the autocorrelation function, it turns out that it is sufficient to solve the discrete Schrödinger equation

$$ \partial_t \varphi_n(t) = -b_{n+1} \varphi_{n+1}(t) + b_n \varphi_{n-1}(t), \quad n = 0, 1, 2, \dots, \quad (1.1) $$

subject to the conditions

$$ \varphi_n(0) = \delta_{n,0}, \quad \varphi_{-1}(t) = 0. \quad (1.2) $$

Equation (1.1) describes a fictitious single particle on a semi-infinite discrete interval, which is described by the wave function $\varphi_n(t)$. At $t = 0$ it is in the origin and propagates in time. The set of positive coefficients $b_n \ge 0$ in eq. (1.1) are known as Lanczos coefficients and play a major role. The solution of eq. (1.1) is a set of real functions $\varphi_n(t)$ for $n = 0, 1, 2, \dots$. It satisfies the normalization condition $\sum_{n=0}^{\infty} \varphi_n(t)^2 = 1$. At $t = 0$, this is obvious due to the initial condition (1.2). At $t > 0$, differentiating the normalization condition it follows $\sum_{n=0}^{\infty} \varphi_n(t) \partial_t \varphi_n(t) = 0$, which is satisfied due eq. (1.1). Once eq. (1.1) is solved, the autocorrelation function is simply given by

$$ C(t) = \varphi_0(t). \quad (1.3) $$

We note that $C(0) = 1$.

Lanczos coefficients $b_n$ of eq. (1.1) depend on the Hamiltonian $H$ and the operator $\mathcal{O}$. They can be obtained by the orthogonalization of the sequence of operators $\mathcal{L}^n\mathcal{O}$ for $n = 0, 1, 2, \dots$, and represent the norms of the orthogonal set of operators. Here $\mathcal{L}$ denotes the Liouville superoperator that acts on operators by the commutator as $\mathcal{L}\mathcal{O} = [H, \mathcal{O}]$. Therefore, in order to obtain $\mathcal{L}^n\mathcal{O}$ for a given $n$, nested commutators should be evaluated. The computational complexity of this task grows exponentially with $n$, and therefore, in practice, only a finite number $n_{\max}$ of Lanczos coefficients is typically available. As a consequence, eq. (1.1) cannot be solved. While the truncation of this equation at $n = n_{\max}$ provides an accurate approximation for $C(t)$ at short times, it breaks down at longer times. This was a severe limiting factor of the recursion method for decades [7].

The crucial step to resolve the above stalemate was made in ref. [4] (see also a precursor work [9–13]), where a universal operator growth hypothesis has been put forward. It states that for a generic (in particular, nonintegrable) many-body system, Lanczos coefficients grow asymptotically linearly with $n$ (with an additional logarithmic correction in one-dimensional systems – a case not addressed in the present paper):

$$ b_n = \alpha n + o(n), \quad n \to \infty. \quad (1.4) $$
The universal operator growth hypothesis (1.4) has been subsequently confirmed explicitly for various many-body models [14–23]. We note that the asymptotic behavior (1.4) is typically slowed down for integrable systems [4]. The universal operator growth hypothesis implies the exponential growth of the Krylov complexity [4]

$$K(t) = \sum_{n=0}^{\infty} n \varphi_n(t)^2, \qquad (1.5)$$

which is regarded to be a measure for the operator growth [24, 25]. The latter is an important and ubiquitous phenomenon that appears in a variety of contexts ranging from quantum optics [26] and quantum networks [27] to cosmology [28], black hole physics [29, 30], holography [31–34], quantum billiards [35] and conformal field theories [26, 30, 31, 36, 37].

The universal operator growth hypothesis (1.4) signals that the recursion method might be augmented by replacing the unknown Lanczos coefficients $b_n$ for $n > n_{max}$, by their extrapolated counterparts $b_n^{extr} = \alpha n$, where $\alpha$ is found by fitting the known Lanczos coefficients. This procedure can admit a perturbative guise as follows. One first introduces an unperturbed Schrödinger equation of the form (1.1) with the coefficients exactly linear in $n$,

$$b_n^0 = \alpha n, \quad n = 0, 1, 2, \dots \qquad (1.6)$$

The actual Schrödinger equation is obtained by perturbing each coefficient according to $b_n = b_n^0 + \delta b_n$ with $\delta b_n = b_n - \alpha n$, cf., ref. [38]. The universal operator growth hypothesis (1.4) guarantees that the perturbation $\delta b_n$ vanishes for large $n$. Note that a possibly large perturbation at a few first sites of the semi-infinite chain can be addressed separately [4, 39–41]. Importantly, the unperturbed dynamics governed by the Lanczos coefficients (1.6) can be solved exactly [4], leading to $C(t) = \mathrm{sech}(\alpha t)$.

While a perturbative scheme along these lines appears feasible, its actual implementation is far from being straightforward [4, 17]. One particularly serious issue is that the subleading terms in the asymptotic expansion of $b_n$ (those hidden in $o(n)$ in eq. (1.4)) can have a strong impact on $\varphi_0(t)$, to the extent that the qualitative behavior of the autocorrelation function is altered [4, 17, 32, 36, 38, 42–48]. It is therefore highly desirable to have a large set of sequences that satisfy the universal operator growth hypothesis (1.4), encompass various subleading terms, and lead to exactly solvable Schrödinger equation (1.1). We will refer to such sequences satisfying the latter requirement as *exactly solvable Lanczos sequences*. A wise choice of a suitable exactly solvable Lanczos sequence as a zero-order approximation could considerably improve the perturbative scheme for a particular many-body model.

The sequence (1.6) is in fact a particular case of a more general exactly solvable sequence given by [4]

$$b_n = \sqrt{n(n-1+\eta)}, \qquad (1.7)$$

with

$$\varphi_n(t) = \sqrt{\frac{(\eta)_n}{n!} \tanh^n(t) \cosh^n(t)}. \qquad (1.8)$$
Here $(\eta)_n = \Gamma(\eta + n)/\Gamma(\eta)$ is the Pochhammer symbol, and $\eta \ge 1$ is a free parameter. By varying the parameter $\eta$, one can control the first subleading term in eq. (1.7) at large $n$. In eq. (1.7) and in the following we have omitted the superscript from $b_n^0$ and set $\alpha = 1$ for simplicity, which amounts to an appropriate rescaling of time. Indeed for the Lanczos coefficients $\tilde{b}_n = \alpha b_n$, the solution of eq. (1.1) is given by $\varphi_n(\alpha t)$, where $\varphi_n(t)$ is defined in eq. (1.8). Many other examples of exactly solvable sequences are studied in ref. [8], where some are in accordance with the universal operator growth hypothesis (1.4) and some are not.

In the present paper we introduce and analyze several families of presumably unknown exactly solvable Lanczos sequences consistent with the universal operator growth hypothesis. We demonstrate that these multiparameter families enable fine tuning of the subleading terms and lead to autocorrelation functions with a diverse qualitative behavior. The remaining part of the paper is organized as follows. In section 2 we provide a brief recap of the recursion method and its relation to the orthogonal polynomials. A family of models based on the continuous Hahn polynomials is introduced and explored in section 3. In section 4 we introduce another family of models whose distinctive feature is the sign alteration in the subleading terms. In section 5 we explore the Lanczos coefficients for the correlation functions with nonzero stationary late-time values. In section 6 some broader implications of the obtained rigorous results are discussed. Some technical details about the continuous Hahn polynomials are given in appendix A.

## 2 Recursion method and orthogonal polynomials

Equation (1.1) can be solved by employing the method of orthogonal polynomials [5, 6, 8, 49]. In order to achieve that, we introduce the system of polynomials satisfying the standard recurrence relation of the form

$$x\pi_n(x) = b_{n+1}\pi_{n+1}(x) + b_n\pi_{n-1}(x), \qquad (2.1)$$

where $b_n$ is a given set of the coefficients. The initial members of the system are $\pi_{-1} = 0$ and $\pi_0 = 1$. Then the formal solution of eq. (1.1) can be expressed as

$$\varphi_n(t) = i^n \pi_n(i \partial_t) C(t), \quad \varphi_0(t) = C(t). \qquad (2.2)$$

Thanks to the Favard theorem [50], the polynomials $\pi_n$ obeying the three-term recurrence relation (2.1) are orthogonal with respect to some weight (orthogonality measure) $\rho(x)$. We thus have

$$\int_{-\infty}^{\infty} dx \rho(x) \pi_n(x) \pi_m(x) = \delta_{n,m}, \quad \int_{-\infty}^{\infty} dx \rho(x) = 1. \qquad (2.3)$$

Requiring that the initial condition $\varphi_n(0) = \delta_{n,0}$ is satisfied for the solution (2.2) one can deduce a connection between $C(t)$ and the weight $\rho(x)$. Introducing the Fourier transform by $C(t) = \int_{-\infty}^{\infty} dx e^{-itx} \hat{C}(x)$, eq. (2.2) becomes

$$\varphi_n(t) = i^n \int_{-\infty}^{\infty} dx e^{-itx} \pi_n(x) \hat{C}(x). \qquad (2.4)$$
Comparing the initial condition (2.4) at $t = 0$ with the orthogonality condition (2.3) taken at $m = 0$, we infer that the weight corresponds to the Fourier transform of the autocorrelation function, $\rho(x) = \hat{C}(x)$. Therefore we obtain

$$C(t) = \int_{-\infty}^{\infty} dx e^{-itx} \rho(x). \qquad (2.5)$$

Equation (2.5) completes the solution of eq. (1.1). Therefore, any set of orthogonal polynomials that satisfies the three-term recurrence relation (2.1) and obeys the orthogonality condition (2.3), provides one solution for the semi-infinite discrete Schrödinger equation (1.1). Since the sets of orthogonal polynomials are much widely studied than eq. (1.1), the latter knowledge is useful to explore eq. (1.1) and its consequences. As a simple example let us consider the classical Hermite polynomials $H_n(x)$. After proper rescaling, we obtain the polynomials $\pi_n(x) = H_n(x/\sqrt{2})/\sqrt{2^n n!}$ that satisfy eq. (2.1) with $b_n = \sqrt{n}$ (note that this does not conform with the universal operator growth hypothesis (1.4)). Their weight is $\rho(x) = e^{-x^2/2}/\sqrt{2\pi}$. Using eq. (2.4) we then find the corresponding wave functions, $\varphi_n(t) = t^n e^{-t^2/2}/\sqrt{n!}$. The latter can also be checked by a direct substitution into eq. (1.1).

Note that the wave functions $\varphi_n(t)$ together with the orthogonal polynomials $\pi_n(x)$ allow for a remarkable representation of the evolution operator. Indeed, let us consider an operator function

$$F(t, \hat{A}) = \sum_{n=0}^{\infty} i^n \varphi_n(t) \pi_n(\hat{A}). \qquad (2.6)$$

Using the initial condition $\varphi_{-1} = \pi_{-1} = 0$ it is easy to show that $F(0, \hat{A}) = 1$ and $\partial_t F(t, \hat{A}) = i\hat{A}F(t, \hat{A})$. This is the first-order linear differential equation that has a unique solution

$$F(t, \hat{A}) = e^{it\hat{A}}. \qquad (2.7)$$

Here $\hat{A}$ is an arbitrary time-independent operator. In the special case when $\hat{A}$ is the Liouville superoperator we thus obtain

$$e^{it\mathcal{L}} = \sum_{n=0}^{\infty} i^n \varphi_n(t) \pi_n(\mathcal{L}). \qquad (2.8)$$

Equation (2.8) is the expansion of the evolution superoperator in terms of the polynomials, see for example refs. [51–55]. In fact, eq. (2.8) can serve as a starting point for a systematic perturbative expansion, whose convergence can be made uniform in time by choosing a set of $\varphi_n(t)$ with the long-time asymptotics fitting that of the actual autocorrelation function. This remark goes beyond the scope of the recursion method, highlighting the importance of orthogonal polynomials in a broader context of quantum dynamics.
# 3 A family of models based on the continuous Hahn polynomials

## 3.1 General model

Equations (1.7) and (1.8) provide an example [4] of the one-parameter family of the solutions of eq. (1.1) based on the Meixner-Pollaczek polynomials<fnref n="1" /> [58]. Even though the polynomials are very involved, the wave functions (1.8) have a relatively simple form.

Here we generalize the latter one-parameter solution to a two-parameter solution of (1.1) with the linear growth of the corresponding Lanczos coefficients. Instead of Meixner-Pollaczek polynomials, our solution is based on the continuous Hahn polynomials [58]. It has the Lanczos coefficients given by

$$b_n^2 = \frac{4n(n+2a-1)(n+2b-1)(n+2a+2b-2)}{(2n+2a+2b-3)(2n+2a+2b-1)}, \quad (3.1)$$

and the wave functions

$$\varphi_n(t) = \mathcal{F}_n \frac{\tanh^n(t)}{\cosh^{4a}(t)} F\left(\begin{array}{l} a-b+\frac{1}{2}, 2a+n \\ a+b+n+\frac{1}{2} \end{array}; \tanh^2(t)\right). \quad (3.2)$$

Here the prefactor is given by<fnref n="2" />

$$\mathcal{F}_n = \frac{1}{n!} \left( \prod_{k=1}^{n} b_k \right) = \sqrt{\frac{(2a)_n (2b)_n (2a+2b-1)_n}{n! (a+b-1/2)_n (a+b+1/2)_n}}, \quad (3.3)$$

and

$$F\left(\begin{array}{c} a, b \\ c \end{array}; x\right) = \frac{\Gamma(c)}{\Gamma(a)\Gamma(b)} \sum_{k=0}^{\infty} \frac{\Gamma(a+k)\Gamma(b+k)}{\Gamma(c+k)} \frac{x^k}{k!} \quad (3.4)$$

is the Gauss hypergeometric function. The solution (3.2) depends on two parameters a and b that can be either both positive real numbers

$$a > 0, \quad b > 0, \quad (3.5)$$

or a pair of complex conjugate numbers with a positive real part,

$$a = r + i\omega, \quad b = r - i\omega, \quad r > 0. \quad (3.6)$$

The wave functions (3.2) have several equivalent representations that are discussed in appendix A. Here we only list one such representation where the symmetry to the interchange of a and b is obvious,

$$\varphi_n(t) = \mathcal{F}_n \frac{\tanh^n(t)}{\cosh^{2a+2b-1}(t)} F\left(\begin{array}{l} a-b+\frac{1}{2}, b-a+\frac{1}{2} \\ a+b+n+\frac{1}{2} \end{array}; -\sinh^2(t)\right). \quad (3.7)$$

<footnote n="1">¹While the Meixner-Pollaczek polynomials depend on two parameters, here we have in mind their one-parametric subclass with the weight given by the square of the gamma function.</footnote>

<footnote n="2">²Note that a seeming similarity between $b_n$ and b is accidental as they denote different objects.</footnote>
The form (3.7) manifestly demonstrates that $\varphi_n(t)$ is real for all admissible values of $a$ and $b$ as all the summands in eq. (3.4) are real for both choices (3.5) and (3.6). Moreover, $\varphi_n(t)$ is positive at any $t$ for real and positive $a$ and $b$. We note that the weight that corresponds to the exactly solvable sequence (3.1) is given by

$$ \rho(x) = \frac{\Gamma(2a + 2b)}{8\pi\Gamma(2a)\Gamma(2b)\Gamma(a + b)^2} |\Gamma(a + ix/4)\Gamma(b + ix/4)|^2. \quad (3.8) $$

Via eq. (2.5) it determines $C(t) \equiv \varphi_0(t)$.

## 3.2 Particular cases

The exactly solvable sequence (3.1) with the solution (3.2) has two parameters and in special cases reduces to more elementary solutions that will be studied in the following. This simplification happens for the values of $a$ and $b$ for which the hypergeometric functions appearing in eqs. (3.2) or (3.7) reduce to a simpler form, or for $a$ and $b$ for which the Fourier transform (2.5) of the weight (3.8) simplifies. For real parameters, the latter is achieved in the cases: (i) $a$ and $b$ are both positive integers, (ii) $a$ and $b$ are both odd half-integers, (iii) one parameter between $a$ and $b$ is a positive integer and the other is an odd half-integer, and (iv) the difference $b-a$ is an odd half-integer. There we can use the relations

$$ |\Gamma(1 + k + ix/4)|^2 = \frac{\pi x/4}{\sinh(\pi x/4)} \prod_{j=1}^{k} (j^2 + x^2/4), \quad k \in \mathbb{N}, \quad (3.9) $$

$$ |\Gamma(1/2 + k + ix/4)|^2 = \frac{\pi}{\cosh(\pi x/4)} \prod_{j=1}^{k} ((j - 1/2)^2 + x^2/4), \quad k \in \mathbb{N}, \quad (3.10) $$

$$ \Gamma(z)\Gamma(z + 1/2) = \sqrt{\pi} 2^{1-2z} \Gamma(2z), \quad (3.11) $$

in conjunction with $\Gamma(z+1) = z\Gamma(z)$ for a complex $z$ with a positive real part.

Let us work out several explicit results. Consider the case $b = a + k + 1/2$ with a positive integer $k$. At $k=0$, for the choice $a = \eta/4$, using, e.g, the Ramanujan formula for the Fourier transform of the square of the Gamma function we obtain

$$ C(t) = \frac{1}{\cosh^n(t)}, \quad b_n = \sqrt{n(n-1+\eta)}. \quad (3.12) $$

Moreover, in this case the hypergeometric function in eq. (3.2) becomes unity and we obtain the wave functions of eq. (1.8). Equation (3.12) corresponds to the solution based on the Meixner-Pollaczek polynomials, which was previously obtained by different means in ref. [4].

For $k=1$ we have $b = a + 3/2$ and

$$ C(t) = \frac{a + \cosh^2(t)}{(a + 1) \cosh^{4a+2}(t)}, \quad b_n^2 = \frac{n(2a + n - 1)(2a + n + 2)(4a + n + 1)}{(2a + n)(2a + n + 1)}. \quad (3.13) $$

Another elementary expression is obtained for the case $a = b = 1/2$. It is given by

$$ C(t) = \frac{2t}{\sinh(2t)}, \quad b_n^2 = \frac{4n^4}{4n^2 - 1}. \quad (3.14) $$
We also notice that for $a = b = 1/4$ the result can be expressed as

$$C(t) = \frac{2K(\tanh^2(t))}{\pi \cosh(t)}, \quad b_n = n - \frac{1}{2}, \qquad (3.15)$$

where K denotes the complete elliptic integral of the first kind. Notice an extremely simple form of $b_n$ in this case.

For the complex conjugated series, a notable example corresponds to $a = b^* = (1 + i\omega)/2$. We find

$$C(t) = \frac{\sin(2\omega t)}{\omega \sinh(2t)}, \quad b_n = \frac{4n^2(n^2 + \omega^2)}{4n^2 - 1}. \qquad (3.16)$$

Therefore, the autocorrelation function shows the damped oscillations. At $\omega \to 0$, eq. (3.16) reduces to the result (3.14).

### 3.3 Asymptotic behavior

Let us study the asymptotic behavior of the autocorrelation function $C(t)$ at late times, $t \to \infty$. In the case of real parameters $b > a > 0$, $C(t)$ can be obtained from eq. (3.2) by setting $\tanh(t)$ to 1 in the hypergeometric function. Using the property<fnref n="3" />

$$F\left(\begin{matrix} a, b \\ c \end{matrix}; 1\right) = \frac{\Gamma(c)\Gamma(c-a-b)}{\Gamma(c-a)\Gamma(c-b)}, \quad \text{Re}(c-a-b) > 0, \qquad (3.17)$$

we obtain

$$C(t) \approx \frac{\Gamma(2a + 2b)\Gamma(b - a)}{\Gamma(a + b)\Gamma(2b)} e^{-4at}, \quad b > a > 0. \qquad (3.18)$$

The case $a > b > 0$ follows from the symmetry as $C(t)$ is invariant to the exchange of $a$ and $b$. Finally the case $a = b$ is more difficult and can be obtained using the expansion<fnref n="4" />

$$F\left(\begin{matrix} a, b \\ a+b \end{matrix}; z\right) = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \left(\ln \frac{1}{1-z} + \frac{\Gamma'(a)}{\Gamma(a)} + \frac{\Gamma'(b)}{\Gamma(b)} + 2\gamma_E\right) (1+O(1-z)), \qquad (3.19)$$

where $z \to 1$ and $\gamma_E$ is the Euler constant. This yields

$$C(t) \approx \frac{4\Gamma(4a)}{[\Gamma(2a)]^2} te^{-4at}, \qquad (3.20)$$

which is consistent with eqs. (3.14) and (3.15). Equations (3.18) and (3.20) reveal that the autocorrelation function for the model (3.1) decays as $C(t) \sim e^{-4\min(a,b)t+\delta_{a,b}\ln t}$. Note the logarithmic correction in the case $a = b$, which makes the decay slower.

In the complex-conjugate case $a = b^*$, we transform eq. (3.2) by making use of the expansion<fnref n="5" />

$$F\left(\begin{matrix} a, b \\ c \end{matrix}; z\right) = \left[ \frac{\Gamma(c)\Gamma(a+b-c)}{\Gamma(a)\Gamma(b)}(1-z)^{c-a-b} + \frac{\Gamma(c)\Gamma(c-a-b)}{\Gamma(c-a)\Gamma(c-b)} \right] (1+O(1-z)), \\ c-a-b \notin \mathbb{Z} \qquad (3.21)$$

<footnote n="3">³https://functions.wolfram.com/07.23.03.0002.01</footnote>

<footnote n="4">⁴https://functions.wolfram.com/07.23.06.0014.01</footnote>

<footnote n="5">⁵https://functions.wolfram.com/07.23.06.0008.01</footnote>
where $z \to 1$. This gives

$$C(t) \approx \frac{2^{4r}}{\sqrt{\pi}} \Gamma\left(2r + \frac{1}{2}\right) \mathrm{Re}\left(\frac{\Gamma(2i\omega)e^{4i\omega t}}{\Gamma(2r + 2i\omega)}\right) e^{-4rt}, \quad (3.22)$$

where the notation of eq. (3.6) is used. Note that at $\omega \to 0$, eq. (3.22) reduces to eq. (3.20), as it must be the case.

We eventually note that eq. (3.22) can also be expressed as

$$C(t) \approx \frac{\Gamma(2a + 2b)}{\Gamma(a + b)} \left( \frac{\Gamma(a - b)e^{-4bt}}{\Gamma(2a)} + \frac{\Gamma(b - a)e^{-4at}}{\Gamma(2b)} \right). \quad (3.23)$$

Equation (3.23) is quite general and applies to both, complex and real cases. It reduces to the special cases (3.18), (3.20), and (3.22) if the corresponding limits are taken.

From the asymptotic expansions (3.18) and (3.20) we note that the decay of $C(t)$ is monotonic for real parameters $a$ and $b$. This should be contrasted to the case of complex parameters (3.22), where the damped oscillations occur. The typical behavior of the autocorrelation function is shown in figure 1.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Two typical behaviors of the autocorrelation function $C(t) = \varphi_0(t)$ given by eq. (3.2). Depending on the choice of parameters indicated in the plots, the autocorrelation function can either feature damped oscillations (left plot) or damping without oscillations (right plot). We also see that the case $a = b = 3/4$ shows slower relaxation from the one of $a = 3/4, b = 5/4$ due to an additional factor $t$ that multiplies the exponential decay in $C(t)$, see eq. (3.20).</figcaption></figure>

Let us connect the late-time asymptotic behavior of $C(t)$ with the asymptotic properties of the Lanczos coefficients $b_n$. Expansion of $b_n$ of eq. (3.1) at $n \to \infty$ is given by

$$b_n = n + a + b - 1 - \frac{(a + b - 1)^2 + (a - b)^2 - 1/4}{2n} + O(1/n^2). \quad (3.24)$$

We can see that if an asymptotic expansion of a generic $b_n$ with linear growth is denoted as $b_n = n + b^{(1)} + b^{(2)}n^{-1} + o(n^{-1})$, then the solution (3.2) can be used to model a system with the subleading coefficients $b^{(1)} > -1$ and arbitrary $b^{(2)}$. To discriminate the two cases of $\omega = 0$ and $\omega \neq 0$ in terms of $b^{(1)}$ and $b^{(2)}$, one needs to compute the sign of

$$\frac{1}{8} - \frac{(b^{(1)})^2}{2} - b^{(2)}. \quad (3.25)$$
It is positive whenever the parameters $a$ and $b$ are real and $a \neq b$, corresponding to $\omega = 0$. Conversely, the expression (3.25) is negative if $a$ and $b$ are complex-conjugate pairs, which corresponds to $\omega \neq 0$. Finally the expression (3.25) nullifies if $a = b$. In this case there is an additional $t$ factor in $C(t)$, see eq. (3.20).

## 3.4 Derivation based on the properties of the hypergeometric functions

Detailed construction of the functions (3.2) [or equivalently (3.7)] that solve eq. (1.1) with the Lanczos coefficients (3.1) and their relation to the continuous Hahn polynomials is given in appendix A. Here we show that eq. (3.7) solves the system (1.1) directly by making use of some basic properties of the hypergeometric functions. Substituting eq. (3.7) into eq. (1.1), after the change of variables we obtain

$$ \left( n + x(2a + 2b - 1) + 2x(1-x) \frac{d}{dx} \right) \Phi_n(x) = n(1-x)\Phi_{n-1}(x) + \frac{b_{n+1}^2}{n+1} x \Phi_{n+1}(x), \quad (3.26) $$

where we have defined

$$ \Phi_n(x) = F \left( \begin{matrix} a - b + \frac{1}{2}, \ b - a + \frac{1}{2} \\ a + b + n + \frac{1}{2} \end{matrix} ; x \right). \quad (3.27) $$

The derivative can be evaluated by means of the property<fnref n="6" />

$$ \left[ c - a - b + (1 - z) \frac{d}{dz} \right] F \left( \begin{matrix} a, b \\ c \end{matrix} ; z \right) = \frac{(c-a)(c-b)}{c} F \left( \begin{matrix} a, b \\ c+1 \end{matrix} ; z \right), \quad (3.28) $$

leading to

$$ (1 - 2x) \Phi_n(x) = (1 - x) \Phi_{n-1}(x) - \frac{4(2a + n)(2b + n)}{4(a + b + n)^2 - 1} x \Phi_{n+1}(x). \quad (3.29) $$

The latter relation holds as one of the contiguous properties for the hypergeometric functions<fnref n="7" />. We have therefore shown that the functions (3.7) [or equivalently (3.2)] are the solution of eq. (1.1).

## 3.5 Krylov complexity

The spread of the operator dynamics can be quantified by the Krylov complexity $K(t)$. It is defined by [4]

$$ K(t) = \sum_{n=0}^{\infty} n \varphi_n(t)^2 = \frac{dK_{\lambda}(t)}{d\lambda} \bigg|_{\lambda=1}, \quad K_{\lambda}(t) \equiv \sum_{n=0}^{\infty} \lambda^n \varphi_n(t)^2. \quad (3.30) $$

Using an integral representation for the hypergeometric function<fnref n="8" /> given by

$$ F \left( \begin{matrix} a - b + \frac{1}{2}, \ 2a + n \\ a + b + n + \frac{1}{2} \end{matrix} ; T \right) = \frac{\Gamma(a + b + n + \frac{1}{2})}{\Gamma(2a + n)\Gamma(b - a + \frac{1}{2})} \int_0^1 du \frac{u^{2a+n-1}}{[(1-u)(1-Tu)]^{a-b+1/2}}, \quad (3.31) $$

<footnote n="6">⁶http://dlmf.nist.gov/15.5.E21</footnote>

<footnote n="7">⁷http://dlmf.nist.gov/15.5.E18</footnote>

<footnote n="8">⁸https://dlmf.nist.gov/15.6.E1</footnote>
where Re(b − a) > − 1 / 2 is assumed, we can represent

$$
\begin{equation}
\begin{split}
\mathcal{K}_{\lambda}(t) = & \frac{(1-T)^{2a+2b} \cos(\pi(a-b)) \Gamma\left(a+b+\frac{1}{2}\right)^2}{\pi \Gamma(2a) \Gamma(2b)} \\
& \times \int_0^1 \int_0^1 \frac{1+T\lambda uv}{(1-T\lambda uv)^{2a+2b}} \frac{\left(\frac{(1-u)(1-uT)}{(1-v)(1-vT)}\right)^{b-a} u^{2a-1} v^{2b-1} du dv}{\sqrt{(1-u)(1-uT)(1-v)(1-vT)}}.
\end{split}
\tag{3.32}
\end{equation}
$$

Here $T = \tanh^2(t)$ and $-1/2 < \text{Re}(b-a) < 1/2$. The normalization condition $\mathcal{K}_1(t) = 1$
enables us to infer the integral

$$
\begin{align}
& \int_0^1 \int_0^1 \frac{1+Tuv}{(1-Tuv)^{2a+2b}} \frac{\left(\frac{(1-u)(1-uT)}{(1-v)(1-vT)}\right)^{b-a} u^{2a-1}v^{2b-1} du dv \nonumber \\
& \qquad = \frac{\pi \Gamma(2a) \Gamma(2b)}{(1-T)^{2a+2b} \cos(\pi(a-b)) \Gamma\left(a+b+\frac{1}{2}\right)^2}. \tag{3.33}
\end{align}
$$

The complexity can be obtained by differentiating eq. (3.32) and reads

$$
\begin{equation}
\begin{split}
\mathcal{K}(t) = & \frac{T(1-T)^{2a+2b} \cos(\pi(a-b)) \Gamma\left(a+c+\frac{1}{2}\right)^2}{\pi \Gamma(2a) \Gamma(2c)} \\
& \times \int_0^1 \int_0^1 \frac{2(a+b)(1+Tuv) + (1-Tuv)}{(1-Tuv)^{2a+2b+1}} \frac{\left(\frac{(1-u)(1-uT)}{(1-v)(1-vT)}\right)^{b-a} u^{2a} v^{2b} du dv}{\sqrt{(1-u)(1-uT)(1-v)(1-vT)}}.
\end{split}
\tag{3.34}
\end{equation}
$$

This formula reduces to $\mathcal{K}(t) = \eta \sinh^2(t)$ in the special case (1.7) [59].$^9$

Note that eq. (3.32) allows one to also address the generalized complexity [56] $K^{(\delta)} = (\lambda \partial_{\lambda})^{\delta} K_{\lambda}|_{\lambda=1}$ for positive integers $\delta$. At the same time, another widespread measure of complexity, the Krylov entropy [57], can not be obtained in an analogous way.

Let us compute the Krylov complexity (3.34) at late times $t \to +\infty$ for which it suffices to consider the leading order as $T \to 1$. In this case, it is obvious that the main contribution to the integrals comes from the regions $u \sim 1$ and $v \sim 1$ and therefore in the leading order we can approximate

$$
\begin{equation}
\begin{split}
\mathcal{K}(t) \approx & \frac{T(1-T)^{2a+2b} \cos(\pi(a-b)) \Gamma\left(a+b+\frac{1}{2}\right)^2}{\pi \Gamma(2a) \Gamma(2b)} 2(a+b) \\
& \times \int_0^1 \int_0^1 \frac{1+Tuv}{(1-Tuv)^{2a+2b+1}} \frac{\left(\frac{(1-u)(1-uT)}{(1-v)(1-vT)}\right)^{b-a} du dv}{\sqrt{(1-u)(1-uT)(1-v)(1-vT)}}.
\end{split}
\tag{3.35}
\end{equation}
$$

$$
\lim_{b \to a+1/2} \cos(\pi(b-a)) \int_0^1 \frac{f(v)}{(1-v)^{b-a+1/2}} dv = \pi f(1).
$$

<footnote>$^9$Notice that this case corresponds to $b = a + 1/2$, so one has to understand the complexity using the
limit</footnote>
Within the same accuracy, we can insert the appropriate powers of $u$ and $v$, to get the integral identical to eq. (3.33) up to shifting $a$ and $b$ by $1/4$. In this way as $t \to +\infty$ we obtain

$$ \mathcal{K}(t) \approx k_{\infty} e^{2t}, \qquad k_{\infty} = \frac{\Gamma(2a + \frac{1}{2}) \Gamma(2b + \frac{1}{2}) \Gamma(a + b + \frac{1}{2})^2}{2\Gamma(2a)\Gamma(2b)\Gamma(a+b)\Gamma(a+b+1)}. \qquad (3.36) $$

The exponential growth of the complexity is consistent with the universal operator growth [4]. Although the result (3.36) is derived under restricted values for $a$ and $b$, it turns out that the latter restriction is irrelevant for eq. (3.36) that remains valid for all admissible $a$ and $b$. We have confirmed this by a numerical analysis.

# 4 A family of models with odd-even alterations in the Lanczos sequence

In the previous section we have encountered examples of the oscillating behavior in the autocorrelation function. Now we will show that there exists another way to model such behavior of the autocorrelation function. It is based on an odd-even alteration in the subleading terms of $b_n$. Let us consider a very simple $C(t)$ of the form

$$ C(t) = \frac{\cos(\omega t)}{\cosh(t)}. \qquad (4.1) $$

We will demonstrate that the corresponding Lanczos coefficients that lead to eq. (4.1) are given by

$$ b_n = \begin{cases} n, & n \in 2\mathbb{Z} \\ \sqrt{n^2 + \omega^2}, & n \in 2\mathbb{Z} + 1 \end{cases}. \qquad (4.2) $$

Note that now we are now solving the inverse problem for eq. (1.1). For a given $\varphi_0(t) \equiv C(t)$ we seek the coefficients $b_n$ such that eq. (1.1) is satisfied subject to the initial condition (1.2). In principle, the solution can be obtained as follows. Once $C(t)$ is known, the weight for $\pi$-polynomials can be found via the Fourier transform (2.5). Knowing the weight, one can construct the orthogonalized set of $\pi$-polynomials that should satisfy the three-term recurrence relation (2.1) that in turn determines $b_n$'s. We will, however, approach this problem differently elaborating upon the method of moments [50].

It is convenient to introduce monic polynomials $P_n \equiv \pi_n \prod_{k=1}^n b_k$, so instead of the relation (2.1) we have

$$ xP_n(x) = P_{n+1}(x) + b_n^2 P_{n-1}(x). \qquad (4.3) $$

These polynomials are orthogonal with respect to the measure $\rho(x)$,

$$ \int_{-\infty}^{\infty} dx \rho(x) P_n(x) P_m(x) = h_n \delta_{n,m}, \quad h_n > 0. \qquad (4.4) $$

Let us recall how $\rho(x)$ (or $C(t)$) and $b_n^2$ are connected. Integrating both sides of eq. (4.3) multiplied by $\rho(x)P_{n-1}(x)$ leads to the connection $h_n = b_n^2 h_{n-1}$. On the other hand, we can define the partition function $Z_n$, which takes the form of a Hankel determinant

$$ Z_n = \frac{1}{n!} \int_{-\infty}^{\infty} dx_1 \rho(x_1) \int_{-\infty}^{\infty} dx_2 \rho(x_2) \cdots \int_{-\infty}^{\infty} dx_n \rho(x_n) \det_{1 \le i,j \le n} (x_i^{j-1})^2 = \det_{0 \le i,j \le n-1} \mu_{i+j}. \qquad (4.5) $$
Here the moments are given by

$$ \mu_k = \int_{-\infty}^{\infty} dx \rho(x) x^k = (i \partial_t)^k C(t) \Big|_{t=0}, \qquad (4.6) $$

see eq. (2.5) for the second equality. By making use of the identity

$$ \det_{1 \le i, j \le n} (x_i^{j-1}) = \det_{1 \le i, j \le n} (P_{j-1}(x_i)) \qquad (4.7) $$

that holds for monic polynomials (we can manipulate with the columns of the matrix on the left-hand side and form polynomials without changing the value of the determinant), we can express $Z_n$ also as

$$ Z_n = \frac{1}{n!} \int_{-\infty}^{\infty} dx_1 \rho(x_1) \int_{-\infty}^{\infty} dx_2 \rho(x_2) \cdots \int_{-\infty}^{\infty} dx_n \rho(x_n) \det_{1 \le i, j \le n} (P_{j-1}(x_i))^2 = \prod_{j=0}^{n-1} h_j, \quad (4.8) $$

where we have used the orthogonality (4.4). Therefore we obtain

$$ b_n^2 = \frac{Z_{n+1}Z_{n-1}}{Z_n^2}. \qquad (4.9) $$

Equation (4.9) should be understood as the connection between $b_n$'s and the moments $\mu_k$ encoded into the Hankel determinant (4.5).

Further simplification to the problem can be achieved if the measure is symmetric, $\rho(x) = \rho(-x)$. In this case following ref. [60] we can express

$$ Z_{2n} = A_n B_n, \qquad Z_{2n+1} = A_{n+1} B_n, \qquad (4.10) $$

where the *partial* partition functions are defined by

$$ A_n = \det_{0 \le i,j \le n-1} \mu_{2i+2j}, \qquad B_n = \det_{0 \le i,j \le n-1} \mu_{2i+2j+2}. \qquad (4.11) $$

It then follows

$$ b_{2n}^2 = \frac{A_{n+1}B_{n-1}}{A_nB_n}, \qquad b_{2n+1}^2 = \frac{A_nB_{n+1}}{A_{n+1}B_n}. \qquad (4.12) $$

Interestingly, the moments (4.6) corresponding to $C(t)$ of eq. (4.1) can be expressed in terms of the Euler polynomials $E_n(x)$, which are defined as

$$ \frac{2e^{xt}}{e^t + 1} = \sum_{n=0}^{\infty} E_n(x) \frac{t^n}{n!}. \qquad (4.13) $$

Indeed, we find $\mu_{2k+1} = 0$ and

$$ \mu_{2k} = (-4)^k E_{2k} \left( \frac{1+i\omega}{2} \right). \qquad (4.14) $$

Therefore,

$$ \det_{0 \le i,j \le n-1} \mu_{2i+2j+2\eta} = 4^{n(n-1)} (-4)^{n\eta} \det_{0 \le i,j \le n-1} E_{2i+2j+2\eta} \left( \frac{1+i\omega}{2} \right). \qquad (4.15) $$
The expressions for the type of Hankel determinants appearing in eq. (4.15) for $\eta = 0$ and $\eta = 1$ are available in the literature [61]. In particular, using corollary 5.2 of ref. [61] we obtain

$$ \det_{0 \le i,j \le n-1} \mu_{2i+2j+2\eta} = \left( (-4)^{\eta} E_{2\eta} \left( \frac{1+i\omega}{2} \right) \right)^n \times \prod_{\ell=1}^{n-1} c_{\eta,\ell}^{n-\ell}, \quad (4.16) $$

with $c_{\eta,\ell} = 4\ell^2[\omega^2 + (2\ell + 2\eta - 1)^2]$. Equation (4.16) enables us to conclude

$$ \frac{A_{k+1}}{A_k} = \prod_{\ell=1}^{k} c_{0,\ell}, \quad \frac{B_{k+1}}{B_k} = (1+\omega^2) \prod_{\ell=1}^{k} c_{1,\ell}. \quad (4.17) $$

We therefore obtain

$$ b_{2k}^2 = \frac{A_{k+1}B_{k-1}}{A_k B_k} = \frac{c_{0,1}}{1+\omega^2} \prod_{\ell=1}^{k-1} \frac{c_{0,\ell+1}}{c_{1,\ell}} = 4 \prod_{\ell=1}^{k-1} \frac{(\ell+1)^2}{\ell^2} = (2k)^2, \quad (4.18) $$

$$ b_{2k+1}^2 = \frac{A_k B_{k+1}}{A_{k+1} B_k} = (1+\omega^2) \prod_{\ell=1}^{k} \frac{c_{1,\ell}}{c_{0,\ell}} = (1+\omega^2) \prod_{\ell=1}^{k} \frac{\omega^2 + (2\ell+1)^2}{\omega^2 + (2\ell-1)^2} = \omega^2 + (2k+1)^2. \quad (4.19) $$

This finishes the proof of eq. (4.2).

# 5 Correlation functions relaxing to nonzero stationary values

## 5.1 Lanczos coefficients from stationary value

In all of the previously considered examples, the autocorrelation function decays to zero at late times. However, in general this is not the case – the late-time stationary value of the autocorrelation function can be nonzero. Here we explore the implications this incurs to the Lanczos coefficients.

Let us consider an autocorrelation function $C(t)$ that decays to zero. We study its deformation according to

$$ C^{(\kappa)}(t) = \kappa + (1 - \kappa)C(t), \quad 0 < \kappa < 1, \quad (5.1) $$

such that

$$ C^{(\kappa)}(t) \xrightarrow{t \to \infty} \kappa. \quad (5.2) $$

We aim to find the changes in the Lanczos coefficients produced by the $\kappa$-deformation (5.1) of the autocorrelation function.

Note that positive $\kappa$ is consistent with physical behavior of the autocorrelation function. Indeed,

$$ C^{(\kappa)}(t) \equiv \text{Tr}\{\mathcal{O}(t)\mathcal{O}\}/\text{Tr}\{\mathcal{O}^2\} = \sum_{E,E'} \langle E|\mathcal{O}|E'\rangle \langle E'|\mathcal{O}|E\rangle e^{-it(E'-E)}/\text{Tr}\{\mathcal{O}^2\}, \quad (5.3) $$

where $E, E'$ and $|E\rangle, |E'\rangle$ are the eigenenergies and the eigenvectors, respectively. If there are degeneracies in the spectrum, we choose the eigenvectors in such a way that $\mathcal{O}$ is
diagonal in invariant subspaces. All oscillating exponents are canceled with the help of the identity $\kappa = \lim_{t \to \infty} t^{-1} \int_0^t dt' C^{(\kappa)}(t')$. Note that the latter limit exists even in those pathological cases when the limit (5.2) does not exist, in which cases it can be regarded as a definition of $\kappa$. Combining the above equalities, we get $\kappa = \sum_E |\langle E | O | E \rangle|^2 / \mathrm{Tr}\{O^2\} \ge 0$.

Using eq. (4.6) one can immediately conclude that the moments corresponding to $C^{(\kappa)}(t)$ are given by

$$ \mu_0^{(\kappa)} = \mu_0 = 1, \quad \mu_{2k}^{(\kappa)} = (1-\kappa)\mu_{2k}, \quad k > 0, \qquad (5.4) $$

where $\mu_{2k}$ are the moments corresponding to $C(t)$. The corresponding deformations for the partial partition functions (4.11) are given by

$$ B_n^{(\kappa)} = (1-\kappa)^n B_n, \quad A_n^{(\kappa)} = (1-\kappa)^n A_n + \kappa(1-\kappa)^{n-1}C_{n-1}, \qquad (5.5) $$

where

$$ C_n = \det_{0 \le i,j \le n-1} \mu_{2i+2j+4}. \qquad (5.6) $$

Note that the determinant (5.6) has not been encountered before. As there is no obvious way to relate $C_n$ to $A_n$ and $B_n$ of eq. (4.11), it is thus not clear how to proceed further from this point following the recipe of the previous section.

In order to circumvent the latter obstacle, we use the moment-generating function instead, $G(z) = \sum_{k=0}^{\infty} \frac{\mu_k}{z^{k+1}}$. Since we assumed that $C(t)$ is symmetric and thus zero odd moments, $\mu_{2k+1} = 0$, we have

$$ G(z) = \frac{1}{z} + \sum_{k=1}^{\infty} \frac{\mu_{2k}}{z^{2k+1}}. \qquad (5.7) $$

A property of the generating function is that it encodes all the Lanczos coefficients $b_n$ once it is presented as a continued fraction [4, 50],

$$ G(z) = \frac{1}{z - \frac{b_1^2}{z - \frac{b_2^2}{z - \frac{b_3^2}{z - \dots}}}}. \qquad (5.8) $$

We now consider the $\kappa$-deformed case. Taking into account the changes of the moments described by eq. (5.4), we easily find that the deformation of the generating function is given by

$$ G^{(\kappa)}(z) = \frac{\kappa}{z} + (1-\kappa)G(z). \qquad (5.9) $$

Accounting for the continued fraction representation of $G^{(\kappa)}(z)$ similar to eq. (5.8), we can calculate the corresponding Lanczos coefficients. The final result reads

$$ \left( \frac{b_{2n}^{(\kappa)}}{b_{2n}} \right)^2 = \frac{1 + \frac{\kappa}{1-\kappa} \delta_{n+1}}{1 + \frac{\kappa}{1-\kappa} \delta_n}, \quad \left( \frac{b_{2n+1}^{(\kappa)}}{b_{2n+1}} \right)^2 = \frac{1 + \frac{\kappa}{1-\kappa} \delta_n}{1 + \frac{\kappa}{1-\kappa} \delta_{n+1}}, \qquad (5.10) $$
where $\delta_0 = 0$, $\delta_1 = 1$, $\delta_2 = 1 + b_1^2/b_2^2$, $\delta_3 = \delta_2 + b_1^2 b_3^2/b_2^2 b_4^2$, etc. The general expression is given by

$$ \delta_n = 1 + \sum_{k=1}^{n-1} y_k, \quad y_k = \prod_{j=1}^{k} \left( \frac{b_{2j-1}}{b_{2j}} \right)^2, \qquad (5.11) $$

where $n = 2, 3, 4, \dots$. Equations (5.10) and (5.11) determine the Lanczos coefficients for the $\kappa$-deformed case characterized by the autocorrelation function (5.1).

It is instructive to derive an asymptotic form of eq. (5.10) in the case the Lanczos coefficients satisfy the universal operator growth hypothesis (1.4). As long as $b_n = n+O(1)$, we obtain $y_k = c_0/k$ at the leading order, where the coefficient $c_0$ depends on the first subleading term in $b_n$ and thus $\delta_n = c_0 \ln n + O(1)$. In this way we obtain

$$ b_{2n}^{(\kappa)} = b_{2n} + \frac{1}{\ln(g(\kappa)n)} + \dots, \qquad (5.12) $$

$$ b_{2n+1}^{(\kappa)} = b_{2n+1} - \frac{1}{\ln(g(\kappa)n)} + \dots, \qquad (5.13) $$

where $g(\kappa)$ depends on the subleading terms and the ellipsis denotes further subleading terms. One can see that the nonzero stationary value of the autocorrelation function is ensured by the odd-even alterations in the Lanczos coefficients that decay as the inverse logarithm of $n$, which is slower than any power law. While it has been known that a nonzero asymptotic value of the correlation function is associated with an alternating term in the Lanczos coefficients [17, 32, 38, 42–45], the precise scaling of the alternating term has not been derived. Note that the inverse logarithmic scaling has been considered among other possibilities in refs. [38, 45].

It is interesting that the *shifted* determinant $C_n$ in eq. (5.6) can be expressed via $A_n$ of eq. (4.11). Indeed, taking into account eqs. (5.5) and (4.12) we obtain

$$ \left( \frac{b_{2n}^{(\kappa)}}{b_{2n}} \right)^2 = \frac{1 + \frac{\kappa}{1-\kappa} \frac{C_n}{A_{n+1}}}{1 + \frac{\kappa}{1-\kappa} \frac{C_{n-1}}{A_n}}. \qquad (5.14) $$

Comparing this with eq. (5.10) we conclude that

$$ C_n = A_{n+1}\delta_{n+1}. \qquad (5.15) $$

Let us study an example. For $b_n$ of eq. (4.2) from the previous section, we obtain

$$ y_k = \frac{\Gamma\left(k + \frac{1}{2} - \frac{i\omega}{2}\right) \Gamma\left(k + \frac{1}{2} + \frac{i\omega}{2}\right)}{(k!)^2 \Gamma\left(\frac{1}{2} - \frac{i\omega}{2}\right) \Gamma\left(\frac{1}{2} + \frac{i\omega}{2}\right)}. \qquad (5.16) $$

Accounting for the determinant (4.15), eq. (5.15) then leads to

$$ \frac{\det_{0 \le i,j \le n-1} E_{2i+2j+4}(z)}{\det_{0 \le i,j \le n} E_{2i+2j}(z)} = 1 + \frac{\sin(\pi z)}{\pi} \sum_{k=1}^{n} \frac{\Gamma(k+1-z)\Gamma(k+z)}{(k!)^2}, \qquad (5.17) $$

where we have introduced $z = (1 + i\omega)/2$. Note, however, that eq. (5.17) is valid at any complex number $z$. We should also note that the two Hankel determinants of Euler
polynomials in the left-hand side of eq. (5.17) have different dimensions. The one known in the literature [61] is given by

$$ \det_{0 \le i,j \le n} E_{2i+2j}(z) = (-1)^{\frac{n(n+1)}{2}} \prod_{\ell=0}^{n} ((\ell+1)^2(z+\ell)(z-\ell-1))^{n-\ell}, \quad (5.18) $$

and the other one, $\det_{0 \le i,j \le n-1} E_{2i+2j+4}(z)$, is calculated here in eq. (5.17).

The problem studied in this section has yet another interpretation. We assumed that the coefficients $b_n$ correspond to the autocorrelation function $C(t)$, which is related to the polynomials that satisfy the three-term recurrence relation (4.3) and the orthogonalization condition (4.4) with $\rho(x)$ as the weight. The obtained coefficients $b_n^{(\kappa)}$ also correspond to a set of orthogonal polynomials with the three-term recurrence relation and the orthogonalization condition, but with respect to the weight

$$ \rho^{(\kappa)}(x) = (1 - \kappa)\rho(x) + \kappa\delta(x). \quad (5.19) $$

We have therefore found the set of coefficients $b_n^{(\kappa)}$ for the weight (5.19) knowing the set of coefficients $b_n$ for the weight $\rho(x)$.

## 5.2 The stationary value from the Lanczos coefficients

Let us consider the reversed problem that consists of finding the stationary, late-time, value $\kappa$ of the autocorrelation function $C^{(\kappa)}(t)$ provided the Lanczos coefficient $b_n^{(\kappa)}$ are known. To achieve that let us consider

$$ y_k^{(\kappa)} = \prod_{j=1}^{k} \left( \frac{b_{2j-1}^{(\kappa)}}{b_{2j}^{(\kappa)}} \right)^2 . \quad (5.20) $$

Using eq. (5.10) we find the connection

$$ y_k = (1 - \kappa) \left( 1 + \frac{\kappa}{1 - \kappa} \delta_k \right) \left( 1 + \frac{\kappa}{1 - \kappa} \delta_{k+1} \right) y_k^{(\kappa)}. \quad (5.21) $$

On the other hand, the left-hand side can be expressed as $y_k = \delta_{k+1} - \delta_k$, as follows from eq. (5.11). We thus obtain

$$ y_k^{(\kappa)} = \frac{1}{\kappa} \left( \frac{1}{1 + \frac{\kappa}{1-\kappa} \delta_k} - \frac{1}{1 + \frac{\kappa}{1-\kappa} \delta_{k+1}} \right). \quad (5.22) $$

For $\delta_k$ that diverges logarithmically with $k$ (or any other $\delta_k$ that tends to infinity as $k \to \infty$), we therefore obtain

$$ \kappa = 1/\delta_{\infty}^{(\kappa)} = \left( 1 + \sum_{k=1}^{\infty} y_k^{(\kappa)} \right)^{-1}. \quad (5.23) $$

In the case of a divergent sum in the right-hand side of eq. (5.23), we obtain $\kappa = 0$. This is consistent with the scaling $y_k \sim 1/k$ obtained previously for the undeformed $C(t)$ that decays to zero.
We find it instructive to also outline a wrong derivation of the stationary value from the Lanczos coefficients. Equation (1.1) involves the sum rule

$$ \sum_{n=0}^{\infty} \varphi_n(t)^2 = 1. \qquad (5.24) $$

If the convergence of $\varphi_n(t)$ to the corresponding stationary values $\varphi_n(\infty)$ (where $\varphi_0(\infty) \equiv \kappa$) were uniform, then the exchange of the limit $t \to \infty$ and the sum would give the sum rule $\sum_{n=0}^{\infty} \varphi_n(\infty)^2 = 1$. At the same time, eq. (1.1) implies $\varphi_{2l+1}(\infty) = 0$ and $\varphi_{2n}(\infty) = (b_{2n-1}^{(\kappa)}/b_{2n}^{(\kappa)})^2 \varphi_{2n-2}(\infty)$, which leads to $\varphi_{2n}^2 = \varphi_0^2 y_n^{(\kappa)}$. Substituting this into the sum rule leads to $\kappa = (1 + \sum_{k=1}^{\infty} y_k^{(\kappa)})^{-1/2}$, which differs from the correct expression (5.23).

This apparent paradox is resolved by noting that the convergence of $\varphi_n(t)$ is not, in general, uniform: at any fixed $t$ an infinite number of $\varphi_n(t)$ with sufficiently large $n$ are not anything close to their stationary values $\varphi_n(\infty)$. Therefore, in general, we cannot exchange the limit $t \to \infty$ and the infinite summation (that contains another limit). One can, however, derive a relaxed version of the sum rule by taking the limit of infinite time in the inequality that involves a sum of finite number of terms, $\sum_{n=0}^{p} \varphi_n(t)^2 \le 1$, where $p$ is finite. This way one obtains

$$ \sum_{n=0}^{\infty} \varphi_n(\infty)^2 \le 1. \qquad (5.25) $$

One can easily verify that this inequality does not lead to any contradiction with eq. (5.23).

The sum rule (5.24) is widely used [8, 45, 57]. The above apparent paradox and its resolution highlights that it should be treated with caution at asymptotically large times.

# 6 Discussion

Our findings highlight the pivotal role of the subleading terms of the Lanczos coefficients $b_n$ at large $n$. In the exactly solvable two-parameter family (3.1) one can separately control two subleading terms proportional to $O(1)$ and $O(1/n)$ by tuning the parameters $a$ and $b$. Asymptotic analysis of the exact solution (see eqs. (3.18) and (3.22)) shows that even the most rough feature of the correlation function – its decay exponent – depends on both subleading terms. Furthermore, subleading terms determine whether the correlation function features damped oscillations or damping without oscillations at large times, as illustrated in Fig. 1. In contrast, the Krylov complexity turns out to be largely insensitive to the subleading terms, see eq. (3.36).

The second family we studied, see eq. (4.2), features an alternating subleading term $O((-1)^n/n)$. This term provides a different pathway to damped oscillations in the auto-correlation function.

Finally, we have shown how Lanczos coefficients should be modified to obtain a deformation of the correlation function with a nonzero stationary value.
Importantly, the strong effect of subleading terms calls for caution when using techniques that rely on approximation and/or extrapolation of Lanczos coefficients. This includes various extrapolated versions of the recursion method [4, 7, 17, 21, 62], approximating Lanczos coefficients by sampling [18] and continuum approximation [4, 8, 44, 45, 57].

The latter technique deserves a separate remark. In the continuum approximation, one treats the Lanczos coefficients $b_n = b(n)$ as a continuous function of $n$ [4, 8, 57]. In the presence of odd-even alterations, one treats $b_{2n-1}$ and $b_{2n}$ as two different continuous functions [45]. By performing a gradient expansion of $b(n)$, one turns a system (1.1) of ordinary equations to a single partial differential equation on the function $\varphi(n, t)$. Usually, subleading terms of the gradient expansion of $b(n)$ are disregarded [4, 8, 45, 57]. Our findings show that this approximation can be unacceptably rough and miss even the basic features of the autocorrelation function. For example, applying the results of the continuum approximation of ref. [45] to the exactly solvable model (4.2), one gets $C(t) \propto e^{-t}$ at large times. This estimate misses the oscillating prefactor present in the actual large-time approximation $C(t) \simeq 2 \cos(\omega t) e^{-t}$, see eq. (4.1).

## A Continuous Hahn polynomials

The continuous Hahn polynomials are defined by the relation [58, 63]

$$p_n(x) = i^n \frac{(a+c)_n (a+d)_n}{n!} {}_3F_2 \left( \begin{array}{c} -n, \ n+a+b+c+d-1, \ a+ix \\ a+c, \ a+d \end{array} ; 1 \right). \quad (A.1)$$

They depend on four parameters, $a, b, c$, and $d$. We consider the parameters that obey

$$\mathrm{Re}(a, b, c, d) > 0, \quad c = a^*, \quad d = b^*. \qquad (A.2)$$

In this case the orthogonality relation of polynomials is given by the integral over the real axis of the form

$$\begin{aligned} & \frac{1}{2\pi} \int_{-\infty}^{\infty} dx \Gamma(a+ix) \Gamma(b+ix) \Gamma(c-ix) \Gamma(d-ix) p_m(x) p_n(x) \\ & = \frac{\Gamma(n+a+c)\Gamma(n+a+d)\Gamma(n+b+c)\Gamma(n+b+d)}{(2n+a+b+c+d-1)\Gamma(n+a+b+c+d-1)n!} \delta_{m,n}. \end{aligned} \qquad (A.3)$$

Using the Barnes integral

$$\frac{1}{2\pi} \int_{-\infty}^{\infty} dx \Gamma(a+ix) \Gamma(b+ix) \Gamma(c-ix) \Gamma(d-ix) = \frac{\Gamma(a+c)\Gamma(a+d)\Gamma(b+c)\Gamma(b+d)}{\Gamma(a+b+c+d)}, \qquad (A.4)$$

valid at $\mathrm{Re}(a, b, c, d) > 0$, we can introduce the weight as

$$w(x) = \frac{\Gamma(a+b+a^*+b^*)}{2\pi\Gamma(a+a^*)\Gamma(a+b^*)\Gamma(a^*+b)\Gamma(b+b^*)} |\Gamma(a+ix)\Gamma(b+ix)|^2, \qquad (A.5)$$
such that it is normalized, $\int_{-\infty}^{\infty} dxw(x) = 1$. The normalization of the polynomials then takes the form

$$ \int_{-\infty}^{\infty} dxw(x)p_m(x)p_n(x) = \frac{(a+c)_n(a+d)_n(b+c)_n(b+d)_n}{(a+b+c+d)_n} \frac{n+a+b+c+d-1}{n!(2n+a+b+c+d-1)} \delta_{m,n}. \quad (A.6) $$

The recurrence relation is

$$ p_{n+1}(x) = \frac{B_n}{A_n}[-x + i(A_n + C_n + a)]p_n(x) + \frac{B_{n-1}B_nC_n}{A_n}p_{n-1}(x), \quad (A.7) $$

where

$$ A_n = \frac{(n+a+b+c+d-1)(n+a+c)(n+a+d)}{(2n+a+b+c+d-1)(2n+a+b+c+d)}, \quad (A.8) $$

$$ B_n = \frac{(n+a+c)(n+a+d)}{n+1}, \quad (A.9) $$

$$ C_n = \frac{n(n+b+c-1)(n+b+d-1)}{(2n+a+b+c+d-2)(2n+a+b+c+d-1)}. \quad (A.10) $$

## A.1 Fourier transform

Consider the Fourier transform involving $p_n(x)$ of the form

$$ W_n(t) = \int_{-\infty}^{\infty} dx e^{-ixt} w(x) p_n(x). \quad (A.11) $$

Let us consider some general properties of this relation. (i) From the orthogonality of the polynomials it follows $W_n(0) = \delta_{n,0}$. (ii) From the definition it follows $W_{-1}(t) = 0$ and $W_0(t) = \int_{-\infty}^{\infty} dx e^{-ixt} w(x)$. Once all $W_j(t)$ are known for $j \le n$, the term $W_{n+1}(t)$ can be found from the recurrence relation (A.7) and the equality

$$ \int_{-\infty}^{\infty} dx e^{-ixt} w(x) x p_n(x) = i \frac{\partial W_n(t)}{\partial t}. \quad (A.12) $$

It yields

$$ W_{n+1}(t) = i \frac{B_n}{A_n} \left[ a + A_n + C_n - \frac{\partial}{\partial t} \right] W_n(t) + \frac{B_{n-1}B_nC_n}{A_n} W_{n-1}(t). \quad (A.13) $$

Equation (A.13) is a convenient way to obtain recursively the expressions for $W_{n+1}(t)$ from the preceding two.

## A.2 Transformation to the form of eq. (1.1)

Let us transform eq. (A.13) to the form

$$ \frac{\partial \phi_n(t)}{\partial t} = -\tilde{b}_{n+1} \phi_{n+1}(t) + \tilde{a}_n \phi_n(t) + \tilde{b}_n \phi_{n-1}(t). \quad (A.14) $$
Introducing

$$W_n(t) = (-i)^n \alpha_n \phi_n(t), \qquad (A.15)$$

we obtain the coefficients

$$\tilde{a}_n = a + A_n + C_n, \quad \tilde{b}_n = B_{n-1}C_n \frac{\alpha_{n-1}}{\alpha_n}, \qquad (A.16)$$

provided

$$\left( \frac{\alpha_{n+1}}{\alpha_n} \right)^2 = -B_n^2 \frac{C_{n+1}}{A_n} \qquad (A.17)$$

is satisfied for all nonnegative integers $n$. Since $A_n < 0$ and $B_n, C_{n+1} > 0$ at $n \ge 0$, the solution can be taken as

$$\alpha_0 = 1, \quad \alpha_n = \prod_{j=0}^{n-1} B_j \sqrt{\frac{C_{j+1}}{-A_j}}, \quad n \ge 1. \qquad (A.18)$$

Then we find

$$\tilde{b}_n = \sqrt{-A_{n-1}C_n}. \qquad (A.19)$$

The functions $\phi_n(t)$ can be generated using the recurrence relation

$$\phi_{n+1}(t) = \frac{1}{\tilde{b}_{n+1}} \left( \tilde{a}_n - \frac{\partial}{\partial t} \right) \phi_n(t) + \frac{\tilde{b}_n}{\tilde{b}_{n+1}} \phi_{n-1}(t), \quad n \ge 0. \qquad (A.20)$$

Here we should use

$$\phi_{-1}(t) = 0, \quad \phi_0(t) = \int_{-\infty}^{\infty} dx e^{-itx} w(x). \qquad (A.21)$$

### A.3 The case $\tilde{a}_n = 0$

In the following we consider the case $\tilde{a}_n = 0$, which occurs if $A_n + C_n = -a$. For real $a, b, c$, and $d$ the admissible set is the one with

$$a = c > 0, \quad b = d > 0. \qquad (A.22)$$

Another possibility with $\tilde{a}_n = 0$ involves the complex parameters and we set $c = a^*, d = b^*$. Direct inspection shows that $\tilde{a}_n = 0$ occurs if $\mathrm{Im}(a+b) = 0$ and $\mathrm{Re}(a-b) \times \mathrm{Re}(a+b-1) = 0$. Therefore two cases arise. One is

$$a = d = r + i\omega, \quad b = c = r - i\omega, \quad r > 0, \quad \omega \neq 0. \qquad (A.23)$$

The other is

$$a = r + i\omega, \quad b = 1 - r - i\omega, \quad c = r - i\omega, \quad d = 1 - r + i\omega, \quad 0 < r < 1. \qquad (A.24)$$
### A.4 Evaluation of $\phi_n(t)$

Let us evaluate $\phi_n(t)$ for $\tilde{a}_n = 0$. We will use the integral representation of continuous Hahn polynomials [64]

$$p_n(z) = \frac{i^n}{2^{a+b^*-1}} \frac{\Gamma(a+b^*+n)}{\Gamma(a+iz)\Gamma(b^*-iz)} \times \int_{-\infty}^{\infty} dx e^{-i2xz} (1-\tanh x)^a (1+\tanh x)^{b^*} P_n^{(a+a^*-1,b+b^*-1)}(\tanh x), \quad (A.25)$$

where $P_n^{(\alpha,\beta)}(x)$ are the Jacobi polynomials. The latter are defined by

$$P_n^{(\alpha,\beta)}(x) = \frac{\Gamma(\alpha + n + 1)}{n! \Gamma(\alpha + \beta + n + 1)} \sum_{k=0}^{n} (-1)^k \binom{n}{k} \frac{\Gamma(\alpha + \beta + n + k + 1)}{\Gamma(\alpha + k + 1)} \left(\frac{1-x}{2}\right)^k . \quad (A.26)$$

We consider $\mathrm{Re}(\alpha, \beta) > -1$, which is fulfilled for $\mathrm{Re}(a, b) > 0$ that is assumed in eq. (A.25). Using

$$\int_{-\infty}^{\infty} dz e^{-ixz} \Gamma(b+iz) \Gamma(a^*-iz) = \frac{2\pi \Gamma(a^*+b)}{(1-e^{-x})^b (1+e^x)^{a^*}}, \quad (A.27)$$

we obtain

$$W_n(t) = i^n (a+b^*)_n \frac{\Gamma(a+b+a^*+b^*)}{\Gamma(a+a^*)\Gamma(b+b^*)} e^{bt} \\ \times \int_0^1 dx \frac{(1-x)^{a+a^*-1} x^{b+b^*-1}}{[1-(1-e^t)x]^{a^*+b}} P_n^{(a+a^*-1,b+b^*-1)}(2x-1). \quad (A.28)$$

Using eq. (A.26) and the integral representation<fnref n="10" />

$$F\left(\begin{matrix} a, & b \\ c & ; z \end{matrix}\right) = \frac{\Gamma(c)}{\Gamma(b)\Gamma(c-b)} \int_{0}^{1} dx \frac{x^{b-1}(1-x)^{c-b-1}}{(1-zx)^a}, \quad \mathrm{Re}(c) > \mathrm{Re}(b) > 0, \quad (A.29)$$

we obtain

$$W_n = \frac{i^n}{n!} (a+b^*)_n (a+a^*)_n e^{bt} \\ \times \sum_{k=0}^{n} (-1)^k \binom{n}{k} \frac{(n-1+a+b+a^*+b^*)_k}{(a+b+a^*+b^*)_k} F\left(\begin{matrix} a^*+b, & b+b^* \\ a+b+a^*+b^*+k & 1-e^t \end{matrix}\right). \quad (A.30)$$

The initial term $\phi_0(t)$ of eq. (A.21) is identical to $W_0(t)$, see eq. (A.11). Therefore from eq. (A.30) we obtain

$$\phi_0(t) = e^{bt} F\left(\begin{matrix} a^* + b, & b + b^* \\ a + b + a^* + b^* & 1 - e^t \end{matrix}\right). \quad (A.31)$$

Here $F(a, b; z)$ is the Gauss hypergeometric function defined in eq. (3.4). Equation (A.31) gives complex $\phi_0(t)$ for the case (A.24). On the other hand, we want to study real $\phi_0(t)$,

<footnote n="10">¹⁰https://dlmf.nist.gov/15.6.E1</footnote>
which occurs in the cases (A.22) and (A.23). In the following we thus consider the latter two sets of parameters. Noting that in this case we have

$$ \alpha_n = \sqrt{\frac{(2a)_n (2b)_n [(a+b)_n]^2 (2a+2b+n-1)_n}{n!(2a+2b)_{2n}}}, \quad (A.32) $$

from eqs. (A.15) and (A.30) we obtain

$$ \phi_n(t) = (-1)^n \sqrt{\frac{(2a+2b)_{2n}(2a)_n}{n!(2b)_n(2a+2b-1+n)_n}} e^{bt} \\ \times \sum_{k=0}^{n} (-1)^k \binom{n}{k} \frac{(2a+2b-1+n)_k}{(2a+2b)_k} F\left(\begin{array}{l} a+b, 2b \\ 2a+2b+k \end{array}; 1-e^t\right). \quad (A.33) $$

Equation (A.33) contains a finite sum that we were not able to calculate directly. Instead of that, we found an alternative route to calculate $\phi_n(t)$. The case $n = 0$ of eq. (A.33) can be expressed as<fnref n="11" />

$$ \phi_0(t) = F\left(\begin{array}{l} a, b \\ a+b+\frac{1}{2} \end{array}; -\sinh^2(t/2)\right). \quad (A.34) $$

Using the recurrence relation (A.20), we then calculated $\phi_1(t)$, $\phi_2(t)$, $\phi_3(t)$, etc., and concluded that they satisfy the general form

$$ \phi_n(t) = \left( \prod_{j=1}^{n} \tilde{b}_j \right) \frac{2^n}{n!} [\sinh(t/2)]^n F \left( \begin{array}{l} a+n/2, b+n/2 \\ a+b+\frac{1}{2}+n \end{array} ; -\sinh^2(t/2) \right). \quad (A.35) $$

This is our final expression for $\phi_n(t)$.

The properties of the Gauss hypergeometric function enabled us to transform eq. (A.35) to several other forms. Some of them are given by

$$ \phi_n(t) = \left( \prod_{j=1}^{n} \tilde{b}_j \right) \frac{2^n}{n!} [\sinh(t/2)]^n F \left( \begin{array}{l} 2a+n, 2b+n \\ a+b+\frac{1}{2}+n \end{array} ; -\sinh^2(t/4) \right) \quad (A.36) $$

and

$$ \phi_n(t) = \left( \prod_{j=1}^{n} \tilde{b}_j \right) \frac{1}{n!} e^{bt} (e^t - 1)^n F \left( \begin{array}{l} a+b+n, 2b+n \\ 2a+2b+2n \end{array} ; 1-e^t \right). \quad (A.37) $$

In the main text, two other forms are given. Note that there are several forms of the prefactor that have some similarities with the structure of arguments in the hypergeometric

<footnote n="11">¹¹https://functions.wolfram.com/07.23.17.0101.01</footnote>
functions. Some forms of the prefactor are

$$
\begin{align*}
\prod_{j=1}^{n} \tilde{b}_j &= \left( \frac{n!}{4^{2n}} \frac{(2a)_n (2b)_n (2a + 2b - 1)_n}{(a + b - 1/2)_n (a + b + 1/2)_n} \right)^{1/2} \\
&= \left( \frac{n!(2a)_n (2b)_n (a+b)_n}{2^{2n} (a+b+1/2)_n (2a+2b+n-1)_n} \right)^{1/2} \\
&= \left( \frac{n!(2a)_n (2b)_n [(a+b)_n]^2}{(2a+2b)_{2n} (2a+2b+n-1)_n} \right)^{1/2}. \tag{A.38}
\end{align*}
$$

Using the Clausen identity

$$
\left[ F \left( \begin{array}{c} a, b \\ a+b+\frac{1}{2} \end{array} ; z \right) \right]^2 = {}_3F_2 \left( \begin{array}{c} 2a, 2b, a+b \\ a+b+\frac{1}{2}, 2a+2b \end{array} ; z \right), \quad (A.39)
$$

we note another interesting representation

$$
\phi_n(t) = \left( \frac{(2a)_n (2b)_n (a+b)_n}{n!(a+b+1/2)_n (2a+2b+n-1)_n} \right)^{1/2} [\sinh(t/2)]^n \\
\times \left[ {}_3F_2 \left( \begin{array}{c} 2a+n, 2b+n, a+b+n \\ a+b+n+\frac{1}{2}, 2a+2b+2n \end{array} ; -\sinh^2(t/2) \right) \right]^{1/2} . \quad (A.40)
$$

The proof that $\phi_n(t)$ indeed satisfies eq. (A.20) is given in the main text.

As a side result of previous considerations, a comparison of eq. (A.33) and the expres-
sion (A.37) leads to the identity

$$
\begin{align}
& \sum_{k=0}^{n} (-1)^k \binom{n}{k} \frac{(2a + 2b - 1 + n)_k}{(2a + 2b)_k} F \left( \begin{matrix} a + b, 2b \\ 2a + 2b + k \end{matrix} ; z \right) \nonumber \\
& = \frac{(a + b)_n (2b)_n}{(2a + 2b)_{2n}} z^n F \left( \begin{matrix} a + b + n, 2b + n \\ 2a + 2b + 2n \end{matrix} ; z \right). \tag{A.41}
\end{align}
$$

The sum in the left-hand side of eq. (A.41) can be found in the literature, see the expression 5.3.5.3 in ref. [65]. It is however expressed in terms of a $_3F_2$ hypergeometric function. After comparing with the right-hand side of eq. (A.41) we obtain the equality

$$
\frac{(a+b)_n (2b)_n}{(2a+2b+n)_n} z^n F \left( \begin{matrix} a+b+n, 2b+n \\ 2a+2b+2n \end{matrix} ; z \right) = (1-n)_n \times {}_3F_2 \left( \begin{matrix} 1, a+b, 2b \\ 1-n, 2a+2b+n \end{matrix} ; z \right). \quad (A.42)
$$

The right-hand side of eq. (A.42) is not defined at positive integers $n$, but can be understood
as

$$
(1-n)_n \times {}_3F_2 \left( \begin{matrix} 1, a+b, 2b \\ 1-n, 2a+2b+n \end{matrix} ; z \right) = \sum_{k=n}^{\infty} \frac{(a+b)_k (2b)_k}{(2a+2b+n)_k (k-n)!} \frac{z^k}{(k-n)!} \quad (A.43)
$$

After shifting the index of summation we then obtain the left-hand side of eq. (A.42).
Therefore both derivations of $\phi_n(t)$ are consistent.

Let us verify the sum rule (5.24). Using the representation (A.37) we can perform the summation using the expression 6.7.2.3 from ref. [65]. Note that the sum rule can also be verified using eq. (A.40) and the expression 6.8.1.31 from ref. [65].

In the notation of the main text, we have used the rescaled time such that $\phi_n(t) = \phi_n(4t)$ and $b_n = 4\tilde{b}_n$.
Acknowledgments

OL thanks A. Teretenkov and N. Il'in for useful discussions. Work at Laboratoire de Physique Théorique was supported in part by the EUR grant NanoX ANR-17-EURE-0009 in the framework of the “Programme des Investissements d'Avenir”.

References

[1] H. Liu and J. Sonner, *Holographic systems far from equilibrium: a review*, *Rep. Prog. Phys.* **83** (2019) 016001.

[2] M.C. Bañuls, M.B. Hastings, F. Verstraete and J.I. Cirac, *Matrix product states for dynamical simulation of infinite chains*, *Phys. Rev. Lett.* **102** (2009) 240603.

[3] R. Orús, *Tensor networks for complex quantum systems*, *Nat. Phys. Rev. 1* (2019) 538.

[4] D.E. Parker, X. Cao, A. Avdoshkin, T. Scaffidi and E. Altman, *A Universal Operator Growth Hypothesis*, *Phys. Rev. X* **9** (2019) 041017.

[5] H. Mori, *A continued-fraction representation of the time-correlation functions*, *Prog. Theor. Phys.* **34** (1965) 399.

[6] M. Dupuis, *Moment and continued fraction expansions of time autocorrelation functions*, *Prog. Theor. Phys.* **37** (1967) 502.

[7] V.S. Viswanath and G. Müller, *The Recursion Method: Application to Many-Body Dynamics*, vol. 23, Springer, Berlin (2008), 10.1007/978-3-540-48651-0.

[8] W. Mück and Y. Yang, *Krylov complexity and orthogonal polynomials*, *Nucl. Phys. B* **984** (2022) 115948.

[9] J.-M. Liu and G. Müller, *Infinite-temperature dynamics of the equivalent-neighbor XYZ model*, *Phys. Rev. A* **42** (1990) 5854.

[10] J. Florencio, S. Sen and Z. Cai, *Quantum spin dynamics of the transverse Ising model in two dimensions*, *J. Low Temp. Phys.* **89** (1992) 561.

[11] V. Zobov and A. Lundin, *Second moment of multiple-quantum NMR and a time-dependent growth of the number of multispin correlations in solids*, *J. Exp. Theor. Phys.* **103** (2006) 904.

[12] T.A. Elsayed, B. Hess and B.V. Fine, *Signatures of chaos in time series generated by many-spin systems at high temperatures*, *Phys. Rev. E* **90** (2014) 022910.

[13] G. Bouch, *Complex-time singularity and locality estimates for quantum lattice systems*, *J. Math. Phys.* **56** (2015) 123303.

[14] X. Cao, *A statistical mechanism for operator growth*, *J. Phys. A: Math. Theor.* **54** (2021) 144001.

[15] J.D. Noh, *Operator growth in the transverse-field Ising spin chain with integrability-breaking longitudinal field*, *Phys. Rev. E* **104** (2021) 034112.

[16] R. Heveling, J. Wang and J. Gemmer, *Numerically probing the universal operator growth hypothesis*, *Phys. Rev. E* **106** (2022) 014152.

[17] F. Uskov and O. Lychkovskiy, *Quantum dynamics in one and two dimensions via the recursion method*, *Phys. Rev. B* **109** (2024) L140301.
[18] A. De, U. Borla, X. Cao and S. Gazit, Stochastic sampling of operator growth dynamics, Phys. Rev. B 110 (2024) 155135.

[19] N. Loizeau, J.C. Peacock and D. Sels, *Quantum many-body simulations with PauliStrings.jl*, SciPost Phys. Codebases (2025) 54.

[20] N. Loizeau, J.C. Peacock and D. Sels, Codebase release 1.5 for PauliStrings.jl, SciPost Phys. Codebases (2025) 54.

[21] A. Teretenkov, F. Uskov and O. Lychkovskiy, *Pseudomode expansion of many-body correlation functions*, Phys. Rev. B 111 (2025) 174308.

[22] I. Shirokov, V. Hrushev, F. Uskov, I. Dudinets, I. Ermakov and O. Lychkovskiy, *Quench dynamics via recursion method and dynamical quantum phase transitions*, arXiv:2503.24362 (2025).

[23] I. Ermakov, *Operator growth in many-body systems of higher spins*, arXiv:2504.07833 (2025).

[24] P. Nandy, A.S. Matsoukas-Roubeas, P. Martínez-Azcona, A. Dymarsky and A. del Campo, *Quantum dynamics in Krylov space: Methods and applications*, arXiv:2405.09628 (2024).

[25] S. Baiguera, V. Balasubramanian, P. Caputa, S. Chapman, J. Haferkamp, M. P. Heller, and N. Y. Halpern, “Quantum complexity in gravity, quantum field theory, and quantum information science,” arXiv:2503.10753

[26] P. Caputa, J.M. Magan and D. Patramanis, *Geometry of Krylov complexity*, Phys. Rev. Research 4 (2022) 013041.

[27] J. Kim, J. Murugan, J. Olle and D. Rosa, *Operator delocalization in quantum networks*, Phys. Rev. A 105 (2022) L010201.

[28] K. Adhikari and S. Choudhury, *Cosmological Krylov complexity*, Fortschr. Phys. 70 (2022) 2200126.

[29] S.-K. Jian, B. Swingle and Z.-Y. Xian, *Complexity growth of operators in the SYK model and in JT gravity*, J. High Energy Phys. 2021 (2021) 1.

[30] A. Kar, L. Lamprou, M. Rozali et al., *Random matrix theory for complexity growth and black hole interiors*, J. High Energy Phys. 2022 (2022) 16.

[31] K. Adhikari, S. Choudhury and A. Roy, *Krylov complexity in quantum field theory*, Nucl. Phys. B 993 (2023) 116263.

[32] A. Avdoshkin, A. Dymarsky and M. Smolkin, *Krylov complexity in quantum field theory, and beyond*, J. High Energy Phys. 2024 (2024) 66.

[33] S. Chapman, S. Demulder, D.A. Galante, S.U. Sheorey and O. Shoval, *Krylov complexity and chaos in deformed Sachdev-Ye-Kitaev models*, Phys. Rev. B 111 (2025) 035141.

[34] P. Caputa, H.-S. Jeong, S. Liu, J. F. Pedraza, and L.-C. Qu, “Krylov complexity of density matrix operators,” J. High Energy Phys., 2024, 05, 337, 2024.

[35] H. A. Camargo, V. Jahnke, H.-S. Jeong, K.-Y. Kim, and M. Nishida, “Spectral and Krylov complexity in billiard systems,” Phys. Rev. D, 109, 4, 046017, 2024.

[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, *J. High Energy Phys.* **2021** (2021) 188.

[38] M. Dodelson, *Black holes from chaos*, arXiv:2501.06170 (2025).

[39] L. Banchi and R. Vaia, *Spectral problem for quasi-uniform nearest-neighbor chains*, *J. Math. Phys.* **54** (2013) 043501.

[40] O. Gamayun, O. Lychkovskiy and J.-S. Caux, *Fredholm determinants, full counting statistics and Loschmidt echo for domain wall profiles in one-dimensional free fermionic chains*, *SciPost Phys.* **8** (2020) 036.

[41] M. Ljubotina, S. Sotiriadis and T. Prosen, *Non-equilibrium quantum transport in presence of a defect: the non-interacting case*, *SciPost Phys.* **6** (2019) 004.

[42] V.S. Viswanath, S. Zhang, J. Stolze and G. Müller, *Ordering and fluctuations in the ground state of the one-dimensional and two-dimensional s=1/2 XXZ antiferromagnets: A study of dynamical properties based on the recursion method*, *Phys. Rev. B* **49** (1994) 9702.

[43] D.J. Yates, A.G. Abanov and A. Mitra, *Lifetime of almost strong edge-mode operators in one-dimensional, interacting, symmetry protected topological phases*, *Phys. Rev. Lett.* **124** (2020) 206803.

[44] D.J. Yates, A.G. Abanov and A. Mitra, *Dynamics of almost strong edge modes in spin chains away from integrability*, *Phys. Rev. B* **102** (2020) 195419.

[45] B. Bhattacharjee, X. Cao, P. Nandy and T. Pathak, *Krylov complexity in saddle-dominated scrambling*, *J. High Energy Phys.* **2022** (2022) 1.

[46] H.A. Camargo, V. Jahnke, K.-Y. Kim and M. Nishida, *Krylov complexity in free and interacting scalar field theories with bounded power spectrum*, *J. High Energy Phys.* **2023** (2023) 1.

[47] O. Lunt, T. Kriecherbauer, K. T.-R. McLaughlin, and C. von Keyserlingk, “Emergent random matrix universality in quantum operator dynamics,” arXiv:2504.18311v2

[48] G. Pinna, O. Lunt, and C. von Keyserlingk, “Approximation theory for Green’s functions via the Lanczos algorithm,” arXiv:2505.00089v1

[49] R. Sasaki, “Towards verifications of Krylov complexity,” *Prog. Theor. Exp. Phys.*, **2024**, 063A01.

[50] T.S. Chihara, *An Introduction to Orthogonal Polynomials*, Gordon and Breach, New York (1978).

[51] H. Tal-Ezer and R. Kosloff, *An accurate and efficient scheme for propagating the time dependent Schrödinger equation*, *J. Chem. Phys.* **81** (1984) 3967.

[52] A. Vijay, R.E. Wyatt and G.D. Billing, *Time propagation and spectral filters in quantum dynamics: A Hermite polynomial perspective*, *J. Chem. Phys.* **111** (1999) 10794.

[53] R. Chen and H. Guo, *The Chebyshev propagator for quantum systems*, *Comput. Phys. Commun.* **119** (1999) 19.

[54] A. Weiße, G. Wellein, A. Alvermann and H. Fehske, *The kernel polynomial method*, *Rev. Mod. Phys.* **78** (2006) 275.

[55] R.D. Soares and M. Schirò, *Non-unitary quantum many-body dynamics using the Faber polynomial method*, *SciPost Phys.* **17** (2024) 128.
[56] Z. Y. Fan, “Generalised Krylov complexity,” *arXiv preprint arXiv:2306.16118*, 2023.

[57] J. Barbón, E. Rabinovici, R. Shir and E. Trivella, *On the evolution of operator complexity beyond scrambling*, *J. High Energy Phys.* **2019** (2019) 264.

[58] R. Koekoek, P.A. Lesky and R.F. Swarttouw, *Hypergeometric Orthogonal Polynomials and Their q-Analogues*, Springer Berlin Heidelberg (2010), 10.1007/978-3-642-05014-5.

[59] B. Bhattacharjee, X. Cao, P. Nandy and T. Pathak, *Operator growth in open quantum systems: lessons from the dissipative SYK*, *J. High Energy Phys.* **2023** (2023) 1.

[60] P.A. Clarkson, K. Jordaan and A. Loureiro, *Generalized higher-order Freud weights*, Proc. R. Soc. A **479** (2023) 20220788.

[61] K. Dilcher and L. Jiu, *Orthogonal polynomials and Hankel determinants for certain Bernoulli and Euler polynomials*, J. Math. Anal. Appl. **497** (2021) 124855.

[62] J. Wang, M.H. Lamann, R. Steinigeweg and J. Gemmer, *Diffusion constants from the recursion method*, Phys. Rev. B **110** (2024) 104413.

[63] R. Askey, *Continuous Hahn polynomials*, J. Phys. A: Math. Gen. **18** (1985) L1017.

[64] H.T. Koelink, *On Jacobi and continuous Hahn polynomials*, Proc. Amer. Math. Soc. **124** (1996) 887.

[65] A.P. Prudnikov, Yu.A. Brychkov and O.I. Marichev, *Integrals and Series, Volume 3*, Gordon and Breach, New York (1990).