Complexity Growth in Integrable and Chaotic Models

We use the SYK family of models with $N$ Majorana fermions to study the complexity of time evolution, formulated as the shortest geodesic length on the unitary group manifold between the identity and the time evolution operator, in free, integrable, and chaotic systems. Initially, the shortest geodesic follows the time evolution trajectory, and hence complexity grows linearly in time. We study how this linear growth is eventually truncated by the appearance and accumulation of conjugate points, which signal the presence of shorter geodesics intersecting the time evolution trajectory. By explicitly locating such"shortcuts"through analytical and numerical methods, we demonstrate that: (a) in the free theory, time evolution encounters conjugate points at a polynomial time; consequently complexity growth truncates at $O(\sqrt{N})$, and we find an explicit operator which"fast-forwards"the free $N$-fermion time evolution with this complexity, (b) in a class of interacting integrable theories, the complexity is upper bounded by $O({\rm poly}(N))$, and (c) in chaotic theories, we argue that conjugate points do not occur until exponential times $O(e^N)$, after which it becomes possible to find infinitesimally nearby geodesics which approximate the time evolution operator. Finally, we explore the notion of eigenstate complexity in free, integrable, and chaotic models.


Introduction
Quantum complexity has been proposed as a quantity relevant for understanding non-perturbative phenomena in quantum gravity, such as the growth of wormholes behind horizons [1][2][3][4][5], the structure of spacetime singularities [6], and the possible appearance of firewalls [7] at late times in an evaporating black hole. The challenge in understanding these conjectures is to have a well-defined measure of complexity in the underlying quantum gravity theory, or, equivalently, in its holographic field theory dual, if the latter exists [8,9]. If the conjectures relating complexity to black hole physics are correct, then we expect that maximally chaotic theories with a holographic dual [10] feature linear growth of complexity for a time exponential in the entropy of the system.
One possibility is that the relevant notion we seek is quantum state complexity. Some progress has been made in computing the circuit complexity of constructing states in some simple free field theories on a lattice [11][12][13][14], but defining state complexity in infinite-dimensional Hilbert spaces that appear in the continuum limit is in general difficult. Free systems have also been used by complexity theorists to build intuition about criteria for complexity growth [15]. An alternative notion that might be relevant is the quantum circuit complexity of the time evolution operator. There has been some progress in computing this quantity in the context of quantum chaotic systems like black holes in holography, and there is evidence that it grows linearly for a long time as expected, given appropriate assumptions [16][17][18]. It is also interesting to consider integrable theories, as some of these do admit quantum gravitational descriptions [19], as well as to discern by comparison what aspects of chaos lead to an exponential time scale for complexity growth. In addition, it may be potentially possible to use complexity as an order parameter in families of theories that interpolate between free, integrable and chaotic limits to distinguish between each regime. The purpose of this paper is to further develop methods for computing the complexity of the time evolution operator in the context of a concrete family of models (the SYK q models) which can be parametrically tuned between free, integrable, and chaotic regimes. On general grounds, the complexity of time evolution is expected to grow linearly with time and then plateau at some fixed value, and subsequently undergo Poincaré recurrences back to small values. The questions of how long this linear growth persists and what height the plateau reaches depend sensitively on the theory under consideration, and will be central issues in this work.
A major drawback of complexity, from a physicist's viewpoint, is its high degree of nonuniqueness. Measuring complexity generally requires many choices, such as a choice of gate set, reference state/operator, or tolerance in preparing the final state/operator. Determining a natural set of these choices for computing complexity in quantum gravity is beyond the scope of this work. Furthermore, computer scientists generally think of complexity in terms of small, discrete operations which are composed to create a complex quantum circuit. As physics generally happens in the continuum, it is advantageous to work with a naturally continuous notion of complexity for operators in physical quantum systems. Such a notion was formulated in terms of minimal geodesic lengths on high-dimensional manifolds of operators [20][21][22], and many recent results on complexity make use of this formalism [11-14, 17, 23-34]. 1 In this setting, there is a relatively natural choice which leads to a unique definition of quantum complexity that is equivalent to the quantum circuit definition: the degree of locality of the Hamiltonian defines a set of "easy" operators (operators which are at most as local as the Hamiltonian). Operators which are more non-local than the Hamiltonian are considered "hard". This choice of splitting into easy and hard operators corresponds to a choice of metric (the "complexity metric") on the group manifold of unitary operators, where directions corresponding to easy operators have low weight and directions corresponding to hard operators have weight of order the Hilbert space dimension. 2 In this geometric formalism for complexity, studying complexity growth is related to studying the growth of the distance function from the identity operator in the complexity metric. As globally length-minimizing geodesics are often difficult to find on generic Riemannian manifolds, the strategy employed by [17] was to look for geodesics that were at least initially globally minimizing, and then to search along those geodesics for possible obstructions to global minimality. On a general Riemannian manifold, such obstructions are either local or global: local obstructions, also known as "conjugate points", imply that the geodesic is not a local minimum of the distance function (i.e., it is a saddle point), while global obstructions, or "geodesic loops", imply that the geodesic is not globally minimal. Any complete picture of complexity growth must include an accounting of both local and global obstructions. Locating global geodesic loops (which are not signaled by conjugate points) in a systematic way is computationally intractable, but (as shown in [17]) conjugate points can be more tractable under certain assumptions, and have a significant effect on complexity growth.
Once global minimality of a given geodesic is obstructed, either by a conjugate point or a geodesic loop, we are guaranteed that the growth rate of the distance function will no longer be exactly linear along this geodesic. However, it may still be approximately linear (with 1 An alternative approach to defining complexity draws intuition from path integrals in quantum field theory, and interprets quantum circuits as optimized procedures for performing such path integrals [35][36][37]. This approach builds on the tensor network formulation of holography [38][39][40][41][42][43][44]. For yet another approach to the analysis of complexity growth, this time making use of unitary k-designs and random circuits, see [45][46][47]. 2 There are proposals for complexity which utilize instead the bi-invariant geometry, which treats easy and hard operators on an equal footing [48][49][50]. A proposal which defines the "infinite cost factor" limit has also been explored [32,51]. a smaller growth rate), if we encounter an isolated conjugate point or geodesic loop, since the new geodesics involved in computing the distance may have growing lengths. We expect, however, that the first conjugate point or geodesic loop along a fixed geodesic associated with time evolution will quickly be followed by the end of complexity growth in general, rather than just a reduction in growth rate, possibly due to a rapid accumulation of subsequent conjugate points/loops. This intuition comes partially from the expected behavior for chaotic Hamiltonians, where after meeting the first obstruction to complexity growth, the complexity is expected to quickly plateau [52]. We will see that free and integrable models also reproduce this expectation, with the first conjugate point signaling the end of complexity growth entirely and a transition to a plateau regime in the distance function within an O(1) time afterward.
Since this paper explores a variety of topics using both analytic and numerical techniques, we now provide a road map by summarizing our main results by section. In Sec. 2, we begin with a review of the geometric formalism developed in [17,22] to keep the discussion self-contained. Since conjugate points play an important role in this work, we explain their significance to complexity growth in detail. We also give new sufficient-but-not-necessary criteria for locating conjugate points in terms of more familiar quantities from thermalization and quantum chaos such as adjoint eigen-operators of the Hamiltonian and infinite-temperature thermal two-point functions.
In Sec. 3, we apply these criteria to the free (q = 2) SYK model. Since free models have relatively simple Hamiltonians, their time evolution operators are simple enough that we can locate all conjugate points and even the geodesic loops which take over after some of these conjugate points. In fact, we find a large number of conjugate points (associated to easy operators) which occur at early (i.e., polynomial) times and signal a rapid end to the linear growth of the complexity of time evolution in the free models, followed by a long plateau. The geodesic loops we study are in one-to-one correspondence with these early conjugate points, which demonstrates that these global obstructions to complexity growth, which are otherwise very difficult to locate, can sometimes be found by leveraging the study of the local obstructions, i.e., conjugate points. These effects place a sharp upper bound on the complexity growth of free systems of N fermions, which is O( √ N ) in the plateau regime.
In Sec. 4 we consider a class of interacting-but-integrable deformations of the free SYK model. We first study a subset of conjugate points in perturbation theory in the coupling constant (which controls the deformation), and find that the deformation causes these conjugate points to move to later times. Going beyond perturbation theory, we also identify certain geodesic loops using the structure of the integrable interaction, which bound the complexity of time evolution in this interacting model. These geodesic loops may not be signaled by conjugate points; if so, this feature of complexity growth distinguishes integrable interacting theories from free theories. The bound on complexity in this interacting integrable model predicts a plateau of height order O(N ) which begins at a time significantly later than in the free model (but still at polynomial time) . Some straightforward generalizations of this simple model show plateaus of height O(poly(N )) for any polynomial in N .
In Sec. 5, we study the possibility of finding conjugate points at sub-exponential times in chaotic theories. In [17], it was argued that in chaotic models, "almost all" of the conjugate points occur at exponential times. One might worry that there are a small number of conjugate points which can nevertheless appear at an earlier time; in particular, prime suspects for this are conjugate points for which the Jacobi field involves only local operators. Indeed, these are precisely the type of conjugate points which obstruct complexity growth in the free SYK model at an early time. Using ideas from random matrix theory and the Eigenstate Thermalization Hypothesis (ETH), we show that in chaotic models, such conjugate points cannot occur before exponential time. This strengthens the arguments of [17] that local obstructions to complexity growth in chaotic models do not occur at sub-exponential times.
In Sec. 6, we numerically study conjugate points for various integrable and chaotic SYK Hamiltonians up to N = 8 (i.e., four qubits). We emphasize that this gives us a concrete (albeit numerical) way to locate obstructions to complexity growth for SYK models, which can in principle be extended to larger N . The numerical results show that a class of conjugate points associated to simple operators (i.e., where the Jacobi field mostly involves simple operators) stay at a fixed time scale as we crank up the weighting of the hard directions, while those associated to hard operators (i.e., where the Jacobi field mostly involves hard operators) rapidly shift to late times proportional to the weighting factor (which is taken to be exponential in N ). Together with the results of Sec. 5, this provides further evidence that the complexity in chaotic models does not plateau until exponential times, modulo global obstructions. Our results on the behavior of conjugate points and geodesic loops in the complexity geometry illustrate how rich geometric structure underlies the growth of the complexity of time evolution in free, integrable, and chaotic theories.
In Sec. 7 we revisit and explore the Eigenstate Complexity Hypothesis (ECH) of [17]. We present the ECH matrices calculated from the eigenstates of SYK models with varying degrees of integrability. The off-diagonal matrix entries for free SYK show fluctuations that scale with the size of the system N , while those of the chaotic models are suppressed uniformly for all N , modulo discrete symmetries of the system. The distribution of the off-diagonal elements is discrete for the free system, which echoes the strong reduction of the number of degrees of freedom analyzed in Sec. 3. This reduction already does not occur in interacting systems even if they are integrable. As expected, the off-diagonal distributions of the interacting-integrable and the chaotic systems have continuous support.
Finally, we end with a discussion of interesting points and future work (Sec. 8).

Conjugate points and complexity growth
In this section, we will discuss conjugate points and their effect on complexity growth. We will see that conjugate points can be studied very concretely in terms of more familiar quantities such as Hamiltonian eigenvectors, thermal two-point functions etc. Of course, global geodesic loops (which are not signaled by conjugate points) should ultimately also play an important role in any complete picture of complexity growth, but a systematic study of these appears to be intractable for now.

Conjugate points in the Euler-Arnold formalism
Let U(H) be the group of all unitary operators on a finite dimensional Hilbert space H, 3 and let {T i } be an orthogonal basis for its Lie algebra with respect to the Killing norm. Quantum circuit complexity is polynomially equivalent to a distance function on U(H), with a certain right-invariant metric G ij (the "complexity metric") which weights tangent space directions corresponding to non-local operators heavily [22]. The choice of which operators are to be considered non-local is not unique; a common choice for spin systems or systems with clear notions of site-based locality is to consider as local all operators which are at most k-local (act on at most k sites or k degrees of freedom) for some fixed k that does not scale with N , the total number of degrees of freedom. Then, any operators which are (k + 1)-local or greater are considered nonlocal and are weighted in the complexity metric. The weighting of the "hard" directions is O(e S = dim H) to ensure the polynomial equivalence to circuit complexity. To implement this weighting, we choose a metric that splits the tangent space into easy directions {T α } and hard directions {Tα}, and weights the hard directions in the length functional by a "cost factor" (1 + µ): When the cost factor is µ = 0, all operators are equally weighted. In Nielsen's setup, the cost factor is taken to be µ ∼ e αS for some O(1) coefficient α, but we will let µ be arbitrary throughout.
Quantum circuits in this context are paths on the unitary manifold, and the complexity of a unitary U is measured by the length of a minimal geodesic connecting the identity to U . An efficient formulation of the geodesic equation on Lie groups equipped with right-invariant metrics was given by Arnold and is known as the Euler-Arnold equation [53,54] where G ij is the metric on the Lie algebra defined in (2.1), and f ij k are the structure constants of the Lie algebra. The Euler-Arnold equation determines a velocity vector V (s), which can then be integrated to give the path followed by the geodesic: where P stands for path ordering. We will always parametrize our paths with s ∈ [0, 1].
Understanding the growth of complexity for a family of operators U (t) now essentially reduces to the question of when a minimal geodesic becomes non-minimizing, and subsequently finding the new minimal geodesic. While the latter problem is difficult, there is actually a local (in the space of paths) signature that a geodesic is non-minimizing: conjugate points. 5 Conjugate points, which were the main objects of study in [17], represent deformations of a geodesic which leave the length and the endpoint locations fixed to first order in the deformation parameter. More precisely: Deformations which leave the length (but possibly not the endpoints) fixed to first order in the parameter η above can be represented as vector fields dU dη along the geodesic, and are called Jacobi fields. If we imagine deforming the geodesic along a Jacobi field, we are not guaranteed that the endpoint of the geodesic will remain fixed at leading order in the deformation parameter. If we do find such a Jacobi field with fixed endpoints along some segment of a geodesic, a shorter path between the initial point and a later point along the path can be found by deforming the geodesic along the Jacobi field between the points which are conjugate, and subsequently smoothing out the resulting kink where the deformed and original paths meet (this relies on the endpoint deviation vanishing). This smoothing reduces 4 The Euler-Arnold equation, while not used in the original formulation of geodesic complexity [22], has been previously used in the context of geodesic complexity and holography [17,24,[32][33][34]. 5 Encountering a conjugate point is sufficient, but not necessary, for a geodesic to become non-minimizing.
1 e itH t ⇤ Figure 1: A cartoon of what happens when a geodesic between 1 and e −itH encounters a conjugate point. The green geodesic is initially the locally minimizing geodesic, before it reaches t * where it encounters a conjugate point (the blue point). For t > t * , the green geodesic is no longer locally minimizing, and a different geodesic (shown in red) will be the local minimum. Note that even though the conjugate point indicates this transition, the new geodesic which takes over after t * is not infinitesimally close to the original one (although we can reach it by gradient flow from the original geodesic).
the length at a lower order in the deformation parameter than the deformation's leading order effect on the length. Thus, the question of whether the endpoint of a geodesic segment is conjugate to the initial point is equivalent to whether there exists a Jacobi field along the segment that fixes the endpoints at leading order in the deformation parameter. Importantly for us, the conjugate point is a signature that the original geodesic is no longer locally minimizing, and is in fact a saddle point after that time. It is worth emphasizing also that the new minimal geodesic which takes over may not be infinitesimally near the original one, and can be highly non-trivial. See Fig. 1 for a depiction of a conjugate point on a compact manifold.
The unitary operator studied in both this work and [17] is the time evolution operator e −itH , where H is the system Hamiltonian. At small enough times t, the globally minimizing geodesic between the identity and e −itH solving (2.2) is the "linear geodesic", a specific geodesic with constant velocity V (s) = Ht. Since the linear geodesic is constant in s, the path ordering in (2.3) is trivial, and the path of unitaries is U (s) = e −istH . By perturbing the Euler-Arnold equation with V → V + δV and keeping the O(δV ) terms, we obtain the Jacobi equation; plugging in V (s) = Ht for the original background geodesic around which we are perturbing, we obtain the Jacobi equation specialized to the linear geodesic: where the L, N L subscripts represent projections to the easy and hard subspaces of generators in su(2 N/2 ). As this is a first-order ordinary differential equation, any initial condition δV (0) can be integrated to a solution δV (s). To find conjugate points, [17] defined a super-operator Y µ (where the subscript µ denotes the cost factor) which takes as input a tangent vector at the identity δV (0), produces the corresponding solution of the Jacobi equation δV (s), and then computes the first order deviation of the endpoint e −itH under deformation of the linear geodesic by δV (s): By expanding the path ordered exponential in a Dyson series, the super-operator effectively computes Solving for δV (s) in terms of the initial velocity deformation δV (0) using equations (2.4) and (2.5), we obtain where the L, N L subscripts denote projections to the purely local and purely nonlocal operator subspaces, and {Tα} is a new orthogonal basis of generators for the nonlocal subspace which diagonalizes the super-operator [H, · ] N L with eigenvalues λα. The intuition for this formula, derived in detail in [17], is essentially to sum up the total deviation along the geodesic by translating the Jacobi field back to the identity and integrating. Functionally, it is the first order correction term in a Dyson series expansion of the path ordering (2.3) in the Jacobi field δV , as written in (2.7). The cost factor µ should be taken to be O(e S ) in the complexity geometry. A conjugate point appears when the first order deviation in the endpoint vanishes for some initial tangent vector δV (0). Therefore, time evolution encounters a conjugate point at time t if the super-operator Y µ has a zero mode at time t. In particular, the zero modes must be Hermitian so that they are valid elements of su(2 N/2 ).

General criteria for locating conjugate points
In this section, we give general criteria for locating conjugate points. Our conditions are sufficient for the existence of conjugate points, but not necessary. Their utility lies in the fact that they relate the locations of conjugate points to more familiar properties of quantum systems such as Hamiltonian eigenstates, adjoint eigen-operators, infinite-temperature thermal two-point functions etc. Further, the hypotheses for these criteria are crucially independent of the cost factor µ and the precise form of H (so long as it is at most k-local).
where µ is the cost factor.
Proof : The proof proceeds by evaluating the super-operator Y µ on the given adjoint eigenoperators of H: (i) Notice first that evaluation of Y µ on a purely local operator O involves only the first term in the square brackets in (2.8). The second and third terms do not contribute since they depend only on the nonlocal components δṼα(0), which are all zero for local δV (0) = O (by assumption).
By evaluating matrix elements of the output in the energy eigenbasis or by expanding out the exponential in the first term of (2.8), we conclude that if O is a k-local, adjoint eigen-operator of the Hamiltonian, then O is also an eigen-operator of the super-operator Y µ : The eigenvalue φ(λt) becomes zero at the locations t * = 2π λ Z, and so we have conjugate points at these locations. Of course, to have a conjugate point we must have a Hermitian zero mode of Y µ , and indeed we do after observing that under these conditions we also have (ii) Likewise, the evaluation of Y µ on a purely nonlocal operator O involves only the third term inside the square brackets in (2.8). The first term inside the square brackets does not contribute because it involves a local projection which will vanish for a purely nonlocal δV (0) = O . To see why the second term does not contribute, observe that it involves the commutator [H,Tα] followed by a projection to the local subspace. Since we have assumed [H, O ] = λ O for a purely nonlocal O , we may take a singleTα to lie along the O direction, and set the rest of δṼα(0) to zero. Then, every term of the form δṼα(0)[H,Tα] L vanishes; all but one vanish due to δṼα(0) = 0, and the final term with Tα ∝ O vanishes due to the projection after the commutator. Again evaluating matrix elements in the energy basis or expanding out the exponential in the third term in equation (2.8), we find that if an adjoint eigen-operator O exists such that O lies entirely along the hard directions, then O is also an eigen-operator of the super-operator Y µ : In this case, the eigenvalue becomes zero at the locations t * = 2π(1+µ) λ Z. Again, we have in mind that the zero modes which lead to conjugate points at these times are really the Hermitian combinations of O and O † , where we have a minus sign in the argument of φ for O † .
We will encounter examples of such conjugate points when we discuss the free SYK model in the next section. In fact, all conjugate points at q = 2 belong to either type (i) or (ii) in Claim 1. As another non-trivial example, consider the q = 4 SYK model. Let the gate set be chosen such that 2-local and 4-local operators are treated as easy, while all other operators are treated as hard. 6 Since the Hamiltonian has a fermion-number symmetry, we can label eigenstates with the corresponding ±1 eigenvalue. Any adjoint eigen-operator of H of the form |m n| where |m and |n have opposite fermion number will therefore entirely lie along the hard directions, and will thus give conjugate points at exactly t * = 2π(1+µ) Em−En Z.
Note that in case (ii), the conjugate points appear at late times, provided the cost factor is taken to be large. In the geometric setup, this cost factor is often taken to be exponential in S, and so we see that these late-time conjugate points appear as an obstruction to complexity growth at exponential times, which is the expected time-scale for complexity saturation in 6 Note that this is a different notion of locality than the notion we use in the majority of this work, where instead we pick some constant cutoff k for which all operators that are at most k-local are considered easy. chaotic quantum systems. So, chaotic theories may have conjugate points of the sort predicted by the hypothesis of Claim 1.(ii), as indeed exemplified by the above example of the q = 4 SYK model with the gate set protected by fermion number symmetry. On the other hand, in (i), the location of the conjugate points does not depend on µ; in this case, conjugate points could potentially lead to a short-time obstruction to complexity growth, where by "short-time" we mean a time of order poly(S). Indeed this is precisely what happens in the free SYK model (see Sec. 3). Since chaotic systems (or, more precisely, systems with geometric, holographic duals) are expected to have complexity growth for exponential time, then we expect such conjugate points which are "associated to simple operators" do not occur in chaotic systems before exponential times. In order to probe this further, we re-formulate the existence of such conjugate points as follows: (2.14) where T α and T β are simple (i.e., at most k-local) generators. If M αβ (t) has a zero mode at time t * , then time evolution encounters a conjugate point at t * .
Proof : Let X α be the zero mode of M αβ (t) at time t * . Now consider We evaluate Y µ (δV (0)), and compute the Frobenius norm of the resulting operator: where in the second equality we have used the fact that the chosen δV (0) lies entirely along the easy directions. Since at time t * we have β M αβ (t * )X β = 0, then we conclude that at time t * we must have ||Y µ (δV (0))|| F = 0, (2.17) which consequently implies Y µ (δV (0)) = 0. Thus, we have a conjugate point at t * .
We will henceforth refer to such conjugate points (which correspond to zero modes of M αβ ) as simple or local conjugate points. Note that M αβ (t) is the infinite temperature, thermal two-point function between two time-averaged simple operators. Claim 2 above states that the first time t * at which this matrix develops a zero mode is precisely when the time evolution geodesic e −itH encounters a conjugate point, and thus necessarily stops being a locally minimal geodesic. Conceptually, this relates complexity growth with a more familiar quantity, namely the thermal two-point function. (In Appendix A we write a general expression relating the full super-operator to the infinite-temperature thermal two-point function which may be of interest for future work.) On the practical side, note that M is a much smaller matrix (polynomial in size) as compared to Y µ (which is exponential in size), and thus gives a useful sufficient-but-not-necessary criterion for locating conjugate points. Such conjugate points, should they exist, will be at a time t * which is independent of µ.
We can also give a physical interpretation to the smallest eigenvalue of M αβ (t). Let λ min (t) be the smallest eigenvalue of M αβ (t). From equation (2.16), we have where we minimize with respect to all (non-zero) local operators δV (0). Physically, this means that it is possible to find an infinitesimally nearby curve with a local initial velocity V (0) = Ht + δV (0) (for infinitesimal ) which satisfies the geodesic equation up to O( 2 ), such that the end point displacement from the target unitary e −itH satisfies: where the subscript F stands for Frobenius norm. Thus, λ min is a measure of the error up to which we can approximate time evolution by an infinitesimally nearby geodesic. If λ min is exactly zero for some t * , then we have a conjugate point at that location. We will call λ min the impact parameter since it measures how close a trajectory with local initial velocity V (0) = Ht + δV (0) comes to hitting the exact final unitary e −iHt . We will return to this in Sec. 5, where we will argue that in chaotic models, λ min ∼ O(e S ) for t < e S , but becomes small thereafter. Consequently, local conjugate points, should they exist, cannot appear before exponential time in chaotic theories.

Relevance of conjugate points in complexity growth
In AdS/CFT, several conjectures relate the quantum complexity of the CFT time evolution operator to the growth of a bulk quantity like an extremal volume or action. While there has been progress in understanding the details of such bulk volume or action calculations, a field-theoretic formulation of circuit complexity in infinite-dimensional Hilbert spaces which reduces to the standard notion of quantum complexity in finite dimensions is still incomplete. This has led to the development of toy models for the complexity geometry which are designed to reproduce certain coarse-grained features of distances on the full unitary manifold with a right-invariant complexity metric [55,56]. For example, one such toy model involves a particle moving on a high-genus Riemann surface with metric induced from its universal covering space, the hyperbolic disk H 2 [55].
While such toy models have led to interesting insights into the behavior of holographic complexity, they lack a crucial feature of the finite-dimensional complexity geometry: conjugate points. In the example of the particle moving on a Riemann surface, there are no conjugate points because the sectional curvature of the induced metric is strictly negative.
In an attempt to justify this shortcoming, one might appeal to results of Milnor on sectional curvatures of Lie groups [57], which roughly imply that most sectional curvatures on "complicated enough" Lie groups with right-invariant metrics are negative. Crucially, however, the results of [57] do not imply that all sectional curvatures are negative. In fact, the most important result in [57] for our purposes is the fact that any right-invariant metric on SU (n) for n > 2 is required to either have some strictly positive sectional curvature or else be completely flat. Some of these curvatures were recently computed explicitly for complexity metrics in [58] and were found to be positive.
There is an obvious tension between the lack of conjugate points in the toy models and the fact that in the finite-dimensional complexity geometry (a right-invariant metric on the unitary group), conjugate points are guaranteed to exist and obstruct the complexity growth of time evolution with arbitrary Hamiltonians. 7 This fact was emphasized in the original formulation of complexity geometry [22], and also in its adaptation to the Euler-Arnold formalism [17].
A possible perspective on this tension is to imagine that, in the context of finite-dimensional holographic systems like the SYK model, conjugate points may move off "to infinity" or simply disappear from the relevant minimal geodesic as the cost factor µ is increased, leading to a situation where there are never any conjugate points along the geodesic relevant to complexity. Unfortunately, as was briefly discussed in [17], this is impossible due to two facts: 1) the initial linear growth of time evolution's complexity is captured by the linear geodesic, and 2) the right-invariant complexity metric depends continuously on the cost factor µ. Using these two facts, we will explain in more detail an argument sketched in [17] which demonstrates that conjugate points must exist along the linear geodesic for arbitrary local Hamiltonians at finite distance and cost factor.
We begin by noticing that the case of zero cost factor, µ = 0, corresponds to a bi-invariant 7 See [59] for a simpler Lie group geometry where conjugate points are guaranteed to be the first obstruction to complexity growth. metric on the Lie group. In this case, the exponential maps of the Lie group and Riemannian manifold coincide, which means that all geodesics take the form e −isH for some Hamiltonian H. In the bi-invariant metric, conjugate points are known to exist at finite distance [22]. 8 Since they begin at finite distance, they cannot move "to infinity" since they are zero modes of the super-operator Y µ , and these zero modes depend continuously on µ. If they were to move to infinity at some finite value of µ, there would be a discontinuity in the super-operator before and after this value.
The only other possibility is that the conjugate points could "disappear", which would correspond to a zero mode of the super-operator becoming complex. That is to say, the Jacobi field which gives the conjugate point could pick up a non-Hermitian contribution at some finite value of µ, and in order to have a true conjugate point the Jacobi field must be purely Hermitian. We do not have a guarantee from simple continuity that this cannot happen, since, for example, the same thing happens for the polynomial equation There is no discontinuity in µ on the left hand side but the solutions become complex as µ goes from negative to positive. So too could the Jacobi fields generating the conjugate points become non-Hermitian at some finite value of µ. However, it turns out that this also cannot happen.
To understand why conjugate points cannot disappear, we apply Morse theory on the space of paths. Let Ω(U 1 , U 2 ) be the space of paths on the Lie group between unitary operators U 1 and U 2 . The dimensionality of this space is formally infinite, but this subtlety turns out not to affect any conclusions [60][61][62]. 9 For the complexity of time evolution, the relevant path spaces are Ω t,H ≡ Ω(1, e −itH ). (2.20) That is to say, Ω t,H is the space of all smooth paths γ(s) with γ(0) = 1 and γ(1) = e −itH . For convenience, we parametrize all paths with s ∈ [0, 1]. We can consider a real-valued function on Ω t,H which is often called the energy functional where we have made use of the splitting of the Lie algebra into local and nonlocal directions (labeled by α andα, respectively), the right-invariance of the complexity metric, and also the velocity along the path V (s) ≡ dγ/ds. 8 In particular, they appear at t * = 2π (Em−En) Z for all eigenvalues Em, En of H. 9 The original work of Morse, reviewed by Milnor in section III of [60], relies on finite-dimensional approximations of the full path space, to which Morse's theory is then applied. By contrast, [61,62] prove the same results by working directly in the infinite-dimensional setting.
Critical points of the energy functional E (µ) on Ω t,H are precisely the paths with velocity V (s) which are geodesics between the identity and e −itH . The most important of these for us is the linear geodesic, which is simply the path V (s) = Ht. Since the linear geodesic is independent of µ, the point in Ω t,H to which it corresponds is fixed as µ increases. Call this point L ∈ Ω t,H . The tangent space to L, and more generally to any point γ in the path space, is the space of vector fields δV (s) along γ for which δV (0) = δV (1) = 0. 10 With this notion of tangent space, one can define the Hessian of the energy functional E (µ) evaluated at L, which we will denote E (where the derivatives are taken in the space of paths), keeping all dependence on µ, t, and H implicit.
One can now apply the Morse index theorem on Ω t,H using E (µ) as the Morse function. The Morse index theorem applied to our situation states that the number of negative eigenvalues of E is equal to the number of conjugate points (counted with multiplicity) along the geodesic L, and that E only has a zero eigenvalue if the endpoint e −itH is conjugate to the identity along L [60]. Since E (µ) depends continuously on µ, and L is independent of µ, the eigenvalues of E must also depend continuously on µ. Therefore, the only way we can "lose" a conjugate point along L is for an eigenvalue of E to pass continuously through zero. In other words, the conjugate point must move beyond e −itH along the linear geodesic. This means that conjugate points cannot simply disappear; the only way to get rid of them is to boost the cost factor µ high enough to push them past the endpoint of the geodesic L. So, by taking t large enough (but still finite), we can extend the endpoint of L to always find conjugate points along L at finite distance and cost factor, just as we claimed. This also amounts to a non-perturbative proof that zero modes of Y µ are always Hermitian matrices because if a zero mode were non-Hermitian then the corresponding conjugate point would disappear, but the zero modes are in one-to-one correspondence with the conjugate points.
All of this means that conjugate points are relevant for any complexity calculation which employs complexity geometry and involves the linear geodesic L, and toy models which ignore them are useful but incomplete representations of the total complexity geometry. It would be interesting to find a toy model which can include conjugate points.

Free theories
We now study the growth of complexity in free and integrable models, starting with the quadratic free fermion model, with Hamiltonian where J ij is an anti-symmetric matrix and the sums run from 1 to N . We consider this model as a q = 2 instance of the SYK q family of models [10,63], 11 There, J i 1 ...iq is totally antisymmetric and is drawn from a Gaussian distribution with mean zero and variance parameterized by J , In our context, we consider a particular instance of the model where we have sampled the couplings J ij from such a distribution. The matrix J ij is antisymmetric and therefore can be written as where V is an orthogonal matrix, and D is block-diagonal with antisymmetric blocks: The matrix V is constructed as follows. First, write the usual diagonalization J = U ΣU † . Since J is antisymmetric, the matrix U is unitary and the eigenvalues of J are ±iω p /2, p = 1 . . . N/2. Next, define the unitary matrix out that U Ω is always a real matrix, so we can identify V = U Ω and then J = V DV T . Now we can define new fermion operators which also satisfy the same anti-commutation relations The notion of locality is unchanged by this transformation, since the new fermion operators are linear in the old ones and V is orthogonal. In terms of these new operators, the Hamiltonian becomes Finally, we define the ladder operators with all other anti-commutators vanishing. In terms of these, the Hamiltonian becomes In this Dirac fermion language, there is a new useful basis of the 2 N − 1 operators which span the algebra su(2 N/2 ). To define this basis, we begin by writing a vector of 4 operators With the entries of this vector labeled by indices in the order β p ∈ {0, −, +, 3}, the operator basis is then the set of products over all choices of {β p }, where we discard the identity β 1 = . . . = β N/2 = 0. The Hamiltonian can be written compactly as and J (p) 3 has eigenvalues ±1 in the energy eigenbasis. Thus, the 2 N/2 eigenvalues of H are for every possible choice of the coefficients σ p from {±1}. The natural notion of locality in the Dirac basis, derived by considering an operator with k Majorana operators to be k-local, is to consider J + and J − as 1-local operators but J 3 as a 2-local operator and J 0 as a 0-local operator. Then, the locality of a general product of J βp 's is simply the sum of the individual localities. Since the Hamiltonian is 2-local, then we will take k = 2 in the rest of this section. Free fermion time evolution was also studied in [15]; we will see that geodesic complexity techniques both reproduce the results found there and allow us to uncover new features of free theories.

Conjugate points
We are interested in the complexity of the unitary operator (3.17) First, we study conjugate points for the linear geodesic. Let us look at the super-operator Y µ derived in [17], whose zero modes as a function of t correspond to conjugate point locations.
For free theories, it turns out that every conjugate point corresponds to a local or non-local eigen-operator of ad H .
To understand the free theory, we observe that the adjoint action of the Hamiltonian is already diagonal in the Dirac fermion basis (3.14) and, recalling that β p ∈ {0, +, −, 3}, we can write it as For the 2 N/2 − 1 operators in the basis that involve only J 0 or J 3 , the adjoint eigenvalue is zero. We take 3-local and higher operators to be nonlocal, since the Hamiltonian is quadratic in the Majorana fermions, i.e., k = 2. Since the adjoint eigen-operators of the Hamiltonian split nicely into simple and hard operators, we can obtain all the conjugate points using Claim 1 in Sec. 2.2. The locations of conjugate points associated to local operators are given by Claim 1.(i). They are where p 1 = p 2 , ω p 1 > ω p 2 > 0, and ω p > 0. We may always define all ω p > 0 for the price of introducing a minus sign in the definition of σ 3 p , and we order the ω p so that ω p > ω q for p < q. These families of conjugate points are associated with operators of the forms respectively. These are two-fold degenerate conjugate points; there are corresponding partner operators, such as i(A † p − A p ) for the second operator in (3.20). Similarly, the locations of conjugate points corresponding to the purely nonlocal operators are given by Claim 1.(ii), where we cannot pick the same ω p twice, and all possible combinations of plus and minus signs can occur in the denominators subject to the constraint that the overall result should be positive. The associated operators are respectively

Exact geodesics
As we showed above, the conjugate points associated to nonlocal directions are quite far from the identity due to the cost factor, while those associated with local directions occur at a time of O(poly(S)) since numerical experiments reveal the range of non-zero ω p to be between O(1/N ) and O(1), with a typical spacing of 1/N . 12 Therefore, as one might expect in the free theory, obstructions to complexity growth occur nearly immediately. We would like to go beyond just identifying the location of such obstructions and actually find the new globally length-minimizing geodesics which replace the linear geodesic in the complexity calculation.
A general geometric strategy for finding these new geodesics will be to isolate relevant subalgebras of su(2 N/2 ) where the effect of the conjugate point can be completely understood. While technically there are conjugate points associated with 2-local operators which occur sooner, it is illustrative to begin with the family of points corresponding to a single ladder operator A † p . Again, there is a two-fold degeneracy of these conjugate points which arises due to the A p . To understand the behavior of the conjugate point at t * = π/ω p , we must at least study the algebra generated by A p and A † p . Furthermore, whatever our choice of subalgebra, we must also include the relevant terms in the Hamiltonian, namely the projection of H to our subalgebra. The smallest possible subalgebra that fits our needs is just a copy of su(2), generated by [J While J + and J − are not Hermitian, they are traceless, and we may Hermiticize them by taking linear combinations to obtain valid su (2) generators. This su (2) is essentially a copy of the Pauli algebra, and exponentiates to an SU (2) subgroup within our total manifold SU (2 N/2 ). Since SU (2) is simply S 3 , we have an immediate interpretation of the conjugate points at t * = π/ω p . The path traced by e −itH in SU (2 N/2 ) begins at the north pole of this S 3 , and the conjugate point sits at the south pole.
There is a corresponding algebraic avenue for understanding this result. Because J (p) 3 has eigenvalues ±1 in the energy eigenbasis, and we have [J ] = 0, the time evolution operator splits as The conjugate point occurs at t * = π/ω p because it is precisely at this time that we have the equivalence Thus, we are led to conclude that the new geodesics should be written by simply modifying the coefficients in the Hamiltonian in the appropriate way to trace out the other half of the great circle and to return to the north pole on S 3 at time 2t * .
We must simply reverse the velocity in the J (p) 3 direction and decrease the coefficient appropriately so that it returns to zero as we return to the north pole of the relevant SU (2). The appropriate operation which achieves this, and takes into account the necessary changes when encountering all other conjugate points in the family t * = πZ/ω p , can be neatly written as where we always take the principal branch of the logarithm with a cut on (−∞, 0). This function effectively computes ω p t modulo 2π, with a result in the range (−π, π]. The new velocity, including the effects of all conjugate points associated to 1-local operators, is This velocity is still constant (i.e., s-independent) and purely local, so it is still a geodesic. This family of geodesics was found in [15] as fast-forwarding Hamiltonians for free fermion time evolution. However, as we will now see, there is more structure in the free theory which allows for a constant factor improvement over the above construction.
Before dealing with the conjugate points corresponding to 2-local operators, we note that the family of geodesics just described induces a self-averaging behavior for the complexity in precisely the way which was understood in a toy model developed in [17]. The toy model made use of the single-qubit complexity geometry, which is simply SU (2). An N -member ensemble of single-qubit Hamiltonians was defined, and the ensemble-averaged complexity in that situation behaved in precisely the manner we have just described for the free SYK model with N Majorana fermions. This is a quantitative instance of self-averaging, an effect which is generally difficult to understand analytically.
With that being said, there are effects at times of O(poly(S)) in the N -Majorana theory which are "intrinsically quantum", and do not arise from self-averaging. These are the conjugate points like t * = πZ/(ω p 1 + ω p 2 ) associated with 2-local operators, which have no analog in the ensemble average toy model. We will see explicitly why this is the case by again analyzing these conjugate points from both geometric and algebraic viewpoints. On the geometric side, we search for a subalgebra which includes the operators generating the (two-fold degenerate) conjugate point, any other operators necessary for the subalgebra to close, and the projection of the Hamiltonian to this subspace. It turns out that we can again manage with just a single su(2) subalgebra generated by .

(3.30)
It may seem like the third operator is not actually a projection of the Hamiltonian, since H involves a weighted sum of J (p) 3 's with different coefficients. The point here is that we are projecting H along a particular direction which involves operators from both the p 1 and the p 2 subalgebras discussed in the 1-local case. In other words, we rewrite and project along the J direction. This subalgebra exponentiates to a copy of SU (2), and we again have an interpretation of the conjugate point as the south pole of an S 3 .
The algebraic viewpoint is a bit more instructive in this case as opposed to the 1-local situation. In that case, we observed that iJ had eigenvalues ±i, and so its matrix exponential was 2π-periodic, which led to a conjugate point at π/ω p where ω p was the coefficient in H. However, in general, for sums of different J (p) 3 , the most we can say is that the eigenvalues are integers. Luckily, for the sum of precisely two J (p) 3 , the eigenvalues are ±2 or zero. Therefore, the matrix exponential of i(J ) is actually π-periodic. This explains why the conjugate point sits at t * = π/(ω p 1 + ω p 2 ) as opposed to 2π/(ω p 1 + ω p 2 ).
With an understanding of these conjugate points, we can write the new velocity. Again we simply make the replacement which handles all conjugate points in the family t * = πZ/(ω p 1 + ω p 2 ). When we encounter the conjugate point at t * = π/(ω p 1 − ω p 2 ), the same replacement will occur on the second term on the right hand side above. Unfortunately, as much as we would like to write a single expression which incorporates the changes in the geodesic after all conjugate points associated to 1-and 2-local operators, we cannot accomplish this with our logarithm branch cut trick. We will simply provide a description of the total velocity.
The linear geodesic begins with velocity Ht. As we increase t, the endpoint of the geodesic moves, and we may encounter a conjugate point. To keep track of these changes, we keep a table of coefficients c p (t), and we will periodically update these with t so that the velocity of the globally minimal geodesic is always (before times of order µ) (3.33) Initially, at very small times t, we have c p (t) = ω p t, and these coefficients will always locally increase linearly with ω p t. To know when we should update a particular c p (t), we keep track of three types of quantities: c p 1 (t) + c p 2 (t), c p 1 (t) − c p 2 (t), and c p (t) themselves. We will update c p (t) so that all such quantities are in the range (−π, π]. Whenever one of the first or second type increases beyond π, we rewrite the velocity as in (3.32) and update c p 1 (t) and c p 2 (t). In the case of the first type, their sum is updated to be in the range (−π, π] but their difference is unchanged. For the second type, their difference is updated but their sum is unchanged. Similarly, when one of the third kind increases beyond a multiple of π, we simply update the individual c p (t) to be in the range (−π, π]. Notice that in the first two cases we needed a second linear relationship (keeping one of the sum or difference fixed) in order to update both c p 1 (t) and c p 2 (t).
It is this "quantum" effect which separates the exact N -Majorana free theory from the ensemble average of single-qubit theories. Indeed, the quantum interference effects of the conjugate points related to 2-local operators actually prevent us from reaching the conjugate The ω p /2 which control the growth of the coefficients c p (t) in (3.33) are the positive eigenvalues of the antisymmetric coupling matrix J ij whose entries are independent Gaussian random variables with mean zero and variance points associated with 1-local operators: if we ever had c i (t) = π for some c i (t), we would certainly have some sum or difference of c i (t) equal to π already, unless all the other c j =i (t) are zero, which is quite finely tuned. The fact that we are able to understand all the globally minimizing geodesics before exponential times as a function of t by studying only conjugate points on the linear geodesic is due to the geometric description of all local operator conjugate points as south poles of 3-spheres. A plot of complexity for the free SYK model is shown in Fig. 2. As all O(N ) terms in the diagonalized Hamiltonian are upper-bounded by π due to the local conjugate points and corresponding geodesic loops, there is a hard upper bound on the free complexity of O( √ N ). The conjugate points associated with non-local operators are not relevant for this discussion because they occur at far later times of O(µ ∼ e αS ). Thus, we have essentially determined the full structure of geometric complexity in the free theory at sub-exponential times, up to the existence of geodesic loops which are not signaled by conjugate points.
We can make progress on this front by ruling out at least one simple class of potential geodesic loops which are not signaled by conjugate points. Though we have demonstrated that conjugate points corresponding to nonlocal (3-local and higher) operators occur at times of order µ, and are thus not relevant for complexity growth below such times, we may wonder if a similar algebraic effect as (3.32) can occur for e.g. a sum of three J (p) 3 's even without a conjugate point. It is clear that there is an algebraic relationship which would allow such a replacement: the sum of three or more J (p) 3 's is still integer valued, so the matrix exponential will be at most 2π-periodic. This would be a geodesic loop that occurs without a conjugate point in the free theory. However, this cannot occur, because of the way the coefficients scale. In general, a sum of m J Notice that for m = 1 and m = 2, the right hand side is π, and this is led to our update rules for the c i (t). However, for m = 3 it is 4π, which means the average value of the c i (t) is 4π/3, which is greater than π. We cannot reach this regime, because the c i (t) are all valued in (−π, π] due to effects of the 1-and 2-local operators. For m = 4, the average value is π, but this also cannot occur because (since c i (t) ≤ π) we must have c i (t) = π for all i for the average of them to be π. This violates the conditions placed by the 2-local operators, namely that the sum of any two c i (t) is less than or equal to π. A similar story holds for all m > 4. So, no periodicity effects arise for this number of J 3 's, and indeed there are no conjugate points associated with such effects.
Throughout this discussion, we have assumed that k = 2, or in other words that 3-local and greater operators are considered nonlocal from the perspective of the complexity metric. However, the classification and locations of conjugate points at arbitrary µ that we described in Claims 1.(i) and 1.(ii) in Sec. 2.2, and then applied to the free theory, does not actually depend on this assumption. The reason our analysis cannot be extended to k > 2 is more subtle. Let us consider k = 3 for concreteness. By Claim 1.(i), there is a conjugate point family at πZ/( The next step to understand these conjugate points is to analyze this operator and the Hamiltonian projection from the geometric or algebraic perspective. From the geometric perspective, the situation is significantly more complicated than the 1-and 2-local cases because the relevant subalgebra is no longer su(2). The two Hermitian operators associated to the conjugate point and the Hamiltonian projection do not close under the Lie bracket, and more operators must be added to ensure closure. Moreover, beyond S 3 , none of the higher-dimensional spheres are Lie groups, so the geometric interpretation of the conjugate point will no longer simply be arrival at the south pole of a sphere. The algebraic perspective has an analogous difficulty: the sum of three or more J (p) 3 's can certainly have an eigenvalue of ±1 or ±2, which is less than the multiplicity we would need to explain the appearance of the conjugate point so soon by some periodicity condition on the matrix exponential.
In a certain sense, this result is not surprising. The 3-local and higher operators do not have such simple interpretations because physically they represent "shortcuts through chaos" which generate free time evolution faster than the free system itself. That is, after the linear geodesic (corresponding to time evolution with respect to the free Hamiltonian) is replaced by a new globally minimizing geodesic at a non-local conjugate point, the shorter trajectory along the new global minimizer can be thought of as Hamiltonian evolution with respect to a different, chaotic effective Hamiltonian. These shortcuts would be interesting to understand, as they utilize chaos in a structured way. 13 In other words: "Chaos isn't a pit. Chaos is a ladder." Of course, it could be that the chaotic deformation "wraps around" a submanifold in the same way as the conjugate points we were able to understand above, and leaves us with a globally minimizing velocity that does not actually involve 3-local or more terms. This observation does not change our conclusion that there are special chaotic deformations which allow speedups for free time evolution; it only means these speedups are not optimal.

Summary
We found all conjugate points along the linear geodesic in the complexity metric, and we determined the associated geodesic loops. To find the conjugate points, we determined all eigenvectors of the super-operator Y µ at arbitrary µ. Using this information, we constructed a geodesic (as a function of t) which is globally length-minimizing from the identity to e −iHt , up to the existence of possible geodesic loops which were not associated with any conjugate points. The length, and therefore the complexity, was bounded at O( √ N ).

Integrable theories and deformations
In Sec. 3, we studied obstructions to complexity growth along the linear geodesic associated with time evolution in the free SYK model. In this section, we will study a class of interactingbut-integrable Hamiltonians. To this end, consider adding a quartic interaction H 1 to the free (quadratic) Hamiltonian H 0 which preserves integrability. An example of such an interaction is a term which is quadratic in the J (p) 3 , so that the total Hamiltonian is 3 ] = 0, we may take M ij = M ji to be a symmetric matrix. In order to avoid introducing a nonzero trace, we take M kk = 0. Since H commutes with all the J (i) 3 's, this interaction term preserves integrability. It is important to note that since our full Hamiltonian is now quartic, we will treat k ≤ 4-local operators as easy in the complexity metric.
The analysis of conjugate points for the above integrable Hamiltonian is somewhat involved, and so we will approach it via perturbation theory in the coupling . However, we note that this analysis becomes much simpler if we modify our gate set slightly by allowing ourselves access to one new elementary operation. To see this, observe that the adjoint eigenvectors of the Hamiltonian H are given by where there are p Dirac excitations in |m and q Dirac excitations in |n , and c p,q is a constant. Note that these operators are almost like "local" operators built out of products of individual fermions, except for the inclusion of the projector P 0 in the product above. In principle, we could consider a gate set where operators of the form (4.2) with (p + q) ≤ k are treated as local/simple, while the rest are treated as hard. We can then again use Claim 1 from Sec. 2.2 to compute the locations of all conjugate points in this case. These will be given by for the simple operators, and for the hard operators, where E m are the eigenvalues of H. At any rate, we will not consider this choice of gate set any further in this work, instead focusing on the more standard choice with k ≤ 4-local operators being treated as simple.

Perturbative conjugate points
To begin with, we will study the effect of the interaction term on the location of conjugate points perturbatively in the coupling constant . 14 Since the general perturbative analysis 14 Readers who do not wish to follow the detailed perturbative calculations may skip ahead to the summary at the end of this section, and proceed to Sec. 4.2.
is very complicated, we focus specifically on the conjugate points of H 0 associated with Jacobi fields that are 1-local operators, such as A † p + A p . Recall that the conjugate points associated with these operators correspond to zero modes of the super-operator Y µ which appear at certain times, and in general they are eigen-operators of Y µ (defined using H 0 ) with eigenvalue Since these zero modes (and all others in the free theory) are two-fold degenerate, we must employ degenerate perturbation theory. In fact, since we have expanded our definition of easy operators to include up to 4-local terms, there are additional 3-local operators which have the same eigenvalues under Y µ . These are e.g.
3 for p = q, and these 3-local operators lead to conjugate points at the same times as the 1-local operators above, so the degeneracy is enhanced.
We proceed by perturbing Jacobi fields and conjugate point times in response to the perturbation of the Hamiltonian (4.1), We reproduce here the equations governing the Jacobi equation and the super-operator, in which we will make the above replacements and expand: Subsequently, we will proceed order by order to see the effect of the perturbation on the locations of conjugate points.

Zeroth order
At O( 0 ), the total Hamiltonian is the free Hamiltonian H 0 , so we can pick δV 0 (0) to be any linear combination of the form: where z i and z j =i are complex numbers. We then obtain a corresponding conjugate point family at t Note that in this case δV 0 (s) = δV 0 (0), and lies entirely along the easy directions. Importantly, we assume µ > 0 here: if µ = 0, then there is a much larger degeneracy due to operators of the form (with j 1 = . . . = j p = i) and the above ansatz needs to be modified.

First order
At O( 1 ), the Jacobi equation reads Since H 0 is quadratic, it does not mix between local and non-local directions. Thus, the solutions are From here, we can compute the perturbative terms in the super-operator, where we have assumed [H 0 , H 1 ] = 0 as is the case for the particular integrable deformation (4.1). In order to extract the change in the location of the conjugate points, we set the O( ) term to zero so that zero modes of the super-operator with respect to the free Hamiltonian remain zero modes of the perturbed Hamiltonian: In order to make further progress, we project this equation into the local and non-local directions: Local : Non-local : Here we have again used the fact that H 0 is quadratic, and so does not mix between local and non-local operators. Plugging in the Hamiltonian deformation H 1 and the ansatz for In taking the local projection, we need to be careful because (J (p) 3 ) 2 = 1. So the local projection becomes Now going back to the local constraint: Local : we take its overlap with A i and A i J 3 respectively. This kills the δV L term above, and we get t Here taking i = 1 suffices to show the general structure. The equations then can be written in the matrix form X z = 0, where (4.26) There are only nontrivial solutions when the determinant of X vanishes. The determinant can be written . So we find N/2 − 2 pairs of conjugate points do not move, while the remaining two pairs move to the new locations: The z's which correspond to the non-trivial displacement at first order are given by The z's corresponding to t In addition, we need to also work out δV (1) (0). From (4.19) it follows that δV but the coefficients w i and w j are not determined at this order. A priori, we could have had other operators appearing in this expansion, but their coefficients must be zero by (4.19). On the other hand, we can solve for δV (1) N L from equation (4.20). Note that which suggests the following ansatz for δV (1) N L (0): where the c jk 's are some complex coefficients to be determined. Substituting this into equation (4.20), we find c jj = 0, while for j = k,

Second order
At second order in , the Jacobi equations are given by d ds δV (2) (1 + µ) d ds δV We will only need to know the explicit form of δV (2) L , which is given by In order to study the displacement of conjugate points at second order, we now compute U −1 δU : (4.38) Then, we must substitute t = t (2) * , and extract the second order terms. In doing this, we should be careful to keep in mind that δV also depends on t.
As in the first order case, the displacement of the conjugate points is determined by taking the overlap of this equation with the local directions, in particular with A i and A i J (for j = i). The terms proportional to δV (2) drop out of these overlaps, and so we do not need to explicitly compute δV (2) at this stage.
In order to simplify the computation, we will only track the conjugate points which do not already move at first order, i.e., which have t (1) * = 0. For these points, we have So, we need to compute In addition, we also need Finally, it is easy to check that the displacement projected along A i vanishes for the conjugate points which do not move at first order, if we take j M ij w j = 0. For these, the displacement projected along A i J (j) 3 is given by: Using φ(x) = e ix/2 sin(x/2) x/2 , one can check that the imaginary parts inside the square brackets precisely cancel. Therefore, these equations take the form: where α is the real constant equal to the term in brackets above. Assuming that the M ij are all different, then we can solve these equations for the z j : Finally, we need to impose the constraint j =i M ij z j = 0 assuming w i = 0, which translates to Therefore, the second order displacements of the N/2 − 2 conjugate points are generically nonzero and can be obtained from the zeros of the complex function f (τ ). These zeroes are always real, as can be checked by explicitly substituting τ = x + iy into (4.47). Note that if any two M ij coincide, then we lose a zero, and that zero corresponds to looking for a solution with w i = 0. We will not consider these additional special cases here.

Integrable geodesic loops
While we have focused on perturbation theory for the conjugate point locations, it is also possible to find certain geodesic loops in this model analytically. Recall that in the free theory, we found many geodesic loops, each of which came from a conjugate point associated with an easy operator. In the deformed theory, we do not have an exact analytic handle on conjugate point locations, but the same sorts of loops can occur because the two terms in (4.1) commute. This means that the time evolution operator splits as The loops we found for the operator e −iH 0 t in Sec. 3 also apply here. Furthermore, since the product J 3 . To take these into account, we follow (3.33) and define coefficients d ij (t) with bounded range where we define the π modulus to take values in (−π/2, π/2]. Then, using the global velocities (3.33) for H 0 , a bounded-length path to e −iHt is The complexity is upper-bounded by the length of this path: An instance of this function is shown in Fig. 3. Qualitatively, we may conclude that the complexity reaches a plateau here as well, but with greater height than the free case. The free complexity is upper-bounded by O( √ N ) since there are N/2 coefficients c p (t) with maximum value π, but the integrable perturbation allows for N (N − 1)/2 more terms in the d ij (t), which leads to an upper bound of O(N ). A strict upper bound in this case is in fact where we simply took the upper limits c p = π and d ij = π/2. We have been careful to say upper-bounded in this discussion because we have not exactly located the conjugate points in this model, and there may be some which are closer to the identity than any of the geodesic loops we considered here. Of course, as in the free case, we also do not have analytic control over every geodesic loop. This and other integrable interacting models could furnish interesting examples of geodesic loops in complexity geometry which are not signaled by a conjugate point in a straightforward way. where we require k ≥ 2c so that H c−1 is an easy operator in the complexity metric, and λ is symmetric in all indices and vanishes when i j = i for any j = (so it is strictly 2c-local). Following the same procedure as before, the complexity of e −i(H 0 + H c−1 )t is upper-bounded by Thus, we have a family of integrable models with complexity of time evolution that is upperbounded by a polynomial O(poly(N )) that depends on the order of the interaction c.

Summary
We calculated the first and second order shifts in location of the conjugate points associated with 1-local operators in the free theory under the integrable deformation (4.1). At first order, all but two of the N/2 degenerate conjugate points remain fixed, and the two which move do so by a distance which depends on the perturbation couplings M ij but not on the cost factor µ. At second order, the N/2 − 2 points which did not move at first order begin to move, and are shifted by a distance which is sensitive to µ. As this shift can become large for µ 1, the perturbation theory may break down. We also found geodesic loops which were analogous to certain loops found in the free theory, but for which we did not find associated conjugate points. These represent potential examples of geodesic loops which are not signaled by conjugate points.
The perturbative results suggest that the complexity grows linearly for a long time as the conjugate points we studied move to later times as µ is increased; however, the existence of these geodesic loops shows otherwise. There are also other conjugate points associated to operators of higher locality which may be independent of µ, the existence of which will be suggested by our numerical results in Sec. 6

Impact parameter and local conjugate points in chaotic theories
We now turn to the interesting case of chaotic Hamiltonians. In [17] it was argued that in a chaotic model, the super-operator Y µ takes a simple form in the energy eigen-operator basis: where under appropriate assumptions the Frobenius norm of the correction term · · · was shown to be exponentially small. Thus, the diagonal entries of the super-operator Y µ in the |m n| basis are O(1) for t (1+µ) Em−En . Since the off-diagonal entries are small, we thus expect that the eigenvalues will also be bounded away from zero, and given that µ scales exponentially with S, we conclude that conjugate points do not occur at sub-exponential times. However, there is a caveat: while the off-diagonal elements of Y µ are suppressed, at the same time there are an exponentially large number of such off-diagonal entries. So although "almost all" of the eigenvalues of Y µ will be O(1) for sub-exponential times, we cannot be certain that a small number of zero modes cannot occur. In fact, local conjugate points (see Claim 2 in Sec. 2.2) are prime suspects at sub-exponential times, as their locations do not depend on the cost factor µ. In Claim 2, we re-formulated such conjugate points in terms of zero modes of the positive semi-definite matrix M αβ , which is the matrix of infinite temperature thermal two-point functions between time-averaged simple operators: We will now argue that in chaotic systems, zero modes of M αβ -and hence local conjugate points -can only potentially arise at exponential times. Our strategy will be to show that the minimum eigenvalue λ min (t) of M αβ is exponentially large for t < e S , and becomes small only thereafter. We will refer to λ min as the impact parameter (see Sec. 2.2).
By expanding in the energy eigenbasis and evaluating the integrals, the matrix M αβ can be written as: where |m , |n are energy eigenstates with energies E m , E n . With the above formula for M αβ (t), we can now estimate the time t * at which we expect a zero mode by using intuition from random matrix theory and the Eigenstate Thermalization Hypothesis (ETH). In this context, we assume ETH is satisfied for the k-local operators that we consider easy in the complexity metric. 15 First, notice that for times less than the inverse maximum energy difference 1/(E max − E min ), we have g(t(E m − E n )) ≈ 1. If we were to make this replacement in M αβ , we would find . We can replace these diagonal terms with 2 N/2 δ αβ by rearranging the above equation as: Now our basic strategy will be to argue that for t < e S , (i) the diagonal entries of M αβ are O(e S ), while (ii) the off-diagonal entries of M αβ are O(1). Since the matrix is polynomial in size (as the α, β indices run over simple operators), this then implies that the eigenvalues will all be O(e S ). On the other hand, when t e S , the diagonal entries can become O(1), and thus the impact parameter, i.e., the minimum eigenvalue of M αβ , can become small, and zero modes could potentially arise.
Diagonal elements: In general, the sum over m, n in the second term above for α = β involves a sum of many g functions along with incommensurate complex numbers m|T α |n and n|T β |m . However, the diagonal of M αβ obeys and so the sum of g functions appears here with all strictly non-negative coefficients. At this point, we invoke ETH, which in this context states that (for m = n) | m|T α |n | 2 ∼ 2 −N/2 |r α,mn | 2 , (5.9) where r α,mn is a random matrix with entries of O(1) magnitude whose squared elements |r α,mn | 2 are all roughly equal and O (1). What this means is that the sum m,n m =n (g(t(E m − E n )) − 1) , (5.10) must become O(2 N ) before the diagonal entries M αα can vanish. This will only occur when almost all of the g functions are close to zero, which can only happen when t e S .
More quantitatively, let us try to approximate the timescale at which this occurs. Notice that we can expand the sum above to include m = n, since these terms have g(0) = 1. Then, we must determine when the sum m,n g(t(E m − E n )) becomes small, i.e., O(1). At large N , we can approximate the double sum as a double integral over two copies of the spectral density ρ(E): Strictly speaking, we should use the SYK spectral density for the function ρ(E). 16 However, we expect that our conclusions about conjugate points should apply to other chaotic systems as well. The key feature of the spectral density for q = 4 SYK is that there is an exponential number of states, e S = 2 N/2 , within a polynomial size window −N ≤ E ≤ N . The precise size of the window is not important for the argument, only that it is polynomial in N . Similarly, the relevant information about the exact height of the spectral density is that it is exponential in N . These properties also hold in e.g. a microcanonical ensemble of black hole microstates, where the window is actually O(1) in size with O(e N ) states. Since we are only interested in these very coarse features of the spectral density, we may approximate ρ(E) above by a constant distribution on −N ≤ E ≤ N : Of course, for sufficiently abnormal models, this density will not be a good approximation, but for chaotic SYK or a black hole microstate ensemble it is sufficient. The result of the integrals is where Si(x) ≡ x 0 dz sin(z)/z. The above estimate is generically an underestimate because the ansatz of a constant spectral density gives additional support to pairs of eigenvalues E m and E n which have separation larger than O(e −S ). The most important feature of (5.13) is that the function Si(x) ≈ π for x 1, so M αα is bounded away from zero by roughly Notice that we did not assume anything about the structure of the matrix r α,mn in making this argument. We only needed the entries to be distributed so that the squares |r α,mn | 2 took roughly the same O(1) value for any m and n, but the entries themselves did not need to be independent random variables. This is less than the usual statement about the ETH ensemble, where the variance of any given r α,mn is not only fixed, but the r α,mn themselves are all independent random variables.
Off-diagonal elements: Having understood the rough order of magnitude for the diagonal entries of M αβ , we now turn to the off-diagonal pieces. For these, we have again a sum of g functions from equation (5.7), but now the coefficients in the sum can be negative. We can get some rough intuition for the order of this quantity by again invoking ETH on the local operator matrix elements.
m|T α |n n|T β |m ≈ e −S/2 r α,mn r β,nm , If the r α,mn were drawn from independent Gaussian distributions with mean zero and O(1) variance, it would be straightforward to compute the typical (expectation) value of the above expression. We would simply find zero for the typical value since r α and r β are independent matrices and have mean zero. To ensure that the fluctuations of this quantity are not excessively large, we could also estimate the variance, which involves a calculation of r α,mn r β,nm r α,m n r β,n m in the aforementioned ensemble. Since r α,mn and r β,mn are independent, this four-point function factorizes into a product of two-point functions. ETH would then tell us that these two-point functions r α,mn r α,pq are proportional to δ mp δ nq since the entries of r α,mn are supposed to be independent Gaussian random variables. Going back to (5.7), we thus conclude that the off-diagonal entries of M αβ are always O(1). However, it cannot be precisely correct to employ ETH in this manner for any choice of eigenstates |m and |n because the operators T α have a known spectrum (all eigenvalues are ±1) which greatly differs from the spectrum of a random matrix with independent Gaussian random entries at large N . So, we will need a different sort of ensemble to get a consistent estimate of the mean and variance of M αβ for α = β.
One candidate which is consistent with all constraints on the matrices T α is the Haar ensemble of unitary matrices employed in the following manner. We pick some fixed basis |i P (for instance, the Pauli basis) in which the form of T α is known by construction to be relatively sparse or simple. Then, we assume that the eigenvectors |n E of the chaotic Hamiltonian H can be roughly thought of as a Haar random unitary rotation of this basis via 17 The off-diagonal terms in M αβ are given by 17 A similar ensemble was used to model a microcanonical window of states in quantum gravity in [65], although in that context the ensemble had a physical interpretation as the dual of a gravitational path integral in the spirit of [66]. Here, by contrast, we use the Haar ensemble to extract information about the typical value and variance of certain matrix elements with the understanding that we are really studying the expected behavior of a quantum chaotic system with fixed Hamiltonian, such as a single instance of the SYK model.
We would like to get an estimate for the mean value of the quantity To compute the typical value, we integrate this expression over the Haar ensemble for U by making use of (5.19) The asymptotic forms [67] and exact expressions [68] for such integrals are well known. With an eye toward the sums over m and n in M αβ , we notice that any term with δ nm must vanish in the full expression since m = n. We obtain (writing · H for the Haar expectation) (5.20) The Haar integration has given us the typical value of M αβ in terms of traces of the operators T α and T β . By construction, we have tr T α = 0 and tr(T α T β ) = 0 for α = β, so the typical value in a chaotic Hamiltonian ensemble defined this way is Incidentally, this calculation also shows that the diagonal terms α = β have a Haar average of order e S until exponential times. By setting α = β in the large parentheses of (5.20), we conclude which is consistent with our estimate that relied on the ETH ensemble (5.9), so the conclusions from that discussion concerning M αα match the results of the Haar ensemble.
If the off-diagonal elements of M αβ are all approximately zero for a given chaotic Hamiltonian, the only way a zero mode can arise is by the vanishing of a diagonal element, which we have shown does not occur until t ∼ e S . To be complete, we should also study the variance M 2 αβ H and ensure it is not too large. A small O(1) variance will ensure that fluctuations in the off-diagonal elements are small relative to the diagonal magnitude.
where α = β and ∆ mn ≡ E m − E n . The basic quantity which we would like to integrate against the Haar measure is The relevant Haar integral is On general grounds, the overall result for the Haar expectation of (5.24) will be written in terms of traces or products of traces of the operators T α , T β , T α , and T β . There are three such combinations which can be nonzero: We deal with each of these three case by case.
The first combination in (5.27), the double trace factor yielding 2 N , is produced by certain products of delta functions from the Haar integral The second combination in (5.27) can be formed with a variety of delta function combinations appearing from the Haar integral. Fortunately, because the trace factor is only e S in this case, the only possible dangerous term which may contribute beyond O(1) must take the form O(e −4S )δ iq δ jp δ ks δ r , (5.29) which is the unique term that appears at O(e −4S ) Haar suppression without any additional delta functions which would cancel the sums in (5.23). However, a term of this form does not lead to tr(T 2 α T 2 β ), but instead leads to tr(T α T β ) tr(T α T β ), which vanishes. So, the leading contribution of the Haar integral to the coefficient of the second term in (5.27) is O(e −5S ), and this is enough to absorb the e S trace factor appearing in all e 4S terms of the four sums in (5.23), yielding an at most O(1) contribution to M 2 αβ H .
The third and final combination in (5.27) must also contribute at most an O(1) result to M 2 αβ H , as the same argument concerning the unique form of the possible dangerous term holds in this case as well, since the trace factor is again only O(e S ).
Putting it all together, we have shown that the Haar average of M 2 αβ , assuming the eigenvectors of our chaotic Hamiltonian are related to some simple basis by a Haar-random unitary transformation, is An analogous argument shows that the diagonal variance is similar, . The next-to-leading term from the trace factor generated by (5.29) actually contributes O(1) to the diagonal variance rather than O(e S ), since we will have O(e −5S ) Haar suppression along with at least one delta function to cancel one sum in the analog of (5.23) for α = β. This is because (5.29) is the unique permutation leading to tr(T α T β ) tr(T α T β ), which is the first of only two new non-vanishing trace factors when α = β, without any such delta functions. The permutation leading to the second new pairing tr(T α T β ) tr(T β T α ), where the first T α is multiplied instead with the second T β in (5.23), comes with two delta functions δ nn δ mm at O(e −4S ) and one delta function at O(e −5S ), just as in (5.28), so there are only O(1) contributions due to this trace factor. The other non-vanishing trace factors for α = β are all captured by the three cases in (5.27), and the suppression arguments we made for those when α = β also apply when α = β. Thus, all diagonal variance contributions are O(1) as claimed. Note that the diagonal variance may have some mild dependence on t; here we have only argued that it is O(N 0 ). The numerical structure of M αβ at large t is shown in Fig. 4.
As we discussed above, this estimate is sufficient to argue that there should be no zero mode of M αβ (t) before times t ∼ e S , as the diagonal of the matrix is overwhelmingly large compared to the off-diagonal elements, and in addition, the matrix size scales as poly(S). Thus, for t < e S , the impact parameter will be O(e S ). On the other hand for t > e S , the diagonal elements of M αβ (t) are O(1) and in particular of the same order as the off-diagonal elements; we thus expect the impact parameter to become small (see Fig. 5).

Summary
We have argued in this section that the minimum eigenvalue of M αβ (t) must be O(e S ) for t < e S , but becomes small for t > e S . This implies that local conjugate points in chaotic theories can potentially occur only beyond exponential time. Even if exact zero modes of M αβ do not occur, we expect the minimum eigenvalue λ min of M αβ to become very small after t ∼ e S (see Fig. 5). Physically, this means that for t e S , it is possible to find an infinitesimally nearby curve with a local initial velocity V (0) = Ht + δV (0) which satisfies the geodesic equation up to O( 2 ), such that the end point is very close to the target unitary e −itH : Thus, it becomes possible to approximate time evolution by an infinitesimally nearby geodesic with a small error after exponential time. If the impact parameter λ min is exactly zero for some t * > e S , then we have a conjugate point at that location, and then we can find a shorter geodesic path to e −itH exactly, with no error.
Our arguments in this section were based on ETH and random matrix theory. A fairly similar story was told for the complete super-operator under the Eigenstate Complexity Hypothesis assumption in [17] (also explored in Sec. 7), but there are two key differences here. First, since there are only polynomially many entries in M αβ , we need not worry about the off-diagonal entries "backreacting" on the diagonal to force an unexpected zero mode at early times. Instead, a zero mode can only occur when a significant portion of the diagonal becomes suppressed at the same order as the off-diagonal entries, and this does not occur until t * ∼ e S . Second, the zero modes which arise in this way are actually independent of µ, and thus are fixed obstructions to the complexity growth of even maximally chaotic systems. We speculate further on the implications in Sec. 8.

Numerical analysis of conjugate points
We now present numerical calculations of conjugate point locations for free, interacting integrable, and chaotic SYK models. The general method that we use is to explicitly construct a matrix representation of the super-operator and compute its smallest eigenvalue (i.e., the eigenvalue with the smallest absolute value) using the Arnoldi (iterative) algorithm [69]. This gives us a concrete, albeit numerical, method to study obstructions to complexity growth. We will limit ourselves to systems up to N = 8 (four qubits) for computational feasibility, but in principle this method is not limited to small N .
We first reproduce the general expression for the super-operator from previous sections for reference, (6.1) There are two key observations that allow us to represent the super-operator more efficiently. The first is that the integrand simplifies immensely if we construct the super-operator in the basis of {T α ,Tα} where the T α are a basis for the local subspace and theTα are the basis for the nonlocal subspace that diagonalizes [H, · ] N L with eigenvalues λα. In that case, the sums disappear and we need only consider the first term or the last two terms depending on the column of the matrix representation that we wish to construct. The second observation is that the integral can be done analytically provided that we express the basis {T α ,Tα} in the energy eigenbasis |m n|. Note that this is not the same thing as writing the super-operator in the energy eigenbasis, which would not respect the split into local and nonlocal terms. In the energy basis these operators have coefficients: The c mn , as well as the energy eigenstates |m , their corresponding eigenvalues, and the diagonalization of [H, · ] N L , can all be precomputed before constructing the super-operator. The plots of minimum eigenvalue versus time below with N = 6 take one to two minutes per curve to generate on a four-core desktop while at N = 8 they take one to two hours per curve.

Free SYK
We first recall the key point of Sec. 3 regarding the location of conjugate points in the free model. When H is quadratic, the adjoint action of the Hamiltonian does not mix the local and nonlocal operator subspaces. Consequently, we can diagonalize [H, · ] in the local and nonlocal subspaces independently with corresponding eigenvalues λ α and λα. The super-operator then reduces to the simpler expression In the basis of operators {T α ,Tα} the super-operator is therefore diagonal with eigenvalues given by the coefficients above. Consequently, we have two 18 zero modes whenever So in the free model, every conjugate point is associated with an individual easy or hard operator. The easy conjugate points never move with µ, while the hard ones are occur at exponential times when µ is O(e αS ).
As computed numerically, the minimum eigenvalue of the N = 6 free theory is shown for several values of µ in Fig. 6, where we take k = 2 to match the locality of the Hamiltonian. Conjugate points occur where the minimal eigenvalue of Y µ touches the x axis. The conjugate point locations as displayed in Fig. 6 exactly match the analytic expression in (6.7). One can clearly see the shifting of several conjugate points as µ is increased; for example, the first conjugate point near t = 0.7 at µ = 0 gets shifted to three times its value, near t = 2.1, when µ = 2, and subsequently moves off the right end of the figure for larger values of µ. Most of the curves overlap for µ > 2 since once µ is sufficiently large only the local conjugate points, whose locations are not functions of µ, remain in a bounded-time region.

Integrable and chaotic models
We now compute the locations of conjugate points where we deform the free Hamiltonian as in Sec. 4 by H = H 0 + δH. Since the numerics are not restricted to taking to be perturbative we will consider = 1.0 in all plots in this section to illustrate large effects of each type of interaction. 19 We will consider three different choices of δH with the same base Hamiltonian H 0 considered across all cases at fixed N , integrable 4-body, (6.8) J ijk ψ i ψ j ψ k ψ chaotic 4-body, (6.9) Notably, δH 2 and δH 3 are effectively the maximally chaotic SYK 4 and SYK 3 theories [10] while δH 1 is the integrable interaction from Sec. 4. 20 The results for δH 1 , δH 2 , and δH 3 are displayed below in Figs. 7, 8, and 9, respectively. We compare the plots of minimum eigenvalue versus time for µ = 0 and µ = 10, where we have chosen µ = 10 as a numerically feasible upper 19 One reason to keep the H0 term is that the algebra generated by adH acting on a basis of operators has too trivial a structure at small N when H is chosen to be only a single q-local term, which can cause unwanted numerical coincidences. 20 We drop the numerical prefactor of 1/4 on δH1 that was written in Sec. 4, but draw Mij rather than Mij/4 from the q = 4 SYK distribution with J = 1 in this section, so the numerical results in each section are on equal footing.
bound to illustrate the large µ regime. 21 Clearly, in all three cases, a number of conjugate points remain fixed, corresponding to local eigen-operators of the super-operator (which were the main subject of Sec. 5 for large N chaotic theories), while those corresponding to eigenoperators with some nonlocal component will tend to shoot off with some µ-dependent speed. One can see this by examining the zeros of the below plots, which correspond to zero modes of the super-operator and thus conjugate points. Zeros that remain fixed in location between µ = 0 and µ = 10 are independent of µ and therefore correspond to local eigen-operators, while those that move do not.  point and e −iHt fails to be the globally minimizing geodesic. This is remarkable because the complexity geometry corresponds to the manifolds SU (8) and SU (16) which are already very complicated fiber bundles over spheres. What this illustrates is that even in such highly nontrivial geometries conjugate points can become relevant for complexity growth almost immediately. This fact is also interesting since it implies that interacting qubit Hamiltonians can be fast-forwarded at early times.
(b) There exist a family of conjugate points whose times are independent of µ. These correspond to local eigen-operators of the super-operator, since the expression (2.8) for the super-operator shows that eigen-operators with a nonlocal component will generically have eigenvalues that are functions of µ. These nonlocal eigen-operators clearly correspond to the conjugate points that move to later times as µ is increased. The existence of the local eigen-operators is surprising; as the subspace of local operators is small compared to the space of all operators, one might have expected that the typical eigen-operator generically contained nonlocal pieces.
(c) The size of the nonlocal subspace controls the density of nonlocal conjugate points and possibly also the speed at which they approach later times as µ is increased. This is visible in the greatly increased density of conjugate points at early times in Fig. 12, where the degree of locality is k = q = 3 in contrast to the other two cases that take k = q = 4. It appears that many of these conjugate points rapidly shoot off to late times whereas in the other two cases many of the nonlocal conjugate points appear to level off quickly.
Notably, there is not a large distinction between the integrable interaction δH 1 in Fig. 10 and the 4-body chaotic interaction δH 2 in Fig. 11 with regard to the behavior of the conjugate points. It appears that the degree of the locality of the Hamiltonian is the most significant factor in these small N plots.
We note here that we have chosen k to be the same order as the locality of the Hamiltonian, both so that the linear geodesic V (s) = Ht is explicitly a solution to the geodesic equations as well as so that the energies of the system are an "easy" observable to measure. However, we could have also chosen k to be smaller, since k = 2 is sufficient for the geometric complexity to approach the true quantum circuit complexity [22]. In this case, the nonlocal subspace of operators is substantially enlarged and we expect that at even smaller values of µ the nonlocal conjugate points occur at late times.
Although we take N to be small for computational feasibility, we emphasize that in the large N limit the size of the nonlocal subspace vastly exceeds the size of the local subspace. The size of the nonlocal subspace scales as O(e N ) while the size of the local subspace scales as O(poly(N )) regardless of k. In this limit, there will be many more nonlocal conjugate points than local, in contrast to what we see in the N = 6 plots.
(d) The first conjugate point is rapidly followed by many more conjugate points, both local and nonlocal. This substantiates arguments detailed in Sec. 3 that the occurrence of the first conjugate point is rapidly followed by the end of the linear regime for complexity, even though many conjugate points may be required to reach the plateau regime.
That the local conjugate points occur at times which are independent of µ might seem prima facie to be at tension with the result in [17] that the first conjugate point should not occur until times of order µ ∼ e αS for Hamiltonians respecting the Eigenstate Complexity Hypothesis (to be discussed in Sec. 7). However, there is no real contradiction. Firstly, the results of [17] apply to large N systems, and the present numerical studies are at small N . Secondly, as we increase N , the µ-independent time scale at which such local conjugate points occur cannot be sub-exponential in chaotic systems. Indeed, in Sec. 5, we gave a general argument for chaotic Hamiltonians based on the Eigenstate Thermalization Hypothesis [70] and random matrix theory ideas that this time scale should be exponential, resolving the apparent tension. This argument need not apply to integrable Hamiltonians, which may still encounter conjugate points at early times in the large N limit.
We also note that Claim 2 in Sec. 2.2 shows that in order to identify the locations of local conjugate points, we need not compute the full super-operator. Instead, we can compute a smaller matrix of polynomial size, M αβ (t) = where α, β only index the local operators. The zero modes of this matrix M αβ then identify the times of local conjugate points with substantially increased computational efficiency. In Fig. 13 we have demonstrated this by plotting the minimum eigenvalues of the super-operator and of M αβ . The zero modes of each clearly coincide, though the matrix M αβ also appears to be more numerically stable in the sense that the precision of the numerical zero modes locates them closer to exactly zero. 22

Summary
We computed the zero modes of the super-operator numerically for N = 6, 8 fermions by applying the Arnoldi algorithm to compute the minimum eigenvalue of its matrix representation, for various choices of Hamiltonian including free, integrable interacting, and chaotic SYK models. We demonstrated that a large class of conjugate points corresponding to local eigen-operators of the super-operator remain fixed, while those which have nonlocal compo-

Eigenstate Complexity Hypothesis
From the geometric formulation of quantum computation, we have found evidence that there is a difference between chaotic and integrable models in the location of obstructions to complexity growth. To argue for this behavior of the complexity of time evolution in chaotic theories, [17] found that a certain matrix R mn of operator matrix elements in the energy eigenbasis was relevant: The T α are the easy operators (at most k-local), Tα are the hard operators (at least (k + 1)local), and |m and |n are energy eigenstates. In order for chaotic Hamiltonians to show the expected complexity growth, individual off-diagonal matrix elements (m = n) should be suppressed by R mn ∼ e −2S poly(S)r mn , where S is the log of the dimension of the Hilbert space and r mn are independent random numbers of order 1. In [17], the conjecture that (7.2) should hold for chaotic Hamiltonians was called the Eigenstate Complexity Hypothesis (ECH), and we refer to R in (7.1) as the ECH matrix. Intuitively, (7.2) states that the eigenstates of a chaotic Hamiltonian are sufficiently different that acting with a local operator on an energy eigenstate will not yield a different one even approximately.
To make the polynomial factor in (7.2) a bit more precise, we can compute the ECH matrix numerically in some simple chaotic systems; it is also interesting to compare to the integrable case. We will find that the average R mn (with m = n) is in general suppressed by roughly the fraction s(N, k) of operators that are k-local, s(N, k) = (number of k-local operators) (total number of operators) .

(7.3)
This holds regardless of whether the dynamics are chaotic or integrable, although higher moments of the distribution of R mn show more of a sensitivity to the dynamics.
For SYK, this means poly(S) in (7.2) is roughly N k . If there are extra symmetries in the problem, which is true in the case of q = 4 SYK, some elements R mn could be slightly reinforced, but in general R (the average of all off-diagonal R mn ) should still be approximately equal to (7.3).

Pure SYK
Here we check the ECH criterion in a general SYK q model (3.2) for different choices of N and q. Recall that N > q, and we also take k ≥ q so that the Hamiltonian is a simple operator.
We first verify the integrability/non-integrability of the model using level spacing statistics: the integrable q = 2 model shows a clear exponential distribution (i.e., Poisson distribution for the energies) while for q = 3, 4 the model is chaotic and shows Wigner-Dyson statistics (Fig. 14). These distributions match the expected level spacing distributions from random matrix theory [71].
In the SYK model, the ratio (7.3) is exactly For q = 2, since the Hamiltonian is free, we would expect the ECH criterion to fail, i.e. the off-diagonal entries should not be suppressed by the ratio defined by (7.3). Visualizing R mn for q = 2 in Fig. 15, we see that there are indeed large off-diagonal elements that are comparable to the diagonal itself. By contrast, for the q = 3 model (Fig. 16), the off-diagonal elements are fairly uniformly suppressed with the exception of the anti-diagonal entries.
In the higher q models, the fermion number operator F becomes relevant for interpreting our results.  and odd fermion number) on the eigenvectors of H. The conserved charge for q = 4 leads to additional suppression and enhancement within the ECH matrix. This can be explicitly seen in that there are two dominant colors in the ECH matrix across N = 10, 12, 14, 16 ( Fig. 17) when the eigenvectors are organized into even and odd symmetry sectors. The bright colors (larger magnitude) in the diagonal blocks indicate that two eigenstates within a given symmetry sector will have enhanced overlap with k-local operators. The darker colors (smaller magnitude) in the off-diagonal blocks indicate that only small overlaps can be achieved by acting with local operators between states in different symmetry sectors.
We can get a better quantitative picture of the differences between q = 2, 3, 4 by studying the distribution of all of the off-diagonal entries in R mn (Fig. 18). For q = 4, we do see the two separate peaks corresponding to the two superselection sectors generated by the charge F , but the average of all off-diagonal entries (roughly the point in between the peaks, ≈ 8.8 × 10 −2 ) is still suppressed by the ratio (7.5). For q = 3, we have no additional symmetry, and we see that the off-diagonal values of the matrix R mn are all suppressed by approximately the same value (poly(S)e −2S ), following (7.2).
Though we have claimed the ECH matrix should distinguish between integrable and chaotic theories by the magnitude of off-diagonal elements, one cannot discern the integrability of the model from the averaged value of R mn (Fig. 19). The off-diagonal values of both the free q = 2 SYK model and the chaotic q = 3, 4 models are suppressed on average in Fig. 18 by the ratio (7.5) with system size N = 14 and the degree of locality k = 4 fixed. The key distinction is that in the q = 3 case, the R mn are continuously distributed with a small standard deviation (Fig. 19). This is also true for q = 4 within the same superselection sector (a) Figure 16: The R mn matrices of the q = 3 SYK model with up to k = 3-local operators taken to be easy. As one increases N , the overall R mn values are suppressed according to (7.3). of the charge F . The cluster centered at the larger averaged R mn in Fig. 18 consists of the overlaps of the eigenvectors within the same superselection sector. The more suppressed R mn cluster comes from the overlaps between eigenvectors across different sectors. By contrast, for q = 2, R mn essentially has a delta function distribution at discrete points and a very large The degree of locality is controlled at k = 4. Notice that the q = 3 case has a small spike at less than the majority of the off-diagonal values. The spike comes from a more suppressed anti-diagonal values. standard deviation consequently. The masses at different values suggest multiple symmetries in the q = 2 SYK model, and the discrete support reflects the integrability of the free theory, where a very small number of parameters control the energies and correlations of the model. The overlaps in the q = 2 case with N = 14 have delta functions at five major values (Fig. 18). Notably in the free case, a significant number of the overlap matrix elements R mn are 0. This leads to smaller average value of R mn in Fig. 19. This together with the delta-function at outlying value indicates a large coefficient of variation at large N (Fig. 19). Similar observations have been made regarding the matrix elements of local operators in the energy basis of an integrable system, m|A|n , in the context of the Eigenstate Thermalization Hypothesis [71,72]. In the case of the XXZ-chain, for example, the off-diagonal entries of m|A|n are generically either highly suppressed or very large, whereas the entries in the non-integrable case show a more uniform distribution.
We have focused on the SYK model here, but see Appendix C for a companion study of the ECH matrix for the mixed-field Ising chain of length N = 7, 8, 9, 10, which similarly features regimes of chaos and integrability.

Deformations of free SYK
In addition, we consider adding integrable or chaotic deformations to the free model as in Sec. 4 and Sec. 6. Adding an integrable term (6.8) to the free Hamiltonian with Gaussian random coupling M ij does not change the level spacing distribution, regardless of the coupling strength (Fig. 20). By contrast, adding a 3-local chaotic perturbation (6.10) to the Hamiltonian affects the energy spacings. As the coupling becomes large, the level spacing  Figure 21: Level spacing for free SYK with an chaotic cubic perturbation (6.10) (a) when the coupling is small and (b) when the coupling is large. As the chaotic term becomes comparable in magnitude to the free term, the level spacing statistics shift from exponential to Wigner-Dyson.
However, both the integrable and non-integrable perturbation to the Hamiltonian change the distribution of the off-diagonal entries of R mn . At very small perturbation parameter, the peaks concentrated at R mn values of the free SYK model widen. When the perturbation parameter becomes large, the R mn form distributions centered around new values that reflect the symmetry of the system (Figs. 22 and 23). Since the integrable Hamiltonian commutes with the fermion number operator, its R mn distribution is bimodal corresponding to the two symmetry sectors; the cubic chaotic SYK model has no such symmetry as discussed before, so the distribution is unimodal. Apart from the existence of these symmetry sectors, the integrable and chaotic perturbations do not show qualitative differences in their distributions. Indeed, the R mn distribution shows a bigger contrast between free systems and interacting systems, rather than between integrable and non-integrable systems. Similar observations can be made in the R mn distribution of transverse Ising model (Appendix C), where the free spin chain also has a delta-function like probability distribution. The integrable but interacting spin chain, analogous to the SYK model perturbed by an integrable deformation, instead clusters around values that reflect the symmetry of the system.
The number of conserved quantities in an integrable theory scales extensively with the number of degrees of freedom, so we might have expected that the distribution of R mn in free and interacting integrable theories would be similar. However, free and interacting integrable  theories are distinct in that free Hamiltonians can be diagonalized into a sum of O(N ) independent terms whereas interacting integrable theories cannot in general. Because of this, the 2 N eigenstates of a free theory carry redundant information about the system, and this is why the values of the overlap matrix R mn show delta-function support. This reduction in the number of effective variables required to characterize R mn is already impossible for interacting integrable theories.

Summary
From the simulations of the free, integrable, and chaotic SYK model with varying system sizes, we have observed that the mean of the off-diagonal ECH matrix R mn is suppressed by the ratio (7.3) in every type of system. The free systems have outliers compared to the mean but also a large number of zero entries that "balance out" the large entries. However, the standard deviation in the free case scales linearly with the system size N while it remains approximately constant in the chaotic cases. From the distribution of the R mn , one can distinguish the free systems from interacting integrable and chaotic systems by the discreteness of their support, but the latter two are difficult to differentiate. Many characteristics of the ECH matrix have also been observed in the ETH matrix. The reason why interacting-integrable and chaotic theories have similar R mn distribution and the precise relation between the ECH and ETH matrices will be explored in future works.

Discussion
In this paper, we studied conjugate points and geodesic loops in various SYK models. In the free model, we located all conjugate points and characterized the family of geodesic loops associated to the local conjugate points. This allowed us to exactly compute the complexity, which is bounded by O( √ N ), and to specify the fast-forwarding Hamiltonian at subexponential times. We studied the motion of conjugate points with µ under the addition of integrable or chaotic interactions both analytically and numerically. In the integrable case, we first showed how to set up perturbation theory for conjugate point locations in the strength of the coupling constant controlling interactions. We also described a family of geodesic loops which bound the complexity by O(poly(N )) for the class of integrable systems we considered. We then studied local conjugate points in chaotic theories. We argued based on the statistics of the matrix M αβ (t) that such local conjugate points do not occur in chaotic systems before exponential time, thus strengthening the arguments given in previous work [17]. We then studied the locations of conjugate points non-perturbatively using numerics for SYK models up to N = 8. Finally, we explored the Eigenstate Complexity Hypothesis (introduced in [17]) in free, interacting-integrable and chaotic SYK Hamiltonians. We view these results as demonstrating a hierarchy of complexity growth between free, integrable, and chaotic models, and as a preliminary attempt at describing a more complete picture of the growth of complexity of time evolution, in which conjugate points play an essential role. Of course, global loops (which are not signaled by conjugate points) should also play an important part in this story, and perhaps they could even obstruct complexity growth before conjugate points. A more complete picture of complexity growth must therefore also include these.

Quantum error correction and AdS/CFT
The modern picture of bulk reconstruction in AdS/CFT involves interpreting the bulk-toboundary map as an isometry which encodes bulk "logical" degrees of freedom within the set of CFT "physical" degrees of freedom in an approximate quantum error correcting code (QECC) [73]. This picture sheds light on several subtle issues in bulk reconstruction, such as the fact that a single bulk operator can have multiple distinct boundary reconstructions on different subregions of the boundary. However, just as there is some ambiguity in the definition of quantum complexity, there are some choices to be made in the definition of the QECC. One such choice is the notion of the code subspace, which is usually taken to be (roughly speaking) the subspace of states which correspond to a bounded number of bulk operator insertions around a semiclassical background. The QECC then reconstructs bulk operators within this Hilbert subspace, rather than on the full CFT Hilbert space. This allows the AdS bulk to incorporate, for instance, radial locality in the form of commutation between a bulk local operator and a boundary local operator [74].
It has always been relatively ambiguous what precisely the code subspace ought to be in a given situation. There are known restrictions on, for instance, what fraction of black hole microstates in a single microcanonical window may be included in a code subspace while keeping the error in the approximate bulk reconstruction under control [75]. Relatedly, the choice of simple operators in geodesic complexity is somewhat ambiguous. One way to construct the code subspace from a CFT perspective is to start with some holographic state (say, the vacuum) and act on it with a few, not-too-heavy single-trace operators. The span of such states forms a subspace which one could regard as the code subspace. Taking inspiration from this idea, one could regard as simple operators (from a complexity perspective) the span of such not-too-heavy single-trace operators which generate the code subspace. This ties together the complexity-theoretic notion of locality and the error correction notion of locality.
One speculative way of operationalizing these ideas in the context of conjugate points is the following. Suppose we take a CFT state which corresponds to a small number of light operator insertions in some background state. In the bulk, this creates some particles near the boundary in some semiclassical geometry. Now one considers time evolution on the boundary by the (local) boundary Hamiltonian H. At very late boundary times t, we expect that the corresponding linear geodesic e −iHt encounters a conjugate point. After the conjugate point, a new globally minimizing geodesic takes over, which may correspond to evolution with a nonlocal effective Hamiltonian as we have discussed in this paper. Meanwhile, in the bulk this time evolution is dual to a scattering process between the particles that were created near the boundary. Assuming this scattering process did not create a black hole, the local time evolution on the boundary does not take our initial state out of the code subspace. Since the new minimizing boundary geodesic at late times does not lie in the code subspace, we expect it corresponds to a time evolution in the bulk involving black holes, as it will take the initial state out of the code subspace. This suggests that the late-time out states of the scattering process in the bulk could have been reached more efficiently by a scattering process involving black holes in the bulk.

Remarks about switchback effect
The geometric complexity theory that we studied was constructed to be polynomially equivalent to quantum circuit complexity [22]. This required choosing the cost factor µ to be exponential in the system size, µ ∼ e S . However, there are arguments from AdS/CFT involving the so-called switchback effect [2] which appear to imply that the notion of complexity which is relevant for holography is compatible with a more gradual weighting scheme [52]. 23 The graduated scheme roughly involves taking operators below a locality threshold k to have cost 1, and operators above this locality to have weight equal to the exponential of their degree of locality. So, a K-local operator for K > k would have weight e K .
A potential issue with the graduated scheme is that it allows (log N )-local gates to act with a polynomial cost. In the more conventional formulation of complexity theory, a fixed upper bound k ∼ O(N 0 ) is chosen on the locality of polynomial cost gates. However, if we interpret the graduated scheme as setting k ∼ c log N instead of k ∼ O(N 0 ), it is in fact still polynomially equivalent to the standard scheme where only O(N 0 )-local gates have polynomial cost. To see this, notice that an arbitrary unitary operator acting on c log N qubits has complexity at most roughly e c log N = N c , which is still polynomial in N . Therefore, in a polynomial-length circuit formed using the graduated scheme, we may simply replace any (log N )-local gates with polynomially many O(1)-local gates without changing the fact that the total circuit length is polynomial in N .
In fact, the graduated scheme assigns a cost of N to a (log N )-local gate, which is precisely the same as the maximum cost of a (log N )-qubit unitary operator if we had used O(1)-local gates in the standard scheme. Given this observation, it is not hard to imagine that the graduated and standard schemes are actually related by some O(1) factor rather than only a polynomial. Since there are differences in sectional curvature between the graduated and standard schemes [55], it would be interesting to understand whether these differences really appear at the level of the complexity of time evolution. It may be that they are only related to O(1) prefactors in that quantity, and more significant differences can only be seen for more complicated quantities like the complexity of a precursor [2].

Relation of ECH to ETH
The behavior of the R mn is similar to the behavior of matrix elements of local observables in the energy basis that show up in the Eigenstate Thermalization Hypothesis (ETH) [70,76]. The ETH matrix in integrable systems features off-diagonal elements that are mostly either very small or zero, with a sparse number of large values. By contrast, the ETH matrix in chaotic models has more uniformly suppressed entries. We have demonstrated analogous properties in the ECH matrix R mn for both the integrable and chaotic SYK q models (Fig. 18, 22, 23) as well as the mixed-field Ising models in Appendix C.
Although the statistical properties of these matrices are similar, the physical settings behind the ECH and ETH statements are different. ECH is founded in ideas of circuit complexity whereas ETH is based on ideas in many-body physics and thermalization. Regardless, the mathematical expression of ECH involves the matrix elements of local operators in the energy basis that appear in ETH. Specifically, one might wonder if there is a relation between ECH and ETH applied to the local operator α T α which is the sum of all local operators. One approach that might be fruitful in understanding the statistics of the matrix elements of such an operator is to consider perturbing a Hamiltonian with a perturbation δH = α T α . In this case, it may be possible to understand the statistics by relating the eigenstates of the perturbed and unperturbed Hamiltonian since the corrections to the energies and eigenstates generically involve matrix elements like m|δH|n . A related thought experiment investigating the statistics of matrix elements of such perturbations was carried out in [76].
In Sec. 5 we showed that in a chaotic theory complexity will grow linearly until times O(e N ) provided we assume that the energy eigenstates of the chaotic theory are essentially Haar random rotations of a fixed "simple" basis. There is a related expectation in ETH, i.e., that all the eigenstates of a thermalizing Hamiltonian "look" thermal [70]. Indeed, Deutsch showed that for a real, symmetric Hamiltonian, a thermalizing perturbation leads to energy eigenstates that are Gaussian random linear combinations of the unperturbed eigenstates [76]. This averaging suppresses the variance of observables by factors of e S just like in the Haar ensemble we proposed for chaotic systems. In [77] entropic suppression of variance was also described for typical, random states in a quantum microcanonical ensemble. In Sec. 5 we are using similar reasoning to argue that typical energy eigenstates of a chaotic theory will be random combinations of a fixed "simple" basis, and so variances will be suppressed via averaging.

General integrable systems
In this paper we have considered an interacting integrable deformation which is a quadratic function of the local operators J (i) 3 in the diagonalized free theory. Consequently, the structure of the operator dynamics as governed by [H, · ] is simplified, which allowed us to obtain analytic results in perturbation theory in Sec. 4. General integrable systems can look much more complicated. For instance, in App. C we study the ECH matrix for the (integrable) transverse-field Ising model, which appears to be a nontrivially interacting lattice spin model. This model is equivalent to a free fermion model after performing a non-local Jordan-Wigner transformation taking the bosonic spins to fermions. It is natural to choose the k-local operators in the theory to be the bosonic spin operators supported only on contiguous size-k regions of the bosonic spin lattice for the purposes of computing complexity. However, the Jordan-Wigner transformation will not respect this split into local and nonlocal operators. Consequently, it is plausible that more general integrable theories behave similar to chaotic systems with respect to complexity, owing to the fact that local operators in the theory are scrambled into the nonlocal sector when the theory is diagonalized.
In full generality, integrable systems in a finite-dimensional phase space can be written in action-angle variables. Much like the dynamics of the harmonic oscillator (a canonical example of an integrable system) consists of "rotation" in phase space, the dynamics of these more general systems also consists of periodic motion in phase space, albeit with a possibly action-dependent frequency. 24 So there is no guarantee that integrable systems appear free in any basis. Infinite-dimensional integrable systems like the Korteweg -de Vries (KdV) 24 More precisely, the "actions" are first integrals of motion. For the one-dimensional harmonic oscillator, the action is proportional to the energy. For the harmonic oscillator, the constant-energy slices foliate phase space by a set of concentric circles, and the dynamics is rigid rotation around these circles. More generally the phase space of an integrable system need only be foliated by topological tori, with the dynamics corresponding to a periodic motion around each torus whose frequency depends on the corresponding action. In fact, the KAM theorem [78,79] guarantees that most of these tori are preserved given small deformations of the Hamiltonian, so we expect results that hold for integrable systems may also hold for perturbatively chaotic systems. system [80] or the solitonic Sine-Gordon system [81] exemplify this fact. Similarly, a broad class of highly interacting quantum integrable systems consists of the lattice spin models that are exactly solvable using the Bethe ansatz [82]. Nevertheless, the Hamiltonian in all these systems is built out of the commuting charge operators, and so we expect our analysis of Sec. 4.2 to be generalizable to such systems. It would be interesting to further explore whether it is possible to make precise analytic statements about complexity in these highly structured models whose notion of locality in the original variables does not align with locality in the variables that simplify the dynamics.

A version of complexity restricted to local circuit modifications
In AdS/CFT, tensor networks have proven useful in gaining intuition about properties of the bulk semiclassical theory [38][39][40][41][42][43][44]. Roughly speaking, the tensor network lives on a tessellation of a bulk Cauchy slice. There is an approximate notion of quantum complexity for tensor networks which corresponds simply to counting the number of tensors in a region of the bulk spacetime, and this supports the suggestion that a quantity like bulk wormhole volume should be dual to quantum complexity in a two-sided black hole [3]. However, once the complexity saturates at its maximum value (polynomial in the entropy for integrable systems and exponential for chaotic systems), the tensor network which grew to foliate the wormhole interior is no longer expected to be the minimal network, just as the quantum circuit which builds the state will become a non-minimal circuit. In the geometric language, the linear geodesic will encounter conjugate points or geodesic loops.
As physics is at least approximately local, for a geometric wormhole interior it would be surprising if there could be large correlated fluctuations of geometry which act in concert to decrease the tensor network size. Such large fluctuations with global changes to tensor network structure would correspond to geodesic loops in the complexity geometry which have little or no relation to the original linear geodesic. This motivates a notion of "local complexity", where only local updates to the tensor network (equivalently, the quantum circuit) which decrease the length are allowed. 25 Of course, if many sequential local updates are made, we can still achieve a large decrease in the size of the tensor network.
In the geometric complexity language, local complexity is computed by starting with the linear geodesic L and flowing downward in the space of paths, where the downward directions all correspond to conjugate points along L. These downward flows will explore the space of paths at least in the neighborhood of L, and will find the geodesic of smallest length which is continuously connected by upward flows in path space to L. To our knowledge, local complexity has not been explored, and may be an interesting alternative to the usual complexity-theoretic definition. We will not explore it in any great detail here, but we will make the following point: local complexity behaves more or less analogously to circuit complexity in chaotic theories like holographic CFTs. There are conjugate points along L which sit approximately at t ∼ e S , which will terminate the linear growth at the expected timescale just as geodesic loops would [17]. Furthermore, the density of conjugate points along the linear geodesic is roughly constant after an initial growth, as can be seen from a simple calculation in the bi-invariant geometry. So, we expect multiple exchanges of dominance between many geodesics induced by all of these conjugate points, which should generate the plateau.
Of course, after we encounter the first conjugate point, the remainder along the linear geodesic are not relevant for complexity growth since there is a new geodesic which computes the complexity. So, in order for this argument to hold, we need a sort of universality among geodesics under the flow from µ = 0 to µ ∼ e S which ensures that, even as the location of the geodesic changes in path space, the conjugate points which were present at µ = 0 are shifted in roughly the same way as occurs for the linear geodesic. That is to say, all geodesics at µ ∼ e S have a constant density of conjugate points around exponential length, just as we expect for the linear geodesic. We have not proven this, but it seems likely from general intuitions about the complexity growth of chaotic Hamiltonians [17].
One difference in these notions is that geometric complexity is always bounded by the diameter of the manifold, and local complexity may slightly violate this bound. However, it is unknown whether anything physical would be associated with such a slight modification of complexity's behavior. If local complexity is really the notion to consider in holography, it will have implications for the complexity of the AdS/CFT dictionary, following the arguments of [83]. This is because the calculation of local complexity is essentially a local optimization problem in path space, which should be polynomially computable in general, unlike circuit complexity which would involve searching the entire path space for potential geodesic loops. 26

Acknowledgments
We thank Steve Shenker for helpful discussion and for motivating us to study the free SYK model, and Pedro Bernardinelli for use of a personal server on which some numerics were performed. VB, MD, and CL are supported in part by the Department of Energy through grant DE-SC0013528. VB is supported in part by the Simons Foundation through the It From Qubit Collaboration (Grant No. 38559) and by the Department of Energy through grant QuantISED DE-SC0020360. VB also thanks the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611, for hospitality while this work was in progress. AK is supported by the Simons Foundation through the It from Qubit Collaboration. MD is supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1845298.

A Two-point functions and the super-operator
In this appendix we will write an interesting expression relating the super-operator to the matrix of thermal two-point correlation functions at infinite temperature, This identity can be derived simply by comparing the Taylor series of both sides. Recalling the Campbell-Baker-Hausdorff identity e xH Oe −xH = e ad xH O, this new identity is useful because it converts a commutator into a derivative that can be integrated by parts.
Using these substitutions and evaluating the integration by parts, after a short computation one can write the super-operator as The last two terms only involve the averaged R matrix, and so have a smooth limit for large t. The first term however involves the exact R matrix, which oscillates wildly at late times. It would be interesting for future work to explore further connections between the super-operator and thermal two-point functions.

B Symmetries of the ECH matrix in the SYK model
In this appendix, we analyze the symmetry properties of the SYK model to explain some of the structure of the ECH matrix found in Sec. 7. We reproduce the SYK Hamiltonian: where J i 1 ,...,iq ∼ N (0, σ 2 ), σ 2 = 2 q−1 (q−1)!J 2 qN q−1 and {ψ i , ψ j } = 2δ ij .
We now explain the curious fact noted in Fig. 17 that the off-diagonal matrix elements R mn of the ECH matrix are exactly symmetric between the even and odd sectors. There is no a priori reason that operator matrix elements in different superselection sectors should be exactly equal. As we will see, it follows from the form of the ECH matrix and the enhanced symmetries above. For the q = 4 model, the Hamiltonian commutes with F . Write |n, ± for the simultaneous eigenstates of H and F . These satisfy H|n, ± = E n |n, ± and F |n, ± = ±|n, ± . Let N = 2, 6 mod 8 (in particular, this includes N = 14). Then the operator T anticommutes with F . In fact, up to a phase, its inverse T −1 is actually just its adjoint T † , since it is a phase times a product of Majorana fermions which all square to 1. This means we must have where we have made the phase explicit and also written T i for any traceless Hermitian generator of su(2 N/2 ). A basis of such generators is given by the 2 N − 1 appropriately Hermiticized products (excluding the identity) of the Majorana fermions (B.2).
The crucial point now is that the time reversal operator exchanges the F superselection sectors due to the anticommutation relation: It may be the case that T permutes the energy levels by π; however, this does not change the conclusion because these permutations will cancel up to an overall phase in the expression T † T i T . Now all that's left is to analyze the matrix element magnitudes: (B.11) Therefore, the plus and minus charge sectors have identical ECH matrix elements. This explains why the two diagonal blocks in Fig. 17 are exactly equal, rather than only approximately equal. By Bott periodicity, this argument should also apply to N = 10, and we numerically verified that for N = 10, 14 the diagonal blocks are the same (up to machine precision) while for N = 12, 16 the diagonal blocks have similar but different values.
C Mixed-field Ising model This model can demonstrate either chaotic or integrable behavior depending on the regime of the parameter space (h, g). This can be observed from the level spacing statistics after exact diagonalization of the model with different choices of (h, g) [85] as shown in Fig. 24. We focus on three cases, two of which are integrable and one of which is chaotic. The first choice takes h = 0 but g = 0. This is the transverse-field Ising model, which is an interacting integrable model. Although it can be mapped to free fermions through the Jordan-Wigner transformation, the physical interpretation of the model is usually through interacting hardcore bosons, so we refer to it as the "interacting integrable" choice. The second model takes g = 0 and h = 0, which can be reduced to the Ising model in the absence of external fields, so this model is effectively free. The third choice takes h = g = 0 and is generally chaotic. We choose (h, g) = (0.5, −1.05) as a particular chaotic point in the space of couplings, (h, g) = (0, −1.05) as an example of interacting integrable couplings, and (h, g) = (0.5, 0) as the non-interacting integrable choice. For the integrable choices of couplings, the level spacings are exponentially distributed. By contrast, in the chaotic case the level spacing distribution is roughly Wigner-Dyson. When (h, g) takes intermediate values between any of these combinations, the level spacing distributions interpolate between the Wigner-Dyson and exponential distributions. This model also has nontrivial momentum sectors which can lead to zero modes if left unfixed; here we generally consider momentum eigenstates in the k = 1 momentum sector. We represent the spin chain operator algebra su(2 N ) with tensor products of Pauli operators σ x i , σ y i , and σ z i where i = 1, . . . , N is the site index. Since the spin-chain with nearest-neighbor interaction still retains a notion of spatial locality, we define simple operators to be operators which are supported on at most k sites (that is, they act as the identity outside a k-site contiguous region). For example, if we consider easy operators to be at most 2-local, σ z 3 σ z 4 is an easy operator, but σ z 3 σ z 6 is not.
We expect that the off-diagonal matrix elements of the ECH matrix R mn will be suppressed by (7.3), which for the Ising model is Below, we compute the ECH matrices for spin chains of length N = 7, 8, 9, 10 with k = 2local operators taken as simple. As in the case of the SYK model, we expect that in general the off-diagonal entries of the ECH matrix are suppressed according to (C.2), and that the variance in the chaotic regime is O(1). We numerically verify these properties in Fig. 26 taking k = 2-local operators to be simple to demonstrate the behavior. And this behavior persists independent of the choice of degree of locality k, as long as it is larger than the degree of interaction of the Hamiltonian.
In Fig. 27 we compare the ECH matrix elements between the different choices of couplings at fixed N = 9. Because the choice of the specific momentum sector has eliminated all the potential symmetries in the system, the ECH matrix entries does not show the same symmetry-sector structure as the q = 4 SYK model. The greater degree of suppression of  : ECH matrices of the mixed-field Ising models of length N = 9 at the free, integrable interacting, and maximally chaotic points, with k = 2-local operators taken as simple.
the off-diagonal matrix elements can be easily seen in the heat map. Another difference between the regimes comes from comparing the free theory to the two interacting theories. The distribution of the overlaps in the free theory is again a discrete distribution with a sizeable mass at 0. On the other hand, the overlap-distribution at the chaotic point and the interacting integrable point look more like a density function, analogous to what happens in the q = 3, 4 SYK models.
We have also investigated the behavior of the distributions of off-diagonal entries of R mn Figure 28: Distributions of ECH matrix elements for the mixed-field Ising model in the integrable, integrable-interacting and chaotic regimes, calculated using spin chain of length N = 9, with k = 2-local operators taken as simple. as N is increased in the various cases (Fig. 28). Because of the large number of vanishing overlaps, the average value of an off-diagonal entry of R mn in the integrable case is smaller than that of the chaotic case ( Fig. 29(a)). However, the variance in the chaotic case remains significantly smaller when N → ∞; it appears to grow very slowly and approximately linearly while the variance in the free theory grows approximately quadratically ( Fig. 29(b)).