Circuit complexity in quantum field theory

Motivated by recent studies of holographic complexity, we examine the question of circuit complexity in quantum field theory. We provide a quantum circuit model for the preparation of Gaussian states, in particular the ground state, in a free scalar field theory for general dimensions. Applying the geometric approach of Nielsen to this quantum circuit model, the complexity of the state becomes the length of the shortest geodesic in the space of circuits. We compare the complexity of the ground state of the free scalar field to the analogous results from holographic complexity, and find some surprising similarities.


Introduction
Recent years have seen exciting progress in understanding the connection between entanglement and geometry [1][2][3][4][5][6]. However, in the context of the AdS/CFT correspondence, our ability to decipher the bulk geometry (or bulk physics, more generally) from information in the boundary CFT remains very incomplete. The challenges are most pronounced if one considers physics behind the horizon of a black hole. Consider for example the eternal AdS black hole, which is dual to the thermofield double (TFD) state [7] |T F D(t L , t R ) = 1 Z β i e −βE i /2 e −iE i (t L +t R ) |i L |i R . (1.1) This describes an entangled state of the two copies of the CFT associated with the asymptotic boundaries (see figure 1), which are joined by a wormhole, i.e., an Einstein-Rosen bridge (ERB), in the bulk [8]. The AdS/CFT correspondence demands that the interior region have an equivalent description in terms of the boundary field theory. But now, in addition to the usual difficulties involved in probing behind the horizon, we have another conundrum: the boundary field theory reaches thermal equilibrium very quickly, on the order of the thermalization time 1/T , while the ERB continues to grow on much longer timescales [9]. Therefore, there must be some quantity in the field theory that corresponds to this fine-grained information -which is evidently not captured by entanglement entropy [10,11] -that continues to evolve long after thermal equilibrium is reached. These considerations led Susskind to introduce holographic complexity as the boundary entity whose growth corresponds to the evolution of the ERB [12][13][14]. In particular, with his collaborators, he developed two new gravitational observables, both of which successfully probe the late-time growth of the ERB. The first of these is referred to as the complex-ity=volume (CV) conjecture, which posits that the complexity of the boundary state is proportional to the volume of a maximal codimension-one bulk surface B that extends to the AdS boundary, and asymptotes to the time slice Σ on which the boundary state is defined [9,12]: where is some length scale associated with the bulk geometry, e.g., the AdS radius or the radius of the black hole. For example, in the eternal AdS black hole, this bulk surface connects the time slices denoted t L and t R on the left and right boundaries through the ERB; see the left panel in figure 1. The second proposal is the complexity=action (CA) conjecture. This identifies the complexity of the boundary state with the gravitational action evaluated on a bulk region known as the Wheeler-DeWitt (WDW) patch [15,16]: One can think of the WDW patch as the causal development of the spacelike surface B picked out by the CV construction. The right panel in figure 1 illustrates the WDW patch for the example of the eternal AdS black hole, where the CFT state is again evaluated on the t L and t R slices of the left and right boundaries, respectively. Both proposals have their merits, as well as certain shortcomings. In any case, they bring to our attention two new classes of interesting gravitational observables which should certainly be studied in further detail. In fact, various aspects of the proposals and these new observables have been examined in a number of recent papers, e.g., [17][18][19][20][21][22][23]. And while both the CV and CA conjectures appear to provide viable candidates for holographic complexity, this research program is still at a very preliminary stage. In particular, one would like to establish a concrete translation of the new observables in the bulk to a specific quantity in the boundary theory, e.g., as was recently found for holographic entanglement entropy [24][25][26]. However, a stumbling block to this endeavor is finding the answer to an even simpler question: what does "complexity" mean in the boundary CFT?
This question is the focus of the present paper. Specifically, our objective is to provide the first steps towards defining circuit complexity in quantum field theory (QFT). A precise understanding of this quantity will not only shed light on the CV and CA proposals, but is also an interesting question deserving of study in its own right. For example, it may also provide new insights into quantum algorithms for the simulation of quantum field theories [27][28][29][30], or more generally into Hamiltonian complexity [31,32], or the efficient description of many-body wave functions [33,34].
In computer science, the notion of computational complexity refers to the minimum number of operations necessary to implement a given task [35,36]. In the present context, the task of interest will be the preparation of a state in the QFT, and we will define the complexity in terms of a quantum circuit model. That is, we will begin with a simple reference state |ψ R , and construct a unitary transformation U that produces the desired target state |ψ T via |ψ T = U |ψ R . (1.4) The unitary U will be constructed from a particular set of simple elementary or universal gates, which can be applied sequentially to the state. When working with such discrete operations, we should also introduce a tolerance ε so that even if we cannot achieve the precise equality above, we may still judge the transformation to be successful when the two states are sufficiently close to one another according to some distance measure, i.e., |ψ T − U |ψ R 2 ≤ ε . (1.5) Of course, there will not be a unique circuit which implements the desired transformation (1.4): generally there will exist infinitely many sequences of gates which produce the same target state. However, the complexity of the state |ψ T may be defined as the minimum number of gates required to produce the transformation (1.4), i.e., the complexity is the number of elementary gates in the optimal or shortest circuit. The challenge then is to identify this optimal circuit from amongst the infinite number of possibilities. Our work takes inspiration from the geometric approach of Nielsen and collaborators [37][38][39], which itself was developed using ideas from the theory of optimal quantum control, e.g., [40][41][42][43]. In Nielsen's case, the question of interest was to find the minimal size quantum circuit required to exactly implement a specified n-qubit unitary operation U (without the use of ancilla qubits). Neilsen approaches this question as the Hamiltonian control problem of finding a time-dependent Hamiltonian H(t) that synthesizes the desired U , where where the Hamiltonian is expanded in terms of generalized Pauli matrices, denoted here as M I , 1 and the P indicates a time ordering such that the Hamiltonian at earlier times is applied to the state first, i.e., the circuit is built from right to left. In [37], the control functions Y I form a (4 n − 1)-dimensional vector space, and can be seen as specifying the tangent vector to a trajectory in the space of unitaries, In this general space, the paths of interest satisfy the boundary conditions U (t = 0) = 1 and U (t = 1) = U . Neilsen's idea is then to define a cost for the various possible paths D(U (t)) = 1 0 dt F U (t),U (t) , (1.8) and to identify the optimal circuit or path by minimizing this functional. In general, the cost function F (U, v) is some local functional of the position U in the space of unitaries and a vector v in the tangent space at this point. Neilsen further argues that for the present problem, a physically reasonable cost function must satisfy a number of desirable features: 1. Continuity: F should be continuous, i.e., F ∈ C 0 . these four properties come very close to defining a class of geometries known as Finsler manifolds. In particular, if we replace the first condition above with 1 . Smoothness: F should be smooth, i.e., F ∈ C ∞ , then eq. (1.8) defines length functional for a Finsler manifold, a particular class of differential manifolds equipped with a quasimetric structure in which the length of any curve is measured by a length functional of the form (1.8), with a Finsler metric F satisfying the four properties enumerated above, see e.g., [44,45]. While the familiar notion of Riemannian manifolds would fall within this definition, Finsler geometry provides a generalization to a broader class of manifolds where the norm on the tangent space is not (generally) induced by a metric tensor. Hence Neilsen has identified the problem of finding an optimal circuit with the problem of finding extremal curves, i.e., geodesics, in a Finsler geometry, and the complexity is then identified with the length of the geodesic. 2 Of course, this still leaves open the question of the precise form of the cost function, and various possibilities are examined in [37]: 3 (1. 9) In the two measures on the right, p I and q I are penalty factors which can be chosen to favour certain directions in the circuit space over others, i.e., to give a higher cost to certain classes of gates. We do not include such factors in most of our analysis, but we return to this issue in section 5. Of course, the F 2 measure yields a standard Riemannian geometry -and in fact, it will be the focus of much of our discussion. The preceding exposition of Nielsen's approach is of course very incomplete, and the interested reader is referred to [37][38][39] for more details. The key feature of this approach is that it enables one to bring the full power of differential geometry to bear on the problem of constructing the optimal quantum circuit, and this provides an objective manner in which to measure the complexity as the length of extremal paths in the geometry. However, at many points our approach will necessarily differ from that of Nielsen since we are studying a different problem, namely complexity in a quantum field theory. The primary purpose of the above presentation was to provide motivation for our geometrical analysis, but we should add that the details of Finsler geometry will not play any role in the following. Rather, a simpler physics-oriented perspective is to view the problem of finding the optimal circuit as a trajectory in the space of all possible circuits, as a classical mechanics problem for the motion of particle governed by the usual Lagrangian in eq. (1.8).
This paper is organized as follows: we begin in section 2 by examining complexity for a simple free scalar field theory. Following the preceding discussion, this requires identifying a simple reference state, introducing a set of elementary gates, and also identifying a family of interesting target states. However, the first step will be to regulate the theory by placing it on a lattice, which reduces the scalar field theory to a family of coupled harmonic oscillators. Hence, as a warm up problem, we consider the case of a single pair of harmonic oscillators. Then, having built up some intuition, we shall geometrize the problem in section 3. The main ideas from Nielsen's approach are implemented here: we represent the circuit as a path-ordered exponential analogous to eq. (1.7), show that our space of circuits forms a representation of GL(2, R), and construct the appropriate (Euclidean) metric. With this in hand, we proceed to find the geodesics, and identify the complexity of the ground state as the geodesic length of the global minimum. In section 4, we return to the field theory problem by generalizing these results to a lattice of coupled oscillators. Given the complexity for the (regulated) field theory, we then ask how our results compare to holographic complexity, and we find some surprising similarities. In section 5, we conduct a preliminary exploration of the effects of introducing penalty factors for nonlocal gates. Finally, we close in section 6 with a brief discussion of our results and directions for future work. Various technical details have been relegated to several appendices: we construct some explicit example circuits using the elementary gates given in section 2 in appendix A, elaborate on some geometrical details in appendix B, derive the normal-mode frequencies for a one-dimensional lattice in appendix C, find a closed-form approximation to the circuit complexity for the d-dimensional lattice in appendix E, and compute an approximation to the optimal circuit in the presence of penalty factors in appendix F.

Complexity for harmonic oscillators
As a first step towards understanding circuit complexity in QFT, we will consider for simplicity a free scalar field in d spacetime dimensions. However, having identified this particular QFT, we must first regulate the theory by placing it on a lattice, 4 which reduces the system to an infinite family of harmonic oscillators. This in turn suggests the much simpler warmup problem of two coupled harmonic oscillators. As it turns out, this simple model retains enough of the structure of the original problem that we will be able to learn several important lessons, which we can then carry over to the problem of circuit complexity in our scalar field theory. As in the general case, to study complexity in the two oscillator problem, we must identify a target state, a reference state, and a suitable family of elementary gates.
We begin with the Hamiltonian of a free scalar field in d spacetime dimensions, (2.1) As mentioned above, our first step is to regulate the theory by placing it on a (square) lattice with lattice spacing δ, in which case the Hamiltonian becomes: wherex i are unit vectors pointing along the spatial directions of the lattice. The resulting theory is essentially a quantum mechanical problem with an infinite family of coupled (one-dimensional) harmonic oscillators. We can make this description manifest by redefining X( n) = δ d/2 φ( n), P ( n) = p( n)/δ d/2 , M = 1/δ, ω = m and Ω = 1/δ, whereupon the Hamiltonian (2.2) takes the familiar form Hence the frequency of the individual masses is given by ω = m, and the inter-mass coupling is given by Ω = 1/δ. Now, the above suggests that we begin with an even simpler warm-up problem, namely, the case of two coupled harmonic oscillators: where x 1 , x 2 label their spatial positions, and we have set M 1 = M 2 = 1 for simplicity. Of course, to solve this system, one simply rewrites the Hamiltonian in terms of the normal modes, 5) 4 Our experience with holographic complexity suggests that we will not be able to sensibly define complexity in a QFT without a UV regulator in place [19].
where 5x This recasts the problem as that of two decoupled simple harmonic oscillators, and hence it is now straightforward to solve for the eigenstates and eigen-energies of the Hamiltonian. For example, we can write the ground-state wave function as the product of the ground-state wave functions for the two individual oscillators: where the normalization has been chosen such that d 2 x |ψ 0 | 2 = 1. We may also express this wave function in terms of the physical positions of the two masses: where We note in passing that our notation for the wave function in eq. (2.8) is slightly more general than necessary; however, these Gaussian wave functions constitute an interesting family of target states for the present exercise. 6 The next step is to identify a simple reference state. Motivated by discussions of holographic complexity [12][13][14], as well as cMERA [49], we choose a reference state where the two masses are unentangled, namely a factorized Gaussian state, (2.10) For the time being, we will simply leave ω 0 as a free parameter which characterizes our reference state. We shall examine specific choices of this frequency in section 4.1.
Having chosen our reference and target states, it remains to identify a simple set of unitary gates with which to construct the desired unitary U , which implements ψ T = U ψ R . The natural operators appearing in the quantum mechanics problem of the two coupled oscillators are the positions x 1 , x 2 and the momenta p 1 = −i∂ 1 , p 2 = −i∂ 2 , which satisfy the canonical commutation relations [x a , p b ] = i δ ab . We can use these operators to build an interesting set of elementary gates for our problem: (2.11) 5 When working in the normal-mode basis, we denote variables (e.g., positions, frequencies), with a tilde to clearly distinguish from the physical basis. The utility of this convention will become apparent later. 6 For example, Gaussian states play an important role in quantum optics, and much of our analysis is closely related to ideas developed in the quantum information literature for this purpose, e.g., [46][47][48].
where x 0 and p 0 are c-number constants. A key point is that we have introduced an infinitesimal parameter 1 into the exponent of each one of these operators. This ensures that the action of any one of these gates only produces a small change on the wave function. The action of each of these gates can be understood with the following general examples: When working with position-space wave functions, the momentum shift produced by K 1 (or K 2 ) amounts to introducing a small plane-wave component in the wave function, as illustrated in (2.12). We refer to Q 11 and Q 22 as scaling gates, for the obvious reason that these operators scale the corresponding coordinate by a small amount. Note that they also introduce an overall normalization factor, which ensures that the norm of the wave function is preserved. The operators Q 21 and Q 12 mix the positions of the two masses, thereby increasing (or decreasing) the entanglement between the two oscillators; hence we refer to these as the entangling gates. The scaling and entangling gates will play a key role in the circuits we construct below. Of course, one could extend the ensemble of gates introduced in eq. (2.11) with operators like Furthermore, one could also introduce gates with even higher powers of x's and p's in the exponent. However, we know that the collection of gates in eq. (2.11) is sufficient to implement the unitary transformation from the specified reference state (2.10) to the desired target state (2.8). Hence for simplicity, we shall work within this subset of all possible unitary gates. A circuit then consists of a sequence of these gates, whose action on ψ R produces the desired state ψ T . For example, consider the following circuit: (2.14) Here, Q 11 acts first, and by acting with the appropriate number of times α 1 , we will increase the reference frequency ω 0 appearing in front of x 2 1 in eq. (2.10) to the desired frequency ω 1 appearing in eq. (2.8). Similarly, the number of times that the Q 21 and Q 22 gates are required to appear in the circuit, namely α 2 and α 3 , are uniquely fixed by the desired ω 2 and β in the target state. The details of the corresponding calculations are given in appendix A, and the final result is We then define the circuit depth as the total number of gates in the circuit. In the above example, we have simply Note the use of the absolute values in the first line. At a pragmatic level, this is required because α 2 is negative in this particular example, i.e., β < 0. But this means that we are giving an equal complexity cost for the inverse gates Q −1 ij as for the original gates Q ij , i.e., we count the appearance of Q −1 ij as one gate in a circuit. We refer to the result in eq. (2.16) as the circuit depth of the particular circuit U given in eq. (2.14). But we must distinguish this from the complexity of the target state ψ T , which is the minimum number of gates required to produce the desired transformation. In other words, the complexity is the circuit depth of the optimal circuit. At present, we have no reason to believe that the simple circuit proposed in eq. (2.14) is the optimal circuit, and in fact, our calculations below will show that it is not.
We can describe the general form of the result in eq. (2.16) as being an overall factor of 1/ , and a coefficient determined by the various physical parameters characterizing the target and reference states. More generally, the circuit depth might be given by an expansion in , beginning with a 1/ term followed by a finite term and then potentially terms involving positive powers of . However, since 1, determining the complexity essentially requires finding the circuit which minimizes the coefficient of the leading 1/ term. For further discussion and additional examples, the interested reader may turn to appendix A.
In the next section, we apply Neilsen's approach of geometrizing the circuit complexity to find the optimal circuit. Before leaving present example however, for comparison to later results it is convenient to express the circuit depth in eq. (2.16) in terms of the normal-mode frequencies using eq. (2.9). This substitution yields

Geometrizing complexity
In the introduction, we discussed Neilsen's approach [37][38][39] of geometrizing the problem of finding the optimal circuit. We now wish to apply this geometric approach to the problem of finding the optimal preparation of the ground-state of two coupled harmonic oscillators. Our first step is to represent the circuit U as a path-ordered exponential, This structure replaces the representation of the circuits as products of the discrete gates in eq. (2.11). The connection with these gates comes about since we choose the operators O I appearing in the exponential to be precisely those appearing in the scaling and entangling gates introduced previously; that is, we write Our notation in eq. (3.1) is that the sum over I runs over the pairs ab, i.e., I ∈ {11, 12, 21, 22}. Hence in the path-ordered exponential, we can think of s as parametrizing a (continuous) product of gates, and the functions Y I (s) as indicating whether the I'th type of gate is turned on or off in this sequence (analogous to the control functions in Nielsen's time-dependent Hamiltonian (1.6)). In the integral appearing in the exponent, the differential ds plays a role analogous to that of the infinitesimal parameter . Finally, the path-ordering symbol indicates that we build the circuit from right to left, i.e., the operators at smaller values of s act on the wave function before those at larger values of s. Furthermore, with this framework, we consider a particular circuit as being constructed by following a particular trajectory, specified by Y I (s), through the space of unitary circuits. Hence we begin with U (s = 0) = 1, and have the family of unitaries Eq. (3.1) then specifies the final unitary at the end-point s = 1, which corresponds to the desired circuit that generates the target state, i.e., U fin = U (s = 1) with ψ T = U fin ψ R . From this perspective, Y I (s) specifies the velocity vector tangent to this trajectory, in a manner in which we will make precise below. In more geometric language which may be familiar from general relativity, we would say that the Y I (s) are the components of the velocity in a particular frame basis, rather than in a coordinate basis. As in the example in section 2 above, the circuit depth is determined by counting the total number of gates appearing in the full sequence comprising the circuit, cf. eq. (2.16). For our path-ordered exponential (3.1), the analogous expression becomes 7 (3.4) 7 Actually this expression (3.4) is the continuum limit of the cost function D(U ) = |αi|. Including the extra factor of in the sum eliminates the 1/ factor, so the circuit depth remains finite in the limit → 0. This cost function corresponds to the F 1 metric in the notation of [37] -see eq. (1.9). Our goal of finding the optimal circuit then amounts to finding the functions Y I (s) which yield the desired unitary U fin while minimizing this cost function. However, having also identified Y I (s) as the velocity along the trajectories U (s), we can use our physical intuition to think of this as a classical mechanics problem where we aim to find the extremal trajectory given a particular set of boundary conditions and the somewhat unusual Lagrangian in eq. (3.4).
A mentioned in the introduction, we can also make other choices for the cost function, and the analysis will go through essentially unchanged. Hence in order to develop the present problem most easily, we shall consider the F 2 or F q metric in eq. (1.9). That is, we replace eq. (3.4) with This expression should be familiar as the action of a particle moving in a curved space, and hence the optimal path corresponds to a geodesic in the corresponding (Riemannian) geometry. As we mentioned above, Y I (s) are the components of the velocity in a particular frame, for which the metric G IJ then defines the inner product. In our examples, G IJ is taken to be a purely constant (and usually diagonal) matrix. We will begin by studying the simple Euclidean metric G IJ = δ IJ , which corresponds to the F 2 metric above. With this choice, motion in every direction in the space of unitaries is assigned the same cost, i.e., the cost of each type of gate is the same. However, our notation is sufficiently general to allow for the assignment of penalty factors for particular gates, as in the F q metric. We shall return to this possibility in section 5.
To proceed further, we must find a prescription to explicitly identify the functions Y I (s). Given eq. (3.3), it is straightforward to show that However, this expression is not particularly useful as it stands. In Neilsen's construction [37][38][39], one works with unitary matrices acting on qubits, rather than operators acting on wave functions. Hence the components of the velocity analagous to eq. (3.6) can be isolated by simply tracing over the corresponding matrix generators. This procedure does not immediately lend itself to eq. (3.6), so in order to make progress, we shall re-express our problem in terms of matrices. Recall that we reduced the problem to evaluating the complexity of the ground state (2.8) of two coupled harmonic oscillators, starting from a factorized Gaussian reference state (2.10). That is, we begin and end with a Gaussian wave function; furthermore, it is straightforward to show that the scaling and entangling operators preserve the general Gaussian form of the wave function, i.e., all of the intermediate wave functions take a form analogous to eq. (2.8). Therefore, since we're only working with Gaussian states, we may think of the space of states as the space of (positive) quadratic forms. In other words, the states under consideration are all of the form and thus we may think of the relevant space of states as the three-dimensional space of 2×2 positive symmetric matrices A, with A ab = A ba , detA > 0, and A 11 , A 22 > 0. 8 In particular, the reference and target states become, respectively, where ω 1 , ω 2 and β are given by eq. (2.9). We now translate the scaling and entangling gates to this matrix representation. That is, we build a representation of these operators as 2×2 matrices which act on the symmetric matrices A. In particular, one finds that the gate matrices act as where (3.10) In this notation, [M ab ] cd is a 2×2 matrix, where c and d denote row and column indices, respectively. 9 Explicitly, we shall denote the basis of generators M I as With this new matrix formulation of our problem, we readily observe that the action of the gates Q ij -or more generally, circuits constructed from Q ij -on the vector (x 1 , x 2 ) T produces a vector whose elements are linear combinations of x 1 and x 2 . Furthermore, since the gates are invertible, this is precisely the definition of the group of transformations GL(2, R). 10 Thus our circuits form a representation of GL(2, R), i.e., the U (s) are trajectories in the space of GL(2, R) transformations. Now, in this matrix formulation, the path-ordered exponentials in eq. (3.1) are replaced by where M I are the generators given in eq. (3.11). The advantage of this formulation is that eq. (3.6) becomes These positivity constraints ensure that both eigenvalues of A ab are positive. 9 A quick way to construct these matrices is to consider the action of O ab on the column vector (x1, x2) T , and then build the matrix M T ab which yields the same result. One can verify that the commutators of the M ab match those of the O ab . Note that, while the action of the Q ab in eq. (3.2) leaves the wave functions properly normalized at each step, we lose track of this normalization when working with the A ab . 10 Note that one can also see the emergence of this group by observing that the algebra of the original operator generators O ab in eq. (3.2) close to form the algebra gl(2, R).
That is, we now have a simple expression which yields the components of the velocity vector Y I (s). Before we can utilize this expression however, we must explicitly construct a parametrization of the GL(2, R) transformations. We proceed with this task in the next subsection, but first let us make a few comments. Our task will be to find the shortest geodesic in some right-invariant metric on GL(2, R) that connects the initial and final states, A R and A T , as in eq. (3.12). We emphasize shortest geodesic because in fact, we will find that there is a continuous family of geodesics connecting the desired states. This non-uniqueness arises because our space of circuits is four-dimensional (since dim (GL(2, R)) = 4) whereas our space of states is only three-dimensional (since the 2 × 2 matrices A ij are symmetric). As a result of this mismatch, we should expect to find a one-parameter family of geodesics U (s) which yield the desired transformation A T = U (s = 1) A R U T (s = 1). However, as we have explained, the complexity is defined as the cost of the minimal or optimal circuit that obtains the specified target state. Hence this one-parameter family of solutions is merely the set of all possible circuits within this class. To find the optimal circuit, we simply need to find the geodesic within this family with the shortest length (3.5).
Since our ultimate aim will be to return to free scalar field theory, we note in passing that the notation introduced in the last two subsections generalizes very easily from two coupled oscillators to N coupled oscillators. We would then build a right-invariant metric on GL(N, R). Furthermore, note that the dimension of the space of circuits becomes N 2 , while the dimension of the space of Gaussian states or quadratic forms is only N (N + 1)/2. Hence the non-uniqueness involved in finding the most efficient circuit U (s) which produces the desired transformation grows quickly. We shall discuss the extension to a lattice of oscillators in section 4.

Geodesics on circuit space
To proceed with constructing the desired geodesics, we must choose an explicit parametrization of a general element U ∈ GL(2, R) = R × SL(2, R). Let us first consider U ∈ SL(2, R), which can be written as We recognize the constraint imposing det U = 1 as the embedding of (Lorentzian) AdS 3 in R 2,2 . Indeed, the appearance of AdS 3 could have been anticipated since the latter is the universal cover of SL(2, R). Our familiarity with this embedding then motivates the following choice of coordinates: x 0 = cos τ cosh ρ , x 1 = sin τ cosh ρ , x 2 = cos θ sinh ρ , x 3 = sin θ sinh ρ , (3.15) where τ, ρ and θ are the usual time, radius, and angle, respectively, of global coordinates on AdS 3 . We can easily extend this parametrization to U ∈ GL(2, R) = R × SL(2, R) by introducing an additional coordinate to parameterize the determinant of U , i.e., Hence we extend eq. (3.15) to x 0 = e y cos τ cosh ρ , x 1 = e y sin τ cosh ρ , x 2 = e y cos θ sinh ρ , x 3 = e y sin θ sinh ρ , (3.17) where, as before, τ, ρ, θ are coordinates on the SL(2, R) subgroup, and y parametrizes the R fibre. With these coordinates, we can express a general U ∈ GL(2, R) as U = e y cos τ cosh ρ − sin θ sinh ρ − sin τ cosh ρ + cos θ sinh ρ sin τ cosh ρ + cos θ sinh ρ cos τ cosh ρ + sin θ sinh ρ . (3.18) We are now equipped to construct the geometry implicit in the cost function (3.5), where the velocity components are given by eq. (3.13). As mentioned above, we begin by choosing G IJ = δ IJ , which assigns an equal cost or weight to every gate. This choice then defines the following right-invariant metric: (3.19) For later use, it is also convenient to express this in the form ds 2 = 2dy 2 + 2dρ 2 + 2dx 2 + 2 cosh(4ρ) dz 2 − 4 cosh(2ρ) dx dz , (3.20) where we have defined the pseudo-lightcone coordinates Note that our metric (3.19) is Euclidean, as is appropriate for defining a cost function, and so does not contain the (Lorentzian) AdS 3 geometry noted above. Indeed, a Lorentzian signature would not be suitable for the problem at hand, since certain directions would then carry negative or zero cost. We discuss the relation between our geometry and that of AdS 3 in appendix B. With the geometry in hand, we now wish to find the geodesics, and thereby the optimal circuit. Inspecting the metric (3.19), we can see three obvious Killing vectors: ∂ y , ∂ τ , ∂ θ . However, the metric is right-invariant by construction, meaning eq. (3.19) remains unchanged if we right-multiply U (s) by a constant GL(2, R) transformation. Therefore there must be one Killing vector for each generator of GL(2, R), namely, four. 11 In fact, it turns out that choosing G IJ = δ IJ results in an extra "accidental" symmetry, and so the metric above has a total of five Killing vectors. These Killing vectors (k I ) i ∂ i are explicitly constructed in appendix B, and are given in eqs. (B.8) and (B.9).
Of course, the existence of five Killing vectors implies an equal number of conserved momenta, c I ≡ (k I ) i g ijẋ j , which we will use to solve for the geodesics. Given the Killing vectors in eqs. (B.8) and (B.9), it is straightforward to evaluate the corresponding conserved quantities: where the dot denotes differentiation with respect to some affine parameter s along the geodesic. We are free to choose this parameter such that the normalization of the tangent vector is constrained to be constant, i.e., In a GR calculation, we would typically choose the normalization (for a spatial geodesic) to be +1, but this choice would leave the final value of s at the end of the circuit undetermined. However, recall that our notation for the path-ordered exponentials above is such that the circuits run over 0 ≤ s ≤ 1, cf. (3.1). Hence we shall scale the affine parameter s to lie in this range. The normalization constant k then gives the length of the geodesic, i.e., the depth of the corresponding circuit, since from eq. (3.5) we have The minimum value of k is then the depth of the optimal circuit, and by extension, the complexity of the target state ψ T . Next, we must establish the boundary conditions for our geodesics. The geodesics (and paths in the circuit geometry in general) are described by x(s) = {τ (s), ρ(s), θ(s), y(s)}. Now, our initial condition is that U = 1 at s = 0, and by comparing with the parametrization in eq. (3.18), we find that all coordinates except θ are initially zero, i.e., Note that the fact that θ = θ 0 is undetermined is not surprising since this is an angular coordinate, but the geodesic starts at the origin ρ = 0. Hence the freedom to specify θ 0 is the freedom that the geodesic leave the origin in any direction. In part, this freedom reflects the fact that we do not expect the boundary conditions to uniquely fix the geodesic, but to instead give rise to a one-parameter family thereof-see the discussion at the end of the previous subsection. Now, the end-point of the geodesic is determined by A T = U (s = 1) A R U T (s = 1), as in eq. (3.12), where the quadratic forms for the reference and target states are given in eq. (3.8). Substituting the initial state A R = ω 0 1 and the explicit representation of the unitaries (3.18) into this relation, we have (3.26) where the subscript 1 denotes the value of the coordinate at s = 1, e.g., y 1 = y(s = 1). Comparing the entries of the matrix on the right-hand side with those of A T in eq. (3.8), we arrive at the following boundary conditions for the end of the geodesic: β/ω 0 = e 2y 1 cos(θ 1 + τ 1 ) sinh(2ρ 1 ) .

(3.27)
Implicitly, these constraints allow us to identify the final coordinates x(s = 1) for the geodesics corresponding to circuits which produce the desired transformation. Explicitly, we may solve this system to obtain However, there is an obvious ambiguity here since θ 1 and τ 1 appear only in the combination θ 1 + τ 1 . Since only this linear combination is fixed by eq. (3.28), we have a one-parameter family of final boundary conditions-the linear combination θ 1 − τ 1 remains unspecified. Naïvely, this might lead one to suspect a two-parameter family of allowed solutions, since the initial conditions left θ 0 unfixed as well. But this is not the case: rather, the geodesic equations of motion relate the freedom in the boundary conditions at s = 0 and s = 1, and the freedom in the initial and final conditions combine to yield the one-parameter family of geodesics anticipated above. This situation is illustrated in figure 2, which shows a one-parameter family of solutions beginning at the origin and ending on the spiral given by θ + τ = θ 1 + τ 1 and radius ρ = ρ 1 . To determine the complexity of the final state A T , we must find the minimum length geodesic within this family, and thereby the optimal circuit. Having specified the boundary conditions, we proceed to solve for the geodesics by examining the conserved momenta (3.22). The first of these gives the simplest constraint: c 1 = 2ẏ. Integrating with respect to the affine parameter s then yields: y(s) = c 1 s/2 + y 0 . In this case the undetermined coefficients are easily fixed by the boundary conditions, y(s = 1) = y 1 and y(s = 0) = 0, hence: c 1 = 2 y 1 and y 0 = 0 =⇒ y(s) = y 1 s .
(3.29) Figure 2: Sketch of the one-parameter family of geodesics. The vertical axis is τ , the horizontal plane is described by the radius ρ and the azimuthal angle θ, and the y direction is suppressed. The circuits which produce the transformation from AR to AT are described by geodesics running from the origin to the blue spiral at θ + τ = θ1 + τ1 and ρ = ρ1 (shown here for the special case θ1 + τ1 = π, which appears in eq. (3.48) below). The black curves represent (non-minimal) geodesics within the one-parameter family of solutions with different values of θ0. The minimum geodesic corresponds to the green line in the τ = 0 plane with ∆θ = 0 (i.e., θ0 = θ1), whose length is given by eq. (3.38).
Next, we consider c 4 and c 5 . These two constraints may be solved to obtaiṅ We then observe thatθ diverges at the origin ρ = 0 unless c 4 = −c 5 , which we must therefore impose in order to be compatible with the initial conditions. Implicitly, we are setting the angular momentum, i.e., the conserved momentum associated with the Killing vector ∂ θ , to zero, which is characteristic of geodesics passing through the (radial) origin ρ = 0. With this condition, the θ equation can be trivially integrated to yield θ = c 5 s + θ 0 , where we have already imposed θ(s = 0) = θ 0 . Imposing the final boundary condition then yields Furthermore, the above allows us to simplify theτ equation tȯ Now, combining our expressions forθ andτ with the constraints c 2 and c 3 in eq. (3.22), we find a relatively simple equation forρ: In principle, we should now solve for the general solutions of eqs. (3.32) and (3.33) subject to the boundary conditions in eqs. (3.25) and (3.28). While it is possible to carry out this exercise, the final solutions are not particularly illuminating. 12 Instead, let us point out the particularly simple solution that arises for ∆θ = 0. In this case the expressions forτ andρ reduce toτ  Note that an explicit path-ordering is not needed in the second expression since it is simply the exponential of a fixed matrix. 13 The circuit depth of U 0 , i.e., the length of the geodesic, is given by eqs. (3.23) and (3.24), which for this simple solution yields Again, in principle, we should determine all of the other geodesics satisfying the appropriate boundary conditions, and compare their respective circuit depths to D(U 0 ) in order to determine the minimum. However, we shall instead provide a more indirect but less technically challenging proof that this simple straight-line solution is in fact the shortest possible geodesic, and hence that it describes the optimal circuit. 12 The general solution for eq. (3.33) is given by where c 2 = c 2 2 + c 2 3 is fixed by substituting the boundary condition ρ = ρ1 at s = 1 into this equation. Furthermore, given this result, it is possible to integrate eq. (3.32) to obtain τ (s); one finds Substituting τ = τ1 at s = 1 into this expression fixes τ1 in terms of ∆θ and c 2 . Combining this result with the boundary condition for θ1 + τ1 in eq. (3.28), we can then determine θ1. In turn, θ0 is now fixed since we know ∆θ and θ1. 13 For this simple case, it is straightforward to identify the exponential form in the second line of eq. (3.37) given the expression appearing in the first. In general however, one would apply eq. (3.13) to identify the components of Y I (s) and then substitute these into eq. (3.12).
To prove that the straight-line solution above is the geodesic whose length is the (global) minimum, we recall from eq. (3.23) that the length of any geodesic is given by the normalization constant k. Now into this expression, we substitute our general solutions for y(s) and θ(s) from eqs. (3.29) and (3.31), respectively, as well as the expression forτ from eq. (3.32), whereupon we find This equation holds point-by-point along any geodesic satisfying the appropriate boundary conditions, but what we would like to argue (without explicitly solving for ρ(s)) is that k 2 is minimized by choosing ∆θ = 0.
To begin, consider motion in the ρ-direction along any of our geodesics. The average velocity is given by and hence we may conclude that 1 0 dsρ 2 ≥ ρ 2 1 , and that this inequality is only saturated whenρ = ρ 1 along the entire geodesic. Now, examining the coefficient of ∆θ 2 in eqn. (3.39), we have where the lower inequality is only saturated at ρ = 0, and the upper inequality is saturated at ρ → ∞. 14 Given that all of our geodesics must start at ρ = 0 and end at ρ = ρ 1 , upon averaging over any of these geodesics, we find Finally, let us average eq. (3.39) over any of our geodesics: Comparing this to eq. (3.38), we have established the inequality k ≥ D(U 0 ). Furthermore, our argument has established that this inequality can only be saturated withρ(s) = ρ 1 (i.e., ρ = ρ 1 s) and ∆θ = 0 (i.e., θ(s) = θ 0 , which implies τ (s) = 0 via eq. (3.32)). We have therefore proved that the simple straight-line geodesic indeed constitutes the global minimum for our cost function (3.24), and hence that eq. (3.38) is in fact the complexity of the Gaussian wave function in this framework: As an exercise, we can compare the above result for D(U 0 ) in eq. (3.38), which was evaluated using eq. (3.24), with the result found by evaluating eq. (3.5). In this case, we must identify the components Y I (s), which is easily done by examining the exponential expression in eq. (3.37): 15 Since these components are all constant, the integral over s in eq. (3.5) is trivial, and the circuit depth (with G IJ = δ IJ ) reduces to in agreement with eq. (3.38).

Normal-mode subspace
To properly interpret the complexity, we must re-express our result (3.38) in terms of the physical parameters of the two coupled oscillators (2.4), as well as the frequency ω 0 in the reference state (2.9). However, one finds that the complexity is most elegantly described in terms of the normal-mode frequenciesω + andω − given in eqs. (2.5) and (2.6). Using eq. (2.9), the final boundary conditions (3.28) simplify to Substituting these expressions for y 1 and ρ 1 into eq. (3.38) then yields the complexity of the ground state, At this point, let us also note that the boundary condition θ 1 + τ 1 = π (along with ∆θ = 0 and τ (s) = 0) implies that the initial angle is θ 0 = π. This straight-line geodesic is illustrated by the green line in figure 2. The corresponding circuit (3.37) simplifies to with y 1 and ρ 1 given by eq. (3.48). 15 Recall that our GL(2, R) generators are given in eq. (3.11).
The simple and elegant form (3.49) of the complexity in terms of the normal-mode frequencies suggests that we should investigate the optimal circuit (3.37) in terms of the normal modes. The relationship between the physical positions of the masses and the normal-mode coordinates was given in eq. (2.6), but we can understand this change of coordinates in terms of a simple rotation. In particular, we can perform the coordinate transformation via the orthogonal rotation matrix R, 16 Introducing the short-hand notation x = (x 1 , x 2 ) T andx = (x + ,x − ) T , the transformation (3.51) may be concisely writtenx = R x, and the inverse transformation becomes x = R Tx . Of course, we can also use this transformation to re-express the target Gaussian wave function in terms of the normal-mode coordinates, whereÃ T denotes the quadratic form describing the ground state in the normal-mode space. Explicitly performing this rotation, one finds That is, the target state becomes a factorized Gaussian in the normal-mode basis, cf. eq. (2.5). Of course, this decoupling was the essential point of introducing the normal-mode coordinates in the first place. Furthermore, if we apply this transformation to the reference state in eq. (3.8), we see that it retains its simple form, i.e., That is, the reference state remains a factorized Gaussian when written in terms of the normal modes. Now, given the action of the gates and circuits on the quadratic forms, cf. eq. (3.12), we can transform our minimal circuit (3.50) to act in the normal-mode space: (3.55) 16 Our transformation matrix certainly satisfies R R T = R T R = 1. However, with the conventions adopted above, we note that detR = −1 and as a result, we actually have that as a numerical matrix R is symmetric, as shown with the eq. (3.51). However, we still distinguish R and R T in the following because R provides a mapping from the physical positions to the normal coordinates, while R −1 = R T provides the inverse mapping. In other words, the columns of R are labeled 1,2 while the rows are labeled +,-and vice versa for R T .
This transformation effects a remarkable simplification of the circuit (3.50) tõ where in the second line we have used eq. (3.48).
The important lesson learned here is as follows: from the perspective of the normal modes, both the target state and the reference state are factorized Gaussians, as shown in eqs. (3.53) and (3.54). The optimal circuitŨ 0 (s) then simply acts in a diagonal fashion to "amplify" each of the diagonal entries in the corresponding quadratic forms, taking ω 0 toω ± in a simple linear manner. It is rather intuitive that this should be the optimal way to prepareÃ T from A R , since if any off-diagonal entries (i.e., entanglement) were introduced along the circuit, they would simply have to be removed by the time the trajectory reaches its end-point. This feature of the optimal circuit will greatly simplify our considerations of a lattice of coupled oscillators in the next section.
Before turning to this generalization however, we wish to emphasize that the original circuit (3.37) is performing the same operation of amplifying the normal modes-this is simply a matter of re-expressing U 0 in an alternative basis of generators. To properly clarify this, we need to introduce some additional notation. In the above, we adopted a tilde to denote various quantities in the normal-modes basis. 17 We also introduced the index notation I = {11, 22, 12, 21} to label the components of the velocity Y I (s) and the generators M I . Here we would like to combine these two conventions to introduce a new index labelĨ = {++, +−, −+, −−} to denote the same objects with components acting in the normal-mode basis. Thus the natural basis of generatorsMĨ with which to construct the circuits acting on the states described in the normal-mode basis arẽ As numerical matrices, theseMĨ are of course identical to the M I given in eq. (3.11), but the two sets of generators act in different spaces. Via the transformation (3.51), we can also transform these generators to act on the states in the original position basis, i.e., MĨ = R TMĨ R: The action of these generators can be read off from the indices, e.g., M ++ scales the x + coordinate or amplifies the corresponding normal mode. Of course, we could also transform the original generators M I in eq. (3.11) withM I = R M I R T to construct the corresponding normal-mode basis. For example,M 11 would still scale the x 1 coordinate but would act on states in the normal-mode basis, i.e., it acts on Gaussian wave functions written in terms of x ± . With this new notation in hand, we would like to express our optimal circuit U 0 in terms of the generators MĨ . It is easily shown, either by examining eq. (3.50) directly or by transforming the expression in eq. (3.56) with U 0 (s) = R TŨ 0 (s) R, that the optimal circuit can be expressed as where M ±± are the linear combinations of the original generators given in eq. (3.58). In this form, we again recognize that the optimal circuit is simply amplifying the two normal modes, without introducing (and then having to remove) any entanglement between x ± . We can also observe that this simple circuit only involves two commuting generators, M ++ and M −− . Since the generators commute, it is straightforward to show that the geometry of corresponding normal-mode subspace is flat. That is, if we consider general circuits of the form then the corresponding metric becomes 18 Hence we recognize the normal-mode subspace as precisely the (θ, τ ) = (π, 0) plane in our extended geometry (3.20). 19 This perspective also makes clear why the optimal geodesic remains in the normal-mode subspace. Examining the full metric (3.20), it is clear that motion in the θ and τ directions only extends the length of the trajectory. Thus since the start and end points both lie in this plane, there is no advantage to be gained by moving out of the normal-mode subspace. This argument also relies on the fact that g yy and g ρρ in the full metric (3.20) are constants, independent of θ and τ , which precludes the existence of "short-cuts" to be found by moving off the normal-mode subspace (we return to this point in section 5). This is another important feature that extends to the case of a lattice of coupled oscillators in the next section.
To close this section, we wish to introduce some additional technology which will prove useful in those that follow. Thus far, we have two particularly useful sets of generators for our gates and circuits, namely, M I and MĨ given in eqs. (3.11) and (3.57), respectively. While these generators all act on states and circuits in the physical basis, M I acts to scale or entangle the physical positions x 1,2 , while MĨ scales or entangles the normal-mode coordinates x ± . The transformation between the two bases is given in eq. (3.58), but we would like to build an explicit transformation matrixR: Note that in the final equality, R is the rotation matrix in eq. (3.51), and we are identifying the indices as follows:Ĩ = (k ) with k, ∈ {+, −}, and J = (ab) with a, b ∈ {1, 2}. 20 This identification is really the origin of the interesting tensor product structure R = R ⊗ R.
The expression in eq. (3.62) indicates that the first (second) R is rotating the first (second) component of the pairs which comprise theĨ and J indices on the two generators. Given this expression, we immediately see that R is also an orthogonal rotation matrix. Hence we can easily invert the transformation between the basis generators via M I = ( R T ) IJ MJ = RJ I MJ . Similarly, this transformation acts on the velocity components as Y I = YJ RJ I . These transformations will prove useful in examining the complexity with cost functions written in different bases. For example, in the present context, we can see that the cost function remains unchanged if we express it directly in the normal-mode basis. We can also transform the metric (3.19) as follows: where we have used the fact that R is an orthogonal matrix. In going from the second to third line, we have used the identity RĨ I δ IJ ( R T ) JJ = δĨJ . Note that the invariance of the metric under this change of basis was already used in evaluating the metric on the normalmode subspace in eq. (3.61). We extend this discussion of changing between the position and normal-mode bases to the case of a linear lattice of N oscillators in appendix D.

A lattice of oscillators
In this section, we wish to return to the original problem of a free scalar field regulated by a lattice, cf. (2.2). That is, we will consider evaluating the complexity of the ground state of a lattice of coupled oscillators (2.3). Drawing on our experience with the two coupled oscillators, this becomes a straightforward calculation. In particular, as we saw above, both the ground state and the reference state are described by factorized Gaussians in the normalmode space. And in this space, the optimal circuit simply amplifies each of the diagonal entries in the corresponding quadratic forms in a linear manner. To simplify the technicalities in the following discussion, we will explicitly consider the case of a one-dimensional lattice, and discuss more general dimensions in the next subsection. Hence, we begin with N oscillators on a one-dimensional circular lattice, with periodic boundary conditions x a+N = x a . 21 As in the two oscillator problem, we have set the masses M a = 1 for simplicity but we should think of the frequencies as being related to the field theory parameters by ω = m and Ω = 1/δ, as in eq. (2.3). The Hamiltonian (4.1) then corresponds to the lattice version of a (one-dimensional) free scalar field on a circle of length L = N δ. Of course, to solve the above system, one simply rewrites the Hamiltonian in terms of the normal modes, where the transformation to the normal-mode basis is achieved by a (discrete) Fourier trans- where k ∈ {0, . . . , N −1}, and we note thatx † k =x N −k . 22 The normal-mode frequenciesω k are defined in terms of the physical frequencies ω and Ω in the Hamiltonian (4.1) as follows: (see appendix C). As desired, eq. (4.2) reduces the problem to N decoupled harmonic oscillators, which enables us to easily write the ground-state wave function as As before, this ground state will be the target state in our complexity computations. While eq. (4.7) will suffice to describe the ground state, in principle, one would also like to express the wave function in terms of the original variables x a in the position basis. This transformation is facilitated using notation introduced in section 3.2. In particular, following eq. (3.51), we write the Fourier transformation (4.3) between the position and normal-mode bases asx = R N x, with (4.8) 22 We can see this result as a combination of two simpler identities:x † k =x −k , which follows from the complex conjugation of eq. (4.3), andx k =x k+N , which follows from the periodicity of the lattice. Note that our convention for the range of k was chosen to match the range of the position labels a, rather than shifting the range of k to run over positive and negative values, i.e., k ∈ {− N/2 + 1, − N/2 + 2, , · · · , N/2 }, which is a more typical convention. Furthermore, for future reference, note that we can define u k ≡ [u k ]a = exp (−2πi k a/N ) as the orthogonal basis of an N -dimensional vector space, satisfying the normalization condition Hence we use the usual definition for the normal-mode momentã where µ ≡ exp (−2πi/N ). 23 Since R N is a unitary matrix, i.e., R N † R N = 1, the inverse transformation is given by x = R N †x . Now let us adopt the notation of section 3 (and in particular, of eq. (3.7)) to write the target state (4.7) as Using the rotation (4.8), we can write this target state in terms of the physical coordinates, 24 We are now prepared to extend our complexity calculations to this lattice of coupled oscillators. We have already identified the target state as the ground state (4.7). In analogy with eq. (2.10), the reference state will be a factorized Gaussian state, where the individual oscillators are completely unentangled. 25 An important feature of our reference state is that it is invariant under translations on the lattice, i.e., the Gaussian of each oscillator has the same width ω 0 . As a result, it remains a factorized Gaussian when expressed in terms of the normal-mode coordinates: Lastly, we need to consider the elementary gates with which we will build the circuit U that implements the desired transformation ψ T = U ψ R . With the notation introduced in eq. (2.11), the set of gates (particularly the entangling and scaling gates) is easily enlarged for the present problem by simply extending the range of the indices: a, b ∈ {1, 2} −→ a, b ∈ {0, 1, 2, · · · , N − 1}. These discrete gates are then easily extended to the path-ordered exponentials introduced in eqs. 2). In discussing the target and reference states with the notation of eq. (3.7), we also anticipated mapping these exponentials to the matrix formulation introduced 23 As discussed in footnote 16 for the matrix R in eq. (3.51), we distinguish RN from R T N even though the numerical matrix in eq. (4.8) is symmetric. Note that if we write out the transformation to show the indices, we havex k = [RN ] ka xa. That is, the row index of RN has values in the momenta k while the column index has values in the lattice position a. In passing, we also observe that eq. (4.8) reduces to eq. (3.51) for the special case N = 2, for which we have µ = exp (−iπ) = −1. 24 The relationω k =ω N −k ensures that AT is real. 25 Recall our tilde notation to distinguish the normal-mode space from the physical space. In particular, the reference frequency ω0 is independent of the normal-mode frequency with k = 0, i.e., ω0 =ω0! in eqs. (3.9), (3.10) and (3.12) for Gaussian states. In fact, the generators have precisely the form given in eq. (3.10), where again the indices run over the range a, b, c, d ∈ {0, 1, 2, · · · , N − 1}. That is, we now have N 2 generators which are N ×N matrices. This extends the GL(2, R) group found in section 3 to the group GL(N, R) in the present problem. Following the analysis in section 3, we use the analogous F 2 cost function, i.e., 26 Hence the optimal circuit will correspond to a geodesic in the GL(N, R) geometry given by a right-invariant metric, analogous to eq. (3.19). To simplify the discussion of the metric here (and in the next section), we introduce the following notation: However, extending the detailed calculations above to the full N 2 -dimensional geometry would be very involved. In particular, the next step would require finding the analog of eq. (3.18), i.e., a convenient parametrization of a general group element U ∈ GL(N, R), which would naturally involve N 2 coordinates. Thus at this point, we rely on the lessons learned from the case of two coupled oscillators in the previous section.
In particular, there we found that since both the ground state and the reference state are described by factorized Gaussians in the normal-mode basis, the optimal circuit simply acts to amplify each of the diagonal entries in the corresponding quadratic forms in a simple linear manner. We have already noted by way of eqs. (4.9) and (4.12) that the former statement about factorized Gaussians also applies in our lattice problem. Hence it is natural that the most efficient circuit simply amplifies the Gaussian width for each of the normal-mode coordinates, i.e., ω 0 →ω k . In particular, the circuit does not introduce any entanglement between the normal modes at any stage, since this entanglement would have to be removed before arriving at the final target state (4.9). Via eq. (3.55), let us write the optimal circuit acting in the normal-mode basis; we have and thus the straight-line circuitŨ 0 (s) becomes This circuit certainly accomplishes the desired transformation withŨ 0 (s = 1) = exp M 0 , but the intuition from the previous analysis of two coupled oscillators suggests that it is also the optimal circuit. We can add to this intuitive picture as follows: in the discussion around eqs. (3.60) and (4.14), we identified the normal-mode subspace as consisting of those circuits U which only involve the scaling generators for the normal modes. Consequently, it is straightforward to show that the geometry of the normal-mode subspace is flat since these generators all commute with one another. In the present case, the normal-mode subspace becomes a N -dimensional Substituting this expression into eq. (4.14), one finds the following flat Cartesian metric induced on this subspace: Therefore any geodesic within the normal-mode subspace will simply take the form of a straight line. It is then straightforward to show that if we confine the circuit to this normalmode subspace (4.18), the optimal circuit is described by the simple circuit in eq. (4.15), which we write as whereŨ 0 (s) andM 0 are defined in eq. (4.16).
There are actually some subtleties in the preceding argument which make the conclusion somewhat premature. The first is that eq. (4.14) which defines the metric is written in the position basis, whereas eq. (4.18) was implicitly calculated for an expression (4.17) written in the normal-mode basis. That is, in eq. (4.17), we worked withM n-m =ỸĨMĨ with a particular choice ofỸĨ . 28 However, we show in appendix D that this was nonetheless a valid approach since the metric takes precisely the same form when written in terms of the normalmode space. This requires extending the discussion around eq. (3.62) describing the change of bases for the case of two couple oscillators to the analogous transformation for our linear lattice of N oscillators.
Secondly, to properly establish that the optimal circuit follows a straight line in the normal-mode subspace, as in eq. (4.19), we must show that no shorter path can be found by making an excursion outside this subspace. To begin, we note that implicitly we assumed in eq. (4.17) that all of the other coordinates in the GL(N, R) geometry could be set to zero. Recall that in the GL(2, R) metric (3.19), the metric on the normal-mode subspace was 27 In general, the coordinatesỹ k are complex but satisfy the normal-mode "reality condition"ỹ † k =ỹ N −k . 28 Again, the tilde on the indexĨ indicates that it runs over pairs of momentum labels (k ), while the tilde on M indicates that these generators act on Gaussian wave functions written with the normal-mode coordinates completely independent of the other coordinates, i.e., we had ds 2 n-m = 2dy 2 +2dρ 2 irrespective of the values of θ and τ . In particular, recall that the optimal circuit was a straight line in this subspace with θ = π and τ = 0.
We would like to establish a similar result for the present N 2 -dimensional geometry. For simplicity, we will work in the normal-mode space. We proceed by expressing general circuits U using the Iwasawa (or KAN) decomposition of GL(N, R); see for example [50]. This states that anyŨ ∈ GL(N, R) can be uniquely written as the product of three matrices,Ũ = K A N , where K is an orthogonal matrix, A is a diagonal matrix with positive entries, 29 and N is an upper triangular matrix with every diagonal element equal to 1. Clearly, we are interested in the A component as this describes the normal-mode subspace, as in eq. (4.17).
As a warm up exercise, let us consider translatingŨ n-m by some fixed angles and shifts. In particular, we writeŨ = K 0Ũn-m N 0 where only theỹ k inŨ n-m vary (cf. eq. (4.17)) and ask what is the metric on the corresponding subspace. Since N 0 acts on the right, and the metric is right-invariant by construction, it has no effect on the geometry. Our experience in changing bases in appendix D allows us the eliminate the K 0 rotation as well: following eq. (D.7), we write the differentials dỸĨ = tr dŨŨM † I as Now in the last factor, K 0 acts by a similarity transformation on the generators which effectively produces a change of basis. However, using the special form of the generators (D.3), it is straightforward to show -following a series of steps analogous to those given in eqs. (D.4) or (D.6) -that It follows that K 0 is an orthogonal matrix since K 0 is orthogonal, and hence this rotation of the generator basis leaves the metric unchanged. Therefore we find that the induced metric on this subspace is which again yields the simple answer given in eq. (4.18). This result establishes that there are indeed no short-cuts to be found by running the circuit through the angle and shift directions. That is, the circuit must run fromỹ k = 0 toỹ k = 1 2 log(ω k /ω 0 ) as in eq. (4.16). Eq. (4.22) further establishes that there will be a fixed distance or cost associated with this displacement, irrespective of the orientation of the normal-mode subspace in the full geometry, that is, irrespective of the angles and shifts chosen 29 We denote this diagonal matrix with the traditional A, but it should not be confused with the quadratic forms specifying the Gaussian states, cf. eq. (3.7). Similarly, K here should not be confused with the gates producing a momentum shift in eq. (2.11), nor should N be confused with the total number of oscillators. We trust that these distinctions will be clear from context. in K 0 and N 0 . Since the full geometry is Euclidean, moving in these "orientation directions" will only add to the distance. Thus the best strategy is to fix the shifts and angles at the beginning of the circuit (to zero, as required by U (s = 0) = 1) and then move only in the normal-mode directions.
This argument is still not quite sufficient to establish that the simple straight-line circuit is a geodesic in the full N 2 -dimensional geometry. In particular, non-vanishing off-diagonal terms in the metric which mixỹ k with the other coordinates would force the geodesic to move away from the normal-mode subspace in the additional angle and shift directions. But since evaluating the full metric would require a rather lengthy and involved calculation, we instead consider small deviations of the circuits around the subspace specified byŨ = K 0Ũn-m N 0 , i.e., we extend our initial ansatz to allow small excursions in the K and N directions, where θĨ , ηĨ 1. Here, the (small) change in K only involves the (antisymmetric) rotation generators while the (small) change in N only involves the shift generators The rotation generators are, of course, a linear combination of the original generators given in eq. (D.3), and hence are not orthogonal to the shift generators in the sense that Of course, all of these generators are orthogonal to the diagonal generators appearing inŨ n-m , which will become the key point momentarily. With our extended circuits (4.23), we now evaluate the differentials dỸĨ = tr dŨŨM † I on the normal-mode subspace, i.e., at θĨ = 0 = ηĨ , where K 0 is the orthogonal matrix given in eq. (4.21) and dỸĨ n-m are the differentials along the normal-mode directions identified in eq. (4.22). As before, the rotation of the differentials by K 0 can be ignored since this transformation leaves δ IJ in the metric unchanged. Next we observe that the only non-vanishing components of dỸJ n-m are along the diagonal directions, i.e.,J = (kk). It is then easy to show that the other two differentials are orthogonal to these. Given the explicit form of the rotation generators in eq. (4.24), it is clear that the second term only contributes in the off-diagonal directions, i.e.,J = (k ) with k = . Similarly, one can show that the same is true of the third term via eqs. (4.17) and (4.25), The key point then is that dỸJ n-m are orthogonal to the other two differentials in eq. (4.27).
In fact, it is straightforward to show that the full metric on the normal-mode subspace (i.e., θĨ = 0 = ηĨ ) becomes Hence there are no off-diagonal terms in the metric, which would drive the geodesic away from the normal-mode subspace. We may therefore conclude that the optimal circuit indeed takes the form of the simple straight-line circuit in eq. (4.16). Thus, for the cost function (4.13), the complexity for our lattice of oscillators is obtained by simply summing up the circuit elements in the normal-mode basis. Using eqs. (4.16) and (4.17), we find where the normal-mode frequencies are given in eq. (4.6). Recall that in our lattice regularization (2.2), we had ω = m and Ω = 1/δ, and so we can express the complexity (4.30) in terms of the field theory parameters viã Furthermore, we can replace N = L/δ where L is the total length of the one-dimensional lattice of oscillators. Of course, ω 0 remains the (as yet unspecified) frequency which specifies the Gaussian reference state (4.11). The entire discussion in this section is easily extended (albeit with a somewhat tedious extension of the notation) to the evaluation of the complexity of a (d−1)-dimensional spatial lattice of N d−1 oscillators, and the final result is where k i are the components of the momentum vector k = (k 1 , k 2 , · · · , k d−1 ), and the normalmode frequencies are given byω The linear size of each spatial direction here is L = N δ, and so the total (spatial) volume of the system is V Hence the total number of oscillators can be expressed as which will prove useful below.

Comparison with holography
Eq. (4.32) gives our result for the complexity of the ground state of a free scalar field in d spacetime dimensions. We would now like to compare this result with the analogous results arising from the proposals for holographic complexity discussed in the introduction. Of course, we must note that we are trying to compare complexities for disparate QFTs, i.e., a free theory with a single degree of freedom in the present case versus a strongly coupled theory with a large number of degrees of freedom in holography. Hence there is no a priori reason to expect that the results should agree in the two cases. Nevertheless, we will find that with certain choices, our QFT calculations share a number of qualitative features with holographic complexity. We can interpret these similarities as providing guidance towards understanding the cost function that underlies the holographic complexity conjectures. Examining eq. (4.32), we see that the expression under the square root essentially involves an integration over the spatial momenta. Our experience with QFT thus suggests that the result will be dominated by the UV modes, i.e., by modes withω k ∼ 1/δ. Hence as an approximation which allows us to identify the leading contribution to the complexity, we may replace all of theω k with 1/δ in eq. (4.32) to obtain 30 where we have used eq. (4.34) to re-express the leading power of N in terms of V /δ d−1 .
The leading UV divergence in holographic complexity for both the CA and CV proposals was studied in some detail in [19]. Hence we can compare our QFT result (4.35) with the analogous results for holographic complexity; denoting the latter collectively as C holo , these were found to take the form C holo ∼ V δ d−1 . (4.36) Thus we see that the leading terms in the QFT and holographic complexities differ by the power of 1/2 appearing in eq. (4.35). However, the origin of this square root is clear: it is simply the overall square root appearing in eq. (4.32), which in turn arises from our use of the F 2 cost function in eq. (4.13). Now, there is nothing wrong with the result in eq. (4.35) per se, but it does suggest that if our QFT complexity is to emulate the leading behaviour found in holographic complexity, then we should make an alternative choice for the cost function. 31 For example, the cost function defined by the F 1 measure in eq. (1.9), which involves the first power of a single sum over the modes, would produce the desired behaviour. More generally, a natural family of cost functions which would reproduce the divergence in eq. (4.36) is where the natural choice would be that κ is a positive integer, but any positive real value (with κ ≥ 1) will suffice for most of the following discussion. Note that we have defined the cost function here with the sum running over the normal-mode basis-we return to this point below in section D.1.
Here, we should note that only κ = 1 (equivalently, the F 1 cost function), satisfies the condition of positive homogeneity as described in the introduction; that is, for general κ > 1, doubling the amplitude of YĨ does not double the cost. This issue can be described as saying that only the case κ = 1 yields a reparametrization-invariant cost function, i.e., that replacing s →ŝ(s) leaves the cost unchanged. However, we may proceed with the physics intuition that we can think of D κ as different kinds of actions describing the motion of a particle in the space of circuits.
Now it is relatively straightforward to show that in fact the straight-line circuit in eq. (4.16) minimizes all of these cost functions. This circuit only acts with scaling gates on the various normal modes. It is clear that if the circuit were to make an excursion away from the normal-mode subspace, i.e., if the path also moved in the entangling directions, this would only turn on new components of the velocity YĨ and thereby increase the cost of the circuit. To establish that the straight-line path is favoured by the general κ cost function, let us consider a more general trajectory in the normal-mode subspace,Ũ 1 (s) = exp M 1 (s) , whereM  ds |∂ s f k (s)| κ ≥ 1 (for κ ≥ 1), and that the inequality is saturated if and only if ∂ s f k (s) = 1, i.e., f k (s) = s. 32 That is, the general κ cost function (4.37) is minimized by the straight-line circuitŨ 0 (s) = exp M 0 s . Furthermore, working with the UV approximationω k = 1/δ, the leading contribution to the complexity then becomes (4.41) An interesting feature of this result is that in limit δ → 0, this contribution appears to diverge faster than the power law 1/δ d−1 in the first factor.
This last observation, however, depends on the choice of ω 0 which defines our reference state (4.11), which we have hitherto left unspecified. Considering this choice, there seem to be a number of reasonable options. First, ω 0 could be associated with some ultraviolet frequency at the lattice scale. For example, ω 0 = e −σ /δ where e −σ provides a numerical scale that ensures ω 0 >ω k for all k. In this case, the leading contribution in eq. (4.41) reduces to With this choice, the extra logarithmic factors in the δ → 0 divergence have been eliminated and we are only left with the 1/δ d−1 factor. However, this choice also entails the interesting feature that the (subleading) infrared contributions to the complexity will involve the UV cut-off scale. That is, the full sum over momenta in eq. (4.37) includes summing over the infrared modes, i.e., modes withω k ∼ m. These infrared contributions will take the form An alternative choice would be to associate ω 0 with some infrared scale, i.e., ω 0 1/δ. One might choose this scale to be a physical scale in the problem, such as the mass m or the volume V , but this would tie the reference state to the properties of the QFT. 33 This appears problematic if we wish to compare the complexities of states in different theories, e.g., with different masses-see below. Hence it seems that we are instead led to choose some arbitrary IR scale to define the reference frequency, which then becomes a part of our definition of the complexity of QFT states. In this sense, the appearance of ω 0 here is not very different from the appearance of the arbitrary numerical factor σ in the complexity with the previous UV choice. Of course, if ω 0 is a fixed IR frequency, the additional logarithmic factor in the complexity (4.41) survives and contributes to the leading divergence in the limit δ → 0.
With the above in mind, we find surprising similarities when we compare our result with the CA proposal for holographic complexity. The leading divergence appearing in the latter (1.3) takes the form [19] where δ is the short-distance cut-off scale in the boundary CFT, L AdS is the AdS curvature scale of the bulk spacetime, and α is an arbitrary (dimensionless) coefficient which fixes the normalization of the null normals on the boundary of the WDW patch. Since C A is a quantity which is to be defined in the boundary CFT, it should not depend on the bulk AdS scale. However, we can eliminate this factor with the freedom in choosing α, i.e., we set α = ω 0 L AdS where ω 0 is some arbitrary frequency. In this case, eq. (4.44) reduces to 33 Furthermore, if we choose ω0 ∼ V −1/(d−1) , the complexity becomes superextensive.
While this choice eliminates the AdS scale, the holographic complexity still depends on the choice of ω 0 , just as in our QFT result (4.41). Furthermore, all of the issues discussed above with respect to this choice also appear in the case of the CA conjecture. In particular, we emphasize that with the choice ω 0 = e −σ /δ, the UV cut-off appears in infrared contributions to the holographic complexity arising from joint terms deep in the bulk [19,23]. Whereas in [19] this ambiguity and the associated issues seemed problematic for holographic complexity, here can view them as a natural feature of complexity for QFTs. We can go further in comparing our QFT results with holography. In particular, in order for the leading divergence in eq. (4.41) to match the holographic result (4.45) more closely, we should choose κ = 1 in eq. (4.37), i.e., the F 1 cost function. Of course, this reasoning only applies when the reference frequency is chosen in the IR and the logarithmic factor modifies the form of the leading divergence. In the case where ω 0 is set by the cut-off scale, the leading divergence is a simple power law and the exponent κ only modifies the overall numerical prefactor, which we have not specified here. However, κ also appears in the IR contribution in eq. (4.43). Again the analogous contributions in holography would be linear in log(δ) because of the form of the corresponding boundary terms in the gravitational action [17]. Hence this reasoning again favours the choice κ = 1. That said, given the aforementioned disparity between the field theories that we are comparing, it is not clear how much weight to give this observation.
One can also look at the form of subleading corrections to the leading divergence. As discussed in appendix E, the first subleading correction comes from the mass and has the form V m 2 /δ d−3 . This form could be anticipated by simple dimensional analysis, and analogous results can be found in holographic complexity as well [52]. Specifically, if the boundary CFT is perturbed by a relevant operator of dimension ∆, the corresponding coupling λ will have dimension d − ∆, and the first subleading correction to the holographic complexity then takes the form V λ 2 /δ 2∆−(d+1) . For the CV conjecture, these calculations follow in close parallel with the analogous calculations of corrections to holographic entanglement entropy induced by relevant operators [53]. However, we expect that analogous results will appear for the CA conjecture. Such corrections to holographic complexity were considered in [19] as arising from placing the boundary theory in a curved spacetime or from evaluating the state of a curved time slice. It may be interesting to understand how to extend our present QFT calculations of complexity to incorporate such situations.

Penalty factors
In eq. (3.5), the cost function was written with a general metric G IJ , which allows us the freedom to include penalty factors to weight certain directions or classes of gates more heavily than others. This is particularly relevant for the lattice of oscillators representing the regulated scalar field. In the previous section, our circuits implicitly included entangling gates Q ab , which coupled points on the lattice which were arbitrarily far apart. However, if we want complexity to be a physical attribute of a QFT, then we would expect it to reflect the notion of locality. That is, gates which couple far-separated points should be more expensive -i.e., incur a higher cost in the geometric distance function -than those which couple nearest neighbours.
To gain some experience with this idea, we return to the problem of two coupled oscillators and we introduce a penalty factor weighting the entangling gates, which act on two oscillators (or sites), more heavily than the scaling gates, which act on a single oscillator (or site). Specifically, we may penalize the "off-diagonal" directions by choosing with a > 1. As a result, our original metric (3.19) is replaced by the following more complicated metric: Of course, this geometry reduces to eq. (3.19) upon setting a = 1. The metric has a slightly simpler expression in terms of the pseudo-lightcone coordinates (3.21), where θ = x + z, τ = x − z: which reduces to eq. (3.20) when a = 1. Although these coordinates somewhat obscure our physical intuition for the geometry, they are computationally much simpler. Therefore we will work with the metric in the form (5.3) for most of the following. As in the unpenalized case, this metric enjoys the four Killing vectors (k I ) i given in eq. (B.8). 34 For the metric (5.3), the associated conserved quantitiesĉ I = (k I ) i g ijẋ i are 34 Again, these arise from the right-invariance of the expression in the first line of eq. (5.2). The fifth "accidental" Killing vector ∂x in eq. (B.9) no longer gives rise to a symmetry for the penalized metric, as is clear from eq. (5.3).
(5.6) Substituting eq. (5.6) back into the expressions for derivatives thus renders them well-behaved at the origin, as required by the initial condition ρ(s = 0) = 0, but since their forms are not appreciably simpler we shall not write them out here. Now with the metric (5.3), the normalization of the tangent vector k 2 = g ijẋ iẋj becomes In principle, one can substitute the expressions in eq. (5.5) forρ,ẋ, andż into this expression to obtain an explicit formula for the geodesic length, with no derivatives. Unfortunately, the resulting expression appears quite intractable, and a general solution remains beyond our reach. However, we are ultimately only interested in the optimal trajectory. Given our experience with the unpenalized metric, one might reasonably conjecture that the global minimum is again obtained with the simple straight-line circuit (3.50), and as we now verify, this trajectory remains a geodesic in the penalized geometry (5.3). Recall that the first constraint in eq. (5.4) yielded the desired behaviour for y, i.e., y(s) = y 1 s, as in eq. (3.29). The straightline solution also had τ and θ fixed with τ (s) = 0 and θ(s) = π. This then implies that x and z are fixed with x(s) = π/2 and z(s) = π/2. Combining the latter with eq. (5.6) then yieldsĉ 2 = 0. Substituting these values of x and z into the last two expressions in eq. (5.5), we obtainẋ = −ĉ and we arrive at the desired solution: ρ(s) = ρ 1 s. Having shown that this simple trajectory remains a geodesic in the penalized geometry (5.3), we substitute this geodesic into eq. (5.7) to obtain where we have introduced the label k 0 to denote the geodesic length of the straight-line circuit to avoid confusion with other lengths considered below. Note that eq. (5.10) is the natural generalization of eq. (3.38) to the case with a > 1. However, it turns out that this is not the minimum geodesic for the penalized metric: shorter trajectories can be found. In particular, examining the geometry (5.3) more closely, we see that g ρρ = 2 a 2 − a 2 − 1 sin 2 (2x) (5.11) depends on the x coordinate. This contrasts with the unpenalized metric (3.19) (or the lattice metric (4.29)) where we found that the geometry of the normal-mode subspace (i.e., the metric for the y and ρ directions) was independent of the other coordinates. The latter property was essential to showing that the straight-line circuit was indeed the optimal trajectory. Examining eq. (5.11), it is clear that there should be short-cuts for the motion along ρ if we move away from the normal-mode subspace, i.e., away from x = π/2. 35 For example, we might consider the following simple path consisting of two segments: a) 0 ≤s ≤ 1 : y = 0 , ρ = ρ 1s , x = π/4 , z = π/4 ; (5.12) b) 1 ≤s ≤ 2 : y = y 1 (s − 1) , ρ = ρ 1 , x = π 4s , z = π/4 , wheres provides some arbitrary parametrization of the path. This segmented path is not a geodesic, but does connect the initial point at the origin to the desired end-point at y = y 1 , ρ = ρ 1 and x = π/2. The first segment moves only in the ρ direction at the optimal value of x, and then the second segment moves uniformly in both the x and y directions to arrive at the required end-point. The total length of this path is where we use the subscript s to denote "segmented", in contrast with k 0 above. Of course, the relation between k 0 and k s now depends on the details of the various parameters y 1 , ρ 1 , and a. It is natural that all three of these coefficients are large, in which case one generally finds k 0 > k s . However, to simplify the analysis and illustrate this result, let us consider the regime where the penalty factor is the largest constant, i.e., a ρ 1 , y 1 and ρ 1 , y 1 1. Then we may approximate the two lengths with k s π 2 √ 2 a + √ 2ρ 1 + 2 √ 2 π y 2 1 a 2 + · · · , (5.14) Thus in the penalized geometry (5.3), the length of the segmented path is much shorter than the straight-line geodesic in this regime. Again, the segmented path is not a geodesic and so cannot describe the optimal path, but we shall find that it gives a remarkably good approximation to the optimal geodesic. We would like to find the optimal geodesic but as noted below eq. (5.7), obtaining the general solution seems out of reach. However we can make progress with a simplifying assumption: if we examine the penalized metric (5.3), we see that as the radius ρ increases, the most rapidly growing component of the metric is g zz ∼ a 2 e 4ρ (for generic x). This suggests that motion in the z direction will quickly be suppressed as the geodesics move out from the 35 As an amusing observation, let us add that it is also clear that there are no such short-cuts if a < 1. That is, if we weight the scaling gates more heavily than the entangling gates, then the straight-line circuit (3.50) will remain the optimal circuit. origin. Therefore we simplify our problem by considering trajectories confined to a constant-z submanifold, for which the relevant metric is given by ds 2 = 2dy 2 + 2 a 2 − a 2 − 1 sin 2 (2x) dρ 2 + 2a 2 dx 2 . (5.15) We will return to justify our assumption of no (or little) motion in the z direction below. Obviously, eq. (5.15) is a much simpler geometry, and the analysis of the geodesics becomes much more tractable. We leave the details of solving for the resulting geodesic to appendix F and only refer to certain key results in the comparison below. Our expectation is that the new geodesic is the optimal trajectory, at least for large a, but we must add that we have not provided an irrefutable proof of this result. Furthermore, using the results of appendix F, we also explicitly show below that the segmented path (5.12) provides a good approximation of this optimal geodesic in the regime where the penalty factor is large. Our analysis in appendix F suggests we use the following quantity to more conveniently compare the lengths of the various paths: For the straight-line geodesic, we have simplȳ where as above, the expansion in the second line assumes a ρ 1 , y 1 . Now for the optimal geodesic, we have from eq. (F.22) Comparing these results, we see that in this large a regime,k for the optimal geodesic is much smaller thank 0 for the straight-line geodesic, and extremely close tok for the segmented path. Furthermore, we note that bothk 0 andk are completely independent of y 1 , while it only appears ink s at order y 2 1 /a 2 . In figure 3, we plotk,k 0 , andk s as functions of a for fixed values of the variable (see eq. (F.17)), as well as for various values of y 1 in the case ofk s . Recall the definition (5.16) which shows that these quantities are giving us direct information about the length of the corresponding paths. Hence one clearly sees in the figure that the new optimal geodesic is shorter than the straight-line circuit for all values of a. In the right panel, we also see that the length of the segmented path quickly approaches the length of the optimal geodesic for large values of the penalty factor, in agreement with eqs. (5.18) and (5.19). In fact, if we take the difference of these two equations in the large a limit, we find Given that in this limit, bothk andk s are diverging, it is impressive to find the simple O(1) difference shown above. Figure 4 examines this difference in more detail numerically. Of course, the above results support the conjecture that the new geodesic represents the optimal geodesic and hence yields the shortest possible distance between the origin and the end-point. Since the segmented path (5.12) is not itself a geodesic, it must have a longer length. However, the impressive agreement in eq. (5.20) seems to indicate that this path is coming very close to the optimal geodesic. We can confirm this very clearly by examining x(s) and ρ(s) numerically. As shown in figure 5, the geodesic essentially has two phases: the first in which ρ increases uniformly with fixed x = π/4 and the second in whichρ → 0 and x increases uniformly from π/4 to π/2. As shown in the figure, these two distinct phases are separated by an abrupt but smooth transition, which becomes particularly obvious for larger a. We might note that the growth in y is uniform throughout the entire span 0 ≤ s ≤ 1. However, for large a, the transition occurs for small s (see further comments below) and so y is growing primarily in the second phase where x increases. Hence we can see very explicitly that the behaviour of the optimal geodesic is indeed very similar to that of the segmented path (5.12) for large a. In both plots, the dashed curves correspond to a = 2, while for the solid curves we have set a = 50. One sees that the amount of "time" for which the circuit remains on the constant x = π/4 segment is inversely proportional to the strength of the penalty factor. For illustrative purposes, we have normalized ρ|s=1 to 2 in the left plot, and normalized bothẋ|s=1 andρ|s=1 to 1 in the right.
Let us examine the behaviour of the transition point in more detail. For computational purposes, we define this as the value of s = s trans at which the two (normalized) curves foṙ x andρ cross in figure 5, i.e.,ẋ(s)/ẋ| max =ρ(s)/ρ| max at s = s trans . In the limit 1, this point can be well-approximated by 36 We plot s trans (a) in figure 6. In conjunction with eq. (F.25), one sees that a large penalty factor strongly suppresses the duration of the first phase, and thus the circuit spends most of its "time" -in terms of some fixed total affine parameter -on the second phase. Additionally, one sees that the switchover point appears to go to zero as a → ∞. We can verify this by first approximating h(a) for large a as and then expanding eq. (5.21) in the limit a → ∞, → 0: where O 0 ≈ 1.11502. Comparing this expression with eq. (F.19) for large a, we see that the leading-order term in s trans (a) above may be written as This result has an intuitive explanation: ρ 1 sets the radial distance that must be covered in the first phase (which is large in the → 0 limit while the "angular" change in x remains fixed), while as explained above, a large penalty factor a compels the circuit to complete this motion as quickly as possible. To close this discussion of the penalized geometry, let us reiterate that we have argued that the geodesic confined to the constant-z subspace (5.15) is the optimal geodesic, and hence that its length gives the complexity of the state. That is, when we penalize the entangling gates with eq. (5.1), the complexity of the ground state becomes in the regime a ρ 1 , y 1 .
Recall thats is some arbitrary parameter along the path. Furthermore, given this form, it is interesting to plot this segmented path in the (ρ, θ, τ ) space in order to visually compare it to the straight-line geodesic-see figure 7. This is essentially a comparison of the optimal geodesic in the original geometry (3.19) to that in the new penalized geometry (5.2). One feature that the figure emphasizes is that these two geodesics end at different points along the allowed (blue) spiral at θ + τ = π, ρ = ρ 1 . Using the expression for a general element of GL(2, R) in eq. (3.18), we write the segmented circuit U s (s) as for 0 ≤s ≤ 1 , U b (s) = e y 1 (s−1) cos π 4 (s − 1) e −ρ 1 − sin π 4 (s − 1) e ρ 1 sin π 4 (s − 1) e −ρ 1 cos π 4 (s − 1) e ρ 1 for 1 ≤s ≤ 2 , (5.28) Note that U a (s = 1) = U b (s = 1), as required by continuity along the path. Furthermore, observe that the circuit along the second segment can be re-expressed as where we have defined the rotation matrix R(s) ≡ cos π 4 (s − 1) − sin π 4 (s − 1) sin π 4 (s − 1) cos π 4 (s − 1) .

(5.30)
The interpretation of eq. (5.29) is that upon completing the first segment with U a (s = 1), the circuit performs a rotation (as well as multiplying by the exponential involving y 1 ) along the Figure 7: Sketch of optimal circuits. The vertical axis is τ and the horizontal plane is described by the radius ρ and the azimuthal angle θ, while the y direction is suppressed. The unpenalized minimum (green) goes straight out along (θ, τ ) = (π, 0). The optimal circuit in the penalized geometry is well-approximated by the segmented circuit Us(s) in red: the first segment, Ua(s), goes straight out along (θ, τ ) = (π/2, 0) until reaching ρ1, whereupon the second segment, U b (s), curves upwards to (θ, τ ) = (3π/4, π/4). One sees that, relative to the unpenalized minimum, the segmented circuit arrives at a different but equally valid point along the one-parameter family of allowed end-points given by θ + τ = π (blue), but with a shorter length. The dashed path has an identical length, and is simply the segmented circuit rotated by 180 o about the (θ, τ ) = (π, 0) axis.
second segment until we reach the desired target state. The additional evolution along this second segment is therefore captured entirely by e y 1 (s−1)R (s). In passing, we also note that at the end-point, i.e.,s = 2, the rotation matrix reduces tō which is closely related but distinct from the rotation matrix R defined previously in eq. (3.51). At its end-point, this new circuit becomes Both of these transformations act on the reference state to produce the target state (see eq. (3.8)) as A T = U A R U T . Since the reference state is proportional to the identity, the additional rotation in eq. (5.34) leaves this state invariant, i.e., A R =R(s = 2) A RR T (s = 2), and so both U s (s = 2) and U 0 (s = 1) will produce the same target state, as required.
It is useful to re-express the new circuit (5.28) in the normal-mode space using eq. (3.55), i.e.,Ũ (s) = R U (s) R T , which yields Ũ a (s = 1) for 1 ≤s ≤ 2 , (5.35) where we have expressedŨ b (s) in a form analogous to eq. (5.29). The key observation to note here is that the optimal circuit involves off-diagonal components when expressed in the normal-mode space. That is, in terms of the normal modes, the new circuit is utilizing entangling gates, i.e., gates which entangle (or disentangle) the normal-mode coordinates. Since we know that the reference state (3.53) and the target state (3.54) are unentangled in this space, it must be that along the first segment U a (s), the circuit is introducing entanglement in the state, but then this entanglement is removed along the second segment U b (s) on the second segment. We can see this explicitly by examining the stateÃ along the trajectory. In particular, along the first segment (i.e., for 0 ≤s ≤ 1), we find Here we see the entanglement (i.e., the off-diagonal terms) begins at zero and steadily grows to a maximum ats = 1 at the end of the first segment. Subsequently, along the second segment (i.e., for 1 ≤s ≤ 2), we find = ω 0 e 2y 1 (s−1) cosh 2ρ 1 − sinh 2ρ 1 sin π 2 (s − 1) − sinh 2ρ 1 cos π 2 (s − 1) − sinh 2ρ 1 cos π 2 (s − 1) cosh 2ρ 1 + sinh 2ρ 1 sin π 2 (s − 1) .
Here, the entanglement shrinks steadily back to zero ass runs over this second interval. Recall that y 1 and ρ 1 are given in terms of the normal-mode frequencies in eq. (3.48).
To reiterate our key observation, with eq. (5.1) penalty factors are introduced to increase the cost of the entangling gates in the position space. As a result (for large a), the optimal geodesic is deformed to be close to the segmented path described in eqs. (5.12) and (5.27). However, in the normal-mode space, the new geodesic is driven off of the normal-mode subspace. That is, even though the initial and final states are unentangled when written in terms of the normal modes, the optimal circuit still introduces entanglement (among the normal modes) at intermediate steps along the trajectory. One gains some insight into this behaviour by transforming the penalized metric (5.1) to the normal-mode basis using the orthogonal matrix R defined in eq. (3.62). The new metric then becomes Here we see that in the normal-mode basis, there is no extra cost attributed to the entangling gates relative to the scaling gates. There are also a number of curious negative entries in the off-diagonal components, but this does not fundamentally distinguish the scaling and entangling gates.

Discussion
In this paper, we took the first steps towards defining circuit complexity in quantum field theory. The key idea, due to Nielsen [37], was to endow the space of circuits with an appropriate geometry which allows one to translate the task of finding the optimal circuit into the task of finding the minimum geodesic (with appropriate boundary conditions). We implemented this approach for a simple free scalar field theory. The first step however was to introduce a UV regulator by placing the theory on a lattice, which reduced the scalar field theory to a family of coupled harmonic oscillators. In this context, we were able to construct a interesting set of elementary gates (2.11), in particular, scaling and entangling gates. We also chose our reference state to be a factorized Gaussian state (4.11), whose simplicity lies in the fact that there is no entanglement between different points on the lattice. For the purposes of this preliminary study, we chose the target state to be the ground state (4.10) of the system, which is also a Gaussian state.
To gain some intuition for the problem, we began by studying the simple case of a pair of harmonic oscillators. The fact that both the reference and target states were Gaussian allowed for the simplification that we could translate from an operator language to a matrix language. It was then straightforward to show that with the F 2 cost function, the desired geometry was given by a right-invariant metric (3.19) on GL(2, R). The optimal geodesic was then a simple straight line, which only moved through a flat two-dimensional subspace of the full, more complicated geometry. Translating this geodesic to the optimal circuit, the latter had a particularly simple interpretation in the normal-mode basis, where it only consisted of scaling gates amplifying the individual normal modes. This was a reflection of the fact that the ground state also takes the form of a factorized Gaussian when written in terms of the normal modes. These results for the two coupled oscillators were then extended to the full problem of a lattice of coupled oscillators with relative ease. In particular, we were able to show that the optimal circuit was given by the analogous straight-line geodesic moving in the normal-mode subspace, without constructing the full right-invariant metric on the N 2 -dimensional geometry of GL(N, R).

Comparison with holography:
As discussed in the introduction, a primary motivation for this paper came from recent efforts to understand "holographic complexity," and so it was interesting in section 4.1 to compare our results to those obtained from the holographic proposals. Here we must reiterate the caveat that this comparison involves two very different QFTs, namely a free theory with a single degree of freedom in our scalar field model versus a strongly coupled theory with a large number of degrees of freedom in holography. Hence there is no a priori reason to expect that the results should agree in the two cases. Nevertheless, we found that if the cost function is chosen appropriately, the scalar field complexity exhibits remarkable similarities with holographic complexity. Our tentative interpretation of this concordance is that it provides insight into the implicit cost functions that underly the holographic complexity conjectures.
In particular, the leading divergences (4.36) in both the CV and CA proposals are extensive, i.e., they are proportional to the volume of the time slice on which the boundary state is evaluated, as shown in [19]. While the F 2 cost function gave a result proportional to V 1/2 in the scalar field theory, it is straightforward to construct a family of cost functions in eq. (4.37), all of which yield an extensive complexity for the scalar field theory. 37 With the new cost functions (4.37), the leading contribution also contained a logarithmic factor, which was ambiguous in that it depended on the choice of the frequency ω 0 specifying the reference state (4.11). However, this precisely matched an ambiguity in the holographic complexity [19] found for the CA construction (1.3). In the latter case, the logarithmic factor came from joint terms [17] in the gravitational action, and the ambiguity arose from the freedom to choose the normalization of the null normals on the boundary of the WDW patch. Whereas this ambiguity had originally been seen as problematic for the CA conjecture, our scalar field calculation indicates that it is a perfectly natural feature associated with the freedom in the choice of the reference state that we can anticipate in any definition of complexity for a QFT.
It might then seem mysterious that no such ambiguity arises for the CV conjecture (1.2). However, as explained in section 4.1, the additional logarithmic factor in the leading term (4.41) becomes a simple numerical coefficient if we choose ω 0 = e −σ /δ, and so such a choice may indeed be an integral part of the microscopic rules implicit in the CV construction. Unfortunately, this does not explain the absence of infrared terms of the form given in eq. (4.43) which might be expected with this choice. While it would be premature to conclude that the CV conjecture is incorrect, we might note that there is an alternative proposal in the literature suggesting that the volume of a maximal time slice in the bulk should be dual to the information metric rather than the complexity in the boundary theory [54].
As further noted in section 4.1, our scalar field complexity emulates the CA proposal (1.3) most closely if we choose the F 1 cost function, i.e., κ = 1 in eq. (4.37). Given the aforementioned disparity in the two field theories in question, it is unclear how much weight to give this observation. However, we might add that the F 1 measure is a natural choice since it adheres most closely to the original definition of complexity, which involved simply counting the number of gates in the optimal circuit. Furthermore, let us add that the F 1 cost function will also feature again in the discussion of cMERA networks below.
Of course, it would be interesting if a more precise connection can be found between the holographic and QFT calculations with regards to the ambiguity in the reference state discussed above. At present, it is actually not clear how the reference state enters in the holographic calculations at all, but perhaps one can draw upon the proposal for a statesurface conjecture in [55]. Undoubtedly, making this connection concrete would bring us closer to an explicit translation for the complexity between the bulk and boundary.

Ambiguities and other miscellaneous complaints:
In the introduction, our "definition" of complexity was rather imprecise, as it left open the choice of the reference state |ψ R , the choice of the set of elementary gates which would be used to construct U , and the choice of the tolerance (and measure) in eq. (1.5). Clearly, even though it is easy to set out interesting questions for complexity (e.g., what is the complexity of a particular state in a particular QFT?), the precise value of the complexity will depend on the details of all of these choices. 38 The ambiguity in our reference state, i.e., the choice ω 0 , was already seen to modify the complexity in an interesting way in the preceding discussion.
Here one might recall early discussions of entanglement entropy from a hep-th perspective: the explicit dependence of the leading contributions on the UV cut-off was certainly seen as problematic (or at least, it was by one of the present authors). However, with some experience, we learned to find universal information in the entanglement entropy and to apply it as a useful diagnostic of QFTs in various ways, e.g., [56,57]. In fact given our experience with entanglement entropy, the non-universality of the leading contributions to the complexity (because of the power-law dependence on δ, as in eqs. (4.35) and (4.41)), was assumed to be self-evident in the present discussion and not even commented upon. Analogously, we would advocate that complexity is again a new quantity with what initially seems to be unusual and perhaps undesirable features, but that we must develop our experience to learn how complexity can inform us about interesting physics and universal properties of QFTs and holography. Hence rather than regarding the ambiguities discussed above as a problem per se, they should be collectively considered as a new feature which we must learn to accommodate in working with complexity.
In the context of the present calculations, clearly evaluating the complexity for a single state, e.g., the ground state, will not be particularly informative. Instead we might compare the complexities of different states, and while extending our calculations to the complexity of excited states would allow such a comparison, it is beyond the scope of the present paper. However, we can certainly compare the complexities of the ground states of different scalar theories, in particular theories with different masses. Here our experience with holographic complexity [52] suggests that there may be interesting information that could be extracted from the finite or logarithmic contributions. For example, for even spacetime dimensions (even d), there will be an interesting contribution that is intrinsically independent of the cut-off, although it may depend on the reference frequency ω 0 . This explicitly appears in eq. (E.12), where we are examining the case d = 2 and κ = 1. In this instance, we may isolate this constant in 39 ( where L is the linear (spatial) size of the system. This result is independent of both the short-distance cut-off scale and the reference frequency. Hence it would be interesting to better understand the meaning of the coefficient of a 0 , and whether it carries some universal information about the underlying theories. Of course, it is straightforward to extend this simple example to higher dimensions. It would also be interesting to compare these results to similar calculations for holographic complexity where the boundary CFT is deformed by a relevant operator-see discussion at the end of section 4.1.
Let us also remark here that the reference state (4.11) is an unusual state from the textbook perspective of QFT. Recall that this state was chosen since it has no entanglement between different points (on the lattice). Such a factorized Gaussian is precisely the kind of reference state that appears in the cMERA construction [49]-see the discussion below. However, such an unentangled state is a very unusual state in standard QFT, which for example would typically have a divergent energy density. Of course, the vacuum energy density of the ground state is also divergent, and so to make this statement meaningful, we may evaluate the difference in the energies of |ψ R and the ground state |ψ T : We therefore see that that generically, if ω 0 δ 1 or ωδ 1 (which were advocated to be the natural choices in section 4.1), the "renormalized" energy density of |ψ R diverges as δ → 0. 40 Hence this would not be a state that would be considered to be part of the standard 39 Of course, we are inspired to formulate this quantity by the constructions using entanglement entropy to examine RG flows of three-dimensional QFTs [58,59]. 40 Note that with some fine-tuning of ω0, we could arrange the difference of energy densities to be finite.
Hilbert space that one builds with particle excitations on top of the vacuum. However, one should simply regard this as another unusual feature of complexity. As we have seen both here with our QFT calculations as well as in holographic complexity, the complexity can only be sensibly defined with a finite value of the regulator, in which case the reference state is certainly a sensible state within the associated Hilbert space.
While introducing a UV regulator was an essential step in sensibly defining the complexity in the scalar field theory, let us add that this does not regulate the size of the Hilbert space in the present case. 41 With the lattice regulator, the scalar field theory is reduced to N d−1 normal-mode oscillators, but the Hilbert space of each of these oscillators is infinite! It is an interesting question whether or not an additional regulator should be introduced to render the total number of states finite as well. Otherwise it would seem that even within the UV regulated theory, there will be states of infinite complexity.

Penalty factors and locality:
In section 5, we experimented with the introduction of penalty factors in the case of two coupled oscillators. In particular, we gave a higher cost to the entangling gates than the scaling gates with (5.1). This certainly resulted in a different optimal circuit, but ultimately the circuit could not avoid incorporating the entangling gates, and so the complexity increased to O(a), as shown in eq. (5.26). Perhaps the most interesting lesson to be learned from these calculations is that the introduction of penalty factors (in the position-space cost function) tends to drive the optimal circuit away from the normal-mode subspace, the restriction to which played a central role in the previous analysis of section 4. However, in our simple experiment in section 5, the optimal circuit was still required to introduce entanglement using the entangling gates, and therefore our calculations did not really address the motivation discussed at the beginning of that section. Namely, we expected that penalty factors could be used to introduce a notion of locality in the complexity of the scalar field theory. In particular, our calculations in section 4 included entangling gates Q ab , which coupled points on the lattice that were arbitrarily far apart, all with equal cost. It seems natural that using gates which couple far-separated points should incur a higher cost than using those which couple nearest neighbors.
To gain some insight into this problem, let us return for a moment to the discrete gates in eq. (2.11) and consider a one-dimensional lattice of N coupled oscillators. We will show that our set of entangling gates is over-complete in the sense that any of these gates can be constructed from nearest-neighbour entangling gates. For example, it is a straightforward calculation (using the Baker-Campbell-Hausdorff formula) to show that 3) 41 We thank Edward Witten for making this observation and raising the following question.
In other words, the next-to-nearest neighbour entangling gate Q 13 is equivalent to 4/ε nearestneighbour entangling gates. A simple generalization of the above result is which implies that the next-to-next-to-nearest neighbour entangling gates are equivalent to 8/ε 2 + 2/ε nearest-neighbour gates (8/ε 2 from the use of the Q 13 's and another 2/ε from the Q 34 's). These calculations can be easily generalized to show that nonlocal gates Q a,a+1+n or Q a+1+n,a , which entangle oscillators that are separated by n intermediate sites, can be constructed by the use of nearest-neighbor entangling gates, where c(0) ≡ 1. Thus to leading order, the "cost" of these nonlocal gates in terms of nearest-neighbour gates grows like c(n) ∼ 1/ε n . Following Neilsen's approach [37][38][39], we would not eliminate these nonlocal gates from the elementary gate set, but would instead modify the geometry by introducing (heavy) penalty factors to discourage the geodesics from moving along the corresponding directions. The structure of eq. (6.5) suggests increasing the penalty factors as a power law to match the growth of the nonlocality, i.e., the directions corresponding to I = (a, a + 1 + n) and (a + 1 + n, a) would be assigned a penalty factor a 2n . Note that this does not penalize the nearest-neighbour gates at all, in contrast to our experiment in section 5. Of course, for a periodic chain of oscillators, the maximum penalty factor would be a N −2 and a N −3 for even and odd N , respectively.
We can gain further insight by translating eq. (6.3) into a macroscopic circuit described by a path-ordered exponential (3.12). In particular, consider the following path: where 0 ≤ s ≤ 4, Θ(x) is the Heaviside theta-function, and the Y I are implicitly zero for all other values of I. That is, we turn on the M 23 generator with amplitude α for the interval 0 ≤ s ≤ 1; M 12 is then turned on with amplitude α for 1 ≤ s ≤ 2; next, M 23 is turned on with amplitude −α for 2 ≤ s ≤ 3; and finally M 12 is turned on with amplitude −α for 3 ≤ s ≤ 4.
Note that the precise parametrization of this path is not important. The circuit then yields U 1 (s = 4) = exp α 2 M 13 , following the same calculation that yields eq. (6.3). Hence we could accomplish the same transformation with U 2 (s) = P exp Now let us compare the costs of these two circuits using the F q measure (1.9), where the nearest neighbour gates are assigned cost 1 while the next-to-nearest neighbour gates are assigned cost a. The cost functions are then easily evaluated to be Hence with an appropriate penalty factor, we can suppress the use of the nonlocal gates in favour of the nearest neighbour gates. While it would be interesting to examine the effect of the above scheme of penalty factors in more detail, we leave this for future work.

cMERA:
The AdS/MERA correspondence was the first proposal for a novel connection between holography and tensor networks [60,61]. This proposal suggests that the MERA (Multiscale Entanglement Renormalization Ansatz) tensor network [34,62,63] provides a discrete representation of a time slice of (three-dimensional) AdS space. As illustrated in figure 8, the MERA network consists of unitary operators which, starting from the simple product state |0 ⊗ . . . ⊗ |0 , generate the ground state in d = 2 critical systems. In other words, the MERA network can be thought of as a quantum circuit. The AdS/MERA correspondence was certainly a source of motivation/inspiration for the early discussions of holographic complexity, in particular, of the CV conjecture [9,12]. Furthermore, in these discussions, it was implicitly considered the optimal circuit for the preparation of the CFT ground state.

MERA as a quantum circuit
"time" Entanglement introduced by gates at different "times" (= length scales) quantum circuit Figure 8: Illustration of MERA as a quantum circuit. Starting from the tensor product state |0 |0 . . . |0 (top), the sequential application of entanglers and isometries (colored dots) efficiently generates the ground state (bottom). These operators can be thought of as unitary gates (although not simple elementary gates), and hence the tensor network can be thought of as the quantum circuit that connects the reference (product) and target (ground) states ψ0 and ψ1. There has been some progress towards developing a continuum version of MERA, however, these constructions are limited to describing very simple QFTs [49]-see also [65][66][67].
In particular, one example is the cMERA description of the ground state of a free scalar field. That is, there is a cMERA circuit which more-or-less performs precisely the transformation for which the circuits studied herein were constructed. Hence, our original expectation was that our analysis would find that the optimal circuit was something like a cMERA network. However, we instead found the straight-line circuit described in section 4. The key difference between the two circuits is that the cMERA circuit is organized to systematically introduce entanglement scale-by-scale, i.e., to order the amplification of the normal modes according to their wavelength [67]. However, by almost all of the measures considered in section 4, including the F 2 cost function and the κ cost functions in eq. (4.37), the straight-line circuit is the optimal circuit. The one exception to this rule is the F 1 (or κ = 1) cost function. This last describes an unusual geometry, 42 which is sometimes called the "Manhattan metric." The key feature of this geometry is that the length is the same for all paths as long as they do not back-track at any point, and hence the straight-line circuit and the cMERA circuit have identical costs for this measure.
This question certainly deserves further study. It appears that there are two possible approaches: the first would be to study more exotic cost functions in order to identify those which favour the cMERA circuit. This may be useful since given the AdS/MERA duality, it may provide better insight into the properties of the cost function that appears in holographic complexity. We also mention that this is likely not a straightforward approach since we found that introducing penalty functions (in position space) seems to drive the optimal circuit out of the normal-mode subspace, whereas the cMERA circuit is confined to this particular slice of the full circuit geometry by construction. A second option might be to introduce new physics in the selection of the "optimal" circuit. That is, while the straight-line and cMERA circuits have equivalent costs according to the F 1 cost function, there may be additional physics considerations, e.g., some relation to renormalization group flows, which lead holography to favour a cMERA-like circuit.

Future Directions:
This paper provides only a preliminary investigation towards understanding circuit complexity in quantum field theory. We already mentioned a number of future directions that we expect will be fruitful. Some examples include extending the present calculations to evaluate the complexity of excited states, producing a more concrete connection between the ambiguities arising in our QFT calculations and those in holographic calculations of the complexity, and studying in detail the effect of penalty factors on the complexity and the structure of the optimal circuit for a lattice of oscillators. Other obvious extensions of the present work would include evaluating the complexity in fermionic theories or in interacting QFTs.
In closing, we would like to draw a comparison with entanglement entropy in QFT. Entanglement entropy has a simple textbook definition: first one must construct the reduced density matrix ρ A of the particular subsystem under study, and then one evaluates the von Neumann entropy of this density matrix as S EE = − λ i log λ i , where λ i are the eigenvalues of ρ A . However, much of the progress in understanding the properties and role of entanglement entropy came from the replica trick, introduced by Calabrese and Cardy [68,69]. The latter applies familiar tools (e.g., path integrals) in a novel setting (e.g., the replicated background geometry) to evaluate the entanglement entropy. Returning to complexity, our present approach is to apply a more-or-less standard textbook definition to evaluating the complexity of states in a QFT, which is a useful preliminary step to gain an understanding of the properties of this new quantity. However, we would really like to develop a new approach, analagous to those developed for entanglement entropy above, which again uses familiar QFT techniques in a presumably novel setting to evaluate some quantity like the complexity. In other words, we are asking what is the new calculation of complexity which is the analog of Calabrese and Cardy's replica trick for entanglement entropy. Indeed, it may be that the first steps in this direction have already been taken in [70,71]-see also [72].

A Example circuits
In this section, we analyze the circuit depth of a few discrete circuits, as introduced in section 2. In particular, let us consider the example given in eq. (2.14), where the target wave function is given in eq. (2.8) and the reference wave function, in eq. (2.10). The question then is to determine the exponents α i in this equation, i.e., the number of times each type of gate is applied in the circuit. Intuitively, U 1 is a string of gates running from the right to the left. The first gates applied are Q 11 and these rescale the coefficient of x 2 1 in the exponent of the Gaussian wave function. Next by applying Q 21 , the two oscillators become entangled, and finally the Q 22 gates rescale x 2 to ensure that x 2 2 appears with the correct coefficient. Hence let us begin the quantitative analysis by considering: where Next applying the Q 21 gates yields Note that the x 1 x 2 cross-term will be rescaled in the next step and so we cannot fix any of the coefficients at this step. So finally, we rescale x 2 with the Q 22 gates: where α 2 and α 3 are determined matching the second and third coefficients in the exponent, i.e., Solving the above constraints then yields As a consistency check, note that with these identifications, the normalization factor of the final wave function becomes which is correct for the wave function to have norm one. Of course, this was expected since, as discussed in the main text, the entangling and scaling gates, Q ij and Q ii , preserve the norm when acting on a wave function. Hence the total number of gates in the circuit U 1 in eq. (A.1) is given by where we have assumed here that ω 1 > ω 0 and ω 1 ω 2 − β 2 > ω 2 0 . 43 As in the main text, we refer to D(U 1 ) as the circuit depth, rather than the complexity, since while it counts the total number of gates in the circuit, we have no reason to expect that U 1 is the optimal circuit. Recall that we introduced absolute values in eq. (A.7) in order to give an equal complexity cost for the inverse gates Q −1 ij as for the original gates Q ij , i.e., we count the appearance of Q −1 ij as one gate in a circuit. At a pragmatic level, this is required because α 2 is negative in our example, i.e., β = (ω + − ω − )/2 < 0.
Note that in evaluating the exponents α i in eqs. (A.2) and (A.5), we are implicitly treating them as real numbers. If we insisted on having integer exponents, then we would would need to round these results up or down to the nearest integer. In this case, we would define a measure of success for our transformation and choose the integer exponents to maximize this measure. For example, we could consider the overlap and choose the precise integer values of α i to minimize χ. Of course, using real exponents α i is very much in line with describing the circuits in terms of path-ordered exponentials (3.1). Of course, this discussion is related to the choice of a tolerance ε in eq. (1.5), i.e., rather than minimizing χ, one might demand that χ ≤ ε. Now let us briefly present a few other examples of simple circuits to further familiarize the reader with the concepts discussed here. First, let us consider first applying the entangling gate before either of the scaling gates: Note that for comparison purposes, our numbering of the exponents is such that they are associated with the same gates as appear in eq. (A.1). The calculation proceeds essentially as above and in the end, we must match the coefficients ω 1 = ω 0 e 2 α 1 , ω 2 = 1 + 2α2 2 ω 0 e 2 α 3 , β = α 2 ω 0 e (α 1 +α 3 ) . (A.10) 43 See further comments at the end of this appendix.
Solving for the exponents α i then yields and hence the circuit depth becomes Comparing the results in eqs. (A.2) and (A.5) with those in eq. (A.11), we see that the exponents for the scaling gates are identical, i.e., α 1 =α 1 and α 3 =α 3 , and only the exponent for the entangling gate has changed. Hence the circuit depth is almost identical to (A.7), except that the second term lacks a factor ω 0 /ω 1 . If we assume ω 1 > ω 0 as before, the present circuit will be slightly longer, i.e., D(U 2 ) > D(U 1 ). As a third simple example, let us consider applying the entangling gate after both of the scaling gates: Again we skip over the details of the calculation; we find that we must match the coefficients Solving for the exponents α i then yieldŝ and hence the circuit depth becomes this becomes (A.18), as claimed. Now, acting with the circuit U 4 and matching coefficients as before, we find Expandingᾱ 2 near ≈ 0, we havē (A. 25) We therefore find that the circuit depth for U 4 is (A.26) Here, as above, we assume ω 0 < ω 1 .
In general, we can describe the form of the circuit depth as being an overall factor of 1/ followed by a coefficient determined by the various physical parameters characterizing the target state and the reference state. More generally, the circuit depth might be given by an expansion in , beginning with a 1/ term followed by a finite term and then potentially terms involving positive powers of . However, since 1, determining the complexity essentially requires finding the circuit which minimizes the coefficient of the leading 1/ term. For comparison to the results for the geometric approach, it is useful to express the present results in terms of the normal-mode frequencies using eq. (2.9). If we focus our attention on the first circuit U 1 in eq. (A.1), the exponents given in eqs. (A.2) and (A.5) become (A.27) As was alluded to above, to proceed further we must decide on the value of the reference frequency ω 0 relative to the normal-mode frequencies. Given the discussion in section 4, there are two natural hierarchies to consider: a)ω + <ω − < ω 0 or b) ω 0 <ω + <ω − . 44 Of course, the ordering of the normal-mode frequencies is fixed and we are really only choosing ω 0 here. In particular, in the first (second) hierarchy, ω 0 is an UV (IR) frequency larger (smaller) than any physical frequency in the coupled oscillator problem. Note that in the first case, all three exponents are negative, while in the second case, α 1 , α 3 > 0 and α 2 < 0. Now evaluating D(U 1 ) = |α 1 | + |α 2 | + |α 3 | in these two cases yields: Recall thatω − >ω + from eq. (2.6). We may compare this result with those derived using the geometric approach of section 3. In particular, if we recall the F 1 measure given in eq. (1.9), the complexity would be given by Further recall that we should compare this with the coefficient of the 1/ factor in the discrete calculations (see footnote 7). Hence we see the second contribution in eq. (A.28) precisely matches the complexity above. However, there is an additional positive term added to this term in D(U 1 ) and therefore, we see that (at least by the F 1 measure) U 1 is not the optimal circuit. We can also describe this circuit as a trajectory in the language of the path-ordered exponentials (3.1). In this case, U 1 as given in eq. (A.1) becomes This form makes clear that the circuit consists of three separate 'straight' segments and so U 3 does not correspond to a geodesic path or an optimal circuit.

B Killing vectors and more geometry
Inspecting the metric in eq. (3.19), we can see that three obvious Killing coordinates: y, τ, θ.
When the penalty factors are introduced in section 5, we find that this is reduced to two Killing coordinates, y and z = (θ − τ )/2, in the geometry described by eqs. (5.2) or (5.3). However, by construction all of these metrics are right-invariant, and hence the corresponding geometries must have one Killing vector for each generator (3.11), namely, four 45 Further as we will see below, the structure of these Killing vectors will be completely independent of the particular choice of G IJ appearing in the metric (assuming that it is a constant matrix), but rather is determined by the structure of eq. (3.13). One way to think of a Killing vector k i is that it gives a coordinate transformation where we have assumed that we are working with the orthogonal basis of generators satisfying tr M I M T J = δ IJ , e.g., see eq. (3.11). We now observe that, since the argument of the trace contains two free indices, we may view this object as a 4 × 4 matrix, which we can then invert to obtain 45 We thank Lucas Hackl for discussions on this point.
Hence we have four independent Killing vectors k I = (k I ) i ∂ i . Given our basis of generators in eq. (3.11) and our parametrization of the circuit space in eq. (3.18), we can easily compute (k I ) i . We can then identify the Killing vectors by simply reading off this matrix row-by-row: where we are using the pseudo-lightcone coordinates of eq. (3.21), with θ = x + z, τ = x − z.
One can explicitly verify that these indeed satisfy the Killing equations, for either of the metrics in eqs. (3.20) or (5.3). However, it is clear that eq. (B.6) does not organize the Killing vectors in the simplest way and therefore we define: However, a simple inspection of the first metric (3.20) reveals that ∂ x is also an independent Killing vector. We therefore definek This is an accidental symmetry that emerges with the choice G IJ = δ IJ . 46 However, as noted above, the four Killing vectors in eq. (B.8) apply for any choice of (constant) G IJ . Of course, the existence of the above Killing vectors implies that there are an equal number of conserved momenta or charges which distinguish the geodesics, i.e., c I ≡ (k I ) i g ijẋ j .
We make use of these momenta in solving for the optimal circuits in sections 3.1 and 5.
AdS 3 geometry: In section 3.1, we noted the appearance of a three-dimensional anti-de Sitter geometry in discussing the parametrization of U ∈ GL(2, R) = R × SL(2, R). Of course, it is natural for 46 We thank Lucas Hackl for discussions on the Killing symmetries.
the AdS 3 to appear since it is the universal cover of the SL(2, R) subgroup. Here we would like to show how the AdS 3 geometry can be realized using the formalism introduced in section 3. In particular, we consider the geometry that results with the choice is a less symmetric metric with only rightinvariance. However, the Lorentzian signature is undesirable for the problem of circuit complexity: since pieces of the circuit that correspond to null-geodesics have zero length or zero cost. This would allow the construction of arbitrarily low-complexity circuits simply by deforming the circuit into null segments. 47 while M 1 describes the remaining R fibre in the GL(2, R) group. With this new basis, the Killing vectors, emerging from the right invariance of the metric, naturally appear in the form given in eq. (B.8). One can easily show that working with these new generators, the metrics appearing in eqs. (3.19) and (5.2) are unchanged, i.e., the corresponding G IJ are left unchanged by the rotation R I J introduced above. Further the AdS 3 geometry in eq. (B.11) is now produced with the choice G IJ = η IJ = diag(1, 1, 1, −1).

C Normal-mode frequenciesω k
The derivation of the normal-mode frequencies, e.g., eqs. (4.6) or (4.33), for a lattice of coupled oscillators is straightforward and can be found in a number of different sources, e.g., any elementary condensed matter textbook. For completeness, we briefly review the result (4.6) for the periodic one-dimensional lattice discussed in section 4. Essentially we only need to apply the inverse Fourier transform to re-express the Hamiltonian (4.1) in terms of the normal modes, i.e., eq. (4.2). In particular, we focus on the potential Considering the second term involving the coupling between the oscillators, we find where in going to the third line we applied the normalization condition (4.4), and in the last step, we usedx kx−k =x kx † k . All sums run from 0 to N −1. The Fourier transform of the first term in the potential (C.2) is trivial, and so we find Hence we have identified the desired normal-mode frequencies, as in eq. (4.6), If we were examining a d-dimensional free scalar field, the lattice would be extended to d-1 (spatial) dimensions and the corresponding normal-mode frequencies becomẽ where k i are the components of the momentum vector, i.e., k = (k 1 , k 2 , · · · , k d−1 ). Implicitly, we have assumed here that the lattice is square with periodic boundary conditions in each direction.

D Change of basis
In this appendix, we would like to extend the discussion, around eq. (3.62), describing the change of bases for the case of two coupled oscillators to the analogous transformation for a lattice of oscillators. In particular, we will focus on the case of a one-dimensional lattice of N oscillators, although it is straightforward to extend the discussion to a lattice extending in d-1 (spatial) dimensions. This transformation is particularly relevant in section 4, where we presented a tentative argument with eq. (4.18), that the metric on the normal-mode subspace is flat, i.e., However, we also noted there that this conclusion is somewhat premature since implicitly we applied eq. (4.14), which defines the metric in the position basis, to a calculation with the diagonal circuit (4.17) written in the normal-mode basis. That is, in eq. (4.14), the indices I, J run over pairs of position labels (ab) and implicitly the generators act on Gaussian wave functions written in terms of coordinates x a . In contrast, in eq. (4.17), we would writẽ M n-m =ỸĨMĨ where the tilde on the indexĨ indicates that it runs over pairs of momentum labels (k ) and the tilde on M indicates that these generators act on Gaussian wave functions written in terms of the normal coordinatesx k . That is, in eq. (4.17), where we are restricting our attention to the normal-mode subspace, we consider the diagonal generators:Ỹ k = δ k ỹ k . Hence to show that the result (D.1) is correct, we must take care to translate between the two bases of generators, discussed above. As applied in eqs. Implicitly the normal-mode generatorsMĨ have the form identical to that given in eq. (3.10), i.e., where we have denotedĨ = (k ) with k, are two momentum labels. Similarly p, q are the row and column indices, respectively, of the N × N matrix, which also take values as momentum labels (since the generator acts in the normal-mode space). Now let us combine these two equations to write 49 In going from the second to third line, we used eq.
In going from the second to third line, we used eq. (D.3) and identified I = (ab) andJ = (pq). We can easily invert the transformation in eq. (D.6) to writeMĨ = [ R N ]Ĩ JM J . As our notation indicates, R N is precisely the unitary matrix appearing in eq. (D.4) and hence it also plays an role in transforming the generators acting in the normal-mode space. Hence by using the special structure of the generators in eqs. With these tools in hand, let us consider re-expressing the cost function (4.13) or the metric (4.14) in terms of the normal-mode basis, using eq. (D.2). Here, we show the calculation 49 Note that the complex conjugation appears on the first factor in RN = R * N ⊗ RN because our convention is that written in terms of the normal modes, the Gaussian wave functions involve bothx k andx † k , e.g., the appearance of |x k | 2 in eq. (4.7).
for the metric and transformation of the cost function follows in a similar manner. Beginning with the differential dY I = tr dU U −1 M † I defined in eq. (4.14), we transform the circuit to the normal-mode space using U = R N †Ũ R N , which yields where we have employed eqs. (D.5) and (D.6) in the second and third equalities above. Hence the metric (4.14) transforms as Note that here we are using the fact that R N is a unitary matrix. Hence we have found that the metric takes precisely the same form whether expressed in terms of the oscillator-position space or the normal-mode space. 50 Of course, the same is true of the cost function (4.13), i.e., it can also be written as Note that this transformation is slightly different than that expressed in eq. (3.63) for the metric for two coupled oscillators. There we are still considering the metric remains in the position basis but it is evaluated with a different basis of generators. The same invariance holds here for a lattice of oscillators, as can be seen by applying eq. (D.4) directly on the metric (4.14) to produce Of course, the same change of basis could also be performed with eq. (D.6), when working in the normal-mode space.

D.1 General cost functions
In eq. (4.37), the κ cost functions were defined with a sum over the components of the velocity YĨ in the normal-mode basis. Here we would like to apply the techniques developed above to examine the differences arising from using the original oscillator position basis. That is, we could equally well define cost functions with (D.11) 50 The fact that this transformation preserves the cost function essentially follows from the Plancherel theorem, which states that the Fourier transform preserves the L 2 norm. We thank Adrián Franco-Rubio for a discussion on this point.
In the previous section with the F 2 cost function, we found that this change of basis had no effect on the complexity but here we will find that, in fact, the complexity is not basis independent. As a simple example, let us consider the case of two-coupled oscillators for which the optimal circuit U 0 (s) appears in eq. (3.50) and so the position-basis velocity components become These two factors are written in terms of the normal-mode frequencies in eq. (3.48), but we can re-express these results as Recall thatω − >ω + , but in the following, we also assume thatω ± > ω 0 , which ensures that both y 1 and ρ 1 are positive quantities. Now we evaluate the cost of U 0 using eq. (D.11) for a few values of κ, 51 Hence we see that it is only for κ = 2 that we reproduce the cost found using eq. (4.37) in the normal-mode basis, i.e., D κ (U 0 ) log κ (ω − /ω 0 ) + log κ (ω + /ω 0 ). These differences in the cost can be understood using the approach developed to implement a change in basis for a lattice of oscillators in section 4. In particular, transforming from the position basis to the normal-mode basis can be described in terms of the unitary matrix R N defined in eq. (D.4). Given the definition of the velocity components in eq. (4.13), we then have 15) or this can be inverted as Y I = R T N IJ YJ . Further, as the discussion there shows the quadratic construction δ IJ Y I (Y J ) * = δĨJ YĨ (YJ ) * is invariant under this change of basis. Therefore the cost evaluated with the F 2 or κ = 2 cost functions are invariant as well.
However, this discussion also makes clear that if we include penalty factors, then these quadratic cost functions are no longer invariant. That is, the penalty factors introduce a more 51 In the case that ω0 >ω±, one should replaceω±/ω0 → ω0/ω∓ in these formulae. Note this substitution only really changes the results for odd κ. general metric G IJ , which transforms nontrivially under the change of basis, i.e., where we assumed symmetry of the metric G IJ = G JI .
This also suggests how we should treat the more general κ cost functions. We should generalize eq. (D.11) to allow for general penalty factors by writing where G I 1 I 2 ···Iκ is a symmetric tensor with κ indices. In eq. (4.37), we are implicitly considering simple 'penalty' tensors of the form In general, it is clear that the unitary transformation will not leave these penalty tensors (or more general choices) invariant. This simply reflects the fact that in choosing different gates, we are treating different gates as fundamental and that in general, we expect the results for the complexity to depend on the choice of the elementary gate set. Of course, this does not mean that the complexity must be evaluated in one particular basis. However, if the cost function is fixed with a certain choice of basis, then changing the basis requires that we properly transform the cost function to the new basis. To gain a better understanding of this situation, let us investigate the case of κ = 1 in more detail. As well as the simplicity of this case, recall that this was the cost function favoured in the comparison to holographic complexity.
Let us begin with the case of N = 2, in which case the transformation matrix R = R 2 takes the simply form given in eq. That is, expressing our κ = 1 cost function (D.11) in terms of the normal-mode basis, we are only penalizing the scaling gate associated withx + ! The other (normal-mode) gates can be inserted in the circuit at zero cost. However, we must add that the transformation in eq. (D.19) is slightly naive since in assumes that the absolute values in the cost function (D.11) play no role, i.e., we are assuming that all YĨ ≥ 0 (or all YĨ ≤ 0). However, one finds that depending on the signs of the various velocity components, only one of the normal-mode gates is penalized at a time. For example, with YĨ ≥ 0 forĨ = ++, −− and YĨ ≤ 0 for Ĩ = +−, −+, 52 one finds GĨ = (0, 0, 0, 2), i.e., only the scaling gate associated withx − is penalized. Similar results arise if we begin with the κ cost functions (4.37) expressed in terms of the normal-mode basis and examine their structure in the position basis. In this case for κ = 1, the original and transformed penalty tensor become GĨ = (1, 1, 1, 1 Hence we have the rather curious result that this cost function is only penalizing the scaling gate associated with x 1 , the position of the first oscillator. Of course, we must again remind the reader that eq. (D.21) assumes that the absolute values in the cost function (4.37) play no role. This assumption is more natural in this case, as with a natural choice of ω 0 , we find that all YĨ ≥ 0 (or all YĨ ≤ 0) for the optimal circuit, i.e., all of the scaling components have a definite sign and all components in the entangling directions vanish. Further using eq. (4.37), we might note that the cost of our straight-line circuit is simply again assumingω ± > ω 0 . Here we emphasize that since only two of the velocity components were non-vanishing, i.e., Y ++ and Y −− , we would arrive at the same cost for a family of penalty tensors of the form GĨ = (1, a 2 1 , a 2 2 , 1) (D. 23) In this case, transforming to the position basis as in eq. (D.21) yields At first sight, this result seems yield a more reasonable penalty tensor, in comparison to eq. (D.21). However, upon closer examination, we see that G 12 = −G 21 and hence one of these penalty factors will be negative. That is, the cost of the circuit will be reduced by including more of one type of the entangling gates, in the normal-mode basis! The only resolution of this unsatisfactory situation is to set the two penalty factors equal, i.e., a 1 = a 2 = a with which eq. (D.24) becomes G I = (1 + a 2 , 0, 0, 1 − a 2 ) , (D. 25) which still requires that a 2 ≤ 1 in order that G 22 ≥ 0. The above results are somewhat unsatisfying. In that, a perfectly reasonable penalty tensor in one basis yields an undesirable or even inconsistent (in the case of negative penalty factors) cost function in another basis. We return to this point in the discussion section, however, we should say that some of these issues arise because of we focused on the simple case of κ = 1. For example, if we consider a κ = 2 cost function (D.11) with our penalized metric (5.1), then transforming to the normal-mode basis yields eq. (5.38). While the resulting metric has negative entries, we know that this in itself is not worrisome. Rather one must examine the eigenvalues of the new metric but, of course, these have not been changed in the transformation and all remain positive.
To close this section, let us comment on extending this discussion to a lattice of oscillators. In particular, we observe that the essential features of the complexity noted in section 4.1 using the κ cost functions (4.37) constructed in the normal-mode basis remain unchanged when working with eq. (D.11) in the position basis. For a (d-1)-dimensional spatial lattice of oscillators, the κ = 1 penalty tensor in eq. (D.18) becomes GĨ = (N d−1 , 0, 0, · · · , 0) (D. 26) in terms of the normal-modes. Hence as in eq. (D.20), only the scaling gate of lowest normal mode is penalized but the cost of that single gate has been increased to N d−1 , the total number of oscillators in the lattice. 53 Hence cost for the straight-line circuit becomes Hence the cost is still proportional to V /δ d−1 , as desired to emulate the holographic complexity. This factor is again multiplied by a logarithmic factor whose argument depends on the reference frequency ω 0 . However, since only the lowest eigenfrequencyω k=0 = m appears in the (single) logarithmic factor, the cut-off scale can only appear in this result through the reference frequency ω 0 . Hence δ appears if ω 0 is chosen as a UV frequency, e.g., ω 0 = e −σ /δ, but it does not appear if the reference frequency is chosen as an IR frequency. We observe that this is the opposite to the situation discussed in section 4.1.

E Approximating the complexity
In section 4.1, we compared our result (4.32) for the complexity of the ground state a (d− 1)-dimensional spatial lattice of N d−1 oscillators to the analogous results for holographic complexity. We could easily identify the leading contribution in the limit of large N and a small UV cut-off distance, i.e., mδ 1. In particular, this led us to consider the generalized family of κ cost functions given in eq. (4.37), which yield Now to identify the leading contribution to either eq. (4.32) or (E.1), we made the crude approximation of replacingω k ∼ 1/δ for all momenta. In the following, we would like to avoid this approximation and examine the complexity (E.1) in more detail. Our result (E.10) is still an approximation but it allows us to consider the subleading contributions to eq. (4.41). In particular, we will determine the leading corrections involving the mass. First we substitute the normal-mode frequencies (4.33) into the above expression, for general κ, to find Now it is certainly true that the second term in the argument of the logarithm dominates for most of the terms in the sum over all momenta. But we would like to be more careful in retaining the leading corrections arising from the mass term. To simplify the analysis below, we will assume that the reference frequency is an IR frequency with ω 0 < m. As a first step, we should isolate the infrared contributions which come from terms in the sum where the mass actually dominates or is comparable to the momentum term in the argument. For large N , we can take the usual continuum limit for these contributions with p i = 2πk i /(N δ). The IR contribution in eq. (E.2) becomes We have dropped the absolute value symbol here because we are assuming that ω 0 < m. The cut-off Λ IR in this integral is an IR scale which delineates the boundary of the IR contributions to the momentum sum in eq. (E.2). Implicitly, we are also letting k i and p i range over positive and negative values so that all of the IR contributions come in the vicinity of k = 0 -see footnote 22. Choosing this cut-off with Λ IR ∝ m, the IR contribution takes the general form where the numerical coefficients c a are independent of m and ω 0 , but will depend of the space time dimension d. The leading contribution then takes the form: C IR c κ V m d−1 [log m/ω 0 ] κ . Having isolated the IR contribution, we return to the UV contributions to eq. (E.2). In these remaining terms, we can consider m 2 /ω 2 0 to be a small correction to the argument of the logarithm, and so we make a Taylor series expansion and keep only the first correction in mδ, i.e., where the notation {k i }>IR indicates that the summation begins at the IR cutoff, i.e., | k| ≥ Λ IR N δ/(2π). The leading term will produce the expected leading contribution in eq. (4.41) with some overall numerical coefficient which depends on κ and d. Of course, there will also be a subleading dependence on our IR cut-off Λ IR ∝ m.
To go further, we focus on the case κ = 1, which simplifies the calculations slightly and also was the case that we found best emulated the holographic complexity. The following analysis is essentially unchanged for larger values of κ. Substituting κ = 1 into eq. (E.5) yields Since the first term is independent of k, the corresponding sum over the UV modes yields a factor of N d−1 (up to corrections proportional to the Λ d−1 IR ). Summing over the second term is more complicated but numerical fits for a range of N and d suggest that the sum takes the form 1 4 where the a i are fixed numerical coefficients. Note that the constant term a 0 appears of odd or even d, and there may also be a logarithmic correction (i.e., log N ) for odd d. Now turning to the second sum in eq. (E.6), we again carried out numerical fits to show 1 4 Collecting the results, our approximation to the total complexity for κ = 1 (and assuming ω 0 < m) becomes To make a comparison to holographic complexity, as in section 4.1, we substitute N d−1 = V /δ d−1 and we also introduce L = V 1/(d−1) as the linear size of the lattice. The complexity (E.10) then becomes +V m d−1 (c 1 log m/ω 0 + c 0 ) . (E.11) Hence we find the expected leading term, which corresponds to the result in eq. (4.41) with κ = 1. We also find a subleading term proportional to V /δ d−1 , as is found in holographic complexity with CA approach [19]. Next we would highlight the correction proportional m 2 V /δ d−3 , for which analogous results can again be found in holographic calculations -see section 4.1 for further comments. It is interesting that we also see corrections, e.g., of the form V /(L 2 δ d−3 ). Of course, this term is far more suppressed that the previous one, but it also involves a fractional power of the volume, i.e., V d−3 d−1 . Such fractional powers would never arise in holographic complexity.
To make the above formulae more concrete, consider the case of a one-dimensional lattice (d = 2), in which case eq. (E.11) reduces to where we have replaced V = L to emphasize that the volume is only a linear length here.

F Optimal geodesic for penalized geometry
We would like to find the optimal geodesic in the penalized geometry (5.3) but as commented below eq. (5.7), finding the general solution for geodesics satisfying the desired boundary conditions seems out of reach. Recall that we were able to show that the simple straight-line geodesic describing the optimal circuit (3.50) in the unpenalized geometry remains a geodesic in our new penalized geometry. However, it was also easy to show that the segmented path described by eq. (5.12) was shorter than this geodesic when the penalty factor was large, i.e., a ρ 1 , y 1 -see eq. (5.14).
To make progress towards finding the optimal geodesic in the new geometry, we make a simplifying assumption. To begin, we examine the penalized metric (5.3) and observe that as the radius ρ increases, the fastest growing component of the metric is g zz ∼ a 2 e 4ρ (for generic x but this component still grows as g zz ∼ e 4ρ for x = π/2). This suggests that motion in the z direction will quickly be suppressed as the geodesics move out from the origin. Therefore we simplify our problem by considering motions on a constant-z submanifold ds 2 = 2dy 2 + 2 a 2 − a 2 − 1 sin 2 (2x) dρ 2 + 2a 2 dx 2 . (F.1) The particular value of z in question will be fixed below by the boundary condition that z = x at s = 0. We will return to justify our assumption of no (or little) motion in the z direction at the end of the appendix. Working with this simpler geometry (F.1), the analysis of the geodesics becomes much more tractable. First, we observe that both ∂ y and ∂ ρ are now Killing vectors, for which the associated conserved quantities are 2c 1 ≡ 2ẏ , 2ac 2 ≡ 2 a 2 − a 2 − 1 sin 2 (2x) ρ , where the factors of 2 and a were chosen to simplify expressions below. Further we have usedc i to avoid confusion with theĉ i in eq. (5.4). Again the first constraint gives the usual solution (3.29) for y, i.e., y(s) = y 1 s withc 1 = y 1 . The second constraint yieldṡ ρ = ac 2 a 2 − (a 2 − 1) sin 2 (2x) .
The normalization of the tangent vector then becomes k 2 = 2ẏ 2 + 2 a 2 − a 2 − 1 sin 2 (2x) ρ 2 + 2a 2ẋ2 where g(x) ≡ 1 sinh −1 (a cot(2x)) , (F. 12) and F is the incomplete elliptic integral of the first kind, Note that in eq. (F.11), we have fixed the integration constant with the boundary condition: ρ(x 1 = π/2) = ρ 1 . Now, unfortunately (F.7) cannot be inverted to find an analytical expression for x(s), which we could then use to obtain ρ(s) via eq. (F.11). However, we can still study the behaviour of these geodesics numerically. Before doing so, it remains to relate the parameters c 2 and k (ork) to the boundary values ρ 1 , and x 0 (as well as x 1 = π/2). Let us first examine the parameter range for which we obtain a real result. It turns out that F in eq. (F.11) is always complex; hence in order for ρ(x) to be real, the coefficient must also be imaginary, which requiresk 2 >c 2 2 =⇒ k 2 > 2 y 2 1 +c 2 2 . (F.14) Turning now to s(x) in eq. (F.7), the elliptic integral Π in this case is always real and the coefficient will also be real in precisely the same regime (F.14). Therefore, this is the only restriction on our parameters required to ensure a real result. We have fixed the integration constants in both s(x) and ρ(x) via the boundary conditions at the end-point of the geodesic that s(x = π/2) = 1 and ρ(x = π/2) = ρ 1 . For the optimal geodesic, we choose the boundary condition x(s = 0) = π/4, which minimizes the cost of motion in the ρ direction -recall eq. (5.11). 54 However we must be careful in evaluating eqs. (F.7) and (F.11) at this value of x with the limits: where F and K are the complete elliptic integrals of the first and third kind, respectively, defined with K(z) = F (π/2 | z) , Π(n | m) = Π (n; π/2 | m) . (F.16) The parametersc 2 andk must be chosen so that both these limits vanish, since initially we must have s = 0 and ρ 0 = 0. In principle we have two equations and two unknowns, but in practice the elliptic integrals are intractable. Fortunately, for our purposes a general solution 54 Alternatively, we could choose x(s = 0) = 3π/4 but the resulting trajectory is simply of a copy of the geodesic found here rotated 180 o around the (θ, τ ) = (π, 0) axis -see figure 7.
where in the second approximation we have performed an expansion in the limit a → ∞.
Additionally, it will be interesting to compare these geodesics against the segmented trajectory described in eq. (5.12). The length of this path is given by eq. (5.13) and so using eq. (F.5), we define 2k 2 s = k 2 s − 2y 2 1 , k s = 1 4 π 2 a 2 + 8ρ 1 2ρ 1 + π 2 a 2 + (4y 1 ) 2 (F. 23) were in the second line we have replaced ρ 1 using eq. (F.18). Note that unlikek andk 0 in eq. (F.20), the parameter y 1 still appears in this expression -although this contribution is suppressed for a y 1 . Again however, the second line above is more suited to numerics than physical inspection; to compare with (F.22), we shall expand with a ρ 1 , y 1 (as well as a, ρ 1 , y 1 1 and where we assume ρ 1 and y 1 to be roughly the same order of magnitude We mentioned above that the segmented path constitutes a remarkably good approximation to the geodesic. Comparingk in (F.22) andk s in (F.24), one can see evidence for this claim in that the leading-order behaviours are precisely the same; deviations arise only in the subleading terms, which are increasingly negligible for large values of a. We discuss this point further in the main text -see eq. (5.20). We also explicitly confirm that the two paths are very close to one another in the large a regime by examining x(s) and ρ(s) numerically, e.g., see figure 5.
In closing the appendix, we remind the reader that to make progress we confined our attention to motion in a constant-z subspace, with the simpler metric (F.1). Hence we should go back to examine whether or not this was a reasonable assumption. In particular, we wish to argue that, at least in the limit a 1, the particular class of geodesics with x 0 = π/4 and x 1 = π/2 obtained for the constant-z subspace are a good approximation to the corresponding geodesics in the full geometry (5.3). Intuitively, we motivated this restriction by the observation that movement in the z-direction is relatively costly. We can quantify this by considering the behaviour ofż, given in eq. (5.5). Recall that τ 0 = 0, and hence z 0 = x 0 = π/4. Then the finiteness condition (5.6) requires that we setĉ 3 = 0. 55 Now along the initial segment, where x = π/4, the derivatives (5.5) reduce tȯ .
Therefore, in the large a limit that we are considering, motion in both the x-and z-directions is highly suppressed, while only motion along ρ is inexpensive. Along the second segment, where we rotate around to x = π/2, bothẋ andż pick up terms which are O(1) in a, butż is still exponentially suppressed in ρ 1 relative toẋ. (Meanwhileρ decreases sharply to 0 on this segment, in the limit a 1.) Thus the geodesics in the full spacetime (5.3) can indeed be approximated by those in the constant-z subspace (F.1), at least in the limits that we are considering.