Quantum Complexity of Time Evolution with Chaotic Hamiltonians

We study the quantum complexity of time evolution in large-$N$ chaotic systems, with the SYK model as our main example. This complexity is expected to increase linearly for exponential time prior to saturating at its maximum value, and is related to the length of minimal geodesics on the manifold of unitary operators that act on Hilbert space. Using the Euler-Arnold formalism, we demonstrate that there is always a geodesic between the identity and the time evolution operator $e^{-iHt}$ whose length grows linearly with time. This geodesic is minimal until there is an obstruction to its minimality, after which it can fail to be a minimum either locally or globally. We identify a criterion - the Eigenstate Complexity Hypothesis (ECH) - which bounds the overlap between off-diagonal energy eigenstate projectors and the $k$-local operators of the theory, and use it to show that the linear geodesic will at least be a local minimum for exponential time. We show numerically that the large-$N$ SYK model (which is chaotic) satisfies ECH and thus has no local obstructions to linear growth of complexity for exponential time, as expected from holographic duality. In contrast, we also study the case with $N=2$ fermions (which is integrable) and find short-time linear complexity growth followed by oscillations. Our analysis relates complexity to familiar properties of physical theories like their spectra and the structure of energy eigenstates and has implications for the hypothesized computational complexity class separations PSPACE $\nsubseteq$ BQP/poly and PSPACE $\nsubseteq$ BQSUBEXP/subexp, and the"fast-forwarding"of quantum Hamiltonians.


Introduction
In recent years the late time dynamics of general relativity have been examined through various lenses. Two of the most prominent directions in this subject deal with quantities whose classical behavior cannot possibly continue to hold into the asymptotic future due to fundamental quantum-mechanical obstructions. The first is the exponential decay of a CFT two-point function computed using classical gravity in an AdS black hole, which could break down at a time as early as t ∼ S (where S is the entropy of the black hole) as a consequence of the unitarity of quantum mechanics [1]. The second is the linear growth of the Einstein-Rosen bridge in the two-sided eternal AdS wormhole geometry, which led to a conjecture relating bulk volume/action and boundary quantum circuit complexity [2,3]. 1 If this conjecture is correct, then the extrapolation of the gravity result to times beyond t ∼ e S is expected to break down due to quantum effects in a finite-dimensional quantum gravity Hilbert subspace. Various studies of both quantum circuit complexity and correlation function behavior have explored these observations [4,5]. However, in the case of circuit complexity, despite the plethora of analytic results from gravity calculations (see [6][7][8][9][10][11][12][13][14][15] and references therein) assuming the volume/action conjecture, there has been little non-perturbative progress towards a first principles calculation of circuit complexity in CFT. In this paper, we seek to remedy this situation by studying the complexity of time evolution with chaotic Hamiltonians (which are expected to have gravity duals), especially with an eye towards the late-time behavior.
At present, the most accessible method to compute complexity in continuum quantum systems is Nielsen's geometric formulation [16][17][18][19]. In this approach, the circuit complexity of a unitary operator U is the length of the minimal geodesic on the unitary group joining the identity to U . One begins by classifying the Lie algebra of the unitary group into "local" or "easy" directions, represented by operators T α , and "non-local" or "hard" directions Tα. Typically, the local directions will consist of operators with less than k-body interactions, for some k. One then picks a right-invariant metric on the group U with the appropriate cost factors built in, such that motion along hard directions is disincentivized. The geodesic length with such a metric was shown to be polynomially equivalent to the usual notion of circuit complexity, which involves counting elementary unitary gates, provided the cost factors are chosen to scale exponentially with the Hilbert space dimension [18]. Heuristically, one can think of the circuit as a sequence of gates which corresponds to a sequence of geodesic segments on the unitary manifold; the geodesic in the geometric framework is then an everywheresmooth approximation to this piecewise-smooth curve (Fig. 1). In this work, we will be interested in this geodesic notion of complexity. This technique has been applied by various authors to compute complexity in several physical systems [20][21][22][23] (see also [24][25][26][27][28][29][30][31][32][33][34][35] for related work, particularly on time evolution of complexity). However, most applications so far have computed geodesics within a subspace of states or circuits, instead of dealing with the entire unitary group manifold. For instance, much recent work has focused on the subspace of Gaussian states, which are relevant in the context of free quantum field theories. This is because in continuum quantum-mechanical systems, the Hilbert space is often infinite-dimensional and it is difficult to define a tractable algebra of operators which generate the entire unitary group on the Hilbert space. Prior work Figure 1: Schematic of the unitary manifold (gray disk). A geodesic path (black) is depicted from the identity to some target unitary U . The red straight lines represent construction of a circuit using some elementary gates g i , and the final unitary is U = g 3 g 2 g 1 . The geodesic approximates the circuit smoothly by varying a control velocity V (s), analogous to an infinitesimal elementary gate, where s parametrizes the curve.
which attempted to deal with the global structure of the unitary group relied on toy models [4,36] of Lie group geometry. These models were constructed using metrics of strictly negative sectional curvature (or a discretization thereof, in the case of [36]) in order to ensure chaotic behavior of geodesics on the unitary manifold [4]. Here, we approach the problem of circuit complexity by studying aspects of geodesics on the complete group manifold SU (2 N/2 ), which is the unitary group acting on the Hilbert space of N/2 qubits. Our primary motivation is to study complexity growth in chaotic quantum systems as opposed to free field theories. To this end, we will use the (generalized) Sachdev-Ye-Kitaev (SYK) model as a specific example of a chaotic Hamiltonian, although most of our arguments are general and should apply to any chaotic system. where the couplings J i 1 ...iq are drawn at random from a Gaussian distribution with mean zero and variance σ 2 where J is a parameter setting the variance [37]. This model is expected to be chaotic and holographically dual to 2D quantum gravity [38][39][40][41] (see also [42] for a review and additional references). From the SYK perspective, the group SU (2 N/2 ) is the group of unitary operations (modulo an overall phase) acting on the Hilbert space of the N Majorana fermions (with N even) ψ i . Our main tool in studying the complexity in this model will be the Euler-Arnold equation [19,43,44], which was also used in a simpler setting in [45].
t ⇠ e N t ⇠ e e N t C C max Figure 2: The complexity in chaotic systems is conjectured [4] to grow linearly in time until a time of order e N , after which it saturates to (and fluctuates around) its maximum value of C max . At doubly exponential time, the complexity is expected to exhibit recurrences. .
From physical considerations and holographic as well as complexity-theoretic arguments, the complexity in chaotic systems has been conjectured [4] to grow linearly in time until a time of order e N , after which it is expected to saturate to (and fluctuate around) its maximum value of C max ∼ poly(N )e N (see Fig. 2), where by poly(N ) we mean N α for some α ≥ 0. Here N is the number of fermions in the SYK model, but more generally it should be taken to be log of the dimension of the Hilbert space. The motivation of the present work is to better understand the origin of this behavior and the various time scales involved from a field theory perspective, within the geodesic complexity framework. One of our main results will be to establish the existence and local minimality of a geodesic between the identity and e −iHt whose length grows linearly with time t. The existence of such a geodesic only relies on general features such as the Hamiltonian being local (i.e., it should be built from easy generators), and uniformity of the cost factor in the easy directions. However, this is not the whole story -the linear geodesic is not guaranteed to be a local minimum of the distance function (i.e., it could be a saddle point), much less a global minimum. As such, it may not be the relevant geodesic for complexity. In this paper, we will investigate in depth the question of local minimality of the linear geodesic by studying conjugate points along it. Roughly, we say that we have a conjugate point at time t if we can deviate infinitesimally from the linear geodesic at time t = 0 (i.e., deform the initial velocity infinitesimally) and return to it at time t along an infinitesimally nearby curve which satisfies the geodesic equation linearized to first order. The original geodesic stops being minimizing past the first conjugate point (i.e., it is a saddle point thereafter), and so for the physical considerations explained in Fig. 2 to be correct, it is necessary (but not sufficient) that the first conjugate point along the linear geodesic appears at times exponential in N . We will show that this is indeed the case for "sufficiently chaotic" Hamiltonians (such as the SYK model) and for an appropriate choice of the cost factors. Therefore, the linear geodesic is at least locally minimizing for times exponential in N , consistent with the expectations in Fig. 2. Our proof will involve a new criterion on the Hamiltonian from the vantage point of circuit complexity which we will call the eigenstate complexity hypothesis (ECH): Eigenstate Complexity Hypothesis (ECH): Let H be the Hamiltonian with energy eigenstates |m , |n etc., T α be the local generators in the Lie algebra, and Tα be the nonlocal generators. Define We will say that the Hamiltonian and the gate set satisfy the eigenstate complexity hypothesis, if for E m = E n in the large-N limit, where S is the log dimension of the Hilbert space (i.e., N 2 ln 2 for the SYK model) and r mn are O(1) numbers which do not scale with S.
In words, ECH is the condition that off-diagonal eigenstate projectors of the form |m n| which map one energy eigenstate of the Hamiltonian to a different eigenstate should have e −S suppressed overlaps with the easy/local/simple directions in the gate set, or equivalently, such off-diagonal energy eigenstate projectors must necessarily be "complex" (i.e., complicated). 2 For Hamiltonians which satisfy the ECH, the conjugate point analysis simplifies greatly, and the exponential bound on conjugate points can be analytically established. We will provide numerical evidence to show that the SYK model indeed satisfies the ECH.
The rest of this paper is organized as follows: in section 2, we will begin by briefly reviewing the geodesic complexity framework and setting up the Euler-Arnold formalism for studying the complexity of local Hamiltonians in the Lie algebra su(2 N/2 ) for even N . In sections 2.1 and 2.2, we study the simple case of N = 2 where all the geodesics between identity and e −iHt can be worked out and the complexity calculated using analytic and numerical techniques. In section 2.3, we will switch to general N and show the existence of a geodesic whose length grows linearly with time. In section 3, we will explore the local minimality of the linear geodesic by studying conjugate points. We will end with some remarks on late-time saturation, complexity classes, and quantum chaos in the Discussion (section 4).
2 Geometry of SU (2 N/2 ) The Hilbert space of N/2 qubits has a natural tensor factorization (2.1) We wish to study the geometry of the set of (unit-determinant) unitary operators U(H) that acts on this Hilbert space. In this case, this set is In order to study the differential geometry of SU (2 N/2 ) from the quantum computation viewpoint, we must pick a basis for the Lie algebra with some notion of locality, i.e., we should be able to identify some generators in the Lie algebra as local or "simple", and the rest as "complex". In quantum computation, we usually choose some simple unitary operators as the elementary gates to be used in building circuits. On the other hand, in the geodesic framework, it is natural to choose a k-local subspace of the Lie algebra of the unitary group manifold to correspond to "simple directions". We may think of the elementary gates of the quantum computation viewpoint as being exponentials of these simple generators. For general H, there is no guarantee that we can choose a basis for the unitary Lie algebra which respects any sort of locality. Luckily, for the qubit case SU (2 N/2 ), there are a couple of natural ways to proceed. We could pick the "Pauli basis", namely products of Pauli matrices acting on individual qubits, as our basis of generators. However, there is a second choice which is more natural from the point of view of the SYK model: consider the gamma matrices γ a with a ∈ {0, . . . , N − 1} which satisfy the Clifford algebra (with γ † a = γ a ): Now consider distinct ordered products T a 1 ···am = γ a 1 . . . γ am with m ∈ {1, . . . , N } and a p < a q for p < q. We will often denote these operators as simply T i , where i stands for the multiindex a 1 · · · a m . The total number of such ordered products is N m=1 N m = 2 N − 1. This is precisely the dimension of the Lie algebra su(2 N/2 ). It is simple to make such ordered products of gamma matrices Hermitian by inserting appropriate factors of i. We claim this construction is a basis for su(2 N/2 ), and we leave the proof to Appendix A. We can endow the gamma matrix basis with a natural notion of locality as follows: k-local generators of the Lie algebra are simply those involving k or fewer gamma matrices. This is precisely the natural notion of locality in the SYK model -from this point of view, the gamma matrices above correspond to the Majorana fermion operators ψ a in the SYK model.
The basic idea in the geodesic framework is to model circuit complexity [19] in terms of a minimal-length geodesic on SU (2 N/2 ) with respect to a right-invariant metric chosen such that it disincentivizes motion in the directions of nonlocal unitary operators. This corresponds to a choice of gate set in the quantum computation picture, where we allow up to k-local gates (i.e., exponentials of k-local generators in the Lie algebra) in our circuit but do not allow more nonlocal gates. In our context, we want to disincentivize motion in directions which correspond to generators involving products of more than k gamma matrices. Let us begin by constructing such a right-invariant metric. We can use the gamma matrix basis for su(2 N/2 ) to compute the structure constants f ij of the Lie algebra, defined as 3 where recall that the T i = γ a 1 · · · γ am are generators built from products of gamma matrices (or equivalently, products of the SYK fermion operators) labelled by the multi-index i = (a 1 · · · a m ). Using these, we calculate the Cartan-Killing form where h ∨ is the dual Coxeter number) which is a positive-definite 4 bilinear form. In order to build in the notion of simple and hard directions in the Lie algebra, we construct a new 3 When sums are not written explicitly, the Einstein summation convention is adopted. We caution the reader that there will be expressions in which repeated indices appear three times, but this will not cause any ambiguities because the three matching indices will always be summed over together. 4 Some definitions of the Cartan-Killing form instead yield a negative-definite form for compact Lie algebras.
We are only interested in this form up to overall sign and normalization since our only use for it is to define a right-invariant Riemannian metric on SU (2 N/2 ).
positive-definite bilinear form on su(2 N/2 ) where the numbers c i are "cost factors". Then a right-invariant metric g can be defined at an arbitrary point U on SU (2 N/2 ) by simply taking where we have used the group structure to transport the tangent vectors X and Y from U back to the identity and then applied (2.6). The cost factors c i encode the information about our choice of local and nonlocal directions, i.e. our notion of k-locality. We will generally take c i = 1 if the generator T i consists of k or fewer gamma matrices, and c i 1 otherwise. Note that if we chose cost factors c i = 1 for all i, the metric (2.7) would actually be bi-invariant.
Here bi-invariant means that the metric is both left and right invariant. The restriction to right-invariance arises by choosing at least one cost factor to be c i = 1 (or more generally by choosing a symmetric bilinear form for the metric which is not proportional to the identity).
Having chosen our cost factors, the geodesic equation on SU (2 N/2 ) with metric (2.7) is given in terms of the Lie algebra metric (2.6) and structure constants (2.4) by the Euler-Arnold equation [43] where the velocities V i (s) control the unitary path the geodesic follows via and we have made use of the path-ordered exponential to solve the matrix equation for the unitary operator dU ds = −iV i (s)T i U (s). (2.10) Finally, we impose the boundary condition for some target unitary whose circuit complexity we wish to study. This complexity is given by the geodesic length where the minimization is over all geodesics from the identity to U target . Throughout this paper, we will be interested in U target = e −iHt , where H is a suitable k-local Hamiltonian.

Analytics for N = 2
We will mainly be interested in studying the complexity for a large-N chaotic Hamiltonian, with the SYK model as a specific example. However, as a warm up, we will begin with the case of an SYK-like model with N = 2 fermions. The algebra for the N = 2 case is simply the familiar su(2) (see also [48]). There are three generators, built from the Hermitian matrices γ a . (2.13) We can compute the structure constants by using the algebra (2.3). (2.14) We see that, even though we have chosen a slightly unusual basis for the algebra, the structure constants are still f ij k = −2 ij δ k ; this is essentially the usual angular momentum algebra up to a minus sign and a factor of 2. This fact will allow us to solve the Euler-Arnold equation directly. The Cartan-Killing form is given by 6 Let us pick c 1 = c 2 = 1 and c 3 = 1 + µ, where µ is a large suppression factor to discourage motion in the T 3 direction. This corresponds to enforcing k = 1 locality. The equations (2.8) then reduce to 16) and this system can be solved to find the unique solution with integration constants v i .
(2.17) 6 We will always normalize the Cartan-Killing form to δij, regardless of the coefficient obtained by using (2.5).
Thus far we have solved the geodesic equation at the level of the Lie algebra, which allows us to obtain the tangent vector at any point along the geodesic given an initial direction. In fact, we can already compute the complexity of a path connecting U (0) = 1 to U (1) = U target : (2. 18) We see that the integrand is actually independent of s, leading to a simple result. All the information about the path length is contained in the magnitude of the tangent vector at the identity.
We really would like to know the geodesic for fixed boundary conditions U (0) = 1 and U (1) = U target in order to fix the initial tangent vector v i T i . The unitary U (s) along the geodesic path from the identity with tangent vector V i (s)T i is given by the path-ordered exponential (2.9). Now, we want to explicitly evaluate what the final unitary U (1) looks like as a function of the initial velocity v i , and then implement the boundary condition U (1) = e −iHt for some local, Hermitian Hamiltonian, in order to solve for v i . However, solving this would require us to brute-force deal with the path-ordering in (2.9), which is a famously difficult problem and is solved in quantum mechanics (where there is a time-ordering rather than a path-ordering) using perturbation theory. We would like a more nonperturbative approach, and we might hope that one exists since we are only dealing with finite-dimensional matrices rather than the infinite-dimensional Hilbert spaces familiar from other quantum systems like the harmonic oscillator. Indeed, such a nonperturbative method for finite-dimensional matrix equations was found in [49]. We employ their construction here. Given the velocity along the geodesic we wish to solve (2.10) subject to U (0) = 1 without the use of the path-ordering P. Let us define the frequency 20) and the function then the solution is Note that this is a completely coordinate-free description of a path on the unitary manifold SU (2); although SU (2) happens to have a convenient interpretation as S 3 , the higher groups SU (2 N/2 ) are nontrivial fiber bundles, so a coordinate patch-based method is likely difficult to implement.
We can now solve for v i (and hence compute the complexity of time evolution) by implementing the boundary condition U (1) = e −iHt for some Hamiltonian H, which we decompose as H = i J i T i . The time evolution operator can be exactly computed with a simple matrix exponential, yielding (letting (2.23) We can easily see that, if we had chosen all the metric cost factors to be c i = 1 (i.e. taken µ = 0 in (2.22)), the time evolution operator would itself define a geodesic curve. This is because, for bi-invariant metrics, the matrix exponential coincides with the (Riemannian) exponential map. In the next section, we will solve the boundary condition for the velocities v i in terms of the Hamiltonian couplings J i , for each value of t, using numerical techniques. There will be, in general, multiple solutions to any such equation which correspond to different geodesics in SU (2) which begin at the identity and end at e −iHt . We must of course find the one with minimal complexity.

Numerics for N = 2
The matrix equation (2.24) reduces to a system of three independent transcendental equations in the velocities v 1 , v 2 , v 3 . These can be solved numerically by brute force. Choosing the Hamiltonian to consist only of easy (1-local) operators, we work in the special case J 3 = 0. The numerical solution for the complexity (2.18) as a function of time t is presented in Fig. 3.
At early times, we find the expected linear growth of complexity, with slope J = J 2 1 + J 2 2 + J 2 3 . At later times, however, one does not see the plateau predicted from holographic considerations; rather, there is an immediate linear decay of the complexity that may be considered to be a Poincaré recurrence. This behavior may be attributed to the simplicity of the group manifold SU (2) S 3 and may be visualized most easily on the two-sphere (Fig. 4). After a certain maximum distance from the identity, the minimal length path on the sphere switches direction, i.e. the velocities v 1 , v 2 change in sign. This maximum distance is equal to π, the maximum of the complexity in Fig. 3, as follows from results in [50].
The complexity demonstrates an initial linear growth with the slope J = J 2 1 + J 2 2 + J 2 3 , attaining a maximum value of π, followed by linear oscillations (with slopes ±J).
In the higher groups SU (2 N/2 ), we expect a quite nontrivial topology that results in a plateau in the complexity as many different geodesics exchange dominance at late enough times. In the simple case of N = 2 fermions, we obviously do not see this effect, because the geodesic does not have sufficient space on the unitary group to wander around, away from its starting point. However, if we disorder average over the couplings J 1,2 in the Hamiltonian, then this is equivalent to considering an ensemble of systems where after the initial linear growth, we would expect to find "cancellations" between the various oscillating geodesic distances, thus leading to a plateau. In the SYK model, the couplings J are drawn from a Gaussian distribution of zero mean with a variance chosen to simplify the large-N limit, and the disorder average is performed during evaluation of correlation functions [38]. Here we will simply disorder average the complexity directly, taking the couplings in the Hamiltonian J 1 , J 2 to be drawn from Gaussian distributions of mean zero and some variance σ 2 .
The result of this averaging procedure is displayed in Fig. 5. Specifically, we are plotting  the disorder-averaged complexity: where f (t) refers to the triangle wave plotted in Fig. 3 and we have used the fact that J = J 2 1 + J 2 2 is drawn from a Rayleigh distribution when J 1 , J 2 are Gaussian-distributed. We see that the disorder average works beautifully in producing a complexity plateau, and that the complexity continues to grow linearly at early times, albeit with a modified slope π 2 σ. Heuristically, we see that even on as simple a manifold as S 3 , the disorder average causes a kind of interference between many different random samples of the couplings J 1 , J 2 . In terms of the triangle waves seen in Fig. 3, if one imagines a large number of copies of the system with different values of J, after the copy with the maximum value of J hits the first peak the various triangle waves begin to interfere destructively. At any fixed late time, the height of any one copy of the system is uniformly distributed between 0 and π and therefore the average complexity at long times is π/2. 7 It is straightforward to prove the late-time limit using the Fourier expansion of the triangle wave. Namely, the Fourier expansion of the triangle wave with slope J can be written The disorder average can be performed term-by-term following (2.25) assuming the integration can be exchanged with the infinite sum, and the result is x 0 e y 2 dy is Dawson's integral. If we allow ourselves to exchange the longtime limit with the infinite sum and use the identity lim x→∞ xF[x] = 1 2 , one finds from (2.27) that lim t→∞C (t) = π 2 . In the large-N SYK model, we expect that the geometry of SU (2 N/2 ) is sufficiently involved that there is a "self-averaging" effect on complexity, namely that the late-time complexity saturation occurs in a single realization of the SYK model (without disorder averaging). This is further discussed in section 4.1.
It is also of interest to find the time scale for the onset of the plateau. Previous discussions of complexity in holography [4] have noted that in large-N chaotic systems this time scale ought to be exponential in the size of the system / number of qubits. The N = 2 model studied in this section is neither large-N , nor chaotic; nevertheless it is useful to check that our results are compatible with an exponential time scale if N were increased. Since the disorder-averaged complexity scales at early times like σt, the plateau begins approximately upon reaching the asymptotic value at t ∼ Cmax σ . Assuming that C max ∼ poly(N )e N and using the fact that in the SYK model, 1/σ is typically taken to scale polynomially in N , one finds a result t ∼ poly(N )e N that is indeed consistent with the holographic expectation for the onset of the plateau. Figure 6: Disorder-averaged complexity (blue) choosing the magnitude J of the couplings to be Gaussian-distributed rather than Rayleigh-distributed. The peak of Fig. 5 is eliminated and the slope at early times is modified to 2 π σ (green), though the plateau value of π/2 (orange) remains the same.
Lastly, we remark that the existence of a peak in Fig. 5 appears to be a peculiarity of taking J 1 , J 2 to be Gaussian-distributed as our source of disorder; if one chose J to be Gaussian-distributed rather than Rayleigh-distributed, the peak vanishes and the average complexity smoothly approaches a plateau, albeit with slope 2 π σ (Fig. 6).

Linear geodesic for arbitrary N
For general N , the algebra su(2 N/2 ) is quite complicated. We have collected some facts, including a derivation of the structure constants, in Appendix A. However, the most important points for us are the following: firstly, the structure constants in the basis T i (corresponding to ordered products of gamma matrices) are fully antisymmetric by virtue of orthogonality in the trace norm. Secondly, f ij is nonzero if and only if Here we are thinking of the multi-indices i, j · · · as binary numbers; for instance (in the ordering convention of Appendix A), the operator T 3 = iγ 1 γ 2 corresponds to the binary number 00 · · · 0011, T 5 = iγ 1 γ 3 corresponds to the binary number 00 · · · 0101 etc. Further, q i is the number of ones (i.e., the number of fermions) in i, ⊕ stands for the bitwise XOR and ∧ stands for bitwise AND. The (suitably normalized) Cartan-Killing form follows after a short For convenience, we will label operators with k or fewer fermions with undotted Greek indices α, β... and those with more than k fermions with dotted Greek indicesα,β..., where k < N is arbitrary for now. We choose the easy directions, i.e., operators with k or fewer fermions, to have cost factors c α = c and the hard directions, i.e., operators with more than k fermions, to have cost factors cα =c. The Euler-Arnold equation can then be written schematically (since we have not determined overall sign) where we have explicitly written out the sums, and the index i on the left-hand side is not to be summed over. There is an interesting structure to (2.30) that emerges when we split into local and nonlocal directions. We first observe that f αα direction with indexα appears in a local direction β's velocity equation, it also appears in the velocity equation of the local direction α which multiplies it in β's equation. A similar story occurs for fα αβ = −fβ αα for the nonlocal directionsβ andα. The local direction α will occur in both of their velocity equations, appearing with opposite sign. We can introduce antisymmetric matricesṀ α β and Mαβ and rewrite (2.30) as whereṀ (V ) is a matrix with local indices which depends linearly on the nonlocal directions' velocities and M (V ) is a matrix with nonlocal indices that depends linearly on the local directions' velocities. Though this system is tricky to even write at arbitrary N , we can find a simple solution to it within the local subspace using the ansatz: which solves (2.31) becauseṀ = 0. The complexity is then where we have no contribution from the nonlocal directions. This is in accord with our intuitions about quantum circuit construction, where we do not just suppress nonlocal gates but completely disallow them. Since the velocities in (2.32) are constant, the path-ordering in (2.9) is trivial and the unitary path is If we take our target state to be U target = e iHt where H is a k-local Hamiltonian we can solve the boundary condition (2.24) to find one easy solution 9 v α = J α t. (2.36) We will refer to this geodesic as the linear geodesic. Assuming that the linear geodesic is the correct minimum, we find that the complexity (2.33) is The linear growth of the complexity in (2.37) matches expectations from holographic calculations of complexity as well as old observations about complexity growth in the geodesic formalism [4,18]. Our task now is to investigate the validity of the assumption that the linear geodesic is the correct minimum to consider.

Conjugate points and the eigenstate complexity hypothesis
One might wonder where the late-time behavior (i.e., late-time saturation) of complexity is going to appear from the previous discussion. The point, of course, is that the linear geodesic cannot be minimal for all times. After all, SU (2 N/2 ) is a compact manifold, and no geodesic path on a compact manifold can globally minimize the length between the identity and U = e −iHt for all t. In general, there are two ways a geodesic can become non-minimizing in a Riemannian manifold M : 9 One might think that the ambiguity in the logarithm gives multiple solutions here, but this is not the case, because generically the "other solutions" obtained from the log will not be entirely along the easy directions, and so are not admissible. These two conditions can roughly be thought of as local and global obstructions to minimality, respectively. This is because conjugate points along a geodesic segment mean that the segment is a saddle point, 11 not a minimum; the number of conjugate points along the segment is equal to the number of "downward directions". Therefore, conjugate points are an obstruction to a geodesic segment being locally minimizing. On the other hand, the absence of conjugate points but presence of geodesic loops indicates that the geodesic segment is locally minimizing but not globally minimizing. We will address the issue of conjugate points in this section. We will not prove the nonexistence of geodesic loops, but see the Discussion (section 4) for some further comments.
Prior studies of complexity using toy models have largely avoided the question of conjugate points (although see [19], where the importance of conjugate points in circuit complexity was emphasized previously) roughly by assuming all sectional curvatures are negative, so that geodesics originating at the same point generically diverge [4,36]. However, this assumption is worrisome, because it is well-known that any unimodular Lie group with left/right-invariant metric must possess some strictly positive sectional curvature, or else be completely flat [53]. If the sectional curvature cannot be everywhere bounded above by zero, one cannot rule out the existence of conjugate points on general grounds. Therefore, it is crucial to understand conjugate points on the full group manifold in the complexity metric (2.7). Here we will show a lower bound on the distance from the origin to the first conjugate point along the linear geodesic V = Ht.
In order to find conjugate points, we look for a velocity perturbation δV (s), also called a Jacobi field, 12 which obeys a first order differential equation known as the Jacobi equation, with particular boundary conditions which we will state precisely later. In 3.1, we solve the Jacobi equation for the velocity perturbation δV (s). In 3.2, we compute the first order change δU in the target unitary due to a velocity perturbation which obeys the Jacobi equation.
Setting this to zero gives a boundary condition for the Jacobi equation, which corresponds to having a conjugate point. In section 3.3, we will show that with the bi-invariant choice of metric (i.e., with the same cost factors for all generators in the Lie algebra), the linear geodesic has a large number of conjugate points. In particular the first conjugate point appears at t = 2π Emax−E min , where E max/min are the largest and smallest eigenvalues of the Hamiltonian respectively. In section 3.4, we will then return to the right-invariant case with a large cost factor for the hard directions in the set of generators. We will argue that at large-N and for Hamiltonians which satisfy what we will call the eigenstate complexity hypothesis (ECH), the linear geodesic segment from the identity to e −iHt does not have any conjugate points till t ∼ e N , thus showing that the linear geodesic is at least locally minimizing until times exponential in N .

Solving the Jacobi equation
In order to discover a conjugate point, we must deform the base geodesic with a velocity perturbation δV (s) which solves the geodesic equation to first order; this first order equation for δV (s) is called the Jacobi equation. The Jacobi equation in our context is obtained by studying the first order correction to the Euler-Arnold equation (around the original, unperturbed geodesic V = Ht) under a velocity perturbation Ht → Ht + δV (s) (see Fig. 7 for an illustration). We will confine our attention to the case with all the easy cost factors being c = 1 and all the hard cost factors beingc = 1 + µ, but it would be interesting to generalize our analysis to the case where the cost factors along the hard directions vary with the scale of non-locality. We can then write the Jacobi equation as where the subscripts L (local) and N L (nonlocal) denote projection into the local and nonlocal subspaces, i.e.
where T α , Tα are bases for the local and nonlocal subspaces respectively.
Let A L denote the vector space spanned by the local generators in the Lie algebra, and A N L denote the vector space spanned by the non-local generators. In order to solve the Jacobi I e iHt The local solution is

Conjugate points as zero modes
We wish to use (3.6) and (3.8) to determine whether there are conjugate points. This can be done by understanding the first order perturbation to the final unitary U (1) induced by δV . The exact final unitary and first order perturbation are, recalling (2.9), where the δU term is obtained by expanding the path-ordering in a Dyson series and taking the first term, (3.10) We now define a super-operator Y (µ) : δV (0) → U −1 δU (1) which acts by where we have inserted our solution for δV (s) into equation (3.10). A conjugate point, in this formalism, is given by the condition and therefore corresponds to a zero mode of the super-operator (3.11). So, our approach to finding conjugate points will be to study the spectrum of Y (µ) and check for when it develops zero modes. While this super-operator, as it appears in our analysis, is a linear operator on the Lie-algebra su(2 N/2 ), it is convenient to view Y (µ) as acting on the complexification of this vector space, i.e., on sl(2 N/2 , C), and study the spectrum in this complexified space. The reason for doing this is that our Lie algebra is a vector space over a non-algebraically-closed field R, and so the eigenvalues of Y (µ) need not be real, and the eigenvectors need not be real combinations of the elements of the Lie algebra. (This is true for essentially the same reason that solutions to x 2 + 1 = 0 only exist in C even though the equation involves coefficients only in R.) Of course, in order for a true conjugate point to appear for some values of t and µ, the zero mode must be Hermitian and traceless. In other words, it must indeed be a valid element of su(2 N/2 ).

Conjugate points in the bi-invariant case
Solving for the conjugate points at general values of µ is analytically hard. We will only be able to do it approximately in section 3.4 for large-N Hamiltonians which satisfy a certain complexity criterion on their eigenstates. But before doing that, it is useful to look at the much simpler case of µ = 0 where we can obtain all the conjugate points exactly. This is because at µ = 0, where all generators are considered computationally "easy", equation (3.11) simplifies greatly, and we get It is easy to guess the eigenvectors of this super-operator: they are simply the energy eigenstate projectors |m n|, where |m and |n are eigenstates of the Hamiltonian with energies E m and E n respectively. Indeed, we find where the eigenvalue φ mn (t) is given by 14 The eigenvalue φ mn (for E m = E n ) becomes zero at Indeed, at these times, the eigenvalues corresponding to both |m n| and |n m| become zero, and we can construct two Hermitian linear combinations out of these. Therefore, the linear geodesic develops a large number of conjugate points at the times given by equation (3.16), for all the possible choices of E m and E n . The first time t > 0 at which it develops a conjugate point is In the SYK model, the maximum separation is known to be ∆ max ∼ N , and so the linear geodesic stops being minimal after t c ∼ 2π N . However, this model is expected to be chaotic, so how is it that the conjugate points are appearing at a time of O(1/N )? The resolution of course lies in the fact that the bi-invariant metric is not the correct Riemannian metric for complexity. To understand physically relevant conjugate points we need to select a notion of locality for our generators. In other words, we need to choose which operators in the theory are "simple". By choosing the bi-invariant metric on the generators we have allowed arbitrary operators as local, but this is definitely not a physically sensible choice. However, the above calculation emphasizes the importance of conjugate points, and the need to make 14 Notice that diagonal projectors |n n| have constant eigenvalue 1 and therefore cannot lead to conjugate points in the bi-invariant case.
sure that they are absent if we are to establish the minimality of a geodesic. We now turn to the question of what happens to conjugate points for chaotic systems when a suitable notion of locality has been established by turning on cost factors in the complexity metric.

Turning on cost factors
We will turn on a finite cost factor µ, which will separate "easy" and "hard" computational directions, or, more physically, operations that we will consider "local" or "non-local". Our aim is to show that the linear geodesic is locally minimizing for times exponential in N , and so contains no conjugate points till such time. As stated previously, we do not have an exact solution for the spectrum of Y (µ) (although it is possible to calculate this spectrum perturbatively in µ, see Appendix B). However, if the Hamiltonian is sufficiently chaotic, then the situation simplifies greatly. More precisely, if the off-diagonal eigenstate projectors |m n| of the Hamiltonian are "complex", in the sense that their overlaps with the local generators are exponentially suppressed in N , then we can give an approximate formula for the spectrum of the super-operator Y (µ) at finite µ. We will call this criterion the eigenstate complexity hypothesis, or ECH for short: Eigenstate Complexity Hypothesis (ECH): Let H be the Hamiltonian with energy eigenstates |m , |n etc., T α be the local generators in the Lie algebra, and Tα be the nonlocal generators. Define R mn = α | m|T α |n | 2 α | m|T α |n | 2 + α | m|Tα|n | 2 .

(3.18)
We will say that the Hamiltonian and the gate set satisfy the eigenstate complexity hypothesis, if in the large-N limit for E m = E n , R mn = e −2S poly(S) r mn , (3.19) where S is ln dim of the Hilbert space (i.e., S = N 2 ln 2 for the SYK model), poly(S) is some polynomial in S, and r mn are O(1) numbers which do not scale with N . We can equivalently state this as || |m n| L || = O(e −S poly(S)), (3.20) where recall that the subscript L indicates projection to the local/easy subspace in the Lie algebra and the operator norm is defined by ||X|| = [Tr(X † X)] 1/2 .
The physical intuition behind this criterion is that off-diagonal projectors of the form |m n| map the energy eigenstate |n to a different eigenstate |m . For chaotic Hamiltonians, this operation should be complex from the point of view of local generators in the Lie algebra, since we expect these energy eigenstates to differ in their fine-grained microstructure. Another reason to expect ECH is that for sufficiently chaotic Hamiltonians, off-diagonal projectors like |m n| tend to have a uniformly distributed overlap with the generators in the Lie algebra (see Fig. 11 in Appendix C), and since there are exponentially many non-local generators and only polynomially many local generators (assuming k does not scale with N ), the projection of |m n| onto the local directions should be exponentially suppressed in N , as per equation (3.19). The interacting SYK model satisfies the ECH. To demonstrate this, we have shown an array-plot of the matrix R mn for a single realization in the left panel of Fig. 8, for the SYK model at N = 14, J = 1 and k = 3, q = 3. We see that the off-diagonal elements of R mn are indeed exponentially suppressed. Taking where the N -dependent coefficient is precisely the number of easy generators divided by the total number of generators, we have shown the distribution P (r mn ) of all the r mn s (including diagonals) over 100 realizations of the SYK model in the right panel of Fig. 8 (with N = 12 for convenience). The r mn s are distributed with a (sample) mean ofr s = 1, and (sample) standard deviation of σ s = 0.14. 15 We have also checked other values of N and (k, q) (with 15 It is easy to show from the definition of rmn that their mean is one: 1 2 N m,n rmn = 1. The distribution P (rmn) can be roughly approximated by the normal distribution with meanr 1 and standard deviation σ 0.098. A slightly better approximation is provided by Student's t-distribution with the parameters r = 0.994, σ = 0.093 and the number of degrees of freedom ν = 6. q < k) and found similar behavior. One novelty for q even (see Appendix C for further discussion) is that the Hamiltonian has a fermion number symmetry (which additionally is diagonal in the basis involving products of fermions), and this leads to an O(1) splitting of the distribution P (r mn ) into two distributions, corresponding to the off-diagonal projectors which either preserve or reverse the fermion number symmetry.
So far, we have presented some numerical evidence to show that the SYK model satisfies the ECH. More generally, we expect chaotic Hamiltonians to satisfy ECH (provided an appropriate choice is made for the k-local generators) as a consequence of a form of the eigenstate thermalization hypothesis (ETH), which is believed to be true in general chaotic quantum systems [46,47,54] (see [55][56][57] for discussion of ETH in the SYK model). ECH is of course very reminiscent of the ETH. In fact, we can see how the two are related in the SYK model. If we take the generators to be T i ∼ ψ a 1 · · · ψ am , the denominator in the definition (equation (3.18)) of R mn is equal to e S = 2 N/2 ; this just follows from the fact that |m n| has operator norm one, while each of the generators T i has norm e S/2 . Now, if we further assume that the local/easy generators satisfy ETH, then each term in the numerator of R mn is also O(e −S ). Since there are at most polynomially many local/easy generators (assuming k does not scale with N ), we deduce that R mn = O(poly(S)e −2S ), provided the easy generators satisfy ETH. From this perspective, we may view ECH as saying that the easy generators in our choice of the gate set should satisfy ETH, but where our easy generators are k-local, and so involve multi-site operators (not simply 1-local operators). On the other hand, ECH is a logically independent criterion from ETH; it requires that the off-diagonal outer products |m n| have small projection onto the easy/local directions, i.e., that they are complex, or alternatively that they are uniformly distributed in terms of their overlaps with all the e 2S generators of the gate-set. It would be interesting to investigate whether large-N integrable Hamiltonians violate ECH. Certainly, off-diagonal operators of the form |m n| in integrable systems tend to have overlaps with a far smaller subset of the e 2S generators in the gate set (see Fig. 12 in Appendix C.) Since the norm of |m n| is one, this naturally requires the individual overlaps n|T i |m to be larger. However, it may be that a weaker form of ECH is still satisfied in such systems, in which case, a weaker form of our results below may be applicable.
Let us now return to the problem of conjugate points at finite cost factor. If we take the statement (3.20) of ECH as given, then we have where in the last line we have used ECH together with the fact that the Hamiltonian is a linear combination of only polynomially many generators, and so the norm of [H, |m n| L ] N L can at most get a polynomial enhancement over the exponentially suppressed norm of |m n| L . This implies that if we take our initial velocity to be δV (0) = |m n|, then the solution (3.6), (3.8) to the Jacobi equation simplifies substantially Below, we will carefully justify that the corrections to equation (3.23), denoted as · · · above, are exponentially suppressed in N , but for now we will proceed with the main argument. With equation (3.23) in hand, we can evaluate the action of the super-operator Y (µ) on |m n|: where, once again, the correction terms are exponentially suppressed in N , as will be justified below. Performing the s integration, we find up to O(e −S poly(S)) corrections. 16 Therefore, as we crank up the cost factor µ, the conjugate points move off to larger and larger times (see Fig. 9). If we take the cost factor to be where is some small positive number (we really only need µt e S , so at polynomial time we can even take to be infinitesimally smaller than one), then there will be no conjugate points for any time polynomial in S. 17 Indeed, the first conjugate point is now located at (3.27) 16 Each conjugate point is two-fold degenerate at leading order in N , and the exponential corrections may split this two-fold degeneracy. 17 Conversely, if µ is polynomial in S, then there will necessarily be conjugate points at polynomial time. Here we have assumed that ∆ max does not scale exponentially with N . Indeed, for the SYK model ∆ max = O(N ). This shows that the linear geodesic segment from the identity to e −iHt is locally minimizing for times exponential in N . To be precise, we have shown that all the low-lying conjugate points which were present in the bi-invariant case have moved to exponential time upon turning on the cost factor µ = e S . In the bi-invariant case, all the diagonal projectors |m m| are eigenvectors of the super-operator Y (0) with unit eigenvalue and do not correspond to conjugate points. We can argue from continuity that this is still the case when we turn on the cost factor µ: since conjugate points are zero modes of Y (µ) , they cannot simply appear out of nowhere; as we can see in figure 9, they can only move smoothly along the time axis. Therefore, no new conjugate points should appear at finite time with a finite cost factor. This argument can be formalized using Morse theory [52].
We also note here that if there is an off-diagonal projector |m n| which violates ECH "maximally", namely that is has an almost unit overlap with the easy/local directions and a small overlap with the hard/non-local directions, then one similarly show that such a projector corresponds to an approximate eigenvector of Y (µ) with the eigenvalue φ mn (t). Therefore, this situation will lead to conjugate points at the O(1) times t = 2π ∆mn Z, provided ∆ mn is not exponentially small, and thus the linear geodesic would stop being minimizing early on in time evolution. We expect this behavior to be present at small N .
Bounding the correction terms: Now we wish to carefully justify that all the correction terms which were ignored above are indeed exponentially suppressed. To this end, let δV mn (s) be the Jacobi field along the linear geodesic with the initial condition δV mn (0) = |m n|, and define δV mn (s) = c(s)|m n| + δW (s), (3.28) where , and δW (s) is the correction to the leading order result in equation (3.23). We insert this into the Jacobi equation to obtain the differential equations satisfied by δW : δW (0) = 0, (3.29) where the source term S above is given by As long as µ (and t as well) does not scale exponentially with S, or scales mildly such as µ ∼ e S , the source terms have an exponentially suppressed norm by ECH: where S = N 2 ln(2), and we are using the usual operator norm ||X|| 2 = T r(X † X). Actually, at least for polynomial time, we can even take µ ∼ e (1− )S , and the source terms will still be suppressed; beyond this our arguments here will break down. Expressing δW in terms of the basis (T α , Tα) introduced previously and solving the second equation in (3.29), we obtain for the non-local piece of δW : remain O(1) and all other contributions are suppressed. Therefore, no zero modes of the super-operator can develop before at least one φ mn has dropped away from 1, and this does not occur until times parametrically later than t c ∼ e N for 0 < 1.

Discussion
In this paper, we study the quantum circuit complexity of unitary time evolution in qubit systems. Here, complexity measures the minimum amount of "simple" (or k-local) operations needed to build the time evolution operator U (t) = e −iHt . Our main tool is a geometrization of complexity in terms of geodesics on the unitary group manifold [19], which we study using the Euler-Arnold equation [45]. Using this approach, we directly relate complexity growth in a physical theory to its spectral properties, and thus to phenomena like chaos and integrability. We propose the Eigenstate Complexity Hypothesis as a criterion on the energy eigenstates of the theory as a condition for linear complexity growth for exponential times, modulo global obstructions, that would be expected in chaotic dynamics. We apply these ideas to the SYK model. First, for N = 2 fermions where the theory is integrable, we solve exactly and show that complexity grows linearly at initial times but then oscillates. For large-N , where the SYK theory is chaotic, we show numerically that ECH is satisfied, thus giving evidence that complexity grows linearly for exponential time, as predicted by the duality of SYK theory with the physics of black holes [37,38].
Various features of the complexity plot in Fig. 2 can be understood as arising from distinct traits of the underlying quantum system. For example, the appearance of a plateau has nothing to do with a notion of complexity, but rather comes from competition between various geodesics on the unitary group manifold and a self-averaging effect at large-N , which both occur even in the bi-invariant geometry where all operators are considered simple. 18 On the other hand, the location in time of the start of the plateau depends strongly on what we select as "simple" (local) vs. "complex" (nonlocal) operations; for example, if all operators are considered "simple" (corresponding to a bi-invariant metric on the unitary group) the complexity plateau starts at a polynomial time in N , rather than at exponential time when only k-local operators are considered simple. Similar statements apply for the complexity ramp and the length of the ramp, respectively. Additionally, large-N features like the ramp and plateau can be discovered at small N (even N = 2) by utilizing disorderaveraging which appears in, e.g., the SYK model. However, doubly exponential features like the Poincaré recurrences of complexity will be washed out by disorder and so are only present for a single instantiation of the model at large-N . The upshot of all this is that the disorderaveraging commonly employed in studies of the SYK model acts as a sort of crutch which replicates large-N features. These features should properly be interpreted as the effects of self-averaging, which occur even in a single model instance as long as the Hamiltonian is chaotic. Furthermore, the qualitative features of the complexity plot are present without any notion of easy/hard (local/nonlocal) operations, but the particular time scales which appear hinge crucially on the introduction of such a notion (defined, say, through cost factors in the complexity metric on the unitary group).

Late-time saturation
From physical considerations, it is expected that for quantum systems with gravity duals, the complexity will grow linearly until some time exponential in N , after which it will saturate [59]. Conceptually, in any theory with a UV and IR cutoff this saturation will occur because of the finite dimension of the group of unitary operators acting on the Hilbert space. This saturation of complexity is expected to arise in the geometric framework when the linear geodesic on the unitary manifold from the reference operator to the time evolution operator stops being globally minimizing. At this point other geodesics take over. Above, we demonstrated criteria for local minimality of the linear geodesic, i.e. under what conditions we can exclude other geodesics that are deformations of the linear one. However, since the unitary group is compact, globally there may be geodesic loops, and we have not studied these. In the simple case with a bi-invariant metric on the unitary group, i.e. setting the cost µ = 0 to take all operators to be "simple", we can get a glimpse of the relevant physics. In this case, the geodesic equation and boundary conditions are The solutions are given by where |n are the energy eigenstates of the Hamiltonian, and k = (k 1 , · · · , k 2 N/2 ) are integers which sum up to zero because of the traceless condition. Therefore, the complexity is given by where the minimization is over all integer vectors k subject to the constraint that they sum up to zero. In Fig. 10 we show a numerical plot 19 of C(t) for ten different realizations of the SYK model with N = 8, q = 3 and J = 1. In all the cases, complexity grows linearly with time for a while, but then saturates at an O(1) time, t ∼ 2π. The main feature to note here is the saturation for each individual SYK realization, arising from different geodesics (i.e., different integer vectors k) dominating the complexity at late times. In section 2.2, we found similar saturation behavior in the N = 2 case after disorder averaging; at larger N , each individual realization seems to "self-average" to produce a plateau, as is evident from the sum in equation (4.3).
We can think of the minimization problem in equation (4.3) as being roughly equivalent to a particle moving uniformly on a 2 N/2 -dimensional torus T 2 N/2 , starting from some initial point with the velocity (E 1 , · · · , E 2 N/2 ). The complexity is then simply the distance of the particle from its starting point. If the energy eigenvalues are suitably commensurate, then the distance from the starting point will grow linearly with time for some time, but then the particle will return to its origin, and this will result in an oscillating complexity. However, if the energy eigenvalues are incommensurate, then the particle will move away linearly, but will not come back to its origin in a short amount of time. Indeed, it will typically wander around the high-dimensional torus at a fixed average distance from the origin, thus leading to a saturation in the complexity.
It would be interesting to extend this analysis to the physically interesting situation where only k-local operators are taken to be "simple" (µ = 0), with cost 1 + µ = 1 + e S for all other operators. A possible general strategy to make progress is to analyze the appearance of the complexity plateau for theories that satisfy our Eigenstate Complexity Hypothesis.

Quantum computation
As discussed above, to show that the linear solution is the global minimizer for an exponential time we need to globally exclude other geodesics. In the bi-invariant geometry (µ = 0, all operators regarded as "simple"), all geodesics which reach the unitary U = e −iHt from the identity have initial velocity vectors equal to log U . The ambiguity in taking this logarithm gave a family of geodesics indexed by k ∈ Z 2 N/2 with n k n = 0, as explained in the previous section. Now consider some of the operators as "non-local" by turning on a cost factor in the metric for these directions in the unitary group. For large enough N we expect all geodesics that appeared in the bi-invariant analysis other than the linear one will have nonvanishing components along the non-local directions. Thus, when µ = 0 these trajectories should no longer be geodesics. Perturbatively, it is obvious that their length increases with µ, but a complete analysis requires a resummation of the perturbative expansion that accounts for the change in the geodesic trajectory as the metric is changed. The goal should be to demonstrate that, at large enough N , all of the non-trivial geodesic loops (if they still exist when 1 + µ ∼ 2 N ) have greater length than the linear solution for any t ∼ poly(N ), where poly(N ) is a polynomial of any degree.
A precise argument to this effect could be combined with our results to demonstrate a novel complexity class separation. 20 This is due to a theorem of Aaronson and Susskind [61], 20 Hamiltonian simulation has been studied in the context of complexity classes previously in quantum computation [60]. who showed that the complexity class separation PSPACE BQP/poly is true if and only if the time evolution operator e −iHt in general has complexity which grows linearly with t for a time greater than any polynomial in N . 21 Of course, we expect such growth only for chaotic Hamiltonians, and not in integrable systems. Furthermore, following Theorem 2 in [61], if the geodesic loop argument outlined above can be made for exponential times (or, more precisely, for times greater than any subexponential 22 ), then our our results (which show there are no conjugate points up to exponential time t c ∼ e N ) would actually imply the even stronger statement PSPACE BQSUBEXP/subexp. It would be interesting (and necessary for the aforementioned class separations to be established) to see if there is a relationship between our ECH criterion, which is central to the argument for complexity growth, and the complexitytheoretic assumption in [61] where the time evolution step e −iH was taken to implement one step of a reversible computationally-universal classical cellular automaton. 23 When conjugate points exist in our analysis they can be interpreted in terms of "fastforwarding" of the Hamiltonian, and of time evolution regarded as a quantum computation. Fast-forwarding of a Hamiltonian H occurs when time evolution with respect to H for a time t can be simulated on a quantum computer, using a different Hamiltonian, in a time much smaller than t [63]. General Hamiltonian simulation algorithms are well-studied in the quantum computation literature [63][64][65][66][67][68][69]. In particular [63] shows the existence of a family of Hamiltonians (based on Shor's algorithm) where an exponential fast-forwarding does happen. In our language, this means that there there is shorter path from the identity to the operator e −iHt than simply following the linear geodesic on the unitary manifold. The existence of a conjugate point does not signal a parametrically faster algorithm, as in the definition of [63], but the absence of a conjugate point is certainly necessary to rule out such speedups. 21 PSPACE is the class of problems that can be solved given polynomial space and BQP is the class of problems that can be solved in polynomial time by a family of quantum circuits that is constructed by a classical algorithm in polynomial time. Next, BQP/poly is the class of problems that can be solved in polynomial time by a family of quantum circuits, given a polynomial size string of "advice" for each problem size which can be used when constructing the corresponding circuit. The advice string can be different for different problem sizes. Here "BQP" stands for Bounded-error Quantum Polynomial time, where, because quantum computations are effectively probabilistic, we must require that errors occur with a probability less than some bound . Finally, BQSUBEXP is the class of quantum computations that can be done in a subexponential time, t ∈ O(e N α ) (for all α > 0), and BQSUBEXP/subexp is the same class but with subexponential size advice strings. 22 Note that some authors disagree on the definition of the subexponential class, effectively over whether it includes times like 2 N 1/3 (more generally, 2 o(n) ) or whether it only includes times strictly less than 2 N α for all α > 0. Our bound on conjugate points holds for a truly exponential time tc ∼ e N , so the discussion of BQSUBEXP/subexp in the main text holds for whichever definition was used by [61]. In the previous footnote we assumed the weaker SUBEXP = α>0 DTIME(2 N α ). 23 It is an empirical observation in complexity theory that universality in cellular automata is not difficult to achieve; on the contrary, it is quite difficult to avoid [62]. We consider it very likely that a generic chaotic Hamiltonian like the SYK model implements such an automaton via its time evolution steps.
Perhaps there is a connection between the existence of conjugate points (or maybe the failure of ECH) and violations of the computational time-energy uncertainty principle defined in [63] to detect speedups.

Quantum chaos
We have proposed the eigenstate complexity hypothesis (ECH) as a criterion for complexity growth for exponential time, a phenomenon that we should expect in chaotic theories, but not in integrable theories. ECH states, roughly, that the off-diagonal projectors of eigenstates should have exponentially small overlap with k-local operators. We have demonstrated that ECH is indeed satisfied in chaotic systems such as the SYK model. Physically, ECH is satisfied in these cases because a given projector |m n| has nonzero overlap with all the e 2S operators in the e S dimensional Hilbert space, which guarantees that the overlap with a given small set of k-local generators must be small by unitarity (see Fig. 11, Appendix C).
But what about integrable systems, such as the free Ising model H = −J i Z i Z i+1 ? All spin configurations are eigenstates of this model. A spin configuration can be turned into another by the action of local raising and lowering operators at some sites and the action of any number of Zs. Because of this there are ∼ e S generators which will have overlaps with a given eigenstate projector with a few spin flips such as |..., 1, ..., 0, ... ..., 0, ..., 1, ...|. Unitarity then suggests that the overlap of this projector with any given k-local operator (effectively the square root of (3.18)) will be ∼ e −S/2 . This does not satisfy the ECH criterion as we stated it, but suggests there should be more refined criteria separating theories that have, e.g., O(1), O(poly) and various weaker exponential overlaps between the eigenstate projectors and the k-local generators. More generally, we see in Fig. 13 (Appendix C) that when a system has conserved charges which act diagonally on the generators of the gate set, there are superselection sectors in the Hilbert space for the overlaps between the eigenstate projectors and the k-local operators. It is plausible that this phenomenon could be shown to generally lead to ECH violation in integrable models.
All of our results were developed in the context of finite dimensional systems. These could be understood as a discretization of continuum field theories with both an IR and a UV cutoff. It would be interesting to understand how to recover the continuum limit as the cutoffs are removed. and also obey γ † i = γ i (the Majorana condition). We can interpret these objects as 2 N/2 ×2 N/2 Hermitian matrices, and they are precisely the generalized gamma matrices of the Clifford algebra C N (R). The basis for su(2 N/2 ) is constructed by taking products of γ i with appropriate factors of i to ensure Hermiticity. Specifically, we consider ordered products γ i 1 . . . γ in with i 1 < · · · < i n . To be Hermitian, such a product needs a factor of i if n(n − 1)/2 is odd. We can now write the set of generators compactly using the set of binary strings of length N , b ∈ B N . The bits of the string are b = b N . . . b 1 , and let q b be the number of nonzero bits in b. We write and we then have We now show these generators are traceless. By construction, the gamma matrices individually are traceless, so we have tr γ i = 0. Additionally, an even number of them will be traceless since we have anticommutation and cyclicity of the trace: For an odd number, we use the gamma matrix construction in terms of tensor products of Pauli matrices. That is, given a set of N − 2 gamma matrices γ (N −2) a , we may create a set of 24 Note that we are labeling gamma matrices with i, j, etc., whereas in the main text we used a, b, etc. The reason for this is that here we reserve early alphabet letters for the binary form of base-10 integers which form an equally valid labeling of the generators that we employ in calculations.

N γ
There are four cases for tr γ i 1 . . . γ i 2k+1 . First, neither γ (N ) N appear in the product. In that case, we have γ , and since tr A ⊗ B = tr A tr B and tr σ 2k+1 3 = tr σ 3 = 0, we have the required result. The second case is when either γ (N ) appear in the product, but not both. Then the final tensor factor is either σ 1 or σ 2 , and similarly we have tr σ 1 = tr σ 2 = 0 so again the entire trace vanishes. The interesting case is when we have both γ i 2k−1 tr σ 1 σ 2 σ 3 . We now repeat the argument for tr γ i 2k−1 , since these smaller gamma matrices have a similar tensor product structure. Again the only interesting case is when both γ We now turn to linear independence. Notice first that, given generators T a and T b with a = b, we have T a T b = αT c for some α ∈ C and some c ∈ B N \ {0}. Now assume for the sake of contradiction that we have b α b T b = 0 for some constants α b ∈ R. Solve for a specific T a with nonzero coefficient and write α a T a = − b =a α b T b . Multiply both sides by whatever multiple of T a we need to get the identity on the left hand side. We now have However, if we now take the trace of both sides, the left hand side is tr 1 = 2 N/2 but the right hand side is c α c tr T c = 0, so they cannot be proportional. Thus, all the generators must be linearly independent.
Note that it is a basic fact of Lie algebras that the structure constants f ab c are fully antisymmetric since we have chosen a basis in which tr T a T b ∝ δ ab . This can be seen by noticing tr T a [T b , T c ] = tr T a f bc d T d ∝ f bc d δ ad = f bc a , and also by cyclicity of the trace we have therefore we have f bc a = f ab c , (A. 8) and this combined with the fact that f ab c = −f ba c implies full antisymmetry. Now recall q i is the number of nonzero bits in the binary expression of index a, i.e. q a is the number of fermions appearing in generator a. Let a ⊕ b be the bitwise "exclusive-or" of a and b. Let a ∧ b be the bitwise "and" of a and b. Then a lengthy calculation shows Furthermore, the magnitude of any nonzero structure constant is precisely |f ab c | = 2. The exact sign is more difficult to determine with simple calculations, but it will not be so important for our analysis. One may wonder how this description of the structure constants is consistent with the properties we claimed before. For instance, if a ⊕ b = c and q a q b + q a∧b is odd, do we have (by total antisymmetry) b ⊕ c = a and q b q c + q b∧c odd? Indeed we do, by properties of ⊕ and ∧: Another computation shows the parity of q a q b + q a∧b matches that of q b q c + q b∧c .
q a q b + q a∧b ≡ 1 mod 2 This serves as a consistency check on our condition (A.9). Notice that the condition (A.9) acts as a nontrivial selection rule for which velocities can appear in the Euler-Arnold equation (2.8). Let us see how this works. For a local direction a, in order for V b V c to appear on the right hand side of (2.8), we must have (without loss of generality) q b ≤ k and q c > k for some choice of k-locality of the gate set. Then, since q c = q a + q b − 2q a∧b , we must have k > 2q a∧b . For k = 2, this implies q a∧b = 0. Combining with q a q b + q a∧b ≡ 1 mod 2, we see q a q b ≡ 1 mod 2. This is quite a nontrivial condition; naïvely we may have imagined in a k = 2 model that there exist nonzero commutators between generators a and b that had q a = q b = 2 and q a∧b = 0, but this cannot be the case because then q a q b ≡ 0 mod 2. Similarly, we might have expected nonzero commutators for q a = 1 and q b = 2 yielding q c = 3 with q a∧b = 0, but this also does not occur by the selection rules (A.9). The conclusion of this analysis is that, in a k = 2 model, the local velocities evolve via dV α ds = 0.

B Conjugate points in perturbation theory
In the main text it was argued that the bi-invariant metric (with c =c = 1) has conjugate points at t mn = 2π ∆mn Z, where ∆ mn = (E m − E n ) are differences between energy eigenvalues of the Hamiltonian. In this appendix, we will perturbatively track the behavior of these conjugate points when we turn on an infinitesimal cost factorc = 1 + along the heavy directions. Recall that the equation for the Jacobi field (i.e., the Euler-Arnold equation Here the subscripts L and N L stand for projections of the corresponding operators along local and non-local directions respectively. We wish to check whether there exists some initial boundary condition δV (0) and some value of t, such that If so, then the corresponding value of t constitutes a conjugate point along the original linear geodesic U linear (s) = e −istH , at which point the linear geodesic stops being globally minimizing.
At first order in , we find the equations d ds δV Let us now return to the conjugate points t mn = 2π ∆mn Z we had obtained at zeroth order. From equation (B.9), it is clear that their locations have now moved at linear order in , which we can keep track of systematically in perturbation theory: t mn = t ,(0) mn + δ (1) t mn + 2 δ (2) t mn + · · · , (B.10) where t ,(0) mn = 2π ∆mn Z denotes the bi-invariant conjugate points. Substituting this expansion in equation (B.9) and demanding that the result vanish for m = m and n = n, 25 we deduce the shifts in the conjugate points: Note that the numerator on the right hand side involves the projection of δV (0) (0) to the non-local subspace, which makes the above expression somewhat non-trivial to evaluate. Importantly, however, δV (0) (0) at the zeroth order was only determined up to an arbitrary complex number z: δV (0) (0) = z|m n| +z|n m|. (B.12) Requiring that our expression for δ (1) t mn is real determines precisely two possible choices of z, which we may call z ± . (Actually, only the phase of z is determined. The absolute value can be set to one by choice of normalization.) Corresponding to these two fixed numbers z ± , we then have the two conjugate points at their respective locations t ,(0) mn + δ (1) t mn (z ± ), given by (B.11). Therefore, we find that at linear order in , the two-fold "degeneracy" in conjugate points splits. Nevertheless, they continue to exist and we have tracked their locations at O( ) above.
This analysis can be repeated order by order in perturbation theory to determine the location of conjugate points. In the main text, it was shown that for the cost factorc = 1+e S , the conjugate points move to t ∼ t (0) (1 + µ), assuming ECH. We see that the perturbative results derived here are consistent with this. In particular, we reproduce the formula in the main text if we drop the δV

C Some more details on ECH
Here we provide some more numerical evidence for ECH in the SYK model. First, let us consider writing a generic off-diagonal eigenstate projector |m n| in the SYK model in terms of the generators T a = (T α , Tα) consisting of products of fermions: |m n| = 1 2 N/2 a c a T a , c a = n|T a |m , (C.1) where a runs over all the directions, easy and hard. We can get some heuristic understanding of why ECH is true in the SYK model by looking at the distribution of the c a . We see from  the left panel of Fig. 11 that the c a are more or less uniformly distributed over all the e 2S generators. Since R mn is the weight in the easy directions, the uniformity in the distribution of c a implies that R mn will be proportional to the number of easy directions divided by the total number of directions, which is precisely what ECH requires. A related comment is that if we build the distribution of the c a s by pooling together these coefficients for all choices of m and n, then we find a distribution with an exponential tail (see the right panel of Fig.  11). Since the tail is exponential, and the R mn s correspond to normalized sums over the easy coefficients, we expect the distribution over R mn s to be Gaussian in the large-N limit, by the central limit theorem. This is consistent with the distribution in Fig. 8. In Fig. 12 we have shown the distribution of the c a for a typical projector of an integrable Hamiltonian H = ψ 1 ψ 2 . Note that in this case, the overlaps are distributed in a much smaller subset of the generators. Nevertheless, there seem to be overlaps with about e S generators (as opposed to e 2S in the SYK model), suggesting a milder suppression of R mn . We have mainly focused on k = 3, q = 3 in our presentation. But similar results also apply to k = 4, q = 4, with slight modifications. The main novelty is that for q = 4, the Hamiltonian has a fermion-number symmetry. As a consequence, the eigenstates of the Hamiltonian carry an extra quantum number, namely the fermion number which acts diagonally on the generators involving products of fermions. This means that the off-diagonal projectors |m n| are of two types: 1. "Bosonic" or fermion number preserving, 2. "Fermionic" or fermion number reversing. As a consequence of this, the distribution of r mn s in this case splits into two well-localized distributions, see Fig. 13. Since the fermionic projectors cannot have any overlap with the four-fermion operators in the easy part of the Lie algebra, the corresponding r mn s are slightly suppressed (their distribution has moved to the left). On the other hand, since the average is constrained to one, this forces the bosonic r mn s to be slightly enhanced (their distribution has moved to the right). However, these effects are polynomial in N , and do not affect the overall exponential suppression of all the R mn s.