title: Complexity growth of operators in the SYK model and in JT gravity
authors: Shao-Kai Jian Brian Swingle Zhuo-Yu Xian
text: 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  equationsection
  
  
      [ Trsgn𝒥𝒢Ø𝒪𝒞ŁℒℰPf#1⟨#1⟩⟨#|1⟨#1||#⟩1|#1⟩⟨#|1⟩#2⟨ #1|.#2⟩#1#2^#1_#2#1|#1|#1(#1)#1[#1]#1{#1} Re Im sgn mod#1 Tr[#1]#1 Det(#1)
          [ ∂>̄⇒→↑↓2FathbbI kΤFathcalTFathcalEFathcalKFathbbZ#1c#1#1{cc #1. Cliff#1 Pf#1mattergauge#1Eq. (<ref>)#1Fig. <ref>eefftrTr SS̃P̃skjian@umd.eduCondensed Matter Theory Center, Department of Physics, University of Maryland, College Park, Maryland 20742, USAbswingle@umd.eduDepartment of Physics, Brandeis University, Waltham, Massachusetts 02453, USA and Condensed Matter Theory Center, Department of Physics, University of Maryland, College Park, Maryland 20742, USAxianzy@itp.ac.cnCAS Key Laboratory of Theoretical Physics, Institute of Theoretical Physics,
              Chinese Academy of Sciences, Beijing, 100190, China
              The concepts of operator size and computational complexity play important roles in the study of quantum chaos and holographic duality because they help characterize the structure of time-evolving Heisenberg operators. It is particularly important to understand how these microscopically defined measures of complexity are related to notions of complexity defined in terms of a dual holographic geometry, such as complexity-volume (CV) duality. Here we study partially entangled thermal states in the Sachdev-Ye-Kitaev (SYK) model and their dual description in terms of operators inserted in the interior of a black hole in Jackiw-Teitelboim (JT) gravity. We compare a microscopic definition of complexity in the SYK model known as K-complexity to calculations using CV duality in JT gravity and find that both quantities show an exponential-to-linear growth behavior. We also calculate the growth of operator size under time evolution and find connections between size and complexity. While the notion of operator size saturates at the scrambling time, our study suggests that complexity, which is well defined in both quantum systems and gravity theories, can serve as a useful measure of operator evolution at both early and late times.
              
              Complexity growth of operators in the SYK model and in JT gravity
                  Zhuo-Yu Xian
                  December 30, 2023
              =================================================================
              
              
              
              § INTRODUCTION
              
              
              Depending on the timescales of interest, several interrelated concepts have recently been proposed to characterize the presence or absence of chaos in quantum dynamics. The basic expectation in a quantum chaotic system is that states in the Schrödinger picture and operators in the Heisenberg picture become more elaborate as time passes. For the purposes of this work, we wish to compare and contrast two ways to quantify this growth: the notion of size and the notion of complexity. We define these notions in detail below, but in brief, size is a measure of how many degrees of freedom are involved in a state or acted on by an operator while complexity refers to the number of elementary steps of some type needed to prepare a state or implement an operator. These two concepts are certainly interrelated in various ways, for example, a certain minimal complexity is required in order for an operator to have large size. In this work, we study and compare precise versions of these notions in two models: the Sachdev-Ye-Kitaev (SYK) model and 2d Jackiw-Teitelboim (JT) gravity.
              
              It is useful to consider two regimes of time, corresponding to times before or after the system has come to approximate global equilibrium. The crossover time between these two regimes, called the scrambling time, will be defined in detail below. Roughly speaking, it refers to the time after which a small perturbation has spread over the entire system.
              
              Prior to the scrambling time, out-of-time order correlation functions (OTOCs) <cit.> characterize the chaotic growth of Heisenberg operators of the form
              
              	O_β(φ) ≡ e^- 1 - 2φ/πβ H/4 O e^- 1 + 2φ/πβ H/4,  φ=θ+iu,
              
              where H is the Hamiltonian, β is the inverse temperature, θ labels the location of the insertion in the imaginary-time evolution, and u=2π t/β is the real time in the unit of β/2π. For a simple operator O which disturbs only a few degrees of freedom, time evolution causes information about this disturbance to spread over the system whenever O is not conserved, [H,O]≠0, a process known as information scrambling <cit.>. For all-to-all chaotic Hamiltonians, information initially spreads exponentially fast with an exponent called a quantum Lyapunov exponent, until it scrambles over the whole system <cit.>, as illustrated in Fig. <ref>. The number of degrees of freedom affected during this scrambling process is measured by the size n of the Heisenberg operator O_β(φ) <cit.>.
              
              
              It will be convenient to translate the language of operators into the language of states using the Choi-Jamiolkowski mapping. The space of operators acting on a Hilbert space ℋ can be mapped to a state in two copies of the Hilbert space ℋ⊗ℋ by O→|O⟩=O⊗1|0⟩ where |0⟩ is a maximally entangled state in the doubled Hilbert space. It is convenient to fix a Hamiltonian H for one copy of the system (e.g. we fix H for the left system, and H^T for the right system which is also denoted as H without confusions) and take |0⟩ = ∑_n|E_n⟩⊗|E_n⟩ where {|E_n⟩} is a basis of energy eigenstates of H. This mapping is appropriate at infinite temperature; it can be extended to finite temperature using the (unnormalized) thermofield double (TFD) state |1_β⟩ =  ∑_n e^-β E_n/2|E_n⟩⊗|E_n⟩ where 1_β=e^-β H/2 corresponds to O=1 in (<ref>). A general operator of the form (<ref>) is mapped to a so-called partially entangled thermal state (PETS) |O_β(θ+iu)⟩ <cit.>.
              The evolution of the operator O_β(φ) = U(t)^† O_β(θ)U(t) is mapped to the evolution of the state |O_β(φ)⟩≡ U(-t)⊗ U(t)|O_β(θ)⟩, where U(t)=e^-iHt and the two copies evolve in opposite directions in time, e.g. with Hamiltonians H and -H. Without the insertion of O,  the state U(-t)⊗ U(t)|1_β⟩= |1_β⟩ is invariant under time evolution. With the insertion of non-conserved O, the state is no longer invariant and the dynamics can be conveniently diagnosed using correlations between the two copies.
              
              In the context of the SYK model made from N fermions ψ^i obeying {ψ^i, ψ^j}= δ^ij with q-body interactions, the operator growth structure is well understood in the large-q limit <cit.> and the conformal limit <cit.>. Here the notion of size n of an operator is the number of elementary operators—the single Majorana operators ψ^j—contained in that operator. The growth of size can be detected by the decay of correlations between the two systems in the PETS state, which is equivalent to a kind of OTOC. More precisely, the size is n[O_β(φ)]=N/2 - i∑_j⟨O_β(φ)|ψ^j⊗ψ^j|O_β(φ)⟩/Z, where the maximally entangled state |0⟩ is defined by (ψ^j⊗1-i1⊗ψ^j) |0⟩=0, ∀ j = 1,...,N and the normalization factor is Z=O_β(φ)|O_β(φ).
              
              The notion of size growth has also been explored in the context of holography. There the TFD state |1_β⟩ is dual to an eternal black hole in AdS space <cit.>. For θ=±π/2, acting a simple operator O on the TFD state |1_β⟩ corresponds to releasing a particle on the asymptotic boundary. The gravity of the black hole forces the particle to fall into the interior of the bulk and affect the near horizon region, which is the holographic bulk counterpart of the boundary growth of size <cit.>. Operator size has also been conjectured to be dual to the momentum of the particle <cit.>. It is therefore interesting to compare the SYK model and JT gravity model, which are closely related in the conformal limit of low temperatures <cit.>. In particular, the size-momentum relation can be studied using SL(2) generators that function as both generators of spacetime transformations and measures of size <cit.>.
              
              
              
              Now consider the situation after the scrambling time. At these longer times, the operator size has reached its equilibrium value, but operators (and states) are still evolving unitarily in the massive many-body Hilbert space. In particular, the complexity of a Heisenberg operator, a quantity borrowed from quantum information theory, is conjectured to continue to grow with time long after the scrambing time <cit.>. We illustrate the evolution after scrambling using a schematic circuit diagram in Fig. <ref>. The most common type of complexity considered in this context is that of circuit complexity, which is defined as the minimal number of elementary quantum gates that are required to produce a target state or operator from a reference. In a chaotic system, it is believed that the circuit complexity of an initial state will grow linearly with time until a time of order the exponential of the system entropy <cit.>.
              
              
              Various definitions of complexity for unitary operators in quantum mechanics and quantum field theories have been proposed <cit.>.
              In this paper, we consider a different notion of complexity for operators, the K-complexity <cit.>, defined through a Krylov basis that is uniquely determined by the evolution Hamiltonian and the reference operator in question. Using the state-operator mapping, we can define this complexity in terms of states or operators; here we focus on PETS as the reference state. The Krylov basis is obtained as follows. In the operator language, one repeatedly applies the Liouvillian map [H,·] to an operator O_β(θ) to generate a basis of operators. In dual state representation, one generates a sequence of new states starting from a PETS, e.g. |[H, O_β(θ)]⟩ = ( H ⊗1 - 1⊗ H) |O_β(θ)⟩. Note that, in the state language, the evolution generated by H ⊗1 - 1⊗ H is called boost evolution due to the interpretation of this transformation as a boost in AdS/CFT. It is one of the SL(2) generators mentioned above. Given this sequence of states (or operators), one then orthonomalizes them to produce the desired Krylov basis,
              |O_β, j(θ)⟩∝ (H ⊗1 - 1⊗ H)^j |O_β(θ)⟩ + ...,
              
              where j denotes the number of applications of the Liouvillian [The Krylov basis generated by the Liouvillian can be incomplete in the Hilbert space ℋ⊗ℋ.], and ... comes from the Gram-Schmidt orthonormalization process. The target state can then be decomposed into the Krylov basis |Ψ_T(t)⟩ = ∑_j i^j ϕ_j(t) |O_β,j(θ)⟩, where ϕ_j denotes the amplitude at j-th basis and the factor of i^j is for convenience. The K-complexity _K is then defined by declaring that basis element j has complexity j,
              _K|O_β,j(θ)⟩=j|O_β,j(θ)⟩,
              
              so the average K-complexity of the target state |Ψ_T(t)⟩ is
              
              	C_K[|Ψ_T(t)⟩]= ⟨Ψ_T(t)|_K|Ψ_T(t)⟩ = ∑_j j |ϕ_j(t)|^2.
              
              
              K-complexity is somewhat analogous to a circuit complexity definition in which the elementary operation is not a unitary but a Hermitian generator of unitary evolution, the Liouvillian or boost generator H ⊗1 - 1⊗ H. The complexity can be defined as an operator by fixing its eigenvectors and eigenvalues in terms of (<ref>). Then one has a notion of average complexity by expanding general states in the Krylov basis (<ref>). Unlike the circuit complexity, K-complexity is uniquely determined by the reference operator and the Hamiltonian, without the introduction of a set of basic gates and a tolerance parameter. This advantage is crucial when considering its gravitational correspondence, since the dictionaries of operator and Hamiltonian are clear, whereas the artificial basic gates and the tolerance parameter are not. We will give more details of the definition K-complexity in the SYK model in the next section.
              
              
              It was shown in <cit.> that the K-complexity grows exponentially with time before the scrambling time for a variety of chaotic Hamiltonians. At these early times, each application of the Liouvillian will increase the operator size by at most a constant amount. As a result, the Krylov basis has a close relation with the operator size at early times. Indeed, it turns out that the K-complexity bounds any properly defined notion of operator size <cit.>. After the scrambling time, the authors of <cit.> conjectured that the K-complexity of a chaotic system will continue to grow, now linearly with time, until it reaches a value that is exponential in the system size. The argument in <cit.> is based on the eigenstate thermalization hypothesis (ETH) applied to the chaotic Hamiltonian of interest <cit.>. Here we explicitly show that the late-time growth of the K-complexity is bounded by a linear function. We also numerically evaluate the K-complexity growth in the SYK model and indeed find linear growth after the scrambling time. Combined with previous results, our results demonstrate that K-complexity shows a exponential-to-linear growth pattern in the SYK model. We conjecture that this holds true for all chaotic systems, namely that after a short “dissipation time” typically set by the local energy scales of the problem, K-complexity shows a universal exponential growth dictated by a Lyapunov exponent which then gives way to a linear growth after the scrambling time. Although we do not see a saturation of the K-complexity at late time due to the limited working precision in the calculation, the late time saturation of K-complexity is seen in a recent preprint <cit.> after we post our results.
              
              To complete the background for our story, we consider the role of complexity in holography. In that context, circuit complexity has received enormous recent attentions for its potential relation to various features of the dual holographic geometry, especially wormholes inside black holes <cit.>. Two conjectures on the duality of the complexity of a state were proposed. The complexity-volume (CV) conjecture states that the complexity is proportional to the volume of the maximal spatial surface (meaning spacetime codimension-one manifold) connecting the boundary of the dual eternal black hole <cit.>. The complexity-action (CA) conjecture states that the complexity is proportional to the action of the Wheeler-DeWitt patch <cit.>. The proposal of holographic complexity inspires new understanding on the complexity in field theories <cit.> and in gravity <cit.>.
              
              In models with partial holographic duals, such as the SYK model, it is interesting to compare microscopic notions of complexity to holographic proposals. For the situation of interest to us, namely PETS states in the SYK model, we consider eternal black holes in JT gravity with a matter field which is dual to the operator O defining the PETS <cit.>. In this model of gravity, the problem of gravitational backreaction of the matter field can be mapped to the motion of particles in a hyperbolic space <cit.>. Focusing on the CV duality for simplicity <cit.>, the holographic complexity becomes proportional to the geodesic distance between the two boundaries of the eternal black hole in JT gravity. Interestingly, this holographic complexity under time evolution also exhibits an exponential-to-linear growth behavior for PETS. Similar behaviors of the complexity appear in the shock wave geometry, which caused by the matter falling into the black hole <cit.>.
              
              Because of the similarities between the microscopic K-complexity in the SYK model and the “coarse grained” complexity defined through the CV conjecture, it is very interesting to compare these two notions of complexity in detail. At early times, both complexities exhibit an exponential growth at a rate set by a quantum Lyapunov exponent. In JT gravity, what appears is the maximal Lyapunov exponent 2π/β, with β is the inverse temperature. In the SYK model, one finds a temperature dependent Lyapunov exponent in the large-q approximation, and in the conformal limit, this exponent approaches the same 2π/β exponent as in the gravity. We find that the K-complexity precisely matches the holographic complexity up to an unimportant constant before the scrambling time. Similarly, both complexities grow linearly after the scrambling time but with different slopes. The slope in the SYK model is set by a microscopic energy scale while the slope in JT gravity is set by the temperature. Nevertheless, both slopes are extensive in the system size and if we consider the rate of complexity growth ratio, i.e.,
              
                  R(t) = d log C(t)/d t,
              
              where C refers to the complexity, then the two notions of complexity turn out to have the same rate up to times that are exponential in system entropy.
              
              To summarize, both size and complexity provide useful windows into quantum chaotic dynamics depending on the timescale of interest. Still, it would be convenient if there is a single quantity that could capture the relevant physics of both quantities. This is not a implausible request since as the size of an operator grows, the complexity is increasing. Indeed, it has been proposed that the growth of operator size is proportional to the rate of increase of its complexity <cit.>. Actually, three seemingly distinct quantities in holography—operator size, complexity, and radial momentum—are proposed to be closely related to each other <cit.>, as schematically shown by the following equation,
              n/β̃∼ P ∼d C(t)/d t,
              
              where 1/β̃ represents an energy scale that sets appropriate units, and n and P refer to the operator size and the momentum, respectively. To verify this relation, we also carry out a calculation of operator size of the PETS in the SYK model and the corresponding SL(2) charge in JT gravity. The results show agreement with both K-complexity and holographic complexity up to the scrambling time. In this sense, the complexity serves as a useful quantity that can capture the dynamics of a simple Heisenberg operator in chaotic systems at both early and late times. While we focused on chaotic systems here, it would be interesting to explore the notion of complexity in integrable systems as well <cit.>.
              
              The rest of this paper is organized as follows. In Sec. <ref>, we explore the dynamics of the K-complexity in the SYK model. We first review the definition of K-complexity, in which the Lanczos coefficient plays an important role. The dynamics of the K-complexity is mapped to a particle moving in one-dimensional lattice made up of the Krylov basis. Then a proof concerning the late-time linear growth of K-complexity is given using a bound on the Lanczos coefficient. We also evaluate the K-complexity in the SYK model, and show that it exhibits an exponential-to-linear growth. In Sec. <ref>, we calculate the holographic complexity of the PETS in JT gravity. The insertion of a simple operator causes a perturbation to the TFD state, which can be mapped to an insertion of a particle moving in the hyperbolic space. The backreaction from the operator insertion is then easily captured at the Schwarzian limit. Then the holographic complexity is measured by the geodesic connecting two asymptotic boundaries of the AdS_2 spacetime. The dynamics of the microscopic K-complexity and the holographic complexity share many similarities, and in certain aspect the K-complexity is a microscopic candidate of the holographic complexity. In Sec. <ref>, the relation between the operator size and the complexity growth rate is considered. We calculate the size of the PETS in both the SYK model and JT gravity. In particular the size is linearly related to the SL(2) charges of AdS_2 spacetime. We also verify that the growth rate of both K-complexity and holographic complexity of the PETS is given by its size in the Lyapunov regime. As a result, the notion of complexity is able to characterize the dynamics of Heisenberg operators in chaotic systems at both short and long times. In Appendix <ref>, we review the finite temperature generalization of K-complexity and generalize the result for non-zero inserting angles θ 0. In Appendix <ref>, we summarize various coordinate systems of AdS_2. In Appendix <ref>, we obtain the generating function of size operator at generic inserting angles. In Appendix <ref>, we discuss the scrambling time of OTOC from the Schwarizan dynamics.
              
              
              
              
              § K-COMPLEXITY IN THE SYK MODEL
              
              
               §.§ Review of K-complexity
              
              
              This section reviews K-complexity associated with the Krylov basis in the SYK model. Unlike an a priori basis, the Krylov basis is uniquely determined by the evolution Hamiltonian and the initial state. As a result, the K-complexity is a natural notion capturing the intrinsic dynamics of the evolution operator, without the ambiguity of choosing an operator basis or elementary gates. We take the inner product in operator space A to be defined by the inner product in ℋ⊗ℋ, namely O|O'=[O^† O'], so the norm is || O || = O|O^1/2.
              
              We work at infinite temperature in this section, which means that the PETS is actually |O_β = 0(θ)⟩, where θ is now irrelevant since the thermal circle is a point with no size. Hence, for notational simplicity, we neglect the subscript β and the angle variable θ, and denote the time evolved PETS by |O_β = 0(φ)⟩ = |O(t)⟩. It is convenient to normalize this operator such that ||O|| = 1. The Heisenberg evolution of an operator is generated by the Liouvillian Ł = [H,·],  i.e., |O(t)⟩ = e^i tŁ|O⟩. The Krylov basis is defined through the Liouvillian superoperator:
                    |O_0⟩ = |O⟩,    b_0=0, 
                    |O_n⟩ = b_n^-1|A_n⟩,   |A_n⟩ = Ł|O_n-1⟩  - b_n-1|O_n-2⟩,    b_n = ||A_n||,     n ≥ 1.
              
              The second line is merely carrying out a Gram-Schmidt procedure on the states |A_n⟩ to produce the states |O_n ⟩. The iteration stops once the Liouvillian fails to generate a linearly independent state. The set of states generated typically span a space of dimension K that is of order K ∼ℋ^2, but they do not always form a complete basis, K ≤ℋ^2. In particular, this happens when the Hamiltonian has conserved charges. For example, the SYK Hamiltonian preserves fermion parity, so the Krylov basis spans the even (odd) fermion parity subspace if one starts with an even (odd) parity reference state.
              
              In terms of the Krylov basis, the Liouvillian superoperator is a simple tridiagonal matrix,
              
              	L_mn = ⟨ O_m| Ł |O_n ⟩ = δ_n,m-1 b_n+1 + δ_n,m+1 b_n.
              
              The coefficients b_n are also known as Lanczos coefficients.
              
              The Heisenberg operator corresponding to O can be decomposed in the Krylov basis, i.e., |O(t)⟩ = ∑_n=0^m i^n ϕ_n(t) |O_n⟩, with ϕ_n real. Unitary evolution implies ∑_n=0^K |ϕ_n|^2 = 1, so ϕ_n can be understood as a wavefunction for the Heisenberg operator in the Krylov basis. The effective Schrödinger equation obeyed by the wavefunciton is
              ∂_t ϕ_n(t) = b_n ϕ_n-1(t) - b_n+1ϕ_n+1(t),   ϕ_0(0) = 1.
              
              This Schrödinger equation (<ref>) effectively describes a quantum particle moving in one dimensional chain in which each lattice site corresponds to an element of the Krylov basis. We note that this mapping from operators to quantum particles shares similar ideas with the mapping from unitary operators to points in a complexity geometry <cit.>.
              
              Now, the K-complexity is defined as a linear operator _K which is diagonal in the Krylov basis and which simply counts the basis elements <cit.>,
              _K|O_n⟩=n|O_n⟩.
              
              Thus, the average K-complexity of the Heisenberg operator O(t) is the average position of the particle moving in the chain, i.e.,
              
              	C_K[O(t)]= ⟨O(t)|_K|O(t)⟩ = ∑_n=0^K n |ϕ_n(t)|^2.
              
              
              The dynamics of K-complexity is governed by the Lanczos coefficient through the Schrödinger equation (<ref>).
              Following Ref. <cit.>, we can build intuition by considering a continuum limit of the Schrödinger equation obtained by introducing a short-range cutoff ϵ with x= ϵ n. Expanding (<ref>) and keeping the lowest-order term in ϵ, we get a continuous version of the Schrödinger equation,
              ∂_t ϕ(x,t) = - v(x) ∂_x ϕ(x,t) - 1/2 v'(x) ϕ(x,t),
              
              where the position-dependent velocity v(x) = 2ϵ b_n captures the information from the Lanczos coefficient. Using a coordinate transformation defined by dy = dx/v(x), the wavefunction changes to ψ(y,t)= √(v(x))ϕ(x,t), and the Schrödinger equation becomes a solvable wave equation,
              
              	(∂_t + ∂_y) ψ(y,t) = 0,   ψ(y,t) = ψ_i(y-t),
              
              where ψ_i(y) is the initial wavefunction at t=0. Note that ψ can be also understood as a wavefunction with normalization 1= 1/ϵ∫ dy |ψ(y,t)|^2.
              
              The average K-complexity of ψ is then
              
              	C_K(t) = ∑_n n |ϕ_n(t)|^2 = 1/ϵ^2∫ dx x |ϕ(x,t)|^2 =  1/ϵ^2∫ dy x v(x) |ϕ(x,t)|^2 = 1/ϵ^2∫ dy x(y) |ψ(y,t)|^2,
              
              where, again, x(y) is determined by the coordinator transformation dy = dx/v(x).
              
              Given a localized initial condition corresponding to the reference state, |ψ_i(y)|^2 = ϵδ(y), the average K-complexity is
              
              	C_K(t) = 1/ϵ∫ dy x(y) δ(y-t) = x(t)/ϵ,
              
              which is fully determined by x(y) or, equivalently, by the velocity through dx/v(x) = dy.
              
              For Lanczos coefficients given by b_n = α n^δ, the velocity is v(x) = 2αϵ (x/ϵ)^δ. Hence, the K-complexity grows as
              
              	C_K(t) ∼ e^2α t,      δ =1
              
              					(2α t)^1/(1-δ),      δ <1
              	.
              
              In particular, for δ = 1, the average K-complexity grows exponentially, while for δ=0, it grows linearly. Ref. <cit.> showed that the Lanczos coefficients are bounded by a linear function when n ≪ N, where N is the system size, implying that the average K-complexity grows at most exponentially up to the scrambling time, α t ∼log N. In the next section, we show that the Lanczos coefficients are bounded by a constant when n ≫ N, and consequently, that the average K-complexity can grow no faster than linearly in time at late times.
              
              
              
              
               §.§ Dynamics of K-complexity in chaotic systems
              
              
              To bound the Lanczos coefficients, it is useful to consider moments of the Liouvillian superoperator,
              μ_2n≡⟨O_0|Ł^2n|O_0⟩,
              
              which are closely related to the following Green function or auto-correlation function,
              
                  G(t)    =   [O^†(0) O(t)]/[O^† O]= ⟨O_0| e^i t Ł|O_0⟩ = ∑_n(it)^2n/(2n)!⟨O_0|Ł^2n|O_0⟩ = ∑_n(it)^2n/(2n)!μ_2n, 
              μ_2n   =   ∫dω/2πω^2n G(ω),      G(ω) = ∫ dt e^i ω t G(t).
              
              Note that the Green function is normalized such that G(t=0) = 1.
              
              Knowing the moments, one can get the Lanczos coefficients using an explicit relation between the two (e.g., see Appendix A of <cit.>). The relation between Lanczos coefficient and the momemt can be obtained also by using saddle-point approximation <cit.>. If the Lanczos coefficients have the smooth form b_n = α n^δ as n →∞, the moment is dominated by
              μ_2n   =    e^2n ∫ dt [H(t) + δlog(2f(t)) ]×α^2n n^2δ n, 
              
                  H(t)    =    - 1+f'(t)/2log[1+f'(t)] - 1-f'(t)/2log[1-f'(t)],   
              where f(t), 0≤ t ≤ 1 is a function uniquely determined by f(0)=f(1)=0, f(t)>0, |f(t)|<1 and f”(t)/1-f'(t)^2 + δ/f(t) = 0. For δ=1, f(t) = sinπ t/π which leads to μ_2n = ( 4α n/π e)^2n <cit.>. For a generic δ>1/2, f(t) can be solved by inverse hypergeometric function. But in any case the prefactor of (<ref>) is given by e^constant× n, so the the moments have the asymptotic behavior
              μ_2n≈α^2n e^2δ n log n + o[n],    n ≫ 1.
              
              where o[n] denotes terms that are at least in the same order of n. While this calculation gives some intuition of the relation between the Lanczos coefficient and the moment, in the following, what we are going to use is the following rigorous bound <cit.>,
              ∏_k=1^n b_k^2 ≤μ_2n≤ C_n max_{b_k} (b_k^2n),
              
              where C_n = (2n)!/n!(n+1)! is the Catalan number. 
              
              Now we show that the Lanczos coefficients are bounded by a constant for n ≫ N. Take the Hilbert space to consist of N Majoranas (with N an even integer) and consider an all-to-all q-body Hamiltonian H = ∑_x h_x such as the SYK model. Each term in the Hamiltonian is taken to be bounded, ||h_x || ≤. We consider the moments and the Lanczos coefficients generated from a simple operator O (for example, a single Majorana operator in the SYK model). Defining l_x = [h_x, ·], the n-th power of the Liouvillian is
              Ł^n |O⟩ = ∑_x_1,...,x_n l_x_n l_x_n-1 ... l_x_1|O⟩.
              
              Each application of l_x_k increases the size of the operator by at most q, so the largest size of l_x_k... l_x_1|O⟩ is k q + 1. Here, the size of a given operator refers to the number of elementary operators, such as a single Majorana operator ψ in the SYK model, contained in that operator. In order to have a nonzero term l_x_k+1 l_x_k ... l_x_1|O⟩, h_x_k+1 and l_x_k ... l_x_1|O⟩ should have a nonvanishing overlap. For each nonvanishing term l_x_k... l_x_1|O⟩, applying Ł will lead to at most 2 (k q +1)N^q-1 nonvanishing terms. As a result, the total number of nonvanishing terms of type Ł^n |O⟩ is bounded,
              ∏_k=1^n 2 ((k-1) q + 1) N^q-1 < (2 q N^q-1)^n n!.
              
              Importantly, the number of nonvanishing terms increases as a factorial of n. Moreover, each individual term l_x_n ... l_x_1|O⟩ is bounded by || l_x_n ... l_x_1|O⟩ || ≤ (2)^n. So the moment is bounded by
              μ_2n = || Ł^n O||^2  ≤ (2)^2n (2 q N^q-1)^2n (n!)^2 < (4 N^q)^2n (n!)^2.
              
              Thus, according to the bound between the moments and the Lanczos coefficients (<ref>),
              ∏_k=1^n b_k^2 ≤μ_2n < (4 N^q)^2n (n!)^2,
              
              the Lanczos coefficients can grow asymptotically 1 ≪ k < N/q at most linearly, i.e., b_k ∝ k^δ, δ≤ 1. This is also consistent with the saddle-point calculation (<ref>).
              
              However, when the size of Ł^n |O⟩ is greater than N, which occurs when n ≥ N/q, we can improve the bound as follows. If l_x_k ... l_x_1|O⟩ has size N, applying Ł can only lead to at most 2 N^q nonvanishing terms. As a result, when n ≥ N/q, the total number of nonvanishing terms in l_x_k ... l_x_1|O⟩ is
              
              	(2 q N^q-1)^N/q (N/q)! (2 N^q)^n-q/N < (N/q)! (2 N^q)^n.
              
              In contrast to the situation when n<N/q, the number of nonvanishing terms now increases at most exponentially with respect to n. The moments are bounded by
              μ_2n≤ (2)^2n ((N/q)!)^2 (2 N^q)^2n = ((N/q)!)^2 (4 N^q  )^2n.
              
              The Lanczos coefficients are bounded by the relation
              
              	b_N/q+1^2 ... b_n^2 ≤μ_2n/b_1^2 ... b_N/q^2 = ((N/q)!)^2 (4 N^q  )^2(N/q)/b_1^2 ... b_N/q^2 (4 N^q  )^2(n-N/q),
              
              which implies that δ = 0 for n ≫ N/q. Combining above results, we have
              
              	b_n ≤λ_L/2 n,         1 ≪ n ≪ N/q 
              λ_C/2 N,         N/q ≪ n ≪ 2^N
              		,
              
              where λ_L is the Lyapunov exponent, and λ_C is a constant independent of n. The factor N in the second line is to capture the system size dependence of b_n at n ≫ N/q, such that λ_C is independent of the system size (see the following). The Lanczos coefficients of n are thus bounded by a linear function in n followed by a plateau that is independent of n. 
              
              Now as mentioned previously, the late-time plateau of Lanczos coefficients was first discussed in <cit.> based on the ETH conjecture. Here, we provided an explicit proof of this plateau behavior of the Lanczos coefficients, which strengthens the results of <cit.>. As we now review, the ETH conjecture is still useful to give an estimate of the plateau value of the Lanczos coefficients <cit.>. We continue to work in a Hilbert space of N Majorana fermions, so the total dimension is 2^N/2.
              Using the Lehmann representation in the energy eigenbasis |E_a ⟩, the Green function and moments read
              
                  G(ω)    =   1/[O^† O]∑_ab 2πδ(ω - (E_a - E_b)) |O_ab|^2,      O_ab≡⟨ E_a | O | E_b ⟩
              μ_2n   =   ∫dω/2πω^2n G(ω) = 1/[O^† O]∑_a,b (E_a - E_b)^2n |O_ab|^2.
              
              According to ETH, the matrix elements O_ab can be approximated by a random matrix to high accuracy in the thermodynamic limit, i.e., O_ab = A(E_a, E_a) δ_ab + A(E_a,E_b) 2^-N/4 R_ab, where A(E_a, E_b) is a smooth function of energies, R_ab denote a random matrix with zero mean and unit variance. If we assume A(E_a, E_b) = A(0,0) F(E_a - E_b) is a function of the energy difference only, then we have
              μ_2n =  2^-N/2∑_a,b (E_a - E_b)^2n|A(E_a, E_b)|^2/∑_c |A(E_c, E_c)|^2 = 2^-N∑_a,b (E_a - E_b)^2n |F(E_a - E_b)|^2 ≈ (N )^2n,
              
              where we have implicitly averaged over the random matrix R_ab. The moment μ_2n is dominated by the largest energy difference between two many-body energy eigenvalues at large n. This implies b_n ≈ N at n ≫ N, namely, the plateau value of Lanczos coefficient is proportional to the system size N.
              
              
              The linear-to-plateau behavior of Lanczos coefficients in turn implies that, in a chaotic system, the average K-complexity of a simple Heisenberg operator as a function of time exhibits an exponential-to-linear growth,
              
              	C_K(t) ≈
              		e^λ_L t,       t_d ≪ t ≪ t_* 
              λ_C N t,       t_* ≪  t
              		
              where λ_L and λ_C are constants, and t_d = λ_L^-1, t_* = λ_L^-1log N/q are the dissipation time and the scrambling time, respectively.
              
              
              
              
               §.§ K-complexity growth of operators in the SYK model
              
              
              We now consider the example of SYK in detail. Our goal is to demonstrate the expectations (<ref>) explicitly. Once again, the SYK Hamiltonian is defined as
              
              	H = i^q/2/q!∑_j_1,..., j_q J_j_1,...,j_qψ^j_1...ψ^j_q,   J_j_1,...,j_q^2 = (q-1)!J^2/N^q-1 = 2^q-1  (q-1)! ^2/q N^q-1.
              
              where the Majorana fermions satisfy ψ^†_j = ψ_j, and {ψ_i, ψ_j} = δ_ij. Exponential growth of K-complexity at early times in the SYK model has been obtained analytically in the large-q limit and numerically by solving the Schwinger-Dyson equation (see Appendix B in <cit.>). At large q, the operator wave function of a single Majorana fermion at time t is
              ϕ_n(t) = √(2/n q)tanh^n  t,      n ≥  1,
              
              leading to the exponential growth of K-complexity at early times,
              
              	C_K(t) = ∑_n=1^∞ n |ϕ_n(t)|^2 =  1/q (cosh 2 t-1),
              
              where the summation over the basis can be extended to infinity because we work at finite time with N →∞. The exponential growth exponent is 2, consistent with the Lyapunov exponent at infinite temperature.
              
              In getting (<ref>), the large-N Wightman correlation function is used, but this only works for the K-complexity before the scrambling time. After the scrambling time, we expect a linear growth of K-complexity. To verify this conjecture, we choose O = √(2)ψ_1, and calculate the Lanczos coefficients numerically. For practical purposes, we can truncate the Krylov space at some n_max≫ 1 and still capture the dynamics for a finite time related to n_max. In Fig. <ref>, we plot the Lanczos coefficients for different choices of N, with all showing plateau behavior in the regime n> N/q. Moreover, the value of the plateau is proportional to the system size N, which is shown in Fig. <ref>. (At exponentially large n ∼ e^N, the Lanczos coefficient decreases and eventually vanishes because the dimension of Krylov basis is bounded by that of the operator Hilbert space <cit.>.) The K-complexity and its time derivative are calculated by solving the Schrödinger equation (<ref>) for q=4 and N=30, with the results shown in Fig. <ref> and <ref>. At late times, t ≫ t_* = λ_L^-1log N/q, the K-complexity grows linearly as expected. Due to the small number of qubits we simulate, N/2=15, the scrambling time is quite small, so the early time exponential regime is not manifest.
              
              
              
              So far, we considered only the infinite temperature case. One way to generalize K-complexity to finite temperature is to consider the corresponding PETS |O_β⟩≡|O_β(θ=0)⟩ at temperature 1/β. We focus on θ = 0, and the case of nonzero inserting angles θ 0 is considered in Appendix <ref>. This generalization is essentially equivalent to a change of definition of the inner product,
              ⟨ O_1 | O_2 ⟩_β = [e^- β H/2 O_1^† e^- β H/2 O_2 ].
              
              Ref. <cit.> obtained an analytic result at large q prior to the scrambling time. We briefly review their analytic results at short times (see Appendix <ref>, and also the Appendix B in <cit.> for more details), and present a numerical evaluation valid at late times.
              
              
              The moments are now related to the finite temperature Wightman correlation function,
              
                  G(t) = [ρ^1/2O ρ^1/2O(t)]/[ρ^1/2O ρ^1/2O],     ρ = e^-β H,     μ_2n = i^2nd^2n/dt^2n G(t)|_t=0,
              
              Using the large-q Wightman correlation function at finite temperature β >0 (Appendix <ref>), we get the generalized K-complexity at early times,
              
              	C_K(t) = ∑_n=1^∞ n |ϕ_n(t)|^2 =  1/q (cosh 2α t-1),
              	    α = cosαβ/2→,      β≪ 1 
              π/β ,      β≫ 1 .
              
              The exponential growth rate is given by 2α, which is equal to the Lyapunov exponent λ_L = 2α at large q <cit.>.
              
              The analytic wavefunction from the large-q Wightman correlation function also allows us to compute moments of _K operator. To do that, we can introduce the generating function of K-complexity, i.e.,
              e^μ_K   =    1+ 4/qlogα t + ∑_n=1^∞ e^μ n2/n qtanh^2nα t = (1+(1-e^μ) sinh^2 α t )^-2/q.
              
              The n-th moment is obtained by taking n-th derivative of the generating function, but we restrict ourselves to the average K-complexity. It would be interesting to explore the holographic duality of such an generating function in the future.
              
              We also obtain the Lanczos coefficient at large n ≫ N, as shown in Fig. <ref>, which shows a similar linear-to-plateau pattern. The slope of the linear function in the Lanczos coefficient gets smaller at lower temperature, because the dynamics is slower, reflecting the decrease of the Lyapunov exponent λ_L = 2α with temperature. The ETH estimate applied to the operator O_β leads to (O_β)_ab = A(E_a, E_a)  δ_ab + A(E_a,E_b) 2^-N/2 R_ab with A(E_a,E_b)=A(0,0)F(E_a-E_b)e^-β(E_a+E_b)/2. Similar to (<ref>), μ_2n∼ (Nℰ)^2n for n≫ N. As a result the plateau value of the Lanczos coefficient remains unaffected. The time derivative of K-complexity growth for various temperatures is plotted in Fig. <ref>, where the slope of the late-time linear growth is independent from the temperature, and the early-time exponential growth region expands due to the decrease of the Lyapunov exponent at finite temperatures.
              
              
              
              
              
              
              § HOLOGRAPHIC COMPLEXITY IN THE JT GRAVITY
              
              
               §.§ Partially entangled thermal state in JT gravity
              
              
              We now carry out corresponding calculations in JT gravity using complexity-volume duality. First, recall the setup of JT gravity. For a holographic conformal field theory, the bulk representation of the PETS state is a half disk with an operator inserted on its boundary <cit.>. We assume that the dual bulk theory is JT gravity <cit.> with a free matter field χ (the bulk field dual to operator O) coupling only to the metric,
              
              I    =    I_ bdy[g, ϕ] + I_M[g,χ],
              
              I_ bdy[g, ϕ]    =    -ϕ_0/16π G_N∫√(g) R+2∫_∂ K-1/16π G_N∫ d^2x ϕ√(g)(R+2)
              +2∫_∂ dx√(h)ϕ_b K,
              where G_N is the Newton's constant. g and h refer to the determinant of metric g_μν and of induced metric. R, K, ϕ, and ϕ_b denote the scalar curvature, the extrinsic curvature, the dilaton field and its value on the boundary, respectively. We require that the constant ϕ_0≫ϕ. 
              The boundary conditions are
              
                  h = 1/ϵ^2dτ^2,     ϕ_b = ϕ_r/ϵ,
              
              where ϕ_r is chosen to be a constant and τ is the imaginary time of the boundary theory.
              
              The first term of (<ref>) is purely topological and gives the residual entropy S_0=ϕ_0/4G_N. The second term in the bulk gives R+2=0 after the dilaton field ϕ is integrated out. So the metric in the bulk is localized to EAdS_2, ds^2 = (dτ̃^2 + dz^2)/z^2 (in this section, we work in Euclidean signature). Thus, the remaining dynamics of the metric is on the boundary and is governed by the last term of (<ref>). This can be effectively reduced to the reparametrization of the boundary time <cit.>τ̃=f(τ),   z=ϵ f'(τ)+o[ϵ^3],
              
              which automatically satisfy the boundary condition of metric in (<ref>). The time reparametrization field f(τ) is governed by the Schwarzian action <cit.>,
              
              I_ bdy[f] 
                 =   -Q∫ dτ√(h) (K-1)=-Q2π-  ∫ dx^2 √(g)R/2-∫ dτ√(h) 
              = -Q(2π+A -L)
                 ≈    - ϵ Q ∫_0^β dτ  Sch(f(τ), τ), 
                    where  
              Q= ϕ_b/8 π G_N,  
              L = β/ϵ,   Sch(f(τ), τ) = - 1/2(f”/f')^2 +( f”/f')'.
              
              We have used the Gauss-Bonnet theorem and R=-2 in the first line. A is the area enclosed by the boundary. L is the length of the boundary.
              
              We assume that the inserted operator O is a single trace operator with scaling dimension Δ and the dual matter field χ vanishes in vacuum. We only consider the case of Δ>0 which agrees with the situation in the SYK model. The dimensionless inner product of the PETS becomes <cit.>      ϵ^2Δ O_β (θ) | O_β (θ)/1_β|1_β
                = ϵ^2Δ O(τ) Oτ'_β 
                = ϵ^2Δδ^2/δχ_r(τ) δχ_r(τ')∫ Df e^-I_ bdy[f]- I_M^[f,χ]
                 =   ∫ Df (  ϵ^2 f'(τ) f'(τ')/[f(τ)-f(τ')]^2)^Δ e^- I_ bdy[f]≈∫ Df 2cosh D(X,X')^-Δ e^-I_bdy[f], 
                 τ   = β/4(1+ 2θ/π),    τ' = -β/4(1 - 2θ/π),
              
              where X(X') denotes the point on the boundary with the boundary time τ(τ'), and D(X,X') denotes the geodesic distance between the two points X and X'. In getting the above equation, we have used the effective action of the matter field <cit.> and the approximation of geodesic distance under the reparameterization (<ref>)
              
              I_M^[f,χ]=-∫ dτ̃_1dτ̃_2χ̃_r(τ̃_1)χ̃_r(τ̃_2)/|τ̃_1-τ̃_2|^2Δ
              ,  χ_r(τ)=f'(τ)^1-Δχ̃_r(f(τ))
              
              cosh D(X(τ_1),X(τ_2))
              =(τ̃_1-τ̃_2)^2+z_1^2+z_2^2/2z_1z_2
              =f(τ_1)-f(τ_2)^2/2ϵ^2 f'(τ_1)f'(τ_2)+ o[ϵ^0].
               
              where we have normalized the prefactor of the correlation function in the Poincare coordinate, and χ_r(τ) is the source of O(τ).
              
              Combining the gravity part and the matter part, we obtain the effective action up to some constants
              
                  I = -Q(A-L) + μ L_μ,  μ=Δ,     L_μ≈ D(X,X').
              
              Equivalently, the first part describes a massive particle with charge Q moving in a hyperbolic space, where A is the area enclosed by the world line and L is the length of the world line <cit.>. We call it the boundary particle because it locates the boundary. The second part describes a neutral particle, where L_μ and μ are the world line length and the mass of the inserted particle, respectively. We call it the inserted particle because it describes the operator insertion. So the degree of freedom are the trajectories of the boundary particles and the inserted particle, where the total length of the world line of the boundary particles L is fixed to be β/ϵ.
              
              We also need to connect the world lines between the boundary particle, and the inserted particle according to the inserted position θ. The dimensionless parameters in the problem are Q,L,μ,θ, where Q,L,μ are measured in units of the AdS radius.
              
              Throughout we assume a low-energy limit and a classical limit Q≫ L≫ 1. To minimize the action, the classical solution tends to have a large A and a small L_μ. One can have an intuitive picture of the solution in Fig. <ref>. We will discuss quantitatively below. 
              
              
              
              
              To prepare the PETS, it is convenient to use an embeddeding space, -(Y^-1)^2+(Y^0)^2+(Y^1)^2 = -1 with the metric ds^2 = -(dY^-1)^2+(dY^0)^2+(dY^1)^2. Other coordinate systems are summarized in Appendix <ref>. 
              
              The action (<ref>) is invariant under the an overall SL(2) transformation on the left boundary particle, the right boundary particle, and the inserted particle. The SL(2) charges of the three particles are denoted as Z_L^a,Z_R^a,Z_μ^a respectively.
              Their trajectories Y_L^a,Y_R^a,Y_μ^a  are determined by their own SL(2) charges via <cit.>
              	Z_R · Y_R = -Q,    Z_L· Y_L = Q,    Z_μ· Y_μ=0.
              
              Actually, the overall SL(2) transformation is a gauge redundancy due to the killing symmetries of the background AdS_2. So the total charges vanish
              
                  Z_R^a+Z_L^a+Z_μ^a = 0,
              
              which controls the interactions between the three particles.
              In the rest frame of the inserted particle, we make the following ansatz for the SL(2) charges,
              
              	Z_i^a = Q/cosh r_i (s_icoshρ_i, 0, -sinhρ_i),  
              	Z_μ^a = (0,0,-μ),   s_L,R=-1,1,   i=L,R.
              
              where r_i and ρ_i are constants determined by the equation of motion.
              From (<ref>), the trajectories of the two boundary particles are
              
                  Y_i(φ) = e^-s_i ρ_i T_2 (cosh r_i, sinh r_i sinφ, sinh r_i cosφ)^T,
                     T_2 = ( cccc 0     0     1 
               0     0     0 
               1     0     0 ),
              
              which are parameterized by the angle φ in the Rindler coordinates.
              
              Schematic trajectories for these particles are shown in Fig. <ref>: the two segments are world lines of the left and right boundary particles with centers at H_L and H_R, respectively. The red line is the inserted particle. The segments are associated with angles θ_i which satisfy
              tanhρ_i = - tanh r_i cosθ_i/2,    i = L, R,
              
              because the three trajectories join at X and X'. The length of the world lines of two boundary particles is fixed by
              θ_i sinh r_i = L/2 1+s_i 2θ/π,
              
              which also enforces the total length to be L. With Eqs. (<ref>), (<ref>) and (<ref>), we can solve for ρ_i,r_i,θ_i in terms of the parameters Q,L,μ,θ.
              
              The dilaton field can be determined from its equation of motion. At point Y, one has
              ϕ/8π G_N = 
                  Z_L · Y,     Y^1<0 
              
                  -Z_R · Y,     Y^1>0
                  .
              
              Figure <ref> shows the configuration of the dilaton field for the PETS in global coordinates for AdS_2. When the centers H_L and H_R are on the both sides of the trajectory of the inserted particle in Euclidean AdS_2, as shown in Fig. <ref>, the dilaton field reaches extremal values at the two points. We call them horizons H_L, H_R although only the one with the smaller value of dilaton is the true horizon <cit.>. Each of the horizons extends along a light cone in Lorentzian AdS_2, as shown in Fig. <ref>. We call the left(right)-going light cone of the left horizon as the left-outer(inner) horizon, and the right(left)-going light cone of the right horizon as the right-outer(inner) horizon.
              
              
              Using the conservation law (<ref>) and the constraints (<ref>,<ref>), the SL(2) charges of the boundary particles can be obtained numerically.
              One can also can get analytic solutions perturbatively in μ L/Q≪ 1,
              
              	r_i = r - 2+(π + 2θ s_i)tanθ/4π^2μ L/Q,
              
              where sinh r = L/2π is the unperturbed radius.
              
              However, at θ=π/2, the expansion (<ref>) breaks down. To get a meaningful result, one can introduce a regularization, i.e., θ = π/2 - δ, and expand first in δ→ 0 and then in μ L^2/Q≪ 1. The result is
              
                  r_R =r  ,      r_L = r -1/4π^2μ L^2/Q.
              
              In this case, the expansion parameter is μ L^2/Q which differs from the parameter μ L/Q at generic θ. This reflects the non-commutativity of the two expansions in the δ = 0 limit.
              
              
              
               §.§ Holographic complexity growth of the Heisenberg operator
              
              
              We define the holographic complexity C_V(t) of the Heisenberg operator as the holographic complexity of the corresponding PETS. In two-dimensional spacetime, the CV conjecture <cit.> states that the holographic complexity is proportional to the geodesic distance D between the boundary points at times t_L and t_R in Lorentzian signature. It is approximated by <cit.>
              	C_V(t_L,t_R) ≈ϕ_0/G_N D Y_Lπ-i u_L, Y_Ri u_R,  
              	u_i=2π t_i/β_i,  β_i=π+2θ s_i/θ_iβ,
              
              where ϕ_0 dominates the cross-section. The geodesic distance D(Y_1,Y_2) between points Y_1 and Y_2 can be evaluated by the inner product in the embeddeding space,
              cosh(D( Y_1, Y_2))=- Y_1 · Y_2. Finally, letting C̃_V=G_N/ϕ_0 C_V, we have
              coshC̃_V(t_L,t_R)
                 =   cosh(ρ _L+ρ _R) cosh r_L cosh r_R +sinh r_L sinh r_R cosh u_L cosh u_R 
                     - sinh(ρ _L+ρ _R) sinh r_L cosh r_R  cosh u_L + cosh r_L sinh r_R cosh u_R 
                     + sinh r_L sinh r_R sinh u_L sinh u_R.
              
              
              Consider first the case with -t_L=t_R=t. The choice of opposing directions of time evolution corresponds to the Heisenberg evolution of operators, U(t)^† O_β(θ) U(t).
              Without the inserted particle, μ=0, the bulk is unperturbed and recovers the Rinder patch of AdS_2, ρ_L + ρ_R  = 0, r_L = r_R = r and -u_L = u_R = 2π t/β.
              The complexity is C_V[1_β]=2ϕ_0r/G_N, which is independent of time because of the boost symmetry of the TFD state.
              
              In the light operator limit, μ L ≪ Q, and considering |θ|≠π/2, the perturbed solution (<ref>) can be used to get the geodesic length,
              coshC̃_V(-t,t)    =   e^2r/2( 1 + πθcosh u - (2 + 2 θtanθ) /2π^2μ L/Q),       t_d ≪ t ≪ t_∗
              1/2 (θ/8πμ L/Qe^r+u)^2  ,       t≫ t_* ,  
              	     u=2π/βt
              
              where r ∼log L, and the dissipation time and the scrambling time are, respectively,
              
              t_d=β/2π,   t_*=β/2πlog8π Qcosθ/μ L.
              
              At late times, we only need consider the leading time dependence in (<ref>). The complexity grows exponentially at early times and linearly at late times,
              C̃_V(-t,t) ≈
              		2r + θ/2πμ L/Qcosh u,       t_d ≪ t ≪ t_∗
              
              	    2log(θ/8πμ L/Q) + 2r + 2u,         t≫ t_*,
              		
              The early time exponential growth has Lyapunov exponent λ_L = 2π/β, and at late time the linear growth rate of C̃_V is λ_C= 4π/β.
              
              In the θ = π/2 case, we use instead the perturbed solution (<ref>), and expand the geodesic length in the limit δ→ 0 and μ L^2 ≪ Q,
              coshC̃_V(-t,t)    =   e^2r/2( 1 + cosh u -1/4π^2μ L^2/Q),        t_d ≪ t ≪ t_∗
              1/2(1/2π1/8 πμ L^2/Q e^r+u)^2,       t≫ t_*
              	    ,     u=2π/βt,
              
              leading to the following complexity dynamics,
              C̃_V(-t,t) ≈  2r + 1/4π^2μ L^2/Qcosh u,       t_d ≪ t ≪ t_∗
              
                  2log( 1/2π1/ 8πμ L^2/Q) + 2r  +
                  2u,        t≫ t_*
                  .
              
              In particular, the qualitative behavior is not modified by the change in perturbation parameter and the Lyapunov exponent and late time exponent are still given by λ_L = 2π/β and λ_C= 4π/β, respectively. It is also interesting to note that if one regularizes the infinity at π/2 by π/2→L/2π, then (<ref>) can include the case θ→π/2 as well.
              
              
              
              These analytic results are also consistent with a numerical evaluation of the complexity as shown in Fig. <ref>. The numerical solution indeed shows that the complexity grows exponentially at first and linearly after the scrambling time.
              
              
              Now, when H_L, H_R are on both sides of the trajectory of the inserted particle, the geodesic giving the complexity crosses the left/right-inner/outer horizons at points H_L^+, H_L^-, H_R^-, H_R^+ from left to right, as shown in Fig. <ref>. Combining with the two points Y_L, Y_R on the boundaries, they divide the complexity geodesic into five intervals. We plot the contributions of these intervals to the complexity in Fig. <ref> in the light operator limit. The growth of complexity is mainly due to the growth of the length of the geodesic distances inside the black hole interiors D(H_L^+,H_L^-) and D(H_R^-,H_R^+). Note in particular the similarity between Fig. <ref> and the K-complexity in Figs. <ref>, <ref>.
              
              To compare with the SYK model, we identify ϵ = 1/. Then the parameters of the JT gravity and the SYK model are related by
              
              Q= α_S N,   L = β,  μ = 1/q,  ϕ_0/4G_N = s N,
              α_S and s are some numerical constants. At the large q limit, α_S=1/4q^2 and s=1/2log 2<cit.>. Based on (<ref>), if the boundary length segments are subtracted from the holographic complexity, we get the following equality in the conformal limit,
              Δ C_V(t) /L =  α_C/q (cosh2π t/β-1) = α_C C_K(t),   α_C = 2s/πα_S,    t_d < t < t_∗.
              
              where we have set θ = 0 for comparison.
              (Actually we verify in Appendix <ref> that for a nonzero θ, the equality still holds.) 
              Note that the scrambling time for K-complexity,  β/2πlogN/q/α≈ t_∗( 1+ O[(log N)^-1] ), is approximately equal to the scrambling time of the holographic complexity in the large N limit. It is also interesting to compare the rate of complexity growth for both K-complexity and holographic complexity. In the conformal limit, the following equation hold after the scrambling time,
              d logΔ C_V(-t,t)/d t = d log C_K(t)/d t,    t> t_d.
              
              At this level, the K-complexity defined in (<ref>) thus gives a microscopic counterpart of holographic complexity.
              
              
              
              
              If two sides are evolved in the same time direction, then the complexity will grow linearly at first since the complexity of the time evolution operator U will dominate over the simple operator O. The geodesic length and complexity are
              cosh D(t,t)    =   e^2r/2cosh^2 u (1 + πθ u - (2 + 2 θtanθ) /2π^2μ L/Q), 
              C̃_V (t,t)    ≈    2r + 2u,    t ≫ t_d,    u = 2π t/β.
              
              At late times, this growth is equivalent to the late-time linear growth of the complexity in simple Heisenberg operators. This is because the Heisenberg operator becomes complicated after the scrambling time, and the Heisenberg evolutions at two sides become mostly uncorrelated.
              
              Another interesting limit is the heavy operator limit. When μ=2Q, the two segments in Fig. <ref> become two tangent thermal circles now with inverse temperatures β_i=π+2θ s_i/2πβ, respectively. The operator O is already complicated initially, and there is not a large separation between the scrambling time and the dissipation time, as shown in Fig. <ref>. In this limit, the complexity approximately decomposes into the complexities of two wormholes separately,
              C̃_V(t_L,t_R)
                 ≈    2logβ_L/πϵcoshπ t_L/β_L + 2logβ_R/πϵcoshπ t_R/β_R
                 =   C̃_V[e^-1/2β_L+it_L H] + C̃_V[e^-1/2β_R+it_R H]. 
              This form exhibits quadratic growth when t≪ t_d and linearly growth when t≫ t_d.
              The decomposition reflects the fact that the heavy operator O effectively cuts the wormhole into two shorter wormholes while creating a large interior, similar to a Python's lunch geometry <cit.>, as shown in Fig. <ref>. The geodesic length at two wormholes grows with time independently with their own inverse temperature β_i.
              
              
              
              
              § THE SIZE OF THE PARTIALLY ENTANGLED THERMAL STATE
              
              
               §.§ The size from the SYK model
              
              
              In this section we obtain some results for operator size and compare them to the preceding complexity results. The maximally entanglement state |0⟩ defined in the doubled SYK Hilbert space satisfies c^j|0⟩=0, ∀ j where c^j=(ψ_L^j+iψ_R^j)/√(2). The size of an operator ψ in the SYK model can be defined as n[ψ]=⟨ψ|n̂|ψ⟩/ψ|ψ, where the size operator is n̂=∑_j=1^N (c^j)^† c^j= N/2+i∑_jψ^j_Lψ^j_R<cit.>. At large q and early time, the growth of size of operator O_β(θ+iu)=√(2)ψ^1_β(θ+iu) is characterized by
              Δñ(t)=n[O_β(θ+iu)]-n[1_β]/δ_β
              =(π  v/2)  (θ  v) cosh (2π v/βt) -2 tan(π  v/2) (θ  v tan (θ  v)+1)/π  v+2 (π  v/2),
              
              where v=αβ/π, u=2π t/β, and the normalization factor δ_β=2G(β/2) is determined by Δñ(0)=1 when θ=π/2. At late time, the exponential growth will slow and eventually vanish as the size approaches Δñ(∞)=N/2. More generally, one can obtain the generating function of size operator at arbitrary inserting angle θ, and it turns out the generating function of size operator agrees with that of K-complexity <ref>.
              
              Notice that, in the large q limit, the scrambling time obtained from both C_K and C_V is t_*∼λ_L^-1ln(N/q) rather than λ_L^-1ln N. So Δñ(t_*)=N/q rather than its saturation value N/2, which implies that Δñ(t_*) deviates from exponential growth before saturation. We will discuss this point in Appendix <ref>.
              
              
              
              
              
              
               §.§ The size derived from JT gravity
              
              
              
              The size is related to the symmetries of AdS_2 in JT gravity. The SL(2) generators B̃,Ẽ,P̃ in the dual AdS_2 are related to the operators in the two sites SYK model <cit.>B̂   =   β/2πH_R-H_L,
              Ê   =   β/2πH_R+H_L+μ̃n̂ - ⟨1_β| (H_R+H_L+μ̃n̂) |1_β⟩, 
              P̂   =   -i[B̂,Ê],
                    where  μ̃/𝒥=2α_S/Δδ_β2π/β𝒥^2.
              
              For conciseness, we consider that the states are normalized. The normalized change of size can be written as
              Δñ(t)
                 =   1/δ_β⟨ O_β(θ+iu)|n̂|O_β(θ+iu)⟩ - ⟨1_β|n̂|1_β⟩
                 =   1/μ̃δ_β2π/β⟨ O_β(θ+iu)|(Ê-B̂-β/πH_L)|O_β(θ+iu)⟩ + β/π⟨1_β|H_L|1_β⟩
                 =   Δ𝒥/2α_Sβ/2π⟨ O_β(θ+iu)|(Ê-B̂)|O_β(θ+iu)⟩ + ∂_φln.Gβ/21-2θ/π;β1+φ/π|_φ=0,
              
              where G(τ;β)=[e^-β HO^†(τ)O]/[e^-β H].
              Its time derivative is proportional to the momentum
              ∂_t Δñ(t)=Δ𝒥/2α_S⟨O_β(θ+iu)|P̂|O_β(θ+iu)⟩.
              
              According <cit.>, when we insert operators to the thermal circle at π/2-θ+iu and -π/2+θ+iu, the generators in the semiclassical limit are
              ⟨O_β(θ+iu)|B̂,P̂,Ê|O_β(θ+iu)⟩/O_β(θ+iu)|O_β(θ+iu) 
              =Δ/cosθsinθ,sinh u,cosh u
              At large q limit, we obtain
              Δñ(t)= β𝒥/π^2 (-2-πtanθ +πθcosh2π t/β),
              
              which is equal to (<ref>) from the SYK model at β𝒥≫1 limit.
              
              
              
               §.§ Relation between the operator size and the complexity
              
              
              We first discuss the relation between K-complexity and the size operator. At infinite temperature, the time derivative of K-complexity is proportional to the size at early times t_d ≪ t ≪ t_∗. Relating (<ref>) and (<ref>), we find that
              1/2d C_K(t)/dt≈1/qΔñ(t),   β→ 0.
               
              Actually, the generating functions of K-complexity and size operator agree, as shown in Appendix <ref>. For finite temperature and θ=0, the relation will be modified by a temperature dependent factor,
              1/2d C_K(t)/dt≈1/qcos^2 π v/2Δñ(t).
              
              This relation may be extended to late times after the scrambling time.
              The linear growth of K-complexity at late times is proportional to the system size, i.e., 1/2d C_K(t)/dt∝N/2 = Δ n(t) for t ≫ t_∗.
              
              Now we consider the relation between the holographic complexity and the size operator.
              The size of an operators is linearly related to its out of time order correlator (OTOC) with Majorana fermions ψ_i<cit.>. From the geometric interpretation of effective theory <cit.>, we find the following relation between the size in the SYK model and the complexity in JT gravity at the limit q≫1, N≫β𝒥≫1 and the early time,
              π ^2/β𝒥Δñ[ψ^1_β(θ+iu)]/N/2
              =-2-2 θtanθ +πθcosh t
              =2π^2 Q/μ LG_N/ϕ_0(C_V[ψ^1_β(θ+iu)]-C_V[1_β])
              
              which is valid under the dictionary (<ref>). Combining it with the Epidemic relation
              d/dtΔñ = λ_L Δñ at the Lyapunov regime <cit.>, we find
              1/T S_0d C_V/dt = Δñ/N/q,
              
              where entropy S_0=ϕ_0/4G_N= sN and temperature T=1/β.
              
              
              
              § CONCLUSION AND OUTLOOK
              
              
              We calculated the complexity of a Heisenberg operator in both the SYK model and JT gravity. In the SYK model, we used the notion of K-complexity defined through the Krylov basis. In the JT gravity model, we used the CV conjecture to define the complexity. The simplicity of JT gravity allowed us to treat the problem of gravitational back-reaction by mapping it to motions of particles in a rigid hyperbolic space. We found that both complexities show an exponential-to-linear growth behavior. In particular, the two notions of complexity actually match up to a constant before the scrambling time. After the scrambling time, although the characteristic energy scales for the two complexities are different, they both show a linear growth with a slope proportional to system size. We also verified the relation between the complexity growth and the operator size before the scrambling time. The complexity can be used to capture the quantum dynamics at both short and long times.
              
              It is worth noting the temperature dependence of the holographic complexity. For a generic insertion θ±π/2, the complexity is inversely proportional to the temperature, C_V ∝ (β/ϵ) e^2π t/β before the scrambling time, and it is proportional to the temperature in the linear growth regime C_V ∝ (ϕ_0/G_N) 4π t/β after the scrambling time. The former is due to the fact that the inserted operator is perturbing the wormhole with temperature 1/β, while the latter is due to the relation between the Lorentzian time and the Rindler time.
              
              By contrast, the role of temperature is less clear in the context of computational complexity, so it is interesting to attempt to introduce temperature into the definition of computational complexity. For example, the temperature dependence of the holographic complexity means a simple identification of the circuit time and the Lorenzian time is not enough to set the complexities equal up to an overall constant. In our study, we generalized the K-complexity to finite temperature by considering the PETS at temperature 1/β as the reference state. This is essentially the same as a temperature-dependent inner product <cit.>. However, as we showed based on (<ref>), such a generalization does not give a completely consistent identification between the two notions of complexity after the scrambling time. By adjusting the overall normalization of one or the other, one could match the early or late time growth but not both. Hence, it is interesting to consider further refinements that might produce even more harmony between the two notions of complexity. Nevertheless, we emphasize that at the level of the rate of complexity growth, the identification between the notions works perfectly well at early and late times as in (<ref>). This suggests that as far as the temperature dependence is concerned, the rate of complexity growth ratio and the related time scales have simpler holographic interpretation than the absolute value of complexity itself.
              
              A possible reason for the above mismatch is the difference between the reference states of K-complexity and holographic complexity. The reference state of K-complexity is the PETS state |O_β⟩; the reference state of holographic complexity is the maximally entangle state |0⟩ <cit.>. It is worth developing an algorithm of complexity which is independent from the choice of initial operator O. It will also benefit the generalization on the complexity of multi operators, which is related to multi shock waves geometries <cit.>.
              
              The CV proposal used in this paper only depends on the geometry and the dilaton. It is an open question that whether matter fields should have a direct contribution to the holographic complexity besides their indirect contribution via back-reaction on the metric. The answer to this question may be crucial for the complexity of heavy operators, such as (<ref>). From the perspective of complexity-action (CA) conjecture <cit.>, the action of the matter field along the trajectory of the inserted particle can directly contribute. We hope to explore this problem in the future.
              
              It is also interesting to consider higher dimensional generalization of PETS and its gravity dual. For instance, in three dimensions, we may insert an end of world brane behind the horizon of an eternal black hole, corresponding to the geometry worked out in <cit.> in the context of an evaporating black hole. The holographic complexity in this case is then proportional to the volume of a two-dimensional maximal surface connecting the two boundaries. From the viewpoint of evaporating black holes, the geometry of the PETS in our study is effectively dual to a two-dimensional version of the entangled system consisted of the black hole and the auxiliary radiation <cit.>. And the holographic complexity calculated here is the so-called unrestricted complexity for decoding the radiation <cit.>.
              
              
              
              § ACKNOWLEDGEMENT
              
              We thank Leonard Susskind, Anatoly Dymarsky, Matteo Carrega and Yu-Sen An for helpful discussions. S. K. J. and B. S. are supported by the Simons Foundation via the It From Qubit Collaboration. Z. Y. X. is supported in part by the Natural Science Foundation of China under Grants No. 11875053 and No. 12075298 and by the National Postdoctoral Program for Innovative Talents BX20180318, funded by China Postdoctoral Science Foundation.
              
              
              
              
              § K-COMPLEXITY IN THE SYK MODEL AT EARLY TIMES
              
              
              We will summarize the Wightman correlation function and the K-complexity at early times (e.g. see Appendix B in <cit.> for more details).
              Using the large-q approximation, the imaginary time correlation function at temperature β is given by
              
              	2 ⟨𝒯_τψ(τ)  ψ(0) ⟩ = 1+ 2/qlogα/ |cosα (τ - β/2)|,   τ>0,   α = cosαβ/2.
              
              So, the Wightman correlation function is
              
                  G(t) = 1+ 2/qlog1/coshα t.
              
              where we have properly normalize it by G(0) = 1. And accordingly the wavefunction of a simple Majorana fermion is <cit.>ϕ_n(t) =  1+2/qlog1/cosh(α t),    n=0 
              √(2/nq)tanh^n α t,       n ≥ 1. 
              	
              This leads to the exponential growth of K-complexity at early time,
              
              	C_K(t) = ∑_n=1^∞ n |ϕ_n(t)|^2 =  1/q (cosh 2α t-1),   α = cosαβ/2.
              
              
              We also generalize the analytical result of the K-complexity in the SYK model to non-zero inserting angles θ 0. For a nonvanishing θ 0, the diagonal component of the Liouvillian superoperator in Krylov basis is not zero,
              
                  L_mn =  δ_n,m-1 b_n+1 + δ_n,m+1 b_n + a_n δ_n,m,    a_n = ⟨ [H,O_n] | O_n ⟩.
              
              This change will lead to a slightly different Schrödinger equation for the wave function O(t) = ∑_n i^n ϕ_n(t) O_n,
              ∂_t ϕ_n = b_n ϕ_n-1 - b_n+1ϕ_n+1 + i a_n ϕ_n.
              
              The Lanczos coefficient can be mapped to a Toda chain flow and solved by a coupled differential equation <cit.>. A general solution is obtained in <cit.>, and after comparing with the correlation function from the SYK model, these coefficients read
              
              a_n=
              2αtan (v θ)/q,    n=0
              
              2 αtan (v θ) n,    n≥1
              ,   
              b_n=α (v θ) √(2/q),     n=1
              α (vθ) √(n (n-1)) ,     n≥2
              ,   
              v=αβ/π.
              
              It is directly to verify that the following wave function is the solution to the above Schrödinger equation,
              ϕ_n(t)=
              1+2/qlogcos(vθ)/cos(iα t+vθ),    n=0,
              √(2/qn)tanh(α t)/cos(vθ)-isin(vθ)tanh(α t)^n,    n≥1.
              
              
              
              It is suggective to include an angle dependent complexity to the definition of K-complexity for the SYK model, namely, 
              
                  C_K = ∑_n cos(θ/παβ) n |ϕ_n|^2.
              
              With this modification, the K-complexity of the SYK model for nonzero θ is
              
                  C_K = 1/q(θ/παβ)(cosh 2α t-1).
              
              At the conformal limit α→π/β, the K-complexity becomes
              
                  C_K ≈1/q (cosh2π t/β-1),
              
              which agrees exactly with the computation from CV conjecture in JT gravity.
              
              
              
              
              § SUMMARY OF COORDINATE SYSTEMS
              
              We summarize various coordinate systems used in the paper.
              We start with embedding coordinate.
              AdS_2 space can be embedded to
              
              -Y_-1^2 - Y_0^2 + Y_1^2 = -1,      ds^2 = -d Y_-1^2 - dY_0^2 + dY_1^2.
              
              The global coordinate which we use to plot the perturbed AdS_2 spacetime is given by
              
              	Y^-1 = cosν/sinσ,    Y^0 = sinν/sinσ,     Y^1 = σ,      ds^2 = - d ν^2 + d σ^2/sin^2 σ.
              
              The Lorentzian coordinate system is related to the embedded coordinate by
              
              	Y^-1 = z/2[ 1 + 1/z^2 (1-t̃^2)],    Y^0 = t̃/z ,    Y^1 = z/2[ 1- 1/z^2 (1+ t̃^2)],      ds^2 = -dt̃^2 + dz^2/z^2.
              
              Furthermore, a possible Rindler coordinate is
              
                  Y^-1 = cosh r,    Y^0 = sinh r sinhφ,     Y^1 = sinh r coshφ,      ds^2 = d r^2 - sinh^2 ρ dφ^2.
              
              
              § GENERATING FUNCTION FOR SIZE AT GENERIC Θ
              
              
              The generating function of size can be obtained in the large-q limit. A simple generalization of <cit.> gives the generating function of size operator at generic angle θ,
              
                    ⟨O_β(φ)|e^μn̂/δ_β|O_β(φ)⟩/⟨1_β| e^μn̂/δ_β|1_β⟩
              =e^μ2sin ^2γ _μ/-cos(2 γ _μ+ (π/2 - θ ) v _μ) +e^qμcos((π/2 -θ ) v _μ)+(1-e^q μ) cosh(u v _μ)^2/q, 
                    π v_μ = β𝒥sinγ_μ,    
              e^μ qsin(π v_μ/2)= sin(2 γ_μ +π v_μ/2),   δ_β = α/^2/q,  α= cosαβ /2,
                φ=θ+iu. 
              
              At infinite temperature β→ 0, the generating function reduces to
               e^μn̂   =   e^μ/[ 1  + (1 - e^q μ) sinh^2  t ]^2/q,
              
              which agrees exactly with that of K-complexity in (<ref>) if one renormalizes the size by a factor of q since by each step the Liouvillian the size increases a constant amount q. Note that the generating function works prior to the scrambling time since we implement the large-q approximation.
              
              
              
              § ABOUT THE SCRAMBLING TIME
              
              
              In this paper, the scrambling time is defined as the time of the crossover between the exponential growth and the linearly growth of complexities. We will show that this scrambling time also appears in the growth of the size,  i.e. the decay of the OTOC, which should slow down before saturation.
              
              The Schwarizan theory (<ref>) is able to capture the decay of OTOCs at both early time and late time.
              In Ref. <cit.>, the OTOC B_l_2(t̃_1)A_l_1(t̃_2)B_l_2(t_1)A_l_1(t_2) corresponds to the gravitational scattering between the outgoing matter A and the infalling matter B near the horizon of the black hole with initial mass m=π/β. Semi-classically, assuming that the change in the mass of the black hole due to the matter is much smaller than m, and considering the small scaling dimensions l_1,l_2≈0, one find that the time shift of the the outgoing matter A is
              t̃_2-t_2≈1/λ_Lln1+4α C e^λ_L(t_2-t_1-t_R),
              
              where C=ϕ_r/8π G_N=Q/ϵ=α_S N/𝒥, m+α is the mass of the black hole before the matter A goes out, and t_R=(2m)^-1ln(4mC). The exponential time shift slows down when 1∼ 4α C e^λ_L(t_2-t_1-t_R), namely at the scrambling time
              
              t_*∼1/λ_Llnm/α∼1/λ_LlnN/q,
              
              where we find the correspondence of α in the SYK model by matching the energies of PETS on both sides,  i.e.4Cmα=⟨ψ^1_β|H|ψ^1_β⟩-⟨1_β|H|1_β⟩=2𝒥/q. The ln(N/q) dependence in the scrambling time read from the OTOC agrees with the result of the complexities. ] ]
  