Quantum Epidemiology: Operator Growth, Thermal Effects, and SYK

In many-body chaotic systems, the size of an operator generically grows in Heisenberg evolution, which can be measured by certain out-of-time-ordered four-point functions. However, these only provide a coarse probe of the full underlying operator growth structure. In this article we develop a methodology to derive the full growth structure of fermionic systems, that also naturally introduces the effect of finite temperature. We then apply our methodology to the SYK model, which features all-to-all $q$-body interactions. We derive the full operator growth structure in the large $q$ limit at all temperatures. We see that its temperature dependence has a remarkably simple form consistent with the slowing down of scrambling as temperature is decreased. Furthermore, our finite-temperature scrambling results can be modeled by a modified epidemic model, where the thermal state serves as a vaccinated population, thereby slowing the overall rate of infection.


Introduction & Summary
In chaotic quantum many-body systems, operators grow in size as time evolves. For example, in spatially local systems one expects that the extent of an operator O (t) grows as since the new terms generated by taking [H, O (t)] will live on the boundary of the domain of O (t) [1][2][3][4]. Consequently, the extent grows linearly with an effective "speed of light" vt.
Up to exponential error all operators outside the effective light-cone will commute with O (t). This effective speed of light is known as the Lieb-Robinson velocity [1]. This highlights the fact that space can be a derived concept in quantum mechanics, as without the Hamiltonian there may not be a sense in which one piece of the Hilbert space factorization is closer to another. Now we would like to contrast this behavior with that exhibited by q-local systems, where the Hamiltonian couples all the degrees of freedom together in q-body interactions. Consequently, there is no notion of spatial locality, and we accordingly refer to such interactions as coupling together "internal" degrees of freedom. Yet, there remains structure in the evolution of operators in these systems, as we often find that the sizes of operators grow exponentially  [5][6][7][8][9][10][11]. 1 Systems with both spatial locality and a large number of internal degrees of freedom-such as (chaotic) field theories in the large-N limit -display both linear spatial growth and exponential internal size growth [12][13][14][15][16]. The growth of an evolving simple operator W (t) can be probed using another simple operator V , using the (anti-)commutator squared |[W (t) , V ]| 2 or their corresponding out-of-time-order correlator (OTOC) W † (t) V † W (t) V [12,[17][18][19][20][21][22][23][24][25].
In order to develop the coarse-grained profile of operator growth, one must compute many OTOCs. The "chaos bound" [26] obeyed by OTOCs suggests that after an initial dissipation time, they de-correlate no faster than exponentially, with a rate λ L no larger than 2πT where T = 1/β is the temperature. This implies that presence of the thermal state ρ ∝ exp (−βH) slows down the effective growth rate of operators as temperature is decreased.
The Heisenberg evolution of operator O (t) is independent of temperature, so the entire effect of temperature must be contained in the matrix elements of O. Therefore, the natural finite temperature generalization of operator size has remained an open question (one recent proposal is given in [27]).
In this article, we address this issue by characterizing not only the average size of an operator but its entire size distribution. We then can define the effective size distribution of an operator at finite temperature by how it changes the size of the square root of thermal density operator ρ 1/2 (we explain why this is a natural choice in section (3)).
This definition leads to some nontrivial general results, independent of the details of the specific physical system. In particular, we observe that in generic fermion systems, the effective size of a single fermion operator is "thermally renormalized" to a value δ β = G (β/2) smaller than 1, where G (τ ) is the thermal two-point function. The size of the thermal operator ρ 1/2 itself is N 2 (1 − δ β ), determined by the same renormalization factor δ β . To gain a more explicit understanding, we will work in the context of the SYK model [28,29], a q-local Hamiltonian built out of N flavors of Majorana fermions, which saturates the chaos bound at low temperatures [20,29,30] The remainder of the article is organized as follows. We begin in section (2) by building up the notion of operator "size". First, we show that one may expand any operator O (t) in an orthonormal operator basis of the unique products of Majorana flavors. In the doubled theory, the operator basis maps to an orthonormal basis of states in the doubled Hilbert space. We define a "size" operator n in the doubled theory counting the average number of flavors in an operator basis state. We are then able to demonstrate that four-point functions measure the average "size" of an operator. Therefore, the de-correlation of a thermal OTOC is exactly equivalent to the growing average size of the operator ψ 1 (t) ρ 1/2 Tr ρ 1/2 ψ 1 (t) ψ j ψ 1 (t) ρ 1/2 ψ j = 1 − 2 N n ψ 1 (t) ρ 1/2 (1. 3) The average size of the operator ψ 1 (t) ρ 1/2 starts at n ψ 1 ρ 1/2 ≈ n ρ 1/2 = N 2 (1 − δ β ) and then grows sigmoidally in time, eventually saturating (scrambling) at a value of N/2.
Up to this point we have been discussing the average size of the operator ψ (t) ρ 1/2 . In fact, the entire size distribution of this operator has physical significance. Hence, in section (3) we construct generating functions for operator size distributions by inserting a weighting factor exp (−µn). First, we study the size distribution of the thermal operator ρ 1/2 by setting up a generating function Z µ ρ 1/2 , which is similar to a grand canonical partition function. Next, we show that the fractional distance to scrambling for the operator ρ 1/2 is always given by δ β ≡ 1 − n/n * = G (β/2) < 1. Then, we set up the generating function for the size distribution of ψ (t) ρ 1/2 , which we find naturally splits into a product of Z µ ρ 1/2 and a modified two-point function G µ (t). We show that the µ-expansion of G µ (t) determines the growth distribution induced by multiplying ρ 1/2 by ψ (t), and that G µ (t) is simply the two-point function for the original theory with a µ-dependent twisted boundary condition. We conclude the section by noting that on average, the size increase induced by a single fermion is given by the fractional scrambling distance δ β , which leads us to propose that δ β should be interpreted as a thermally renormalized unit of size.
Everything in sections (2) and (3) applies to general fermionic systems. In sections (4) and (5), we apply this methodology to the large-q SYK model. Solving the large-N saddle point equation in the large q limit with our µ-dependent twist, we obtain the full operator growth structure. After a dynamical renormalization of coupling constant J (which after a short amount of time essentially amounts to replacing the coupling with a smaller βdependent constant) and a renormalization of size unit from 1 to δ β , we observe that the full growth structure has the same functional form as the infinite temperature case. The dynamical renormalization of the coupling is the signature of the slowdown of the effective growth rate as temperature is decreased.
We conclude the section by discussing how to understand this finite temperature slowdown of scrambling in an epidemic model, where the thermal factor effectively vaccinates a large subset of the population, thereby slowing down the overall infection rate. We end the paper by discussing implications and future directions in section (6).

Operator Distributions and Two-Sided Wavefunctions
As an operator O (t) evolves in time, it becomes supported along operators of increasing size. This can be inferred from the Heisenberg equation of motionȮ (t) = i [H, O (t)]. Now, in order to properly discuss how much one operator is supported along another, we need an operator inner product. When the Hilbert space is finite-dimensional it is natural to use the Frobenius inner product: We may then expand operators in an orthonormal operator basis, which amounts to inserting a complete set of operators Note that at this point we have set up a Hilbert space of operators. If the original Hilbert space H has dimension L, the operator Hilbert space is H ⊗ H, with dimension L 2 .

Purification
Since H ⊗ H is isomorphic to H ⊗ H, one can always maps each operator to a quantum state in a "doubled" system with Hilbert space dimension L 2 . More explicitly, this mapping is defined by considering two copies of the original physical system, named as L and R, and introducing a maximally entangled state |0 (see Fig. 1(a)). For any maximally entangled state, there is a basis choice of the form {|n L ⊗ |m R } such that For later convenience we have chosen the norm of the state to be 0|0 = L.
Then the operator-to-state mapping is defined by where O L is the operator O acting on the Hilbert space of the left system, as is illustrated in Fig. 1(b). It is easy to verify that the inner product is determined by the Frobenius inner product of the corresponding operators. Our orthonormal basis of operators {Γ I } will thereby serve as an orthonormal basis of states |Γ I . Thus, the problem of understanding how O (t) is distributed across a particular choice of basis operators is equivalent to understanding how the two-sided state |O (t) is distributed across a particular choice of two-sided basis states. Since the choice of maximally entangled state |0 is not unique, the operator-to-state mapping has an ambiguity of U (L) acting on the R system. Since the same transformation is performed to |O and the basis vector |Γ I , all our discussion will be independent from this freedom of basis choice.

Orthonormal Basis of Operators
One of the simplest algebras with an interesting finite dimensional representation is the algebra of N flavors of Majorana fermions, where N is even: Note that this implies that ψ 2 j = 1, which will be convenient for our purposes, unlike the more common convention where {ψ i , ψ j } = δ ij and thus ψ j = 1/2. Such operators are traceless, Hermitian, and unitary. Furthermore, the algebra is invariant under taking any single ψ i → −ψ i , and so the product of any subset of the N fermions is also traceless. Thus, it is easy to construct an orthogonal operator basis by taking unique ordered products of Majorana fermions where the pre-factor has been inserted so that the resultant Γ I matrices are Hermitian. All nontrivial Γ I (with k > 0) are traceless. Since the product Γ I Γ J is also a string of fermions, which is only trivial when I = J (when the two strings are identical and Majorana fermions pairwise cancel), we have Furthermore, the basis operators Γ I have simple algebraic relations, since they either commute or anti-commute according to the relation where |I| is the number of elements in the multi-index I.

Mapping Basis Operators to Basis States
The purification isomorphism is quite simple to realize, as has been discussed by [32,33]. We consider two copies of the original system, which contains 2N Majorana fermions labeled by ψ L j and ψ R j , j = 1, 2, ..., N . We then define a maximally entangled state |0 , We may think of this state as a vacuum (all spins down, all bits set to 0) with regards to a set of entangled complex fermions operators where c † j = (c j ) † . Since state |0 is the ground state of a quadratic Hamiltonian H = j c † j c j , it is straightforward to compute the entanglement entropy [34] and verify that the state is maximally entangled between L and R. As we discussed earlier, the choice of |0 is not unique, but this choice is convenient for our purpose. The basis operators Γ I are mapped to states in the doubled system of 2N Majorana fermions: Therefore each basis operator Γ I is mapped to a particular fermion configuration in the doubled system, with fermions i 1 , i 2 , ..., i k , as is illustrated in Fig. 2(a). Essentially, the identity operator maps to the vacuum and nontrivial operators are mapped to excitations in the doubled theory.

Four-Point Functions Probe Operator Size
At this point, we can discuss the number operator n j ≡ c † j c j , which returns 1 when applied to basis states containing the flavor j and zero otherwise Thus, we see that for a generic operator O, the expectation value of n j returns the percentage of basis operators in O containing flavor j. Furthermore, we note that this expectation value is closely related to a one-sided four-point function (see Fig. 2), since Here we have assumed O to be fermionic. In the first two steps, we simply plugged in the definitions of n j and |O . In the third step, we anti-commuted iψ R j through O L , as right fermionic operators anti-commute with left fermionic operators. Then, we used the definition of |0 (2.8) to replace −iψ R j |0 with ψ L j |0 . Afterwards, we had an expectation value of only left operators for a maximally entangled state, so we traced out the right Hilbert space entirely, leaving us with an infinite temperature four-point function of the left-only system.
The relationship between operator quantities and one-sided correlators is simpler in terms of the anti-commutator squared, since we have where we used (2.12) to replace Tr O † ψ j Oψ j with (2n j − 1). One should note that if O is bosonic, the right-hand side of (2.12) will acquire a minus sign and the anti-commutators in (2.13) will be replaced with commutators. We denote the average value of n j in operator O as We can also define a total number operator (a.k.a. size operator) that returns the number of flavors or size of a basis state with Tr where the anti-commutators are replaced with commutators if O is bosonic. Alternatively, we may flavor average Eq. (2.12) in order to relate the flavor-averaged four-point function to the average size Noting that the number of unique products of k Majoranas Eq. (2.5) goes as N k , we see that the most common size is that of N/2. Indeed, a generic operator should be equally supported across all unique products, leading to a distribution of the form P k = N k /2 N . Thus, a totally scrambled operator has a size n * = N/2, so we see that the flavor-averaged four-point function measures the average fractional distance 1 − 2n/N an operator's size is from this scrambled value. Therefore we define the fractional scrambling distance operator

Operator Size Generating Function
By defining the number operator n, we can now go beyond the average operator size probed by four-point functions. Rather than just the average, we study all moments systematically by introducing a generating function [10] Z By taking derivatives of the generating function we can obtain all moments of n: A more useful expansion is a Taylor expansion in e −µ : in which the coefficients P n [O] is the percentage of terms in O having size n.

Including Temperature
The main goal of the current work is to understand the role of temperature in operator growth. After all, the dynamics of the operator ψ 1 (t) under Heisenberg evolution has no knowledge about temperature. One natural way to introduce temperature is to consider the operator ρ 1/2 where ρ = Z −1 β exp (−βH) is thermal state at inverse temperature β. The purification of ρ 1/2 is the thermofield double (TFD) state |T F D (the other factor of ρ 1/2 from the full ρ is used to make T F D|) where the Hamiltonians H L , H R are required to satisfy the condition (H L − H R )|0 = 0. This state is a natural choice for studying thermodynamic properties, because for each operator O, we can consider the corresponding operator Oρ 1/2 , and its average size will be directly measured by the finite temperature four-point function: Now, we are interested in uncovering the size distribution of a time evolved "thermal" operator O (t) ρ 1/2 . To do so, we will define a generating function for its size moments Z µ O (t) ρ 1/2 . We shall find that such a generating function naturally factorizes into a product of the generating function for the Gibbs state Z µ ρ 1/2 and a "connected" piece G µ [O (t)]. By using this to extract the size distribution of ρ 1/2 from the size distribution of O (t) ρ 1/2 , we find a natural definition for a "thermal" size of O (t).

Thermal State
We begin by studying ρ 1/2 . Taking O = ρ 1/2 in Eq. (2.16), we find the following relation between the thermal two-point function and the size operator This relation tells us that the most de-correlated value of the Euclidean two-point function -G (β/2) -is equal to the fractional distance the operator ρ 1/2 is from being scrambled which implies that the average size of ρ 1/2 is given by In the high temperature limit β → 0, one expects G(β/2) G(0) = 1, since the fermions square to one (2.4), which is consistent with the fact that ρ 1/2 approaches identity and the size shrinks to zero. On the contrary, in the low temperature limit β → ∞, if G(β/2) → 0, the size of ρ 1/2 approaches the scrambled (typical) value N/2. This result is very general since the imaginary time two-point function G(τ ) decays in most physical systems. For example, in all systems with a unique ground state and an excitation gap, G(τ ) decays exponentially at low temperature limit, so that G(β/2) → 0 when β → ∞. In a conformal field theory, G(τ ) decays in power law in the zero temperature limit, which also leads to the same length n β→∞ = N/2.
To learn more than just the average size, we construct the generating function Therefore learning about the operator distribution of ρ 1/2 is equivalent to learning about the fermion number distribution in the thermofield double state.

Thermal Fermion
, we see that the average size of ψ (t) ρ 1/2 is entirely equivalent to an out-of-time-order correlator (a.k.a. OTOC) [17][18][19][20]: Therefore, the statement that the OTOC de-correlates exponentially is equivalent to the statement that the average size of the operator ψ 1 (t) ρ 1/2 grows exponentially. If the OTOC vanishes in long time, that implies that the size of ψ 1 (t)ρ 1/2 reaches the scrambled value n * = N/2.
The size distribution of ψ 1 (t)ρ 1/2 can be uncovered through the generating function The operator e −µn can be viewed as an Euclidean time evolution with time µ and Hamiltonian n, so that the generating function (3.8) is related to the two-point function in a system with time-dependent Euclidean evolution: where T is the Euclidean time ordering symbol and ψ L 1 (τ a,b ) are imaginary time evolved fermion operators. Note that the denominator in (3.9) is, up to a factor of thermal partition function Z β that cancels with the numerator, exactly the size generating function Z µ ρ 1/2 in Eq. (3.6). Therefore, the size generating function Z µ ψ 1 (t)ρ 1/2 naturally factorizes into the product of This equation clarifies that the two-point function G µ (β/4 + + it, β/4 − + it) measures the size change ψ 1 (t) induces upon ρ 1/2 through multiplication. We can see this directly by applying the expansion in Eq. (2.20) to both sides of this equation, in order to obtain the following convolution formula for the size distribution of ψ 1 (t) ρ 1/2 defined by the expansion of G µ (β/4 + + it, β/4 − + it) in powers of e −µ : In this sense, K β m can be viewed as the "growth distribution" caused by applying ψ 1 (t) to the thermal state ρ 1/2 .
Note that the discussion above can be generalized to arbitrary operators. For an arbitrary operator O, as long as we normalize it such that measures the effective size distribution of O when applied to the thermal state.

Twisted Boundary Condition
We are interested in studying the generating function Z µ [O] for O = ρ 1/2 and O = ψ 1 (t)ρ 1/2 . As we discussed earlier, inserting the operator exp (−µn) corresponds to changing the imaginary time evolution. The computation can be simplified by noticing that exp (−µn) is a size distribution Figure 3: Schematic illustration of the size distribution P n ψ 1 (t) ρ 1/2 for the operator ψ 1 (t) ρ 1/2 (black curve) which naturally decomposes into a convolution of a growth distribution K β n [ψ 1 (t)] (blue dashed curve) with the size distribution P n ρ 1/2 of the operator ρ 1/2 (red dashed curve). This is due to the factorization relation of their respective generating functions (3.10).
Gaussian operator, such that its action by conjugation to fermion operators ψ L,R i leads to a simple linear superposition: As a result, inserting the operator-weighting term exp (−µn) is equivalent to twisting the boundary condition of the fermion fields at τ = β/4. It is convenient to "de-purify" the system and return to the single copy of fermion fields, but with a twisted boundary condition. The single field is defined by continuously stitching the left and right fields together: with the requirement of course that ψ (τ + β) = −ψ (τ ). This stitching transforms the purified action for the two fields into the original action for this single field; however, the twist condition must accompany the fields. In conclusion, the two-sided path integral in the presence of the factor exp (−µn (β/4)) equals the original path integral where the fields are twisted according to Therefore, we conclude that calculating the two point function G µ is equivalent to calculating the original two-point function, but with the following twisted boundary conditions We note that while these conditions break time translation invariance, they preserve a set of discrete symmetries. Specifically, if the original Hamiltonian is time-reversal invariant, then G µ (τ 1 , τ 2 ) has reflection symmetry across the lines τ 1 ± τ 2 = nβ/2 for all integers n ∈ Z Thus, we need only to solve for G µ (τ 1 , τ 2 ) in the fundamental domain 0 < τ 1 − τ 2 < β/2 and β/2 < τ 1 + τ 2 < β, as shown in Fig. 4(b) by the union of regions I and II.

Thermally Renormalized Unit of Size
As an interesting application of our formalism, let us note how ψ 1 (t) affects ρ 1/2 by taking t = 0 and consider the change of average size by a single fermion operator ψ 1 .
At infinite temperature β → 0, G (β/2) = 1, which restores the trivial result that ψ 1 increases the size of the density operator (which is proportional to identity operator, with size 0) by 1. At finite temperature, interestingly, the size change induced by a single fermion operator is smaller than 1, and is given by the same imaginary time two-point function as the one that determines the fractional scrambling distance δ β = 1 − n[ρ 1/2 ] N/2 in Eq. (3.4). In general, the size increase induced by ψ i is ∆n β [ψ i ] = G ii (β/2), which may depend on i. The average size increase is exactly Physically, the average size change due to applying a single fermion is generically δ β < 1 at finite β, because in the presence of a nontrivial ρ 1/2 there is a chance that multiplying by ψ 1 decreases the size, as is illustrated in Fig. 5, although the chance of increasing the size is always bigger. The closer the length of ρ 1/2 is to the scrambling value n * = N/2, the smaller is the size increase ∆n β ψ 1 ρ 1/2 . For a fully scrambled operator with n = N/2, δ = 0, multiplying a fermion ψ 1 has equal chance of increasing or decreasing the size, so that the average size stays the same.
It should be emphasized that the discussion above is not restricted to the thermal density operator. For any operator O, we can define the size change and obtain the following identity: The only thing special for the thermal density operator is the relation of ∆n to imaginary time two-point function in a single-copy system. Furthermore, instead of ψ i we can consider a string  from O. We have If we average over all Majorana strings Γ I with the same size k, we obtain This observation suggests that at finite temperature (or more generally, for any density operator ρ), the fractional scrambling distance δ, rather than 1, should be considered as the fundamental unit of size, which is carried by each fermion operator. Indeed, as we will discuss in next section, our calculation in the SYK model in the large q limit suggests universal behavior occurs when size is measured in this renormalized unit.

SYK Model
In this section, we will study the operator size growth in the SYK model [28,29]. This model features q-local interactions with independently random couplings, where each of the couplings is normal distributed . This is because there are O (q) locations at any depth of a given graph to insert another melon; however, there are typically only O (1) locations to thread another melon. Accordingly, we take our coupling J to equal J / √ 2q. Consequently, the O (q) combinatorial enhancement gained for each new melon insertion is canceled by the J 2 = J 2 / (2q) factor accompanying said melon. This q-scaling of the coupling isolates the infinite subset of the planar graphs where the graphs are two copies of a tree that are then glued together (a.k.a. "doubletree" graphs) such as (a). All non-doubletree graphs such as (b) are suppressed in q since they receive factors of J 2 / (2q) for each melon, but do not receive the necessary number of q combinatorial enhancements.
At large N , the two-point function satisfies the saddle-point equations where bracketed terms are Matsubara frequency matrices. One should note that since the fermions square to one, [G 0 ] −1 = −iω/2 rather than −iω.

Large q Approximation
In the language of Feynman diagrams, the Schwinger-Dyson equation (4.2) corresponds to only keeping the leading "melon diagrams" as is shown in Fig. 6. All other diagrams are subleading in large N . In the large q limit, there are two types of diagrams. Those with melons inserted into melons (such as Fig. 6(a)) receive a combinatorial q enhancement, as there are many rungs upon which one may insert (hence the need for a q −1 factor in the self-energy to keep everything finite). In contrast, diagrams where melons are simply threaded together (such as Fig. 6(b)) do not receive this enhancement [35]. Thus, at large q only the former dominate, which corresponds to the following truncation of the Schwinger-Dyson expansion: Combing the equations together and Fourier transforming, one obtains The role played by G 0 in this equation is to require that G → G 0 as τ 12 goes to integer multiples of β. Therefore, if we take G = G 0 e σ/q with σ → 0 at the τ 12 boundaries, we obtain Liouville's equation [30,36] where the field σ is expected to be periodic in both of its arguments, as well as have kinks when τ 12 approaches integer multiples of β. Now in order to find G µ , we will need to solve the above equations with the twisted boundary conditions (3.17). Furthermore, our twisted two-point function G µ also satisfies the reflection conditions (3.18). Thus, we need only solve for G µ in the fundamental domain 0 < τ 1 − τ 2 < β/2 and β/2 < τ 1 + τ 2 < β, as shown in Fig. 4(b) by the union of regions I and II.

The Large-q solution
In the large q limit, each commutator with the Hamiltonian increases the size of operator by ∼ q, so that it is natural to measure the operator size in unit of q. In the generating function, this corresponds to definingμ ≡ qµ (for reasons to be explained in the next subsection we will actually use the slightly smaller variableμ = qδ β µ (4.16)), withμ kept finite in the large q limit. The derivative of the generating function overμ measures the size n in unit of q. If we consider the large q limit withμ being kept finite, and use the large-q ansatz for twisted two-point function G µ (τ 1 , τ 2 ) = G 0 (τ 1 , τ 2 ) e σµ(τ 1 ,τ 2 )/q , (4.6) the boundary condition (3.17) reduces to Thus, to the leading order of 1 q , the two equations for β/4 and 3β/4 decouple. As a reminder, these equations encode the effect of moving a Majorana fermion on the Euclidean circle (i.e. the two TFD half-circles) from one side of the twist operator (3.9) to the other. The factor of exp (−μ/q) = exp (−µ) denotes the fact that such an action changes the size of the doubled state by exactly one whole fermion. However, we shall find that for states of large size such as the thermofield double (i.e. Gibbs state) each Majorana fermion actually increases the total average operator size by a smaller fraction δ β rather than 1, where the size change δ β is smaller when the state's size is larger. As such, we shall find that it will be appropriate to use the variableμ = qµδ β instead ofμ = qµ when taking the large-q limit.
The twisted two-point function in large-q limit can thus be obtained by solving Liouville's equation (4.5) with the µ-dependent boundary conditions (4.7). Here we will skip the tedious details and directly present the solution. When the times are such that the two fermions are on the same side of the twisted boundary, which corresponds to τ 2 > β/4 (region I in Fig.  4(b), we find a seemingly time-translation invariant solution solution However, when the times are such that the two fermions are on opposite sides of the twisted boundary at β/4 and 3β/4, which for our domain amounts to the condition τ 2 < β/4 (region II in Fig. 4(b), the time translation symmetry is explicitly broken Here the parameters α µ and γ µ are functions of βJ and µ, which are determined by the boundary condition G(τ, τ ) = 1 as well as the reflection conditions (3.18) α µ β = βJ sin γ µ , sin α µ β 2 + 2γ µ = e −μ sin α µ β 2 (4.10) In the limit µ → 0, we recover the untwisted two-point function G µ=0 (τ 1 , τ 2 ) = G µ=0 (τ 1 −τ 2 ) = G (τ 1 − τ 2 ) in the whole domain, and the equation for the parameters reduce to the ordinary case [30]: The asymptotic behavior at small values of βJ and large values of βJ respectively are given by

Size renormalization
Before carrying further analysis to the SYK operator growth in next section, we need to discuss an important modification to the two-point function solution due to higher order q effects. If we take τ 1 → β/4 + , τ 2 → β/4 − in Eq. (4.9), we obtain This is the kernel that determines the size change induced by multiplying ψ 1 to ρ 1/2 , which has been discussed in section (3.4). Taking the µ-derivative of G µ , we find However, we also know that the size change is directly determined by the two-point function due to Eq. (3.19): where α ≡ α µ=0 (βJ ) is the smallest positive root of Eq. (4.11). This discrepancy between the two calculations is because δ β → 1 in the large q limit, and the O (q −1 ) difference is neglected in the approximation we made to the boundary condition. The easiest way to resolve this issue and makes a consistent large-q limit is by redefininĝ µ = qµ toμ = qµδ β (4. 16) in Eq. (4.7), which leads to the same substitution in Eqs. (4.8), (4.9), and (4.10). In the following, we will always use this definition ofμ. Physically, this substitution is a consequence of the size renormalization discussed in section (3.4). Each Majorana fermion increase the operator size by δ β rather than 1. Each action of the Hamiltonian increases the operator size by ∼ qδ β . Although in large q limit 1 − δ β is order q −1 , it is important to keep track of this distance, since the same δ β also measures the fractional scrambling distance of ρ 1/2 , as we discussed in section (3.1). The size of ρ 1/2 is which decreases with increasing q, but is always large since we should always take the large N limit before taking large q.

SYK Operator Growth
We are now equipped with everything we need to understand P n ψ 1 (t) ρ 1/2 , the size distribution of ψ 1 (t) ρ 1/2 . According to Eq. (3.10), the generating function Z µ ψ 1 (t)ρ 1/2 for this distribution splits into a product of the thermal state's generating function Z µ ρ 1/2 and G µ (β/4 + + it, β/4 − + it). The latter is simply the twisted two-point function we discussed in the previous subsection with an analytic continuation.

Thermal State
The generating function Z µ ρ 1/2 is the partition function of the system with the insertion exp (−µn (β/4)) divided by that of the original system (see Eq. (3.6)). This quantity can be determined by the twisted two-point function, since one has In theory, we can integrate this equation to obtain Z µ ρ 1/2 . However, many important properties of the distribution can be inferred from just the first and second moment.
The first moment is simply the average size The behavior of n ρ 1/2 for two different q values are plotted in Fig. 7 as a function of βJ . Interestingly, we see that for larger q it takes larger values of βJ to achieve the same average size. Thus, the exponentiation of increasingly heavy SYK Hamiltonians results in relatively lighter thermal states. This is unintuitive, so one might argue that it is simply due to the q-scaling nature of J , but that is only a constant shift in the log scale plot in Fig. 7, which is nowhere large enough to account for the above discrepancy. The true origin of this effect is the power of 2/q in the two-point function. We conclude that this heavy-light relationship is thus a non-trivial consequence of the large-N and large-q limit.
For βJ e q , one expects that the higher order corrections to Liouville's equation (4.5) cannot be neglected, and so one should turn to the conformal approximation [30]. Their expression for the two-point function implies that the average size of ρ 1/2 when N βJ 1 is given by The difference between the large q and low temperature n ρ 1/2 = N 2 1 − G β 2 is captured by the factor c (q). It monotonically increases from c (4) = (2/π) 1/4 ≈ 0.9 when q = 4, and asymptotically approaches 1 when q is large as 1 − 2/q 2 . We expect that q = 4 and large βJ will be where our large q approximation will have the largest error. However, that this error is at worst 10% renews our confidence that the large q approximation captures important analytic features of the large-N SYK model.
The second derivative of ln Z µ determines the width of the distribution: Therefore the width of the distribution σ n ∝ √ N , such that the relative deviation from the average value σ n ρ 1/2 /n ρ 1/2 ∝ N −1/2 is sharply peaked in the large N limit. This is a consequence of large N factorization.

Thermal Fermion
As explicitly discussed in section (3.2), the generating function for the growth distribution K β (3.12) is determined by the twisted two-point function (4.9) where α µ and γ µ depend on µ and βJ through the constraints (4.10).

Average Size
This implies that the average size of the operator ψ (t) ρ 1/2 is given by where δ β = (α/J ) 2/q , and α ≡ α µ=0 (βJ ) is the smallest positive root of Eq. (4.11). We see that the difference in averages sizes of ψ (t) ρ 1/2 and ρ 1/2 is a simple when expressed in the renormalized size unit δ β , which inspires us to define a notion of the "average growth" of ψ 1 (t) as Now, scrambling occurs when the average size of ψ 1 (t) ρ 1/2 given by Eq. (5.6) reaches n * = N/2. This produces a slightly complicated expression for the scrambling time t * ; however, it simplifies dramatically when phrased in terms of the average growth of ψ 1 (t). Manipulating the scrambling time equation n ψ 1 (t * ) ρ 1/2 = N/2, we find that one may equivalently state that scrambling occurs when the average growth of ψ 1 (t) reaches n * = N/2 Figure 8: After either dynamical renormalization (5.12) or time re-parametrization (5.13), the growth distribution K β [ψ 1 (t)] takes the same form as the Heisenberg evolution of the operator ψ 1 (t) (i.e. the infinite temperature size distribution P [ψ 1 (t)]). This distribution is given by Eq. (5.14), and we plot it on a log-log scale. Note that it reaches out towards larger operators exponentially quickly.
This growth is consistent with the known result of large-q Lyapunov exponent [30]. λ L = 2α (5.9)

Full Growth Structure
In the Lyapunov regime, we may expand the generating function of the growth distribution as where δ β = (α/J ) 2/q , and α ≡ α µ=0 (βJ ) is the smallest positive root of Eq. (4.11). Grouping terms by powers of exp (−µ) and using the definition (3.12), we conclude that the growth distribution is given by where we note that (−1) n −2/q n is always positive for integer n. Thus, K β [ψ 1 (t)] ≥ 0 and so we have no negative probabilities in the size distribution of ψ 1 (t) ρ 1/2 , since it is given by P ψ 1 (t) ρ 1/2 = K β [ψ 1 (t)] * P ρ 1/2 as shown in section (3.2).
Interestingly, we see that the growth distribution K β (5.11) has a functional form independent of temperature, which we plot in Fig. (8). We can use either of two methods to expose this phenomenon. One option is to replace the coupling J with the dynamically renormalized couplingJ (t) (plotted in Fig. (9)): The other option is to re-parametrize timẽ Both methods transform the finite temperature growth distribution into that of the Heisenberg evolution of the operator ψ 1 (t) (i.e. the infinite temperature size distribution P 1+qn [ψ 1 (t)]) [11]. For example, usingt gives 14) This temperature-independence is fascinating since the Heisenberg evolution of ψ 1 (t) was obtained in [11] via fully-dressed Feynman graph calculations. In other words, the distribution P [ψ 1 (t)] represents the simple tree graphs such as Fig. 6(a) constructed using the original SYK Hamiltonian. However, we just showed how the growth distribution K β [ψ 1 (t)] can be easily transformed to P [ψ 1 (t)]. Therefore, since P ψ 1 (t) ρ 1/2 = K β [ψ 1 (t)] * P ρ 1/2 and P ρ 1/2 is well-peaked, we are led to the remarkable conclusion the growth dynamics of large-N , large-q SYK model is totally universal. In fact, if one waits an initial period α −1 to enter the Lyapunov regime, then one need simply use the effective size δ β and coupling J = α for the full growth structure of ψ 1 (t) ρ 1/2 to match that of ψ 1 (t).

Finite Temperature Epidemic Model
In this subsection we will discuss the physical interpretation of the SYK operator growth by relating it to an epidemic model. Intuition for operator scrambling behavior has been developed by various authors [6,11,22,37], resulting in an infection picture for operator growth. An operator such as ψ 1 (t) can be expanded in the strings of Majorana fermion Γ I . We consider the fermions already included in the string as "infected". Heisenberg evolution of Γ I generates a term [Γ I , H] which could contain a few more fermions. For example for SYK model with q-body interactions, in the large N limit most of the terms have one fermion replaced by q − 1 other fermions. In order for these q − 1 fermions to be "infected", they must not be already in Γ I . Therefore the infection rate depends on the infectable population.
In the simplest infection model for a population of n * individuals, the rate of infection is proportional to the number of unexposed people times the number of contagious people More generally, in various quantum circuit and Hamiltonian systems, both terms on the right-side of the equation may be raised to various powers or there may even be a sum of such terms, due to the potential multi-body nature of the interaction. For example, in SYK, upon a single commutation with the Hamiltonian, a size 1 operator becomes a size q − 1 operator, so we might expect various powers of q to appear in the above expression. Regardless, in either case sigmoidal behavior will be produced, which is consistent with general expectations of four-point functions.
Let us see just how well such a picture can apply to the SYK model. Taking the derivative of Eq. (5.6) and using Eq. (5.8), we find that during the Lyapunov regime (log N αt 1) d dt n ψ 1 (t) ρ 1/2 ≈ (2J ) 1 − n ρ 1/2 n * q/2 n ψ 1 (t) ρ 1/2 − n ρ 1/2 (5.16) Comparing with the infection equation (5.15), we have the fundamental rate r = 2J as well one of the terms being raised to q/2 due to the q-local nature of the interaction. However, rather than 1 − n ψ (t) ρ 1/2 /n * q/2 , which one may have expected by direct analogy with the infection equation, we have the static term 1 − n ρ 1/2 /n * q/2 = δ q/2 β . During the Lyapunov regime, these two are the same to leading order in N . Lastly, it appears through the final term that of the large population n ψ 1 (t) ρ 1/2 , only the small population n ψ 1 (t) ρ 1/2 − n ρ 1/2 possesses the ability to infect others. Notice that there remains the large population n ρ 1/2 who count as having been exposed, but do not infect others. It is thus natural to view this group as a vaccinated population.
In other words, after waiting for the dynamical renormalization/time re-parametrization to settle down, the physics of the four-point function is well-described by an infection model, with the caveat that only a small population n ψ (t) ρ 1/2 − n ρ 1/2 possess the ability to infect. In this sense, the operator ρ 1/2 vaccinates a finite fraction of the N flavors. Now regardless of whether any particular individual possess the ability to infect, it remains that a large portion of the population has been exposed, and thus the probability for any contagious Figure 10: Illustration of different epidemics at infinite temperature (left panel) and finite temperature (right panel). Green dots represents unexposed individuals, red dots represents contagious individuals, and blue dots represent vaccinated individuals. Finite temperature factors such as ρ 1/2 "use up" some of the available flavors for growth, resulting in collisions like those depicted in Fig. (5). This effect ends up being well-modeled by an epidemic where these flavors or individuals count as having been exposed, but do not spread disease. As a result of this large vaccinated population, it simply more rare for a contagious individual to encounter an unexposed individual, even at the start. Hence the rate of infection -the Lyapunov exponent -slows down, as seen in the right figure. individual to encounter an unexposed individual is decreased. Consequently, the overall rate of infection slows down to as illustrated in Fig. (10).

Discussion
The methodology developed in sections (2) and (3) is very powerful, as it applies to all fermionic systems. Specifically, determining the system's full growth distribution amounts to calculating the twisted (3.17) two-point function G µ followed by inverse transforming in µ (3.12). The large-N saddle point technique and the large-q simplification enabled us to obtain a closed solution in SYK. Even if analytics are too difficult, this analysis can be effectively implemented numerically for many classes of models.
These techniques also allowed us to compute a four-point function, since Eq. (3.2) shows that the average size of the operator ψ (t) ρ 1/2 (5.6) gives the value of a certain four-point function. We can generalize and calculate arbitrary four-point functions by moving the twist (3.17) to other locations. This has the non-trivial consequence that the twisted two-point function solves the ladder kernel [30,31]. In practice solving the former can be substantially easier than solving the latter. As an example, in [38] we use this "twisted" technique to derive an elegant expression for the large-q SYK four-point function at arbitrary coupling and temperature. Like the growth distribution methodology, this new method for calculating four-point functions works for all fermionic systems.
The dynamical renormalization of the coupling (5.12) plays a central role in this work. It will be important to understand this in a deeper and more general context. The success of the modified infection model in capturing the thermal operator growth suggests that the principle underlying the finite temperature slowdown in SYK is competition for Majorana flavors. The presence of various powers of the thermal state exp (−βH) "uses up" some finite fraction of the flavors. Consequently, when we apply a single fermion, there is a fractional probability for it to become absorbed and thus its size is renormalized (3.20) to a value based upon the percentage δ β of "unused" flavors. Now, the renormalized coupling J (t) (5.12) slows down during time-evolution. We believe that this occurs due to the same principle, but have not yet fully understood the mechanism. Our belief is motivated by the empirical observation that the Lyapunov exponent is a power of the percentage δ β of "unused" flavors λ L = lim t t dissipation 2J (t) = 2J (δ β ) q/2 ≡ 2J 1 − n ρ 1/2 n * q/2 (6.1) This kind of sigmoidal operator growth is generic in many-body chaos. However, without Majorana fermions, the manner in which the thermal factors interfere with operator growth must be more complicated, as there is not a bit-like notion of "using up" a flavor. Sigmoidal behavior signals the existence of a competition for some finite resource. For SYK, this resource was flavor; we only have N flavors with which to grow operators, so eventually we will be led to flavor collisions as in Fig. (5). However, flavor competition is only one aspect of competition for a more general resource. The question remains: what do operators compete for during evolution? Is there some sort of "operator entropy"? Perhaps when summed across "all operators" at finite temperature, there is always a fixed amount of total correlation with an initial simple operator due to unitarity. A better understanding of such a resource would give an organization to operator dynamics at different energy scales.
Our results have interesting implications for the holographic dual of the SYK model. This is simplest to understand when we explicitly express our results in terms of the doubled theory. We defined an entangled orthonormal basis for the doubled theory using the eigenstates (2.10) of the size operator (2.14). Taking the state ψ L 1 (t) |T F D 3 , we related its size wave-function squared (i.e. P n ψ (t) ρ 1/2 ≡ n|ψ L 1 (t)|T F D 2 ) to the size wave-function squared for the thermofield double state (i.e. P n ρ 1/2 ≡ | m|T F D | 2 ) n|ψ L 1 (t)|T F D 2 = n m=0 K β n−m [ψ 1 (t)] | m|T F D | 2 (6.2) isolating the time-dependence into the growth distribution K β [ψ 1 (t)] (3.11). Using this, we found the "average growth" of ψ 1 (t) (5.7) at low temperatures to be 1+2 (βJ sinh (πt/β) /π) 2 , which was shown in [41] to exactly match the classical momentum dynamics of a "boundary" particle falling into a near-extremal black hole. That is, the average growth of an SYK fermion exactly matches the average momentum of an infalling particle in a N AdS 2 black hole.
It is a striking result of our analysis that the full size wavefunction squared of the SYK fermion precisely relates to the full momentum wavefunction squared of the infalling particle. The universal form (5.14) of the growth distribution K β [ψ 1 (t)] precisely gives the squared coefficients of the AdS 2 momentum bulk-to-boundary propagator. Exploring this connection will be an important focus of future work 4 .