Operator size distribution in large N quantum mechanics of Majorana fermions

Under the Heisenberg evolution in chaotic quantum systems, initially simple operators evolve into complicated ones and ultimately cover the whole operator space. We study the growth of the operator “size” in this process, which is related to the out-of-time-order correlator (OTOC). We derive the full time evolution of the size distribution in large N quantum mechanics of Majorana fermions. As examples, we apply the formalism to the Brownian SYK model (infinite temperature) and the large q SYK model (finite temperature).


Introduction
Quantum scrambling is a process that local information disperses into the entire system under unitary evolution [1,2].It is vital for the thermalization in isolated quantum systems [3,4].Recent developments of information scrambling have connected different branches of physics, including the quantum theory of gravity [5,6], non-Fermi liquids in condensed matter physics [7] and many-body quantum chaos [8].
In the Heisenberg picture, scrambling manifests itself as the growth of simple operators [9].For a system of N Majorana fermions (with the convention 1 {χ j , χ k } = 2δ jk ), a single operator χ k evolves into a superposition of products of Majorana operators under time evolution χ k (t) = e iHt χ k e −iHt = odd n j 1 <j 2 <...<jn c j 1 j 2 ...jn (t) χ j 1 χ j 2 ...χ jn , (1.1) † guyingfei@gmail.com 1 The convention here differs from the one commonly used in the SYK related literature by a factor of 2.Here we requires the normalization χ 2 j = 1 for the orthonormality of the operator basis {χ j , χ j1 χ j2 , ..., χ j1 χ j2 ...χ jn , ...} where the inner product of two operators X and Y is defined as 2 −N/2 Tr(X † Y ).The prefactor is the inverse of the dimension of the Hilbert space.
where the coefficients c j 1 j 2 ...jn (t) are the (time dependent) amplitudes on the orthonormal operator basis and may be interpreted as the "operator wave function" [10].The length n of the operator string χ j 1 χ j 2 ...χ jn is called the "size" of this basis element.Accordingly, the probability distribution of the size is given as follows, P (n, t) = j 1 <j 2 <...<jn |c j 1 j 2 ...jn (t)| 2 . (1. 2) The conservation of the total probability n P (n, t) = 2 −N/2 Tr(χ k (t) † χ k (t)) = 1 follows immediately from the unitarity of the time evolution.In this measure, the operator growth is reflected as the shifting of the distribution towards a larger size.The notion of operator size is closely related to the out-of-time-order correlator (OTOC) [11,12] in many-body quantum chaos and the momentum of a particle falling towards a black hole horizon in holography [13][14][15].However, it is challenging to obtain the analytic expression for operator size and its distribution.Progress has been made in quantum circuits models [16][17][18][19][20][21], numerical studies and experimental proposals for spin models [9,[22][23][24], and bounds on the operator size in various contexts [25][26][27][28][29].For the large N quantum mechanics like the SYK model [30,31], analytic results were only achieved for the initial growth [10,32,33], leaving the saturation in the late time unexplained.In this paper, we develop a theoretical framework to determine the distribution of size over the full range of time.We apply the framework to two examples: first, the Brownian SYK model where we can compare the analytic result with numerics reported in Ref. [34] and find an agreement to good precision; second, the large-q SYK model at finite temperature.

Relation to OTOC
An operator O in a system of N Majorana fermions (χ j , j = 1, 2..., N ) can be mapped to a state |O in the doubled system consists of N complex fermions c j = (χ j + iψ j )/2, where ψ j labels auxiliary Majorana operators.The mapping can be constructed by acting the operator O on a reference state in the doubled system, e.g., O → |O = O|EPR with the reference state chosen to be a maximally entangled state between the original system χ and the auxiliary system ψ.
There are still plenty of freedom in such a maximally entangled state.For the purpose of this paper, we make a particular choice: we require n j |EPR = 0 (with n j = c † j c j ) for all j, namely, |EPR a zero-occupation state for all N flavors.The merit of this choice is that the size of an operator can now be identified with the total fermion occupation number of the corresponding state [32].To see this, consider a basis element O = χ j 1 χ j 2 ...χ jn in the operator space of length n, then its corresponding state has the following form as χ j = c j + c † j and the operator c j annihilates |EPR by construction.For such state, the operator size is equal to the total particle number n := j n j .Consequently, the size distribution (1.2) of the Heisenberg operator χ k (t) can be written as the following expectation value where Π(n) = δ n, j n j is the projection to the eigenspace with fixed particle number.
In the large N limit, instead of the discrete distribution P (n, t) whose argument n is unbounded, it is more convenient to use the continuous distribution P(s, t) with the continuous (rescaled) variable s ∈ [0, 1] representing the operator size in unit of N : where the projection Π(s) = δ(s − Σ j n j /N ) is a delta function instead of the Kronecker delta.
Here, we have dropped the dagger since χ † k (t) = χ k (t) for Majorana operators.The discrete and the continuous distribution are related at finite N via equation P(s, t) = N P (sN, t), where the prefactor N ensures the normalization 1 0 P(s, t)ds = 1.To compute P(s, t) in the large N limit, it is useful to consider the following generating function S(ν, t) := (2.4) The expectation value on the r.h.s.can be mapped back to a Keldysh correlation function in the original system (i.e.only involves χ fermion) as where T c denotes the contour ordering (gray line) on the double Keldysh contour and the auxiliary fermion fields ψ j (0) have been mapped to iχ j (0) using ψ j |EPR = iχ j |EPR .The discussions up to now are generally applicable to quantum systems of Majorana fermions.Next, we will focus on the interacting systems that exhibit many-body chaos with an exponentially growing OTOC in the early time.This assumption amounts to the following ansatz [35,36] (this is a non-linear generalization of the ansatz in Ref. [37]) where is regarded as the propagator of the scrambling modes.The prefactor C is of order N and therefore suppresses the higher order terms in OTOC at early time.
However, when the center of mass time difference t 1 +t2 −t 3 −t 4 2 is of order the scrambling time t scrambling ∼ ln N , the exponential growing factor is comparable with N .And we need to include the whole series to recover the correct late time behavior, e.g., the saturation of OTOC.Here, the integer "m" counts the number of scrambling modes (see Fig. 1 (a) for an example with m = 3), and the growing exponent κ in the propagator λ is known as the Lyapunov exponent for the many-body quantum system.The retarded and advanced vertex functions (Υ R,m , Υ A,m ) capture the interactions between the scrambling modes (wavy lines) and the fermions (straight lines, external legs): (2.7) Now, for the generating function S(ν, t) of the size distribution expressed as (2.5), its Taylor expansion in parameter ν consists of OTOCs of multiple operators The prefactor e − ν 2 comes from the constant piece in the exponent in (2.5).A typical connected diagram for S(ν, t) is shown in Fig. 1 (b) with n = 2. Now, we sum all such diagrams and obtain the following formula in terms of vertex functions and propagators of scrambling modes, where the time arguments for the vertex functions are zero and have been left out.

General formula
Following the manipulations in Ref. [36], the vertex functions Υ R/A,m can be rewritten as the moments of a distribution 2 h R/A , i.e.Υ R/A,m = +∞ 0 y m h R/A (y)dy.Inserting into (2.9),we obtain a simpler expression for the generating function whose time dependence is only through the propagator λ = C −1 e κt .The function f A in the exponent is the Laplace transform of h A , namely, It has the meaning of a correlation function on the perturbed background characterized by a mean-field parameter y.This function is also relevant in the discussion of quantum teleportation by traversable wormholes [38,39].
The above expression for the generating function (3.1) leads to an integral formula for the size distribution by the inverse Laplace transform In other words, the size distribution is connected to the distribution h R (y) by a change of variable from y to s = (1 − f A (λy))/2.The integral with the delta function results in an alternative expression using inverse function The general form of the size distribution function (3.3) is our main result, it implies the following universal behavior for large N quantum mechanics of Majorana fermions 1.In the early time, i.e. when λ = C −1 e κt 1, we keep to the linear order of λ and therefore the size variable s = 1−f A (λy) 2 ≈ 1 2 λyΥ A,1 grows exponentially with time.More precisely, the expectation value to its first order approximation grows exponentially as the linearized OTOC, consistent with the result in literature, e.g.Ref. [10].
2. The novelty of our general result is that it also captures the saturation at late time!For example, at very late time λ 1, the function f A 1 and we have converge to a delta function with s = 1/2.It means that the initially simple operator is maximally scrambled after a long enough evolution such that a typical component in the linear superposition (1.1) is an operator string of length about half of the total flavors, as it is generated by random decisions of whether or not to include a particular fermion operator.
In the following, we will insert known forms of the functions h R and f A from two SYK-like models to our general formulas, and show the analytical expressions for the size distribution.

Brownian SYK and comparison with numerics
The Hamiltonian of the Brownian SYK model [40] consists of random all to all time-dependent interactions given as follows  As time passes (from smaller λ to larger ones), the distribution shifts towards larger s and finally accumulates at s = 1/2; (b) Comparison between theoretical predictions (solid lines) and numerical results (markers) in the Brownian SYK model with q = 4 and N = 500.We plot the moments (2s) m instead of s m so that they saturate to the same value.
where J j 1 j 2 ...jq (t) are Gaussian variables with zero mean and variance J j 1 j 2 ...jq (t)J j 1 j 2 ...jq (t ) = (q − 1)!Jδ(t − t )/N q−1 .The functions h R and f A are obtained in [36] h In addition, for the propagator λ = C −1 e κt , we have the prefactor C = 2N ∆ 2 and Lyapunov exponent κ = 4(q − 2)J.Putting all these into (3.3),we have3 We plot the size distribution P(s, t) in Fig. 2 (a), showing a shift of weights towards the larger size and final convergence to a delta function at s = 1/2.From (4.3) we can also deduce the moments of the size distribution where U (a, b, z) is the confluent hypergeometric function.This formula generalizes the results obtained in Ref. [42], which are only for the first moment s.
Moreover, the average dynamics of the Brownian SYK model can be efficiently simulated using bosonic collective modes, where the relevant Hilbert space dimension grows polynomially with N [34] (consequently, the authors were able to study OTOCs upto N = 10 6 ).We extend the analysis to numerically evaluate the moments s m for N = 500 and q = 4, see Appendix A for details.We compare them with the predictions from our general results in Fig. 2 (b) and find good agreement to high precision.
5 Finite temperature operator size and the large-q SYK So far, our definition and discussions of the operator size have been limited to the infinite temperature system as we were focusing on the Heisenberg evolution of operators rather than the underlying states/ensembles (e.g. the EPR state that was used in (2.3) is a purification for an infinite temperature density matrix).It is also of great interest to generalize the notion of operator size to finite temperatures.However, there is no consensus yet on how to take into account the renormalization effect from the thermal ensemble.Existing proposals include [32,43].In this paper, we will propose a definition related to the operator "Q" in Ref. [36] (cf.section 4.3), which measures the "magnitude of the TFD deformation".
Recall that in infinite temperature case, we have mapped an operator O to a state in the doubled system |O = O|EPR , and related the operator size of O to the fermion number in |O .In this situation, |EPR (i.e. the infinite temperature TFD) is regarded as the vacuum that the fermion operator c j = (χ j + iψ j )/2 annihilates, and the size counts how many excitations an operator creates from the vacuum.Now, for finite temperatures, the vacuum in the doubled system is identified as the TFD | √ ρ β which purifies the thermal density matrix.We map an operator O to a symmetric4 "excited state" |O = |ρ Here, the fermion operator ξ j is designed to annihilate the TFD state.c β = χ j (iβ/4)χ j (−iβ/4) = (cos vπ 2 ) 2∆ is a normalization factor.We use the following definition [36] with the convention χ j (t + iτ ) = e −Hτ χ j (t)e Hτ for the imaginary arguments.
Next, we apply the definition (5.1) for the finite temperature operator size and calculate the size distribution for a Heisenberg evolved simple fermion operator χ k (t) in the large-q SYK model.The model was introduced by Maldacena and Stanford in Ref. [44] and is analytically tractable at all temperatures.The Hamiltonian is given as follows where J j 1 j 2 ...jq are random Gaussian variables with zero mean and variance J 2 j 1 j 2 ...jq = (q−1)!J 2 2qN q−1 .By definition, N is taken to infinity before q goes to infinity.We keep J fixed in this limit.A more convenient pair of parameters for this model are ∆ = 1/q and v, the latter is related to the coupling through equation v = βJ π cos πv 2 .∆ is the fermion scaling dimension, and v determines the time scale in this model, e.g., the Lyapunov exponent κ = 2πv β .The normalization factor c β = (cos vπ 2 ) 2∆ for this model.The computation for the distribution P(s, t) is similar to the Brownian SYK.The main difference is that there are more independent terms in ξ † j ξ j than its infinite temperature counterpart, see Appendix B for details.The upshot is that the different terms come with different phases comparing to (3.3), and resulting in the following expression for finite temperature size distribution which reduces to the infinite temperature result (3.3) when the Lyapunov exponent is negligible comparing to the temperature, i.e. v = κ/2πT = 0. Inserting the explicit form of the functions h R (y) = y 2∆−1 e −y /Γ(2∆) and f A (x) = (1 + x) −2∆ obtained in Ref. [36], we plot the size distribution with fixed λ = (4N ∆ 2 cos πv 2 ) −1 e 2πv β t and different v in Fig. 3, to emphasize the role of the non-trivial phases e ±iv π 2 .

Summary and discussions
In this paper, we provide a method to determine the size distribution function P(s, t) for large N Majorana systems with many-body chaos, extending the existing literature about initial growth to the full time range.Remarkably, the notion of operator size has been related to the quantum teleportation via traversable wormholes [14,15,38].Indeed, some of our discussions have counterparts in holography.For example, the integral form of the generating function (3.1) is similar to a two point function on the background of a traversable wormhole [45], where variable y corresponds to the momentum and f A (λy) represents a two point function on a shifted background parameterized by y.
Furthermore, the size distribution may be generally regarded as a new type of "order parameter" that refines the OTOC in diagnosing information scrambling.Recently, Refs.[46,47] considered the evolution of operator size in open quantum systems.In particular, [47] showed that the size distribution can be used to exhibit dynamical transitions between the scrambling phase and the dissipative phase for systems coupled to a heat bath.
B Details of the finite temperature operator size and comparison with Qi-Streicher's definition In this section, we present details of the finite temperature calculation with our definition of size (5.1).Again, we start with the generating function The prefactor in front of the expectation value and also in the exponent arises from the normalization in the definition of size (5.1).
To proceed, we represent (B.1) on the double Keldysh contour as where i0 ± = ±iε is an infinitesimal shift on the imaginary time, to indicate the positions of the operators as shown Fig. 4 (a).The four terms in the exponent can be grouped into two classes: (η 3 , η 4 ) = (++), (−−), corresponding to the OTOCs, and (η 3 , η 4 ) = (+−), (−+), corresponding to the normal-order correlators that can be factor out in the large N limit.Therefore, we have 3) Following the same diagrammatic analysis as in the main text, we obtain the following formula are obtained by restricting the general functions h R/A (y, θ), f R/A (x, θ) (derived in Ref. [36], cf.section 6.3) to special case when θ = β/2 (we have further simplified the expressions by eliminating common factor c β ).Inverse Laplace transforming S(ν, t) to P(s, t), we obtain (5.4).For the rest of this section, we would like to compare our definition of the finite temperature operator size with Qi-Streicher's definition proposed in Ref. [32].More specifically, Qi and Streicher considered the "size" of the combination χ k (t) √ ρ β in the sense of the "length" of the operator basis, i.e. measured by j c † j c j as in the infinite temperature case.In this definition, the generating function is given as Now, we apply our method to compute the generating function under Qi-Streicher's definition.There are two key differences: (1) our definition involves two OTOC configurations while Qi-Streicher's definition only involves one; (2) as shown in Fig. 4 (b), the two Keldysh contours are placed in different locations, which invokes more general functions h R/A (y, θ), f R/A (x, θ).Nevertheless, it still yields an explicit formula in the large-q SYK model S(ν, t) = .From the first step, where f A ∈ (0, 1), we can see that the size in Qi-Streicher definition ranges from 1 − c β /2 to 1/2, consist with the result in [32].
Eq. (B.9) also has implications on the relation between Qi-Streicher size and complexity.At short time, we can expand f A in terms of λ 1.This implies the time dependent part of size s is proportional to the null momentum of the particle.In the long-time limit, the result shows the operator size is determined by the correlation between two sides on the perturbed geometry.Assuming the complexity-volume conjecture in [48], this implies in AdS 2 we have − log(1/2 − s) ∝ C, where C is the computational complexity of the quantum state.

Figure 1 :
Figure 1: (a) and (b) are typical diagrams for OTOC and S(ν, t) respectively.Wavy lines represent the scrambling modes with propagator λ = C −1 e κt .At late time t ∼ t scrambling when λ ∼ 1, we need to sum over diagrams with multi scrambling modes (wavy lines).

Figure 2 :
Figure 2: (a) Size distribution of the Brownian SYK model at different times, parameterized by the values of λ = C −1 e κt .As time passes (from smaller λ to larger ones), the distribution shifts towards larger s and finally accumulates at s = 1/2; (b) Comparison between theoretical predictions (solid lines) and numerical results (markers) in the Brownian SYK model with q = 4 and N = 500.We plot the moments (2s) m instead of s m so that they saturate to the same value.
and define the operator size as a new type of fermion number in the state size[O] β = O|n β |O O|O ,

Figure 3 :
Figure 3: Plot of P(s, t) for fixed λ = 100, ∆ = 1/4, and three choices of v, to highlight the role of the phase factor in (5.4).

Figure 4 :
Figure 4: The path-integral contour for the generating function of the operator size distribution: (a) the definition in the main text [36] and (b) Qi-Streicher's definition.