Quantum Complexity of Time Evolution with Chaotic Hamiltonians

TL;DR

This paper investigates the quantum complexity growth in chaotic systems like the SYK model, demonstrating linear growth over exponential times and its relation to geodesics on the unitary manifold, with implications for quantum computational complexity.

Contribution

It introduces the Eigenstate Complexity Hypothesis (ECH) and shows that in large-N chaotic systems, complexity grows linearly for exponential times, supported by numerical evidence.

Findings

Large-N SYK model satisfies ECH and exhibits linear complexity growth.

Integrable N=2 fermion case shows initial linear growth followed by oscillations.

Complexity growth relates to spectral properties and energy eigenstates of physical theories.

Abstract

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.

Equations

$$H=\sum _{i_1<\cdots <i_q}{J}_{i_1\dots i_q}{\psi }_{i_1}\dots {\psi }_{i_q},$$

$${\sigma }^2=\frac{(q-1)!{\mathcal{J}}^2}{N^{q-1}},$$

$$R_{mn}=\frac{{\sum }_{\alpha }|⟨m|T_{\alpha }|n⟩{|}^2}{{\sum }_{\alpha }|⟨m|T_{\alpha }|n⟩{|}^2+{\sum }_{\dot{\alpha }}|⟨m|T_{\dot{\alpha }}|n⟩{|}^2}.$$

$$R_{mn}=e^{-2S}\text{poly}(S)r_{mn},$$

$$\mathcal{H}=\underset{N/2}{\underbrace{{\mathbb{C}}^2\otimes \cdots \otimes {\mathbb{C}}^2}}.$$

$$\mathcal{U}(\mathcal{H})=\mathcal{U}({\mathbb{C}}^{2^{N/2}})=SU(2^{N/2}).$$

$$\{{\gamma }_a,{\gamma }_b\}=2{\delta }_{ab}.$$

$$[T_i,T_j]=i{f_{ij}}^{\ell }{T}_{\ell },$$

$$K_{ij}=-\frac{1}{h^{\vee }}{f_{im}}^{\ell }{f_{j\ell }}^m,$$

$$G_{ij}=\frac{c_i+c_j}{2}{K}_{ij},$$

$$g_U(X,Y)=G(X{U}^{-1},Y{U}^{-1}),$$

$$G_{ij}\frac{d{V}^j}{ds}={f_{ij}}^p{V}^j{G}_{p\ell }{V}^{\ell },$$

$$U(s)=\mathcal{P}\exp \left(-i{\int }_0^{s}d{s}^{\prime }{V}^i(s^{\prime })T_i\right),$$

$$\frac{dU}{ds}=-i{V}^i(s)T_{i}U(s).$$

$$U(1)=U_{\text{target}}$$

$$\mathcal{C}(U_{\text{target}})=\text{min}{\int }_0^{1}ds\sqrt{G_{ij}{V}^i(s)V^j(s)}$$

$$||U_{\text{target}}-\tilde{U}(1)|{|}_{\text{op.}}\le \frac{e^{S/2}}{\sqrt{1+\mu }}\mathcal{C}(U_{\text{target}}).$$

$$\begin{array}{rl}{T}_1& ={\gamma }_1,\\ T_2& ={\gamma }_2,\\ T_3& \equiv T_{12}=i{\gamma }_1{\gamma }_2.\end{array}$$

$$\begin{array}{rl}[T_1,T_2]& =-2i{T}_3,\\ [T_2,T_3]& =-2i{T}_1,\\ [T_3,T_1]& =-2i{T}_2.\end{array}$$

$$K_{ij}={\delta }_{ij}.$$

$$\begin{array}{rl}\frac{d{V}^1}{ds}& =-2\mu V^2{V}^3,\\ \frac{d{V}^2}{ds}& =2\mu V^3{V}^1,\\ (1+\mu )\frac{d{V}^3}{ds}& =0,\end{array}$$

$$\begin{array}{rl}{V}^1(s)& =v^1\cos \left(v^3\mu s\right)-v^2\sin \left(v^3\mu s\right),\\ V^2(s)& =v^2\cos \left(v^3\mu s\right)+v^1\sin \left(v^3\mu s\right),\\ V^3(s)& =v^3/2.\end{array}$$

$$\begin{array}{rl}\mathcal{C}& ={\int }_0^{1}ds\sqrt{(V^1{)}^2+(V^2{)}^2+(1+\mu )(V^3{)}^2}\\ & =\sqrt{(v^1{)}^2+(v^2{)}^2+\frac{1+\mu }{4}(v^3{)}^2}.\end{array}$$

$$V(s)=V^{i}(s)T_{i}=\left(\begin{array}[]{cc}-v^{3}/2&(v^{1}-iv^{2})e^{-iv^{3}\mu s}\\ (v^{1}+iv^{2})e^{iv^{3}\mu s}&v^{3}/2\\ \end{array}\right),$$

$${\omega }^2=(v^1{)}^2+(v^2{)}^2+\frac{1}{4}(1+\mu {)}^2(v^3{)}^2,$$

$$\phi (s)=e^{-\frac{i}{2}{v}^3\mu s},$$

$$U(s)=\left(\begin{array}[]{cc}\varphi(s)\left(\cos\omega s+i\frac{v^{3}(1+\mu)}{2\omega}\sin\omega s\right)&-i\varphi(s)\frac{(v^{1}-iv^{2})}{\omega}\sin\omega s\\ -i\varphi(-s)\frac{(v^{1}+iv^{2})}{\omega}\sin\omega s&\varphi(-s)\left(\cos\omega s-i\frac{v^{3}(1+\mu)}{2\omega}\sin\omega s\right)\\ \end{array}\right).$$

$$e^{-iHt}=\left(\begin{array}[]{cc}\cos Jt+i\frac{J_{3}}{J}\sin Jt&-i\frac{J_{1}-iJ_{2}}{J}\sin Jt\\ -i\frac{J_{1}+iJ_{2}}{J}\sin Jt&\cos Jt-i\frac{J_{3}}{J}\sin Jt\end{array}\right).$$

$$U(1)=e^{-iHt}$$

$$\bar{C}(t)={\int }_0^{\infty }dJ{f}_{△}(t)\frac{J}{{\sigma }^2}{e}^{-\frac{J^2}{2{\sigma }^2}}$$

Full text

1 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\sim S$ (where $S$ is the entropy of the black hole) as a consequence of the unitarity of quantum mechanics Maldacena2003 . 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 Stanford:2014jda ; Brown2016 .111Here “circuit complexity” measures the minimum number of simple, perhaps locally acting, gates necessary to construct a desired state or operator from a fixed reference. If this conjecture is correct, then the extrapolation of the gravity result to times beyond $t\sim 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 Brown2017 ; Cotler2017 . However, in the case of circuit complexity, despite the plethora of analytic results from gravity calculations (see Carmi:2017jqz ; Fu:2018kcp ; Agon:2018zso ; Swingle2018a ; Chapman2018a ; Chapman2018b ; Brown:2018bms ; Caceres2019 ; Caceres:2018blh ; Goto2019 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 Nielsen2003 ; Nielsen2005 ; Nielsen2006 ; Nielsen2007 .222But see Czech:2017ryf ; Caputa:2017yrh for proposed path integral approaches, which have not been shown to be polynomially equivalent to quantum circuit complexity. 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_{\alpha}$, and “non-local” or “hard” directions $T_{\dot{\alpha}}$. 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 Nielsen2006 . 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 everywhere-smooth 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 Myers2017 ; Chapman:2017rqy ; Khan:2018rzm ; Myers2018 (see also Caputa:2017yrh ; Hashimoto:2017fga ; Kim:2017qrq ; Moosa:2017yvt ; Molina-Vilaplana:2018sfn ; Alves:2018qfv ; Magan:2018nmu ; Caputa:2018kdj ; Camargo:2018eof ; Ali:2018fcz ; Belin:2018bpg ; Jiang:2018nzg 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 which attempted to deal with the global structure of the unitary group relied on toy models Brown2017 ; Lin2018 of Lie group geometry. These models were constructed using metrics of strictly negative sectional curvature (or a discretization thereof, in the case of Lin2018 ) in order to ensure chaotic behavior of geodesics on the unitary manifold Brown2017 . 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.

Recall that the SYK model is a quantum-mechanical system comprising $N$ Majorana fermions $\psi_{i}$ with the Hamiltonian

[TABLE]

where the couplings $J_{i_{1}\dots i_{q}}$ are drawn at random from a Gaussian distribution with mean zero and variance $\sigma^{2}$

[TABLE]

where $\mathcal{J}$ is a parameter setting the variance Kitaev2015 . This model is expected to be chaotic and holographically dual to 2D quantum gravity Maldacena2016 ; Polchinski:2016xgd ; Kitaev:2017awl ; Stanford:2017thb (see also Sarosi:2017ykf 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) $\psi_{i}$. Our main tool in studying the complexity in this model will be the Euler-Arnold equation Arnold1966 ; Tao2010 ; Nielsen2007 , which was also used in a simpler setting in Balasubramanian2018 .

From physical considerations and holographic as well as complexity-theoretic arguments, the complexity in chaotic systems has been conjectured Brown2017 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 $\mathcal{C}_{\text{max}}\sim\text{poly}(N)e^{N}$ (see Fig. 2), where by poly$(N)$ we mean $N^{\alpha}$ for some $\alpha\geq 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 no conjugate points appear along the linear geodesic at times sub-exponential in $N$. We will give an argument 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\rangle,|n\rangle$ etc., $T_{\alpha}$ be the local generators in the Lie algebra, and $T_{\dot{\alpha}}$ be the non-local generators. Define

[TABLE]

We will say that the Hamiltonian and the gate set satisfy the eigenstate complexity hypothesis, if for $E_{m}\neq E_{n}$ in the large-$N$ limit,

[TABLE]

where $S$ is the log dimension of the Hilbert space (i.e., $\frac{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\rangle\langle 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).333We discuss the relationship with the well-known Eigenstate Thermalization Hypothesis (ETH) PhysRevA.43.2046 ; 1994PhRvE..50..888S in the main text. For Hamiltonians which satisfy the ECH, the conjugate point analysis simplifies greatly, and the exponential bound on conjugate points can be analytically argued. 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 $\mathfrak{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

[TABLE]

We wish to study the geometry of the set of (unit-determinant) unitary operators $\mathcal{U}(\mathcal{H})$ that acts on this Hilbert space. In this case, this set is

[TABLE]

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 $\mathcal{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 $\gamma_{a}$ with $a\in\{0,\dots,N-1\}$ which satisfy the Clifford algebra (with $\gamma_{a}^{\dagger}=\gamma_{a}$):

[TABLE]

Now consider distinct ordered products $T_{a_{1}\cdots a_{m}}=\gamma_{a_{1}}\dots\gamma_{a_{m}}$ with $m\in\{1,\dots,N\}$ and $a_{p}<a_{q}$ for $p<q$. We will often denote these operators as simply $T_{i}$, where $i$ stands for the multi-index $a_{1}\cdots a_{m}$. The total number of such ordered products is $\sum_{m=1}^{N}{{N}\choose{m}}=2^{N}-1$. This is precisely the dimension of the Lie algebra $\mathfrak{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 $\mathfrak{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 $\psi_{a}$ in the SYK model.

The basic idea in the geodesic framework is to model circuit complexity Nielsen2007 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 $\mathfrak{su}(2^{N/2})$ to compute the structure constants ${f_{ij}}^{\ell}$ of the Lie algebra, defined as444When 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.

[TABLE]

where recall that the $T_{i}=\gamma_{a_{1}}\cdots\gamma_{a_{m}}$ are generators built from products of gamma matrices (or equivalently, products of the SYK fermion operators) labelled by the multi-index $i=(a_{1}\cdots a_{m})$. Using these, we calculate the Cartan-Killing form

[TABLE]

(where $h^{\vee}$ is the dual Coxeter number) which is a positive-definite555Some 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})$. bilinear form. In order to build in the notion of simple and hard directions in the Lie algebra, we construct a new positive-definite bilinear form on $\mathfrak{su}(2^{N/2})$

[TABLE]

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

[TABLE]

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+\mu$ with $\mu\gg 1$ otherwise; we will specify how large $\mu$ has to be shortly. 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}\neq 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 Arnold1966 666The original article is in French, but an English summary can be found in Tao2010 .

[TABLE]

where the velocities $V^{i}(s)$ control the unitary path the geodesic follows via

[TABLE]

and we have made use of the path-ordered exponential to solve the matrix equation for the unitary operator

[TABLE]

Finally, we impose the boundary condition

[TABLE]

for some target unitary whose circuit complexity we wish to study. This complexity is given by the geodesic length

[TABLE]

where the minimization is over all geodesics from the identity to $U_{\text{target}}$. Throughout this paper, we will be interested in $U_{\text{target}}=e^{-iHt}$, where $H$ is a suitable $k$-local Hamiltonian.

Before moving on to calculations, we would like to specify how large the cost factor $\mu$ needs to be. The essential reason for choosing $\mu$ large is that it prevents the geodesic from wandering off in the hard, i.e., non-$k$-local directions, and we can then take such a geodesic to be a reasonable approximation to the true minimal circuit built only out of the allowed $k$-local gates (more precisely, gates of the form $g_{\alpha}=e^{i\epsilon T_{\alpha}}$ for $k$-local $T_{\alpha}$). Now imagine that we start with our minimal geodesic $U(s)$, and then define a new curve $\widetilde{U}(s)$ by simply projecting out from its velocity $V(s)$ all the hard directions. Having done so, this new curve would deviate from the target unitary by some amount, which we would like to be smaller than some fixed error tolerance. It was shown in Nielsen2006 that777In lemma 1 of Nielsen2006 a weaker inequality was proven with the coefficient of the right hand side being $e^{S}$ as opposed to $e^{S/2}$. However, one can do better by using the fact that the operator norm is upper bounded by the Frobenius norm in step 3 of the derivation found in the appendix of Nielsen2006 .

[TABLE]

If we pick $1+\mu=\frac{1}{\epsilon^{2}}e^{S}\mathcal{C}^{2}$, then the right hand side above can be made smaller than $\epsilon$. Therefore, for any target unitary with polynomial complexity888Our arguments should work more generally for any target unitary with sub-exponential complexity. For e.g., when $t$ scales sub-exponentially with $S$, then $e^{-itH}$ necessarily has at best sub-exponential complexity, as will become clear later., we need to take the cost factor to be $\mu\propto e^{S}$, where the proportionality factor may scale at most polynomially in $S$.

2.1 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 $\mathfrak{su}(2)$ (see also Brown:2019whu ). There are three generators, built from the Hermitian matrices $\gamma_{a}$.

[TABLE]

We can compute the structure constants by using the algebra (2.3).

[TABLE]

We see that, even though we have chosen a slightly unusual basis for the algebra, the structure constants are still ${f_{ij}}^{k}=-2\epsilon_{ij\ell}\delta^{\ell 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 by999We will always normalize the Cartan-Killing form to $\delta_{ij}$, regardless of the coefficient obtained by using (2.5).

[TABLE]

Let us pick $c_{1}=c_{2}=1$ and $c_{3}=1+\mu$, where $\mu$ 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

[TABLE]

and this system can be solved to find the unique solution with integration constants $v^{i}$.

[TABLE]

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_{\text{target}}$:

[TABLE]

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_{\text{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 Giscard2015 . We employ their construction here. Given the velocity along the geodesic

[TABLE]

we wish to solve (2.10) subject to $U(0)=1$ without the use of the path-ordering $\mathcal{P}$. Let us define the frequency

[TABLE]

and the function

[TABLE]

then the solution is

[TABLE]

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=\sum_{i}J_{i}T_{i}$. The time evolution operator can be exactly computed with a simple matrix exponential, yielding (letting $J=\sqrt{J_{1}^{2}+J_{2}^{3}+J_{3}^{2}}$)

[TABLE]

We can easily see that, if we had chosen all the metric cost factors to be $c_{i}=1$ (i.e. taken $\mu=0$ in (2.23)), 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

[TABLE]

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.

2.2 Numerics for $N=2$

The matrix equation (2.25) 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.19) 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=\sqrt{J_{1}^{2}+J_{2}^{2}+J_{3}^{2}}$. 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)\simeq 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 $\pi$, the maximum of the complexity in Fig. 3, as follows from results in su2dist .

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 Maldacena2016 . 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 $\sigma^{2}$.

The result of this averaging procedure is displayed in Fig. 5. Specifically, we are plotting the disorder-averaged complexity:

[TABLE]

where $f_{\triangle}(t)$ refers to the triangle wave plotted in Fig. 3 and we have used the fact that $J=\sqrt{J_{1}^{2}+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 $\sqrt{\frac{\pi}{2}}\sigma$. 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 [math] and $\pi$ and therefore the average complexity at long times is $\pi/2$.101010We thank Clélia de Mulatier for discussion on this point. 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

[TABLE]

The disorder average can be performed term-by-term following (2.26) assuming the integration can be exchanged with the infinite sum, and the result is

[TABLE]

where $\mathcal{F}[x]=e^{-x^{2}}\int_{0}^{x}e^{y^{2}}dy$ is Dawson’s integral. If we allow ourselves to exchange the long-time limit with the infinite sum and use the identity $\lim_{x\to\infty}x\mathcal{F}[x]=\frac{1}{2}$, one finds from (2.28) that $\lim_{t\to\infty}\bar{C}(t)=\frac{\pi}{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 Brown2017 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 $\sigma t$, the plateau begins approximately upon reaching the asymptotic value at $t\sim\frac{\mathcal{C}_{\text{max}}}{\sigma}$. Assuming that $\mathcal{C}_{\text{max}}\sim\text{poly}(N)e^{N}$ and using the fact that in the SYK model, $1/\sigma$ is typically taken to scale polynomially in $N$, one finds a result $t\sim\text{poly}(N)e^{N}$ that is indeed consistent with the holographic expectation for the onset of the plateau.

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 $\sqrt{\frac{2}{\pi}}\sigma$ (Fig. 6).

2.3 Linear geodesic for arbitrary $N$

For general $N$, the algebra $\mathfrak{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}}^{\ell}$ is nonzero if and only if

[TABLE]

Here we are thinking of the multi-indices $i,j\cdots$ as binary numbers; for instance (in the ordering convention of appendix A), the operator $T_{3}=i\gamma_{1}\gamma_{2}$ corresponds to the binary number $00\cdots 0011$, $T_{5}=i\gamma_{1}\gamma_{3}$ corresponds to the binary number $00\cdots 0101$ etc. Further, $q_{i}$ is the number of ones (i.e., the number of fermions) in $i$, $\oplus$ stands for the bitwise XOR and $\wedge$ stands for bitwise AND. The (suitably normalized) Cartan-Killing form follows after a short computation111111The appropriate normalization factor in the general case with our choice of generators is $2^{-(N+1)}$.

[TABLE]

For convenience, we will label operators with $k$ or fewer fermions with undotted Greek indices $\alpha,\beta...$ and those with more than $k$ fermions with dotted Greek indices $\dot{\alpha},\dot{\beta}...$, 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_{\alpha}=c$ and the hard directions, i.e., operators with more than $k$ fermions, to have cost factors $c_{\dot{\alpha}}=\bar{c}$. The Euler-Arnold equation can then be written schematically (since we have not determined overall sign)

[TABLE]

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.31) that emerges when we split into local and nonlocal directions. We first observe that ${f_{\alpha\dot{\alpha}}}^{\beta}=-{f_{\beta\dot{\alpha}}}^{\alpha}$. So, if a nonlocal direction with index $\dot{\alpha}$ appears in a local direction $\beta$’s velocity equation, it also appears in the velocity equation of the local direction $\alpha$ which multiplies it in $\beta$’s equation. A similar story occurs for ${f_{\dot{\alpha}\alpha}}^{\dot{\beta}}=-{f_{\dot{\beta}\alpha}}^{\dot{\alpha}}$ for the nonlocal directions $\dot{\beta}$ and $\dot{\alpha}$. The local direction $\alpha$ will occur in both of their velocity equations, appearing with opposite sign. We can introduce antisymmetric matrices ${\dot{M}^{\alpha}}_{\beta}$ and ${M^{\dot{\alpha}}}_{\dot{\beta}}$ and rewrite (2.31) as

[TABLE]

where $\dot{M}(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:

[TABLE]

which solves (2.32) because $\dot{M}=0$. The complexity is then

[TABLE]

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.33) are constant, the path-ordering in (2.9) is trivial and the unitary path is

[TABLE]

If we take our target state to be $U_{\text{target}}=e^{iHt}$ where $H$ is a $k$-local Hamiltonian

[TABLE]

we can solve the boundary condition (2.25) to find one easy solution121212One 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.

[TABLE]

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.34) is

[TABLE]

where in the second equality we have rewritten the coefficient in terms of the energy eigenvalues $E_{m}$ by relating the expressions inside the square roots to $\operatorname{tr}H^{2}$. The linear growth of the complexity in (2.38) matches expectations from holographic calculations of complexity as well as old observations about complexity growth in the geodesic formalism Nielsen2006 ; Brown2017 . Our task now is to investigate the validity of the assumption that the linear geodesic is the correct minimum to consider.

3 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$:

Conjugate points: given a geodesic $U(s):[0,1]\to M$, there exists a variation through curves $U(\eta,s):[-\delta,\delta]\times[0,1]\to M$ such that $U(\eta,s)$ obeys the geodesic equation at first order in $\eta$, $U(0,s)=U(s)$, $U(\eta,0)=1$ and $U(\eta,1)=U(1)+\mathcal{O}(\eta^{2})$.131313See Witten:2019qhl for a recent discussion of conjugate points in general relativity. 2. 2.

Geodesic loops: given a geodesic $U(s):[0,1]\to M$ there is another geodesic $\widetilde{U}(s):[0,1]\to M$ such that $U$ and $\widetilde{U}$ have the same length $L[U]=L[\widetilde{U}]$, $U(0)=\widetilde{U}(0)$, and $U(1)=\widetilde{U}(1)$.

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,141414We mean here a saddle point of the energy functional on the space of paths Morse1934 . 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 Nielsen2007 , 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 Brown2017 ; Lin2018 . However, this assumption is worrisome, because it is well-known that any unimodular Lie group with left- or right-invariant metric must possess some strictly positive sectional curvature, or else be completely flat Milnor1976 . 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$. We will call the time at which the linear geodesic develops this first conjugate point $t_{c}$.

In order to find conjugate points, we look for a velocity perturbation $\delta V(s)$, also called a Jacobi field,151515More precisely, the Jacobi field is the first order deformation of the original geodesic, and $\delta V(s)$ is its derivative pulled back to the identity. We will sometimes loosely refer to $\delta V$ itself as the Jacobi field. which obeys a first order differential equation known as the Jacobi equation, with particular boundary conditions which we will state precisely later. In section 3.1, we solve the Jacobi equation for the velocity perturbation $\delta V(s)$. In section 3.2, we compute the first order change $\delta 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_{c}=\frac{2\pi}{E_{\text{max}}-E_{\text{min}}}$, where $E_{\text{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 for sub-exponential times, and thus the linear geodesic is at least locally minimizing until times exponential in $N$.

3.1 Solving the Jacobi equation

In order to discover a conjugate point, we must deform the base geodesic with a velocity perturbation $\delta V(s)$ which solves the geodesic equation to first order; this first order equation for $\delta 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\to Ht+\delta 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 being $\bar{c}=1+\mu$ (where $\mu\sim e^{S}$ as explained previously), 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

[TABLE]

where the subscripts $L$ (local) and $NL$ (nonlocal) denote projection into the local and nonlocal subspaces, i.e.

[TABLE]

where $T_{\alpha}$, $T_{\dot{\alpha}}$ are bases for the local and nonlocal subspaces respectively.

Let $\mathcal{A}_{L}$ denote the vector space spanned by the local generators in the Lie algebra, and $\mathcal{A}_{NL}$ denote the vector space spanned by the non-local generators. In order to solve the Jacobi equation, note that the second equation in (3.1) involves the super-operator $\mathbf{C}:\mathcal{A}_{NL}\to\mathcal{A}_{NL}$ defined by

[TABLE]

The nonlocal equations can be solved by introducing a new basis $\widetilde{T}_{\dot{\alpha}}$ for the nonlocal subspace such that $\mathbf{C}$ is diagonal161616This is always possible since $\mathbf{C}$ is Hermitian viewed as a matrix acting on $\mathcal{A}_{NL}$, and the spectral theorem of linear algebra states that Hermitian matrices may always be unitarily diagonalized with real eigenvalues.:

[TABLE]

In this basis, we can write $\delta V_{NL}=\sum_{\dot{\alpha}}\delta\widetilde{V}^{\dot{\alpha}}\widetilde{T}_{\dot{\alpha}}$, where we note that the $\delta\widetilde{V}^{\dot{\alpha}}$ are numbers (i.e., the coefficients) while $\widetilde{T}_{\dot{\alpha}}$ are hard/non-local operators in the Lie algebra. The nonlocal equations become (in components)

[TABLE]

with no summation over $\dot{\alpha}$ on the right-hand side. The solution is therefore

[TABLE]

Plugging this into the local equations, we have

[TABLE]

The local solution is

[TABLE]

3.2 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 $\delta V$. The exact final unitary and first order perturbation are, recalling (2.9),

[TABLE]

where the $\delta U$ term is obtained by expanding the path-ordering in a Dyson series and taking the first term,

[TABLE]

We now define a super-operator $\mathbf{Y}_{(\mu)}:\delta V(0)\to U^{-1}\delta U(1)$ which acts by

[TABLE]

where we have inserted our solution for $\delta V(s)$ into equation (3.10). A conjugate point, in this formalism, is given by the condition

[TABLE]

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 $\mathbf{Y}_{(\mu)}$ 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 $\mathfrak{su}(2^{N/2})$, it is convenient to view $\mathbf{Y}_{(\mu)}$ as acting on the complexification of this vector space, i.e., on $\mathfrak{sl}(2^{N/2},\mathbb{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 $\mathbb{R}$, and so the eigenvalues of $\mathbf{Y}_{(\mu)}$ 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 $\mathbb{C}$ even though the equation involves coefficients only in $\mathbb{R}$.) Of course, in order for a true conjugate point to appear for some values of $t$ and $\mu$, the zero mode must be Hermitian and traceless. In other words, it must indeed be a valid element of $\mathfrak{su}(2^{N/2})$.

3.3 Conjugate points in the bi-invariant case

Solving for the conjugate points at general values of $\mu$ 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 $\mu=0$ where we can obtain all the conjugate points exactly. This is because at $\mu=0$, where all generators are considered computationally “easy”, equation (3.11) simplifies greatly, and we get

[TABLE]

It is easy to guess the eigenvectors of this super-operator: they are simply the energy eigenstate projectors $|m\rangle\langle n|$, where $|m\rangle$ and $|n\rangle$ are eigenstates of the Hamiltonian with energies $E_{m}$ and $E_{n}$ respectively. Indeed, we find

[TABLE]

where the eigenvalue $\phi_{mn}(t)$ is given by171717Notice that diagonal projectors $|n\rangle\langle n|$ have constant eigenvalue 1 and therefore cannot lead to conjugate points in the bi-invariant case.

[TABLE]

The eigenvalue $\phi_{mn}$ (for $E_{m}\neq E_{n}$) becomes zero at

[TABLE]

Indeed, at these times, the eigenvalues corresponding to both $|m\rangle\langle n|$ and $|n\rangle\langle 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

[TABLE]

In the SYK model, the maximum separation is known to be $\Delta_{\text{max}}\sim N$, and so the linear geodesic stops being minimal after $t_{c}\sim\frac{2\pi}{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 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.

3.4 Turning on cost factors

We will turn on a finite cost factor $\mu$, 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 $\mathbf{Y}_{(\mu)}$ (although it is possible to calculate this spectrum perturbatively in $\mu$, see appendix B). However, if the Hamiltonian is sufficiently chaotic, then the situation simplifies greatly. More precisely, if the off-diagonal eigenstate projectors $|m\rangle\langle 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 $\mathbf{Y}_{(\mu)}$ at finite $\mu$. 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\rangle,|n\rangle$ etc., $T_{\alpha}$ be the local generators in the Lie algebra, and $T_{\dot{\alpha}}$ be the non-local generators. Define

[TABLE]

We will say that the Hamiltonian and the gate set satisfy the eigenstate complexity hypothesis, if in the large-$N$ limit for $E_{m}\neq E_{n}$,

[TABLE]

where $S$ is $\ln\text{dim}$ of the Hilbert space (i.e., $S=\frac{N}{2}\ln\,2$ for the SYK model), $\text{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

[TABLE]

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||=[\mathrm{Tr}(X^{\dagger}X)]^{1/2}$.

The physical intuition behind this criterion is that off-diagonal projectors of the form $|m\rangle\langle n|$ map the energy eigenstate $|n\rangle$ to a different eigenstate $|m\rangle$. 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\rangle\langle 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\rangle\langle 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$, $\mathcal{J}=1$ and $k=4,q=3$. We see that the off-diagonal elements of $R_{mn}$ are indeed exponentially suppressed. Taking

[TABLE]

where the $N$-dependent coefficient is 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 of $\bar{r}_{s}=1$, and (sample) standard deviation of $\sigma_{s}=0.14$.181818It is easy to show from the definition of $r_{mn}$ that their mean is one: $\frac{1}{2^{N}}\sum_{m,n}r_{mn}=1$. The distribution $P(r_{mn})$ can be roughly approximated by the normal distribution with mean $\bar{r}\simeq 1$ and standard deviation $\sigma\simeq 0.098$. A slightly better approximation is provided by Student’s t-distribution with the parameters $\bar{r}=0.994,\sigma=0.093$ and the number of degrees of freedom $\nu=6$. We have also checked other values of $N$ and $(k,q)$ (with $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 PhysRevA.43.2046 ; 1994PhRvE..50..888S ; DAlessio:2016rwt (see Sonner:2017hxc ; Hunter-Jones:2017raw ; Nayak:2019khe 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}\sim\psi_{a_{1}}\cdots\psi_{a_{m}}$, 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\rangle\langle 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})$.191919More precisely, the ETH suppression to $R_{mn}$ is $e^{-S(\bar{E})}$ where $\bar{E}=(E_{m}+E_{n})/2$, but we expect $e^{S-S(\bar{E})}$ to be polynomial in $S$. Since there are at most polynomially many local/easy generators (assuming $k$ does not scale with $N$), we deduce that $R_{mn}=O(\text{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\rangle\langle 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. We expect that large-$N$ integrable Hamiltonians should violate ECH, but it would be interesting to study this in greater detail. Certainly, off-diagonal operators of the form $|m\rangle\langle 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\rangle\langle n|$ is one, this naturally requires the individual overlaps $\langle n|T_{i}|m\rangle$ to be larger. For instance, in appendix C we show numerical evidence that for $q=2$ (quadratic) SYK-like Hamiltonians, the individual overlaps $\langle n|T_{i}|m\rangle$ are all $O(1)$, and hence $R_{mn}$ only has an $e^{-S}$ suppression, as opposed to $e^{-2S}$ in the chaotic case (or equivalently $||\,|m\rangle\langle n|_{L}||\sim e^{-S/2}$ as opposed to $e^{-S}$). Indeed, our arguments below for complexity growth will crucially rely on this enhanced suppression in chaotic systems.

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

[TABLE]

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 $\left[H,|m\rangle\langle n|_{L}\right]_{NL}$ can at most get a polynomial enhancement over the exponentially suppressed norm of $|m\rangle\langle n|_{L}$. This implies that if we take our initial velocity to be $\delta V(0)=|m\rangle\langle n|$, then the solution (3.6), (3.8) to the Jacobi equation simplifies substantially

[TABLE]

Below, we will carefully justify that the corrections to equation (3.23), denoted as $\cdots$ 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 $\mathbf{Y}_{(\mu)}$ on $|m\rangle\langle n|$:

[TABLE]

where, once again, the correction terms are exponentially suppressed in $N$, as will be justified below. Performing the $s$ integration, we find

[TABLE]

(Note that the function $\phi_{mn}$ was defined in (3.15).) Therefore, under the assumption (3.20), the super-operator $\mathbf{Y}_{(\mu)}$ is a diagonal matrix in the $|m\rangle\langle n|$-basis with the diagonal entries given by $\phi_{mn}\left(\frac{t}{1+\mu}\right)$. The diagonal entry corresponding $|m\rangle\langle n|$ becomes zero at

[TABLE]

up to $O(e^{-S}\text{poly}(S))$ corrections.202020Each conjugate point is two-fold degenerate at leading order in $N$, and the exponentially suppressed corrections may split this two-fold degeneracy. Therefore, as we crank up the cost factor $\mu$, all the diagonal entries become approximately equal to one. Indeed, the first time at which one of these diagonal entries becomes zero moves to later and later time (see Fig. 9) is now approximately located at

[TABLE]

If we take the cost factor to be

[TABLE]

where $\epsilon$ is some small positive number (as will become clear shortly, the above argument works when $\mu t\ll e^{S}$ which is satisfied by this choice at any sub-exponential time), then $\mathbf{Y}$ will be approximately the identity matrix for any time sub-exponential in $S$, and so we do not expect it to have zero modes. Of course, there is an important caveat here – although the corrections to $\mathbf{Y}$ are exponentially small, the matrix in question is exponentially large and so one might worry that the eigenvalues of $\mathbf{Y}$ get corrected at $O(1)$. We will address this issue below.

Here we have assumed that $\Delta_{\text{max}}$ does not scale exponentially with $N$. Indeed, for the SYK model $\Delta_{\text{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 $\mu=e^{(1-\epsilon)S}$. In the bi-invariant case, all the diagonal projectors $|m\rangle\langle m|$ are eigenvectors of the super-operator $\mathbf{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 $\mu$: since conjugate points are zero modes of $\mathbf{Y}_{(\mu)}$, they cannot simply appear out of nowhere; as we can see in Fig. 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 Morse1934 .

We also note here that if there is an off-diagonal projector $|m\rangle\langle 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 can similarly show that such a projector corresponds to an approximate eigenvector of $\mathbf{Y}_{(\mu)}$ with the eigenvalue $\phi_{mn}(t)$. In this situation, we may expect to find conjugate points at the $O(1)$ times $t=\frac{2\pi}{\Delta_{mn}}\mathbb{Z}$, provided $\Delta_{mn}$ is not exponentially small, and if so 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 $\delta V_{mn}(s)$ be the Jacobi field along the linear geodesic with the initial condition $\delta V_{mn}(0)=|m\rangle\langle n|$, and define

[TABLE]

where $c(s)=e^{-\frac{i\mu ts\Delta_{mn}}{1+\mu}}$, and $\delta 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 $\delta W$:

[TABLE]

where the source term $\mathcal{S}$ above is given by

[TABLE]

As long as $\mu t\ll e^{S}$, say for instance $\mu t\sim e^{(1-\epsilon)S}$, then the source terms have an exponentially suppressed norm by ECH:

[TABLE]

where $S=\frac{N}{2}\ln(2)$, and we are using the Frobenius norm $||X||^{2}=Tr(X^{\dagger}X)$. For polynomial times, we can therefore take $\mu\sim e^{(1-\epsilon)S}$, and the source terms will still be suppressed; beyond this value of $\mu$ our arguments here will break down.212121Note that for integrable systems, weaker suppression implies that our argument breaks down at $\mu\sim e^{S/2}$, that is far before the required value of $e^{S}$ for the cost factor. For chaotic systems, we can push the cost factor to $e^{(1-\epsilon)S}$, which is almost the required value. Expressing $\delta W$ in terms of the basis $(T_{\alpha},\widetilde{T}_{\dot{\alpha}})$ introduced previously and solving the second equation in (3.29), we obtain for the non-local piece of $\delta W$:

[TABLE]

Therefore, the norm of this correction term is given by

[TABLE]

where we have used $||\int X||\leq\int||X||$. Repeating the same argument for the local directions, we see that in fact

[TABLE]

Now coming to the action of the super-operator on $|m\rangle\langle n|$, we have the exact statement

[TABLE]

We can bound the norm of the second term above by once again using $||\int X||\leq\int||X||$, together with equation (3.34):

[TABLE]

This completes our justification that the corrections to equation (3.25) are indeed exponentially suppressed in $N$. The upshot of these arguments is that, for the parameter regimes we are interested in, the functions $\phi_{mn}\left(\frac{t}{1+\mu}\right)$ remain close to one and all other contributions are suppressed. Therefore, no zero modes of the super-operator can develop before at least one $\phi_{mn}$ has dropped away from 1, and this does not occur until times exponential in $N$.

We will now address a possible caveat in the above discussion: we have shown that the superoperator Y is an approximately diagonal matrix with the diagonal entries being approximately one at times much smaller than $t_{c}$, and exponentially suppressed off-diagonal entries. So let us write

[TABLE]

where $\Phi$ is the diagonal part and $\delta\textbf{Y}$ is the off-diagonal part. Equation (3.36) shows that the $L_{2}$ norm of any row in $\delta\textbf{Y}$ is bounded by an exponentially small quantity. One might worry that, since there are exponentially many of these rows, they may combine to lead to significant deviations in the eigenvalues of Y compared to $\Phi$.222222We thank Daniel Ranard for emphasizing this point to us. The point, however, is that the constraint on the norm of the individual rows of $\delta\textbf{Y}$ is strong enough that almost all the exponentially many eigenvalues can only receive exponentially small corrections, while only an $O(1)$ number of eigenvalues can be affected significantly. We can see this by estimating the average magnitude of the eigenvalues of $\delta\textbf{Y}$,232323Note that we can focus on the eigenvalues of $\delta\textbf{Y}$ because $\Phi$ is approximately the identity matrix, up to exponentially small diagonal corrections. To be more systematic, we can absorb these corrections inside $\delta\textbf{Y}$, and then make the remainder of the argument. By the Cauchy-Schwarz inequality, the rows of this newly defined $\delta\textbf{Y}$ also obey a bound on their norms. which is exponentially suppressed because of (3.36). Furthermore, the variance in the distribution of the eigenvalues can also similarly be shown to be exponentially small. Thus, almost all the eigenvalues of Y will be unaffected by the correction term $\delta\textbf{Y}$, and thus be bounded away from zero, i.e., almost all the conjugate points (which were present at small $\mu$) will get lifted. It is nevertheless true that this argument does not preclude large corrections to a small number of the eigenvalues, and thus does not completely rule out “accidental” conjugate points; it will be interesting to see if this can be accomplished by using more detailed properties of $\delta\textbf{Y}$.

We emphasize that the potential remaining conjugate points discussed above are “accidental” from the perspective of a random family of Hamiltonians in the following sense: in quantum circuit complexity, we are concerned with families of Hamiltonians and therefore with families of conjugate points. In the bi-invariant analysis, we found that conjugate points were very generic close to the identity, and specifically that any family of random Hamiltonians will have a family of conjugate points in the bi-invariant metric with distances from the identity set by the total spectral range. If the entries of the Hamiltonian have mean zero and unit variance, this conjugate point family is actually moving closer to the identity as we increase $N$. The “accidental” conjugate points above are not generic in this way, and require some fine tuning of the matrix $\delta\textbf{Y}$. Therefore, we do not expect them to exist in families (i.e. for arbitrary $N$), and even if we are unfortunate enough to encounter such a family, we expect that a small perturbation of the Hamiltonians will destroy them.

4 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 Nielsen2007 , which we study using the Euler-Arnold equation Balasubramanian2018 . 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 Kitaev2015 ; Maldacena2016 .

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.242424Recall that a similar dip-ramp-plateau pattern in correlation functions appeared in integrable theories without disorder as studied in Balasubramanian2019 . 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 disorder-averaging 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 disorder-averaging 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).

4.1 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 Susskind2018 . 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 $\mu=0$ to take all operators to be “simple”, we can get a glimpse of the relevant physics.252525See Yang:2018nda for a different interpretation of the physics of the bi-invariant situation. In this case, the geodesic equation and boundary conditions are

[TABLE]

The solutions are given by

[TABLE]

where $|n\rangle$ are the energy eigenstates of the Hamiltonian, and $\vec{k}=(k_{1},\cdots,k_{2^{N/2}})$ are integers which sum up to zero because of the traceless condition. Therefore, the complexity is given by

[TABLE]

where the minimization is over all integer vectors $\vec{k}$ subject to the constraint that they sum up to zero. In Fig. 10 we show a numerical plot262626We used the function NMinimize in Mathematica to make this plot. We cannot guarantee that numerical minimization has converged to the true global minimum, and so these plots should be regarded as upper-bounds on the true complexity. of $\mathcal{C}(t)$ for ten different realizations of the SYK model with $N=8,\,q=3$ and $\mathcal{J}=1$. In all the cases, complexity grows linearly with time for a while, but then saturates at a time of order $t\sim 2\pi/\Delta_{\text{max}}$, which recall is precisely the time when conjugate points appear in the bi-invariant case! Another important feature to note here is the saturation for each individual SYK realization, arising from different geodesics (i.e., different integer vectors $\vec{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},\cdots,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” ($\mu=0$), with cost $1+\mu=1+e^{(1-\epsilon)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.

4.2 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 ($\mu=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 $\vec{k}\in\mathbb{Z}^{2^{N/2}}$ with $\sum_{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 $\mu\neq 0$ these trajectories should no longer be geodesics. Perturbatively, it is obvious that their length increases with $\mu$, 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+\mu\sim 2^{N}$) have greater length than the linear solution for any $t\sim\text{poly}(N)$, where $\text{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.272727Hamiltonian simulation has been studied in the context of complexity classes previously in quantum computation Chuang2018 . This is due to a theorem of Aaronson and Susskind Aaronson2018 , who showed that the complexity class separation PSPACE $\nsubseteq$ 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$.282828PSPACE 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 $\epsilon$. Finally, BQSUBEXP is the class of quantum computations that can be done in a subexponential time, $t\in O(e^{N^{\alpha}})$ (for all $\alpha>0$), and BQSUBEXP/subexp is the same class but with subexponential size advice strings. Of course, we expect such growth only for chaotic Hamiltonians, and not in integrable systems. Furthermore, following Theorem 2 in Aaronson2018 , if the geodesic loop argument outlined above can be made for exponential times (or, more precisely, for times greater than any subexponential292929Note 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^{\alpha}}$ for all $\alpha>0$. Our bound on conjugate points holds for a truly exponential time $t_{c}\sim e^{\epsilon N}$, so the discussion of BQSUBEXP/subexp in the main text holds for whichever definition was used by Aaronson2018 . In the previous footnote we assumed the weaker SUBEXP = $\bigcap_{\alpha>0}$ DTIME$(2^{N^{\alpha}})$.), then our our results (which show there are no conjugate points up to exponential time $t_{c}\sim e^{\epsilon N}$) would actually imply the even stronger statement PSPACE $\nsubseteq$ 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 complexity-theoretic assumption in Aaronson2018 where the time evolution step $e^{-iH}$ was taken to implement one step of a reversible computationally-universal classical cellular automaton.303030It 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 Neumann1966 . We consider it very likely that a generic chaotic Hamiltonian like the SYK model implements such an automaton via its time evolution steps.

When conjugate points exist in our analysis they can be interpreted in terms of “fast-forwarding” 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$ Aharonov2017 . General Hamiltonian simulation algorithms are well-studied in the quantum computation literature Berry2007 ; Childs2011 ; Childs2012 ; Berry2014 ; Berry2015 ; Aharonov2017 ; Aharonov2018 . In particular Aharonov2017 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 Aharonov2017 , but the absence of a conjugate point is certainly necessary to rule out such speedups. 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 Aharonov2017 to detect speedups.

4.3 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\rangle\langle 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\sum_{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 $Z$s. Because of this there are $\sim e^{S}$ generators which will have overlaps with a given eigenstate projector with a few spin flips such as $|...,1,...,0,...\rangle\langle...,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 $\sim 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.

Acknowledgments

We are grateful to Yosi Atia, Adam Bouland, Lampros Lamprou, Clélia de Mulatier, Daniel Ranard, Gábor Sárosi, Edward Witten, Zhenbin Yang, and Ying Zhao for useful discussions. VB, OP, AK and MD were supported in part by the Simons Foundation through the It From Qubit Collaboration (Grant No. 38559), and also by the DOE (Contract No. FG02-05ER-41367 & QuantISED grant DE-SC0020360). MD is supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1845298.

Appendix A Majorana Fermion Basis for $\mathfrak{su}(2^{N/2})$

We begin with a set of Majorana fermion operators $\gamma_{i}$ which obey the commutation relations313131Note 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.

[TABLE]

and also obey $\gamma_{i}^{\dagger}=\gamma_{i}$ (the Majorana condition). We can interpret these objects as $2^{N/2}\times 2^{N/2}$ Hermitian matrices, and they are precisely the generalized gamma matrices of the Clifford algebra $\mathcal{C}\ell_{N}(\mathbb{R})$. The basis for $\mathfrak{su}(2^{N/2})$ is constructed by taking products of $\gamma_{i}$ with appropriate factors of $i$ to ensure Hermiticity. Specifically, we consider ordered products $\gamma_{i_{1}}\dots\gamma_{i_{n}}$ with $i_{1}<\dots<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\in\mathcal{B}_{N}$. The bits of the string are $b=b_{N}\dots b_{1}$, and let $q_{b}$ be the number of nonzero bits in $b$. We write

[TABLE]

and we then have

[TABLE]

We now show these generators are traceless. By construction, the gamma matrices individually are traceless, so we have $\operatorname{tr}\gamma_{i}=0$. Additionally, an even number of them will be traceless since we have anticommutation and cyclicity of the trace:

[TABLE]

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 $\gamma_{a}^{(N-2)}$, we may create a set of $N$ $\gamma_{i}^{(N)}$ by taking

[TABLE]

There are four cases for $\operatorname{tr}\gamma_{i_{1}}\dots\gamma_{i_{2k+1}}$. First, neither $\gamma_{N-1}^{(N)}$ nor $\gamma_{N}^{(N)}$ appear in the product. In that case, we have $\gamma_{i_{1}}^{(N)}\dots\gamma_{i_{2k+1}}^{(N)}=\gamma_{i_{1}}^{(N-2)}\dots\gamma_{i_{2k+1}}^{(N-2)}\otimes\sigma_{3}^{2k+1}$, and since $\operatorname{tr}A\otimes B=\operatorname{tr}A\operatorname{tr}B$ and $\operatorname{tr}\sigma_{3}^{2k+1}=\operatorname{tr}\sigma_{3}=0$, we have the required result. The second case is when either $\gamma_{N-1}^{(N)}$ or $\gamma_{N}^{(N)}$ appear in the product, but not both. Then the final tensor factor is either $\sigma_{1}$ or $\sigma_{2}$, and similarly we have $\operatorname{tr}\sigma_{1}=\operatorname{tr}\sigma_{2}=0$ so again the entire trace vanishes. The interesting case is when we have both $\gamma_{N-1}^{(N)}$ and $\gamma_{N}^{(N)}$ in the product. Then we have $\operatorname{tr}\gamma_{i_{1}}^{(N)}\dots\gamma_{i_{2k+1}}^{(N)}=-\operatorname{tr}\gamma_{i_{1}}^{(N-2)}\dots\gamma_{i_{2k-1}}^{(N-2)}\operatorname{tr}\sigma_{1}\sigma_{2}\sigma_{3}$. We now repeat the argument for $\operatorname{tr}\gamma_{i_{1}}^{(N-2)}\dots\gamma_{i_{2k-1}}^{(N-2)}$, since these smaller gamma matrices have a similar tensor product structure. Again the only interesting case is when both $\gamma_{N-3}^{(N-2)}$ and $\gamma_{N-2}^{(N-2)}$ both appear. Following this chain, we end up with only a single gamma matrix in the first part of the tensor structure, and the trace of any single gamma always vanishes. So we conclude

[TABLE]

Thus, all generators are traceless as desired.

[TABLE]

We now turn to linear independence. Notice first that, given generators $T_{a}$ and $T_{b}$ with $a\neq b$, we have $T_{a}T_{b}=\alpha T_{c}$ for some $\alpha\in\mathbb{C}$ and some $c\in\mathcal{B}_{N}\setminus\{0\}$. Now assume for the sake of contradiction that we have $\sum_{b}\alpha_{b}T_{b}=0$ for some constants $\alpha_{b}\in\mathbb{R}$. Solve for a specific $T_{a}$ with nonzero coefficient and write $\alpha_{a}T_{a}=-\sum_{b\neq a}\alpha_{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 $1\propto\sum_{b\neq a}\alpha_{b}T_{a}T_{b}\propto\sum_{c\neq a}\beta_{c}T_{c}$. However, if we now take the trace of both sides, the left hand side is $\operatorname{tr}1=2^{N/2}$ but the right hand side is $\sum_{c}\alpha_{c}\operatorname{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 $\operatorname{tr}T_{a}T_{b}\propto\delta_{ab}$. This can be seen by noticing $\operatorname{tr}T_{a}[T_{b},T_{c}]=\operatorname{tr}T_{a}{f_{bc}}^{d}T_{d}\propto{f_{bc}}^{d}\delta_{ad}={f_{bc}}^{a}$, and also by cyclicity of the trace we have

[TABLE]

therefore we have

[TABLE]

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\oplus b$ be the bitwise “exclusive-or” of $a$ and $b$. Let $a\land b$ be the bitwise “and” of $a$ and $b$. Then a lengthy calculation shows

[TABLE]

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\oplus b=c$ and $q_{a}q_{b}+q_{a\land b}$ is odd, do we have (by total antisymmetry) $b\oplus c=a$ and $q_{b}q_{c}+q_{b\land c}$ odd? Indeed we do, by properties of $\oplus$ and $\land$:

[TABLE]

Another computation shows the parity of $q_{a}q_{b}+q_{a\land b}$ matches that of $q_{b}q_{c}+q_{b\land c}$.

[TABLE]

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}\leq 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\land b}$, we must have $k>2q_{a\land b}$. For $k=2$, this implies $q_{a\land b}=0$. Combining with $q_{a}q_{b}+q_{a\land b}\equiv 1\mod 2$, we see $q_{a}q_{b}\equiv 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\land b}=0$, but this cannot be the case because then $q_{a}q_{b}\equiv 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\land 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 $\frac{dV^{\alpha}}{ds}=0$.

Appendix B Conjugate points in perturbation theory

In the main text it was argued that the bi-invariant metric (with $c=\bar{c}=1$) has conjugate points at $t_{mn}^{\star}=\frac{2\pi}{\Delta_{mn}}\mathbb{Z}$, where $\Delta_{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 factor $\bar{c}=1+\epsilon$ along the heavy directions. Recall that the equation for the Jacobi field (i.e., the Euler-Arnold equation linearized around the linear geodesic) takes the form

[TABLE]

[TABLE]

Here the subscripts $L$ and $NL$ 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 $\delta V(0)$ and some value of $t$, such that

[TABLE]

If so, then the corresponding value of $t$ constitutes a conjugate point along the original linear geodesic $U_{\text{linear}}(s)=e^{-istH}$, at which point the linear geodesic stops being globally minimizing.

We take $c=1$ and $\bar{c}=1+\epsilon$, and expand

[TABLE]

Expanding equations (B.1) and (B.2) at zeroth order, we find that $\delta V^{(0)}(s)=\delta V^{(0)}(0)$, i.e., it is $s$-independent. Indeed, this is the bi-invariant Jacobi field, which gives the family of conjugate points at

[TABLE]

for any pair of distinct eigenvalues $m$ and $n$ of the Hamiltonian with $\Delta_{mn}=(E_{m}-E_{n})$, and any non-zero complex number $z$. We can set the absolute value of $z$ to one by choice of normalization, but its phase is not determined at this order. This implies a two-fold “degeneracy” in all the bi-invariant conjugate points.

At first order in $\epsilon$, we find the equations

[TABLE]

[TABLE]

which can be solved to obtain

[TABLE]

From equation (B.3), we then find that the final displacement at the present order is given by

[TABLE]

where found it more convenient to write the matrix elements of $U^{-1}\delta U(1)$ in the energy eigenstates $|m^{\prime}\rangle,\,|n^{\prime}\rangle$, and we have defined the function

[TABLE]

Let us now return to the conjugate points $t_{mn}^{\star}=\frac{2\pi}{\Delta_{mn}}\mathbb{Z}$ we had obtained at zeroth order. From equation (B), it is clear that their locations have now moved at linear order in $\epsilon$, which we can keep track of systematically in perturbation theory:

[TABLE]

where $t^{\star,(0)}_{mn}=\frac{2\pi}{\Delta_{mn}}\mathbb{Z}$ denotes the bi-invariant conjugate points. Substituting this expansion in equation (B) and demanding that the result vanish for $m^{\prime}=m$ and $n^{\prime}=n$,323232Of course, we should also demand that the matrix elements of $U^{-1}\delta U(1)$ vanish for $m\neq m^{\prime}$ and $n\neq n^{\prime}$; these constraints partially determine $\delta V^{(1)}(0)$. we deduce the shifts in the conjugate points:

[TABLE]

Note that the numerator on the right hand side involves the projection of $\delta V^{(0)}(0)$ to the non-local subspace, which makes the above expression somewhat non-trivial to evaluate. Importantly, however, $\delta V^{(0)}(0)$ at the zeroth order was only determined up to an arbitrary complex number $z$:

[TABLE]

Requiring that our expression for $\delta^{(1)}t_{mn}^{\star}$ is real determines precisely two possible choices of $z$, which we may call $z_{\pm}$. (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_{\pm}$, we then have the two conjugate points at their respective locations $t^{\star,(0)}_{mn}+\epsilon\delta^{(1)}t^{\star}_{mn}(z_{\pm})$, given by (B.11). Therefore, we find that at linear order in $\epsilon$, the two-fold “degeneracy” in conjugate points splits. Nevertheless, they continue to exist and we have tracked their locations at $O(\epsilon)$ 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 factor $\bar{c}=1+\mu$, the conjugate points move to $t\sim t^{(0)}(1+\mu)$, 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 $\delta V^{(0)}_{L}$ term in equation (B.11).

Appendix 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\rangle\langle n|$ in the SYK model in terms of the generators $T_{a}=(T_{\alpha},T_{\dot{\alpha}})$ consisting of products of fermions:

[TABLE]

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=\psi_{1}\psi_{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\rangle\langle 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.