Circuit complexity in quantum field theory
Ro Jefferson Robert C. Myers
hep-th quant-ph
  
Polarization-Division Multiplexing Based on the Nonlinear Fourier Transform
    Jan-Willem Goossens,1,2,* Mansoor I. Yousefi,2 Yves Jaouën2 and Hartmut Hafermann1
    Received ***date***; Accepted ***date***
======================================================================================





§ INTRODUCTION

Recent years have seen exciting progress in understanding the connection between entanglement and geometry <cit.>. 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 <cit.>
|TFD(t_L,t_R)⟩=1/√(Z_β)∑_ie^-βE_i/2 e^-iE_i(t_L+t_R) |i⟩_L |i⟩_R .
This describes an entangled state of the two copies of the CFT associated with the asymptotic boundaries (see figure <ref>), which are joined by a wormhole, an Einstein-Rosen bridge (ERB), in the bulk <cit.>. 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 <cit.>. 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 <cit.> – 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 <cit.>. 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 complexity=volume (CV) conjecture, which posits that the complexity of the boundary state is proportional to the volume of a maximal codimension-one bulk surface ℬ that extends to the AdS boundary, and asymptotes to the time slice Σ on which the boundary state is defined <cit.>:
𝒞_V(Σ) = max_Σ=∂ℬ [𝒱(ℬ)/G_N   ℓ]   ,

where ℓ is some length scale associated with the bulk geometry, 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 <ref>. 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 <cit.>:
𝒞_A(Σ)=I_WDW/π ħ .

One can think of the WDW patch as the causal development of the spacelike surface ℬ picked out by the CV construction. The right panel in figure <ref> 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, <cit.>. 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, as was recently found for holographic entanglement entropy <cit.>. 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).[We also refer the reader to ref. <cit.> for a recent complementary investigation in this direction.] 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<cit.>, or more generally into Hamiltonian complexity <cit.>, or the efficient description of many-body wave functions <cit.>.

In computer science, the notion of computational complexity refers to the minimum number of operations necessary to implement a given task <cit.>. 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⟩ .
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, ||  |ψ_T⟩-U |ψ_R⟩ ||^2≤ . 
Of course, there will not be a unique circuit which implements the desired transformation circuitDef: 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 (<ref>), 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 <cit.>,[See <cit.> for another application of Nielsen's ideas in holography. We also refer the interested reader to ref. <cit.>, which introduces an interesting connection between quantum algorithms and geodesics on the Fubini-Study metric.] which itself was developed using ideas from the theory of optimal quantum  control, <cit.>. 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, 
U=P exp[∫_0^1dt H(t) ]    where  
H(t)=∑_I Y^I(t) M_I ,
where the Hamiltonian is expanded in terms of generalized Pauli matrices, denoted here as M_I,[Our notation diverges from that of Neilsen, in order to increase the similarity of these equations with our notation in the main text. In particular, note that we have absorbed a factor of -i in M_I so that these are now anti-Hermitian operators.] and the P indicates a time ordering such that the Hamiltonian at earlier times is applied to the state first, the circuit is built from right to left. In <cit.>, 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,
U(t)=P exp[∫_0^t t̃ H(t̃) ] .

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))=∫_0^1t FU(t),U̇(t) ,
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: 
0.8ex

   1.  Continuity: F should be continuous, F∈ C^0.
0.8ex

   2.  Positivity: F(U,v)≥0 with equality if and only if v=0. 
0.8ex

   3.  Positive homogeneity: F(U,λ v)= λ  F(U,v) for any positive real number λ. 
0.8ex

   4.  Triangle inequality: F(U,v+v')≤ F(U,v)+F(U,v') for all tangent vectors v and v'. 
0.8ex

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
0.8ex

   
1'.  Smoothness: F should be smooth, F∈ C^∞,
0.8ex

then eq. costco 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 costco, with a Finsler metric F satisfying the four properties enumerated above, see <cit.>. 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, geodesics, in a Finsler geometry, and the complexity is then identified with the length of the geodesic.[For future reference, when referring to general paths or circuits, we will use “size,” “length,” “cost,” and “depth” interchangeably; however, “complexity” will be reserved for the length of the optimal path or circuit.] 

Of course, this still leaves open the question of the precise form of the cost function, and various possibilities are examined in <cit.>:[The functions F_1 and F_p are not technically Finsler metrics, since both fail to meet the smoothness requirement. However, as explained in <cit.>, they can be approximated arbitrarily well by metrics which are Finsler. This subtlety will not be important for our analysis.] 
F_1(U,Y)=   ∑_I |Y^I| ,        F_p(U,Y)   =   ∑_I p_I |Y^I| ,

F_2(U,Y)=   √(∑_I Y^I^2) ,      
F_q(U,Y)   =   √(∑_I q_IY^I^2) .

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, 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 <ref>. 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 <cit.> 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. costco.

This paper is organized as follows: we begin in section <ref> 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 <ref>. The main ideas from Nielsen's approach are implemented here: we represent the circuit as a path-ordered exponential analogous to eq. row, show that our space of circuits forms a representation of , 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 <ref>, 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 <ref>, we conduct a preliminary exploration of the effects of introducing penalty factors for nonlocal gates. Finally, we close in section <ref> 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 <ref> in appendix <ref>, elaborate on some geometrical details in appendix <ref>, derive the normal-mode frequencies for a one-dimensional lattice in appendix <ref>, find a closed-form approximation to the circuit complexity for the d-dimensional lattice in appendix <ref>, and compute an approximation to the optimal circuit in the presence of penalty factors in appendix <ref>.



§ 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,[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 <cit.>.] which reduces the system to an infinite family of harmonic oscillators. This in turn suggests the much simpler warm-up 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,
H=1/2∫^d-1x[π(x)^2+∇⃗ϕ(x)^2+m^2ϕ(x)^2] .
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:
H=1/2∑_{p()^2/δ^d-1+δ^d-1[1/δ^2∑_iϕ()-ϕ(-x̂_i)^2+m^2ϕ()^2]} ,
where x̂_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()=δ^d/2ϕ(), P()=p()/δ^d/2, M=1/δ, ω=m and Ω=1/δ, whereupon the Hamiltonian eq:Hlattice takes the familiar form
H=∑_{P()^2/2M+1/2 M [ω^2 X()^2+Ω^2∑_iX()-X(-x̂_i)^2]} .
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:
H=1/2[ p_1^2+p_2^2+ω^2x_1^2+x_2^2+Ω^2x_1-x_2^2] ,

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,
H=1/2p̃_+^2+ω̃_+^2x̃_+^2+p̃_-^2+ω̃_-^2x̃_-^2 ,

where[When working in the normal-mode basis, we denote variables (positions, frequencies), with a tilde to clearly distinguish from the physical basis. The utility of this convention will become apparent later. ]
x̃_±≡1/√(2)x_1±x_2 ,    ω̃_+^2=ω^2 ,    ω̃_-^2=ω^2+2Ω^2 .

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:
ψ_0(x̃_+,x̃_-)=ψ_0+(x̃_+)ψ_0-(x̃_-)
=ω̃_+ω̃_-^1/4/√(π) exp[-1/2ω̃_+x̃_+^2+ω̃_-x̃_-^2] ,
where the normalization has been chosen such that ∫ d^2x |ψ_0|^2=1. We may also express this wave function in terms of the physical positions of the two masses:
ψ_0(x_1,x_2)=ω_1ω_2-β^2^1/4/√(π)exp[-ω_1/2x_1^2-ω_2/2x_2^2-βx_1x_2] ,
where
ω_1=ω_2=1/2ω̃_++ω̃_- ,   
β≡1/2ω̃_+-ω̃_-<0 .
We note in passing that our notation for the wave function in eq. eq:targetPhys is slightly more general than necessary; however, these Gaussian wave functions constitute an interesting family of target states for the present exercise.[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, <cit.>.]

The next step is to identify a simple reference state. Motivated by discussions of holographic complexity <cit.>, as well as cMERA <cit.>, we choose a reference state where the two masses are unentangled, namely a factorized Gaussian state,
ψ_R(x_1,x_2)=√(ω_0/π) exp[-ω_0/2x_1^2+x_2^2] .
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 <ref>. 

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:
     H=e^ix_0p_0 ,     J_a=e^ix_0p_a ,    K_a=e^ix_ap_0 ,  

Q_ab   =e^ix_ap_b  (with ab) ,     Q_aa=e^i/2x_ap_a+p_ax_a=e^/2 e^ix_ap_a ,
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:
H ψ(x_1,x_2)   =e^ip_0x_0ψ(x_1,x_2)            (global) phase change

J_1 ψ(x_1,x_2)   =ψ(x_1+x_0,x_2)            shift  x_1  by constant x_0

K_1 ψ(x_1,x_2)   =e^ip_0x_1ψ(x_1,x_2)            shift  p_1  by constant p_0

Q_21 ψ(x_1,x_2)   =ψ(x_1+x_2,x_2)            shift  x_1  by  x_2     (entangling gate)

Q_11 ψ(x_1,x_2)   =e^/2ψe^x_1,x_2           scale  x_1 →e^x_1     (scaling gate)

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 (<ref>). 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. eq:gates with operators like
exp[i p_0/x_0  x_1x_2]   or  exp[i x_0/p_0  p_1^2] .

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. eq:gates is sufficient to implement the unitary transformation from the specified reference state eq:refPhys to the desired target state eq:targetPhys. 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:
ψ_T=Uψ_R≡Q_22^α_3 Q_21^α_2 Q_11^α_1 ψ_R .
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_1^2 in eq. eq:refPhys to the desired frequency ω_1 appearing in eq. eq:targetPhys. 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 <ref>, and the final result is 
      α_1=1/2logω_1/ω_0 ,
    α_2=1/√(ω_0/ω_1)β/√(ω_1ω_2-β^2) ,


                α_3=1/2logω_1ω_2-β^2/ω_0 ω_1 .



We then define the circuit depth as the total number of gates in the circuit. In the above example, we have simply
𝒟(U)   =   |α_1|+|α_2|+|α_3|


   =   1/[1/2 logω_1ω_2-β^2/ω_0^2+√(ω_0/ω_1)|β|/√(ω_1ω_2-β^2)] .

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, β<0. But this means that we are giving an equal complexity cost for the inverse gates Q_ij^-1 as for the original gates Q_ij, we count the appearance of Q_ij^-1 as one gate in a circuit. 

We refer to the result in eq. leach as the circuit depth of the particular circuit U given in eq. eq:eg1. 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. eq:eg1 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. leach 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 <ref>. 

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. leach in terms of the normal-mode frequencies using eq. eq:omega12pm. This substitution yields
𝒟_1=1/[ 1/2log_+/ω_0+1/2log_-/ω_0+_--_+/√(2_+_-)√(ω_0/_++_-) ] . 
Recall that _->_+ from eq. qm3 (and implicitly, we are assuming ω̃_±>ω_0). 





§ GEOMETRIZING COMPLEXITY


In the introduction, we discussed Neilsen's approach <cit.> 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,
U=𝒫 exp∫_0^1s Y^I(s) _I ,    ψ_Tx_1,x_2=Uψ_Rx_1,x_2 .
This structure replaces the representation of the circuits as products of the discrete gates in eq. eq:gates. The connection with these gates comes about since we choose the operators _I appearing in the exponential to be precisely those appearing in the scaling and entangling gates introduced previously; that is, we write 
Q_ab = exp[ _ab]
    with    _ab=(i x_a p_b+1/2 δ_ab) .

Our notation in eq. eq:pathPsi is that the sum over I runs over the pairs ab, 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 (<ref>)). In the integral appearing in the exponent, the differential s 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, 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
U(s)=𝒫 exp∫_0^ss̃ Y^I(s̃) _I .
Eq. eq:pathPsi then specifies the final unitary at the end-point s=1, which corresponds to the desired circuit that generates the target state, 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 <ref> above, the circuit depth is determined by counting the total number of gates appearing in the full sequence comprising the circuit, eq. leach. For our path-ordered exponential eq:pathPsi, the analogous expression becomes[Actually this expression cost1 is the continuum limit of the cost function 𝒟(U)=∑ϵ |α_i|. Including the extra factor of ϵ in the sum eliminates the 1/ϵ factor, so the circuit depth remains finite in the limit ϵ→0.]
𝒟(U)=
∫_0^1s∑_I|Y^I(s)|=
∫_0^1s[|Y^11(s)|+|Y^12(s)|+|Y^21(s)|+|Y^22(s)|] .

This cost function corresponds to the F_1 metric in the notation of <cit.> — see eq. (<ref>). 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. cost1.

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. (<ref>). That is, we replace eq. cost1 with
𝒟(U)=∫_0^1s√(G_IJ Y^I(s) Y^J(s)) .
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, 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 <ref>.

To proceed further, we must find a prescription to explicitly identify the functions Y^I(s). Given eq. path2, it is straightforward to show that
Y^I(s) _I=_sU(s) U^-1(s) .
However, this expression is not particularly useful as it stands. In Neilsen's construction <cit.>, one works with unitary matrices acting on qubits, rather than operators acting on wave functions. Hence the components of the velocity analagous to eq. v12 can be isolated by simply tracing over the corresponding matrix generators. This procedure does not immediately lend itself to eq. v12, 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 (<ref>) of two coupled harmonic oscillators, starting from a factorized Gaussian reference state (<ref>). 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, all of the intermediate wave functions take a form analogous to eq. (<ref>). 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
ψ≃exp[-1/2x_a A_ab x_b] ,

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.[These positivity constraints ensure that both eigenvalues of A_ab are positive.] In particular, the reference and target states become, respectively,
A_R=ω_01 ,    A_T=(

ω_1   β

β   ω_2

) ,

where ω_1, ω_2 and β are given by eq. eq:omega12pm.

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 
A' = Q_ab  A   Q_ab^T ,

where
Q_ab = exp[ M_ab]
    with  [M_ab]_cd = δ_ac δ_bd .

In this notation, [M_ab]_cd is a 2×2 matrix, where c and d denote row and column indices, respectively.[A quick way to construct these matrices is to consider the action of _ab on the column vector x_1,x_2^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 _ab. Note that, while the action of the Q_ab in eq. operate leaves the wave functions properly normalized at each step, we lose track of this normalization when working with the A_ab.  ] Explicitly, we shall denote the basis of generators M_I as

M_11   =[ 1 0; 0 0 ] ,    
M_12=[ 0 1; 0 0 ] ,    


M_21   =[ 0 0; 1 0 ] ,    
M_22=[ 0 0; 0 1 ] .



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 .[Note that one can also see the emergence of this group by observing that the algebra of the original operator generators _ab in eq. operate close to form the algebra gl(2,ℝ).] Thus our circuits form a representation of , the U(s) are trajectories in the space of  transformations. 

Now, in this matrix formulation, the path-ordered exponentials in eq. (<ref>) are replaced by
U(s)=𝒫 exp∫_0^ss̃ Y^I(s̃) M_I ,    with  
A_T=U(s=1) A_R U^T(s=1) ,
where M_I are the generators given in eq. Msimple. The advantage of this formulation is that eq. v12 becomes
Y^I(s) M_I=_sU(s) U^-1(s)
    Y^I(s)=_sU(s) U^-1(s)M^T_I .
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  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  that connects the initial and final states, A_R and A_T, as in eq. pathA. 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=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 cost2.

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,ℝ). 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 <ref>.



 §.§ Geodesics on circuit space
 
To proceed with constructing the desired geodesics, we must choose an explicit parametrization of a general element U∈= ℝ×SL(2,ℝ). Let us first consider U∈SL(2,ℝ), which can be written as
U=[ x_0-x_3 x_2-x_1; x_2+x_1 x_0+x_3 ] ,
    with    x_0^2+x_1^2-x_2^2-x_3^2=1 .
We recognize the constraint imposing detU=1 as the embedding of (Lorentzian) AdS_3 in ℝ^2,2. Indeed, the appearance of AdS_3 could have been anticipated since the latter is the universal cover of SL(2,ℝ). 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ρ , 
where τ, ρ and θ are the usual time, radius, and angle, respectively, of global coordinates on AdS_3. We can easily extend this parametrization to U∈= ℝ×SL(2,ℝ) by introducing an additional coordinate to parameterize the determinant of U, U=[ x_0-x_3 x_2-x_1; x_2+x_1 x_0+x_3 ] ,
    with    x_0^2+x_1^2-x_2^2-x_3^2=e^2y .
Hence we extend eq. slr to
x_0=e^ycosτcoshρ ,   
x_1=e^ysinτcoshρ ,   
x_2=e^ycosθsinhρ ,   
x_3=e^ysinθsinhρ , 
where, as before, τ, ρ, θ are coordinates on the SL(2,ℝ) subgroup, and y parametrizes the ℝ fibre. With these coordinates, we can express a general U∈ as
U
=e^y [  cosτcoshρ-sinθsinhρ -sinτcoshρ+cosθsinhρ;  sinτcoshρ+cosθsinhρ  cosτcoshρ+sinθsinhρ ] .


We are now equipped to construct the geometry implicit in the cost function cost2, where the velocity components are given by eq. (<ref>). 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:
s^2   =δ_IJ   U U^-1  M^T_I U U^-1 M^T_J

   =2y^2+2ρ^2+2cosh(2ρ)cosh^2ρ τ^2+2cosh(2ρ)sinh^2ρ θ^2-2sinh^22ρ τθ .

For later use, it is also convenient to express this in the form
s^2=2y^2+2ρ^2+2x^2+2cosh(4ρ) z^2-4cosh(2ρ)
 x z ,

where we have defined the pseudo-lightcone coordinates
x≡1/2(θ+τ) ,    z≡1/2(θ-τ) .
Note that our metric metric1 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 <ref>.

With the geometry in hand, we now wish to find the geodesics, and thereby the optimal circuit. Inspecting the metric metric1, we can see three obvious Killing vectors: _y, _τ, _θ. However, the metric is right-invariant by construction, meaning eq. metric1 remains unchanged if we right-multiply U(s) by a constant  transformation. Therefore there must be one Killing vector for each generator of , namely, four.[We thank Lucas Hackl for discussions on this point.] 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 <ref>, and are given in eqs. (<ref>) and (<ref>).

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. (<ref>) and (<ref>), it is straightforward to evaluate the corresponding conserved quantities:
c_1   =2 ẏ ,

c_2   =2 sin(θ-τ)ρ̇+cos(θ-τ) [sinh(4ρ)-sinh(2ρ)θ̇-sinh(4ρ)+sinh(2ρ)τ̇] ,

c_3   =2cos(θ-τ)ρ̇-sin(θ-τ) [sinh(4 ρ)-sinh(2 ρ)θ̇-sinh(4 ρ)+sinh(2 ρ)τ̇] ,

c_4   =cosh4ρ-cosh2ρθ̇-cosh4ρ+cosh2ρτ̇ ,

c_5   =(1-cosh2ρ) θ̇+(1+cosh2ρ) τ̇ ,

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, g_ijẋ^iẋ^j
=2ẏ^2+2ρ̇^2+2cosh(2ρ)sinh^2ρ θ̇^2+cosh^2ρ τ̇^2-2sinh^2(2ρ) θ̇ τ̇≡k^2 .
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, eq:pathPsi. Hence we shall scale the affine parameter s to lie in this range. The normalization constant k then gives the length of the geodesic, the depth of the corresponding circuit, since from eq. cost2 we have
D(U)=∫_0^1s√(g_ij ẋ^i ẋ^j)≡k .

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 (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. (<ref>), we find that all coordinates except θ are initially zero, (s=0)={0, 0, θ_0, 0} .
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. pathA, where the quadratic forms for the reference and target states are given in eq. eq:stateMatrix. Substituting the initial state A_R=ω_0 1 and the explicit representation of the unitaries (<ref>) into this relation, we have
A_T=ω_0 UU^T=
ω_0  e^2y_1[ cosh(2ρ_1)-sin(θ_1+τ_1)sinh(2ρ_1)            cos(θ_1+τ_1)sinh(2ρ_1);            cos(θ_1+τ_1)sinh(2ρ_1) cosh(2ρ_1)+sin(θ_1+τ_1)sinh(2ρ_1) ] ,

where the subscript 1 denotes the value of the coordinate at s=1, y_1=y(s=1). Comparing the entries of the matrix on the right-hand side with those of A_T in eq. eq:stateMatrix, we arrive at the following
boundary conditions for the end of the geodesic:
ω_1/ω_0   =e^2y_1[cosh(2ρ_1)-sin(θ_1+τ_1) sinh(2ρ_1)] ,

ω_2/ω_0   =e^2y_1[cosh(2ρ_1)+sin(θ_1+τ_1) sinh(2ρ_1)] ,

β/ω_0   =e^2y_1cos(θ_1+τ_1) sinh(2ρ_1) .

Implicitly, these constraints allow us to identify the final coordinates (s=1) for the geodesics corresponding to circuits which produce the desired transformation. Explicitly, we may solve this system to obtain
e^2y_1=√(ω_1ω_2-β^2)/ω_0 ,  cosh(2ρ_1)=ω_1+ω_2/2√(ω_1ω_2-β^2) ,  tan(θ_1+τ_1)=ω_2-ω_1/2β .

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. fini, 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 <ref>, 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 (<ref>). 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 .

Next, we consider c_4 and c_5. These two constraints may be solved to obtain
τ̇=c_5+c_4-c_5/4cosh^2ρ ,    θ̇=c_5+c_4+c_5/4sinh^2ρ .

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, 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 
c_5=Δθ≡θ_1-θ_0
        θ(s)=Δθ s+θ_0 .

Furthermore, the above allows us to simplify the τ̇ equation to
τ̇=Δθ1-1/2cosh^2ρ .
Now, combining our expressions for θ̇ and τ̇ with the constraints c_2 and c_3 in eq. eq:conserved, we find a relatively simple equation for ρ̇:
ρ̇^2=c_2^2+c_3^2/4-Δθ^2/4 tanh^2ρ .


In principle, we should now solve for the general solutions of eqs. tauDot4 and rhoDot4 subject to the boundary conditions in eqs. eq:coordinit and fini. While it is possible to carry out this exercise, the final solutions are not particularly illuminating.[The general solution for eq. rhoDot4 is given by
sinhρ=c/√(c^2-Δθ^2) sinhs/2√(c^2-Δθ^2) ,

where c^2=c_2^2+c_3^2 is fixed by substituting the boundary condition ρ=ρ_1 at s=1 into this equation. Furthermore, given this result, it is possible to integrate eq. tauDot4 to obtain τ(s); one finds
τ=Δθ s-tan^-1Δθ/√(c^2-Δθ^2)tanhs/2√(c^2-Δθ^2) .
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. fini, we can then determine θ_1. In turn, θ_0 is now fixed since we know Δθ and θ_1. ] Instead, let us point out the particularly simple solution that arises for Δθ=0. In this case the expressions for τ̇ and ρ̇ reduce to
τ̇=0         τ=0 ,

ρ̇=1/2√(c_2^2+c_3^2)       ρ=ρ_1 s ,

which combine with y=y_1 s and θ=θ_0 from eqs. yrun and gonzo to describe a simple “straight-line” geodesic. Substituting this solution into eq. Umatrix, we can write the corresponding circuit as
U_0(s)   =e^y_1s[ coshρ_1s- sinθ_0 sinhρ_1s           cosθ_0 sinhρ_1s;           cosθ_0 sinhρ_1s  coshρ_1s+sinθ_0 sinhρ_1s ]

   =exp[[ 1 0; 0 1 ] y_1 s+[ -sinθ_0  cosθ_0;  cosθ_0  sinθ_0 ] ρ_1 s] .

Note that an explicit path-ordering is not needed in the second expression since it is simply the exponential of a fixed matrix.[For this simple case, it is straightforward to identify the exponential form in the second line of eq. (<ref>) given the expression appearing in the first. In general however, one would apply eq. eq:metricRightInv to identify the components of Y^I(s) and then substitute these into eq. (<ref>).] The circuit depth of U_0, the length of the geodesic, is given by eqs. eq:geodesic and costa4, which for this simple solution yields
D(U_0)=√(2(y_1^2+ρ_1^2)) .

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.

To prove that the straight-line solution above is the geodesic whose length is the (global) minimum, we recall from eq. eq:geodesic 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. yrun and gonzo, respectively, as well as the expression for τ̇ from eq. tauDot4, whereupon we find
k^2
=2y_1^2+2ρ̇^2+1-1/2cosh^2ρ Δθ^2 .
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
∫_0^1s  ρ̇=ρ_1 . 
Additionally, we have
0≤∫_0^1s  (ρ̇-ρ_1)^2=∫_0^1s ρ̇^2 -ρ_1^2  ,

and hence we may conclude that ∫_0^1 s ρ̇^2≥ρ_1^2, and that this inequality is only saturated when ρ̇=ρ_1 along the entire geodesic. Now, examining the coefficient of Δθ^2 in eqn. eq:geodesic3, we have
1/2≤1-1/2cosh^2ρ≤1 ,

where the lower inequality is only saturated at ρ=0, and the upper inequality is saturated at ρ→∞.[In accordance with its interpretation as a radial coordinate, we do not consider negative values of ρ.] Given that all of our geodesics must start at ρ=0 and end at ρ=ρ_1, upon averaging over any of these geodesics, we find
1/2<∫_0^1 s1-1/2cosh^2ρ < 1 .

Finally, let us average eq. eq:geodesic3 over any of our geodesics:
k^2   =2y_1^2+2∫_0^1s ρ̇^2 + Δθ^2 ∫_0^1 s1-1/2cosh^2ρ

   ≥2y_1^2+2ρ_1^2 +Δθ^2/2
≥2(y_1^2+ρ_1^2) .

Comparing this to eq. solver4b, we have established the inequality k≥ D(U_0). Furthermore, our argument has established that this inequality can only be saturated with ρ̇(s)=ρ_1 (ρ=ρ_1 s) and Δθ=0 (θ(s)=θ_0, which implies τ(s)=0 via eq. tauDot4). We have therefore proved that the simple straight-line geodesic indeed constitutes the global minimum for our cost function costa4, and hence that eq. solver4b is in fact the complexity of the Gaussian wave function in this framework:
(A_T)=√(2y_1^2+ρ_1^2) .


As an exercise, we can compare the above result for D(U_0) in eq. solver4b, which was evaluated using eq. costa4, with the result found by evaluating eq. cost2. In this case, we must identify the components Y^I(s), which is easily done by examining the exponential expression in eq. solver4a:[Recall that our  generators are given in eq. Msimple.]
Y^11   =   y_1-ρ_1 sinθ_1 ,  Y^22=y_1+ρ_1 sinθ_1 ,


          Y^12=Y^21=ρ_1 cosθ_1 .

Since these components are all constant, the integral over s in eq. cost2 is trivial, and the circuit depth (with G_IJ=δ_IJ) reduces to
D(U_0)   =   √((Y^11)^2 +(Y^12)^2+(Y^21)^2+(Y^22)^2)


   =   √(
y_1-ρ_1sinθ_1^2+2ρ_1cosθ_1^2
+y_1+ρ_1sinθ_1^2)
=√(2y_1^2+ρ_1^2) ,
in agreement with eq. (<ref>).



 §.§ Normal-mode subspace


To properly interpret the complexity, we must re-express our result solver4b in terms of the physical parameters of the two coupled oscillators qm1, as well as the frequency ω_0 in the reference state eq:omega12pm. However, one finds that the complexity is most elegantly described in terms of the normal-mode frequencies _+ and _- given in eqs. qm2 and qm3. Using eq. eq:omega12pm, the final boundary conditions fini simplify to
y_1=1/4 log_+_-/ω_0^2 ,    ρ_1=1/4 log_-/_+ ,    θ_1+τ_1=π .

Substituting these expressions for y_1 and ρ_1 into eq. solver4b then yields the complexity of the ground state,
(A_T)   =D(U_0)=
1/2√(log^2ω̃_+/ω_0+log^2ω̃_-/ω_0) .
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 <ref>. The corresponding circuit solver4a simplifies to
U_0(s)=e^y_1s[  coshρ_1s -sinhρ_1s; -sinhρ_1s  coshρ_1s ]
=exp[[  y_1 -ρ_1; -ρ_1  y_1 ] s
] ,

with y_1 and ρ_1 given by eq. fin2.

The simple and elegant form eq:CpmPre of the complexity in terms of the normal-mode frequencies suggests that we should investigate the optimal circuit solver4a in terms of the normal modes. The relationship between the physical positions of the masses and the normal-mode coordinates was given in eq. qm3, 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,[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. eq:rot. 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. ] 
R=1/√(2)[  1  1;  1 -1 ]
      
[ x̃_+; x̃_- ]=R[ x_1; x_2 ] .
Introducing the short-hand notation x=  x_1, x_2^T and x̃=x̃_+,x̃_-^T, the transformation (<ref>) may be concisely written x̃ =R  x, and the inverse transformation becomes x =R^T x̃. Of course, we can also use this transformation to re-express the target Gaussian wave function in terms of the normal-mode coordinates,
ψ_T∼exp[-1/2 x^TA_T  x]=exp[-1/2  x̃^T R  A_TR^T x̃]    Ã_T=R A_T R^T ,
where Ã_T denotes the quadratic form describing the ground state in the normal-mode space. Explicitly performing this rotation, one finds
Ã_T=[ ω̃_+    0;    0 ω̃_- ] .

That is, the target state becomes a factorized Gaussian in the normal-mode basis, eq. qm2. 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. eq:stateMatrix, we see that it retains its simple form, Ã_R = R  A_R  R^T =ω_0  1 .

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, eq. pathA, we can transform our minimal circuit solver5 to act in the normal-mode space:
Ũ_0(s) ≡R U_0(s) R^T    where   Ã_T=Ũ_0(s=1)
 Ã_R Ũ^T_0(s=1) .
This transformation effects a remarkable simplification of the circuit solver5 to
Ũ_0(s)   =exp[[ y_1-ρ_1       0;       0 y_1+ρ_1 ] s]

   =exp[[ 1/2 log_+/ω_0             0;             0 1/2 log_-/ω_0 ] s]=[ _+/ω_0^s/2          0;          0 _-/ω_0^s/2 ]  ,

where in the second line we have used eq. fin2.

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. Atilde and Atilde2. 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 Ã_R, since if any off-diagonal entries (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 solver4a 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.[At this point, we wish to alert the reader to a subtle distinction that arises in our notation here: as established in footnote <ref>, we have introduced tilde's to distinguish quantities related to the normal modes from similar quantities in the position basis. Beginning with eq. Atilde, a state, circuit, or generator carrying a tilde acts in the normal-mode space, on wave functions written in terms of normal modes. However, this should be distinguished from the instances described here, where we place the tilde's on the indices. These tilded indices indicate that a normal-mode “basis” may still appear on objects acting in the oscillator position space. For example, above eq. bark, M_Ĩ indicates certain linear combinations of the standard generators Msimple, which still act on wave functions written in terms of x_1,x_2, but in a way that scales or entangles the normal modes.] 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 generators M̃_Ĩ with which to construct the circuits acting on the states described in the normal-mode basis are 
M̃_++=[ 1 0; 0 0 ] ,   M̃_+-=[ 0 1; 0 0 ] ,   M̃_-+=[ 0 0; 1 0 ] ,
   M̃_–=[ 0 0; 0 1 ] .

As numerical matrices, these M̃_Ĩ are of course identical to the M_I given in eq. Msimple, but the two sets of generators act in different spaces. Via the transformation eq:rot, we can also transform these generators to act on the states in the original position basis, M_Ĩ=R^T M̃_Ĩ R:
M_++=   1/2[ 1 1; 1 1 ]        = 1/2(M_11+M_22+M_12+M_21 ) ,

M_+-=   1/2[  1 -1;  1 -1 ]     = 1/2(M_11-M_22-M_12+M_21 ) ,


M_-+=   1/2[  1  1; -1 -1 ]    = 1/2(M_11-M_22+M_12-M_21) ,


M_–=   1/2[  1 -1; -1  1 ] 
   = 1/2(M_11+M_22-M_12-M_21) .
The action of these generators can be read off from the indices, M_++ scales the x_+ coordinate or amplifies the corresponding normal mode. Of course, we could also transform the original generators M_I in eq. Msimple with M̃_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, 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. solver5 directly or by transforming the expression in eq. solver5a with U_0(s)=R^T Ũ_0(s) R, that the optimal circuit can be expressed as
U_0(s)=exp[M_++  (y_1-ρ_1)+M_–  (y_1+ρ_1) s]
 ,

where M_±± are the linear combinations of the original generators given in eq. bark. 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
U(y,ρ)=exp[M_++  (y-ρ)+M_–  (y+ρ)]
 ,

then the corresponding metric becomes[This conclusion is slightly premature, since we have not shown that the metric metric1 is invariant under the change of basis from the original generators Msimple to those in eq. bark, but we shall prove this below in eq. metric1a. Note that we have also used that the new basis of generators still satisfies M_Ĩ M_J̃^T=δ_ĨJ̃.]
s_n-m^2   =δ_ĨJ̃   U U^-1  M^T_Ĩ U U^-1 M^T_J̃

   =(y-ρ)^2+(y+ρ)^2=2y^2+2ρ^2 .

Hence we recognize the normal-mode subspace as precisely the θ,τ=π,0 plane in our extended geometry metric2.[Implictly, we may allow ρ to run over positive and negative values in eq. swer4. Hence this subspace also includes θ,τ=0,0.] This perspective also makes clear why the optimal geodesic remains in the normal-mode subspace. Examining the full metric metric2, 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 metric2 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 <ref>). 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. Msimple and Msimple2, 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. bark, but we would like to build an explicit transformation matrix R̂:
M_Ĩ = R_Ĩ J M_J    where    R_Ĩ J =1/2[  1  1  1  1;  1 -1  1 -1;  1  1 -1 -1;  1 -1 -1  1 ]=R_ka⊗R_ℓb .

Note that in the final equality, R is the rotation matrix in eq. eq:rot, and we are identifying the indices as follows: Ĩ=(kℓ) with k,ℓ∈{+,-}, and J=(ab) with a,b∈{1,2}.[ Recall that as defined in eq. eq:rot, R is the matrix which transforms the `1,2' indices of the oscillator position basis to the `+,–' indices of the normal-mode basis—see footnote <ref>.] This identification is really the origin of the interesting tensor product structure R=R⊗ R. The expression in eq. bark2 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̃M_J̃=R_J̃ I M_J̃. Similarly, this transformation acts on the velocity components as Y^I=Y^J̃R_J̃ 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 metric1 as follows:
s^2   =δ_IJ   U U^-1  M^T_I U U^-1 M^T_J

   =R_Ĩ I  R_J̃ J δ_IJ   U U^-1  M^T_Ĩ U U^-1 M^T_J̃

    =δ_Ĩ J̃   U U^-1  M^T_Ĩ U U^-1 M^T_J̃ ,

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 normal-mode subspace in eq. metric5x. 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 <ref>.



§ A LATTICE OF OSCILLATORS


In this section, we wish to return to the original problem of a free scalar field regulated by a lattice, eq:Hlattice. That is, we will consider evaluating the complexity of the ground state of a lattice of coupled oscillators hqm. 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 normal-mode 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,
H=1/2∑_a=0^N-1[  p_a^2+ω^2  x_a^2+Ω^2x_a-x_a+1^2] ,

with periodic boundary conditions x_a+N=x_a.[Note that for convenience, we have labeled the first oscillator with a=0, rather than a=1, the sum in eq. qm88 runs over a∈{0,1,⋯,N-1}.] 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. hqm. The Hamiltonian qm88 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,
H=1/2∑_k=0^N-1[   |p̃_k|^2+_k^2  |x̃_k|^2 ] ,

where the transformation to the normal-mode basis is achieved by a (discrete) Fourier transform,
x̃_k≡1/√(N)∑_a=0^N-1exp-2πi k/N ax_a .
where k∈{0,…,N-1}, and we note that x̃_k^†=x̃_N-k.[We can see this result as a combination of two simpler identities: x̃_k^†=x̃_-k, which follows from the complex conjugation of eq. eq:Fourier, and x̃_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, 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
u⃗^†_ k·u⃗_k'=∑_a=1^N[u^†_ k]_a [u_k']_a=∑_a=0^N-1exp-2πi(k-k')/N a=N δ_k,k' .
Hence we use the usual definition for the normal-mode momenta
p̃_k≡1/√(N)∑_a=0^N-1exp2πi k/N  ap_a .
Note the change in the sign in the exponential in comparison to eq. eq:Fourier. This definition then produces
the standard commutation relations: [x̃_k,p̃_k']=iδ_kk' and [x̃_k,x̃_k']=0=[p̃_k,p̃_k']. ] 
The normal-mode frequencies _k are defined in terms of the physical frequencies ω and Ω in the Hamiltonian qm88 as follows:
_k^2=ω^2+4Ω^2 sin^2πk/N ,
(see appendix <ref>). As desired, eq. qm288 reduces the problem to N decoupled harmonic oscillators, which enables us to easily write the ground-state wave function as
ψ_0(x̃_0,x̃_1, x̃_2,⋯)
=∏_k=0^N-1 (_k /π)^1/4 exp[-1/2 _k |x̃_k|^2] .
As before, this ground state will be the target state in our complexity computations. 

While eq. targetk 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 <ref>. In particular, following eq. eq:rot, we write the Fourier transformation eq:Fourier between the position and normal-mode bases as x̃ = R_ N  x, with
R_N≡1/√(N)[         1         1         1         …         1;         1         μ       μ^2         …     μ^N-1;         1       μ^2       μ^4         …  μ^2(N-1);         ⋮         ⋮         ⋮         ⋱         ⋮;         1     μ^N-1  μ^2(N-1)         … μ^(N-1)^2;           ]
 ,

where μ≡exp-2π i/N.[As discussed in footnote <ref> for the matrix R in eq. eq:rot, we distinguish  from ^T even though the numerical matrix in eq. logan is symmetric. Note that if we write out the transformation to show the indices, we have x̃_k = [R_ N]_ka  x_a. That is, the row index of  has values in the momenta k while the column index has values in the lattice position a.  In passing, we also observe that eq. logan reduces to eq. (<ref>) for the special case N=2, for which we have μ=exp-iπ=-1.] Since  is a unitary matrix, ^†=1, the inverse transformation is given by x =^† x̃. 

Now let us adopt the notation of section <ref> (and in particular, of eq. gauss) to write the target state targetk as
ψ_T(x̃_k)=∏_k=0^N-1_k/π^1/4 exp[-1/2 x̃^†Ã_T  x̃]    with    Ã_T=diag_0,…,_N-1 .

Using the rotation logan, we can write this target state in terms of the physical coordinates,[The relation _k=_N-k ensures that A_T is real. ]
ψ_T(x_a)=∏_k=0^N-1ω_k/π^1/4 exp[-1/2  x^T A_T   x ]    with    A_T=^†Ã_T  .



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 targetk. In analogy with eq. eq:refPhys, the reference state will be a factorized Gaussian state,
ψ_R(x_a)=ω_0/π^N/4 
exp[-1/2  x^T  A_R   x ]    with    A_R=ω_0 1 .

where the individual oscillators are completely unentangled.[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, ω_0≠_0!]  An important feature of our reference state is that it is invariant under translations on the lattice, 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:
ψ_R(x̃_k)=ω_0/π^N/4 
exp[-1/2  x̃^†Ã_R   x̃ ]    with    Ã_R=  A_R ^†=ω_0 1 .



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. eq:gates, 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. eq:pathPsi and path2, U(s)=𝒫 exp[∫_0^ss̃ Y^I(s̃) _I], where the index I runs over the N^2 values corresponding to pairs (ab), and the operators _I take the same form as in eq. operate. In discussing the target and reference states with the notation of eq. gauss, we also anticipated mapping these exponentials to the matrix formulation introduced in eqs. matrix0, matrix and pathA for Gaussian states. In fact, the generators have precisely the form given in eq. matrix, 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  group found in section <ref> to the group GL(N,ℝ) in the present problem.

Following the analysis in section <ref>, we use the analogous F_2 cost function, i.e.,[In the position basis, there is no need for the complex conjugations appearing in eqs. cost5 or metric5 since all of the relevant quantities are real. However, we are including them here in anticipation that later on, we will transform these formulae to the normal-mode space. These  transformations are accomplished with  in eq. logan, which is a complex unitary matrix.  Hence using, M^† rather than M^T allows us to use precisely the same expressions without change.  Of course, as defined in eq. eq:Fourier, the normal modes are generally complex, but as we commented above, they also satisfy the “reality condition” x̃_k^†=x̃_N-k—which ensures that we have not doubled the number of degrees of freedom.]

𝒟(U)=∫_0^1s√(δ_IJ Y^I(s)Y^J(s)^*) ,     where    Y^I(s)=_sU(s) U^-1(s)M^†_I .

Hence the optimal circuit will correspond to a geodesic in the GL(N,ℝ) geometry given by a right-invariant metric, analogous to eq. metric1. To simplify the discussion of the metric here (and in the next section), we introduce the following notation:
s^2=δ_IJ   Y^IY^J^*     with    Y^I=U U^-1  M^†_I .

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. Umatrix, a convenient parametrization of a general group element U∈GL(N,ℝ), 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. raffle and refk 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, ω_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 raffle. Via eq. rotateU, let us write the optimal circuit acting in the normal-mode basis; we have
U_0(s) = ^† Ũ_0(s)      where    Ã_T=Ũ_0(s=1)
 Ã_R Ũ^†_0(s=1) ,
and thus the straight-line circuit Ũ_0(s) becomes
Ũ_0(s)=exp[M̃_0 s]
    with  
M̃_0=diag(1/2log_0/ω_0, 
1/2log_1/ω_0,
⋯, 1/2log_N-1/ω_0 ) .

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. swer4 and metric5, 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 subspace of U∈GL(N,ℝ) with the form U= Ũ^†, where[In general, the coordinates ỹ_k are complex but satisfy the normal-mode “reality condition” ỹ_k^†=ỹ_N-k.]
Ũ_n-m=exp[M̃_n-m]
    with    M̃_n-m=diag(ỹ_0, ỹ_1, ⋯,  ỹ_N-1 ) .

Substituting this expression into eq. metric5, one finds the following flat Cartesian metric induced on this subspace:
s^2_n-m=|ỹ_0|^2+|ỹ_1|^2+⋯+|ỹ_N-1|^2 .

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 normal-mode subspace metric6, the optimal circuit is described by the simple circuit in eq. nmU, which we write as
U_0(s) = ^† Ũ_0(s)  = exp[^†  M̃_0   s] ,

where Ũ_0(s) and M̃_0 are defined in eq. raffle2. 

There are actually some subtleties in the preceding argument which make the conclusion somewhat premature. The first is that eq. metric5 which defines the metric is written in the position basis, whereas eq. metric6 was implicitly calculated for an expression raffle4 written in the normal-mode basis. That is, in eq. raffle4, we worked with M̃_n-m=Ỹ^ĨM̃_Ĩ with a particular choice of Ỹ^Ĩ.[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 x̃_k.] However, we show in appendix <ref> that this was nonetheless a valid approach since the metric takes precisely the same form when written in terms of the normal-mode space. This requires extending the discussion around eq. bark2 describing the change of bases for the case of two coupled 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. straight9, 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. raffle4 that all of the other coordinates in the  geometry could be set to zero. Recall that in the  metric metric1, the metric on the normal-mode subspace was completely independent of the other coordinates, we had s^2_n-m=2 y^2 +2ρ^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 Ũ using the Iwasawa (or KAN) decomposition of ; see for example <cit.>. This states that any Ũ∈ 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,[We denote this diagonal matrix with the traditional A, but it should not be confused with the quadratic forms specifying the Gaussian states, eq. gauss. Similarly, K here should not be confused with the gates producing a momentum shift in eq. eq:gates, nor should N be confused with the total number of oscillators. We trust that these distinctions will be clear from context.] 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. raffle4.

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 (eq. raffle4) 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 <ref> allows us the eliminate the K_0 rotation as well: following eq. tran1, we write the differentials Ỹ^Ĩ=(Ũ Ũ M̃^†_Ĩ) as
Ỹ^Ĩ=
Ũ_n-m  Ũ_n-m^-1 [K_0^T M̃_Ĩ K_0]^† .
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 lot1, it is straightforward to show – following a series of steps analogous to those given in eqs. lot2 or lot2a – that
K_0^T M̃_Ĩ K_0 = [K_0]_ĨJ̃ M̃_J̃    where    K_0=K_0⊗K_0 .

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
s^2_n-m=δ_ĨJ̃ Ỹ^Ĩ_n-m
( Ỹ^J̃_n-m)^*    where    Ỹ^Ĩ_n-m=Ũ_n-m  Ũ_n-m^-1  M̃_Ĩ^† ,

which again yields the simple answer given in eq. metric6.

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/2log(ω̃_k/ω_0) as in eq. raffle2. Eq. tran6 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 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, we extend our initial ansatz to allow small excursions in the K and N directions,
Ũ=K_0 exp[ M̃^rot_Ĩθ_Ĩ] Ũ_n-mexp[M^shift_Ĩη_Ĩ] N_0 ,

where θ_Ĩ, η_Ĩ≪1. Here, the (small) change in K only involves the
(antisymmetric) rotation generators
[M^rot_kℓ]_pq=δ_kpδ_ℓq- δ_ℓpδ_kq    with    k<ℓ ,
while the (small) change in N only involves the shift generators 
[M^shift_kℓ]_pq=δ_kpδ_ℓq    with    k<ℓ .
The rotation generators are, of course, a linear combination of the original generators given in eq. lot1, and hence are not orthogonal to the shift generators in the sense that
M^rot_kℓ [M^shift_pq]^†=δ_kp δ_ℓq .

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 extend0, we now evaluate the differentials Ỹ^Ĩ=( Ũ Ũ M̃^†_Ĩ) on the normal-mode subspace, at θ_Ĩ=0=η_Ĩ,
Ỹ^Ĩ=[K_0]_ĨJ̃ [ 
Ỹ^J̃_n-m
+M̃^rot_K̃ dθ_K̃ M̃_J̃^†+ Ũ_n-m M^shift_Ĩ  Ũ_n-m^-1 dη_Ĩ M̃_J̃^†]
 ,
where K_0 is the orthogonal matrix given in eq. tran5 and Ỹ^Ĩ_n-m are the differentials along the normal-mode directions identified in eq. tran6. 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 Ỹ^J̃_n-m are along the diagonal directions, 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. rotgen, it is clear that the second term only contributes in the off-diagonal directions, J̃=(kℓ) with kℓ. Similarly, one can show that the same is true of the third term via eqs. raffle4 and shiftgen,
Ũ_n-m M^shift_kℓ Ũ_n-m^-1=e^ỹ_k-ỹ_ℓ
 M^shift_kℓ    with    k<ℓ .

The key point then is that Ỹ^J̃_n-m are orthogonal to the other two differentials in eq. tran7. In fact, it is straightforward to show that the full metric on the normal-mode subspace (θ_Ĩ=0=η_Ĩ) becomes
s^2_n-m=|ỹ_0|^2+|ỹ_1|^2+⋯+|ỹ_N-1|^2
+∑_k<ℓ[(dθ_kℓ)^2 + |dθ_kℓ+e^ỹ_k-ỹ_ℓ dη_kℓ|^2
] .

Hence there are no off-diagonal terms in the metric, which would drive the geodesic away from the normal-mode subspace.[In general, we are asking that there are no source terms in the linearized equations for the θ_Ĩ and η_Ĩ. In turn, this means that we are asking that there are no linear terms in the cost function cost5 when expanding about the straight-line trajectories raffle2. Here we have explicitly shown that no such linear terms arise as off-diagonal terms in the metric, involving the differentials of θ_Ĩ and η_Ĩ. In principle, we should also verify that the metric components g_ỹ_kỹ_k are not varied at linear order in the perturbations. However, our previous analysis shows that the metric on the normal mode subspace is completely independent of the coordinates parametrizing the K and N transformations, which ensures that this class of potential terms linear in θ_Ĩ and η_Ĩ vanishes. Therefore the above discussion is sufficient to ensure that the geodesic equations in the full GL(N,ℝ) geometry have no source terms which would push the straight-line trajectories 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. raffle2. 

Thus, for the cost function (<ref>), the complexity for our lattice of oscillators is obtained by simply summing up the circuit elements in the normal-mode basis.  Using eqs. (<ref>) and (<ref>), we find
=1/2√(∑_k=0^N-1(
logω̃_k/ω_0)^2) ,
where the normal-mode frequencies are given in eq. eq:eigenfreq. Recall that in our lattice regularization eq:Hlattice, we had ω=m and Ω=1/δ, and so we can express the complexity complexityN in terms of the field theory parameters via
_k^2=m^2+4/δ^2 sin^2πk/N .
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 refX.

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
=1/2√(∑_{k_i}=0^N-1(
logω̃_k⃗/ω_0)^2) ,
where k_i are the components of the momentum vector k⃗=(k_1,k_2,⋯,k_d-1), and the normal-mode frequencies are given by
_k⃗^2=m^2+4/δ^2 ∑_i=1^d-1sin^2πk_i/N .
The linear size of each spatial direction here is L=Nδ, and so the total (spatial) volume of the system is V=L^d-1=N^d-1δ^d-1. Hence the total number of oscillators can be expressed as
N^d-1=V/δ^d-1 ,

which will prove useful below.



 §.§ Comparison with holography


Eq. complexityNd 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, 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. complexityNd, 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, 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. complexityNd to obtain[See appendix <ref> for more accurate estimates of the large N behaviour of the complexity complexityNd.]
≈N^d-1/2/2 log1/δ ω_0∼V/δ^d-1^1/2 ,

where we have used eq. vold 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 <cit.>. Hence we can compare our QFT result croot with the analogous results for holographic complexity; denoting the latter collectively as _holo, these were found to take the form 
_holo∼V/δ^d-1 .
Thus we see that the leading terms in the QFT and holographic complexities differ by the power of 1/2 appearing in eq. croot. However, the origin of this square root is clear: it is simply the overall square root appearing in eq. complexityNd, which in turn arises from our use of the F_2 cost function in eq. cost5. Now, there is nothing wrong with the result in eq. croot 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.[An alternative approach <cit.> would be to simply assign each gate the cost N^d-1/2. However, it seems this may be problematic if, we wish to compare complexities for different UV cut-offs. ] For example, the cost function defined by the F_1 measure in eq. eq:Fmetrics, 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. Cholo is
𝒟_κ=∫_0^1s∑
 |Y^Ĩ(s)|^κ ,
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 <ref>.

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, that replacing s→ŝ(s) leaves the cost unchanged. However, we may proceed with the physics intuition that we can think of 𝒟_κ 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. raffle2 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, 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)], where
M̃_1=diag(f_0(s)/2log_0/ω_0 , 
f_1(s)/2log_1/ω_0,
⋯, f_N-1(s)/2log_N-1/ω_0 ) ,

and each of the f_k(s) is an arbitrary function satisfying
f_k(s=0)=0    and    f_k(s=1)=1 .

With this ansatz, the cost function Dalpha evaluates to
𝒟_κ(Ũ_1)=1/2^κ∑_k |log_k/ω_0|^κ×∫_0^1s |∂_s f_k(s)|^κ   .

However, with reasoning along the lines of that in eqs. fun1 and fun2, one can argue that ∫_0^1 s |∂_s f_k(s)|^κ≥1 (for κ≥1), and that the inequality is saturated if and only if ∂_s f_k(s)=1, f_k(s)=s.[An exception to this result arises for κ=1. In this case, the bound is saturated by any functions f_k(s) satisfying ∂_s f_k(s)≥0 everywhere. We return to this point in the discussion section <ref>.] That is, the general κ cost function Dalpha 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
≈V/δ^d-1 
| log1/ω_0 δ|^κ .

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 refX, 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. calpha reduces to
≈σ^κ V/δ^d-1 .

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. Dalpha includes summing over the infrared modes, modes with ω̃_k⃗∼ m. These infrared contributions will take the form
_IR≈-  log^κ(mδ) .



An alternative choice would be to associate ω_0 with some infrared scale, ω_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.[Furthermore, if we choose ω_0∼ V^-1/(d-1), the complexity becomes superextensive.] This appears problematic if we wish to compare the complexities of states in different theories, 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 calpha 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 ccaction takes the form <cit.>
_A∼V/δ^d-1 log(L_AdS/α δ) ,

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 _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 α, we set α = ω_0 L_AdS where ω_0 is some arbitrary frequency. In this case,  eq. leader reduces to
_A∼V/δ^d-1 log(1/ω_0 δ) . 

While this choice eliminates the AdS scale, the holographic complexity still depends on the choice of  ω_0, just as in our QFT result calpha. 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 <cit.>. Whereas in <cit.> 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. calpha to match the holographic result leader2 more closely, we should choose κ=1 in eq. Dalpha, 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 pre-factor, which we have not specified here. However, κ also appears in the IR contribution in eq. calpha3. Again the analogous contributions in holography would be linear in log(δ) because of the form of the corresponding boundary terms in the gravitational action <cit.>. 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 <ref>,  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 <cit.>. 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 <cit.>. However, we expect that analogous results will appear for the CA conjecture. Such corrections to holographic complexity were considered in <cit.> 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. cost2, 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 – 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
G_IJ=diag(1,^2,^2,1) ,

with >1.  As a result, our original metric metric1 is replaced by the following more complicated metric:
s^2   =    G_IJ  U(s) U^-1(s)M^T_I U(s) U^-1(s)M^T_J

   =    2 y^2+2[^2-(^2-1)sin^2(θ+τ)]ρ^2+(^2-1)sin2(θ+τ)sinh(2ρ) ρ (τ-θ)


      +2[^2cosh(2ρ)-(^2-1)cos^2(θ+τ)
sinh^2ρ]cosh^2ρ τ^2


      +2[^2cosh(2ρ)-(^2-1)cos^2(θ+τ)
cosh^2ρ]sinh^2ρ θ^2

      -[2^2-(^2-1)cos^2(θ+τ)]sinh^2(2ρ) τ θ .
Of course, this geometry reduces to eq. (<ref>) upon setting =1. The metric has a slightly simpler expression in terms of the pseudo-lightcone coordinates (<ref>), where θ=x+z, τ=x-z:
s^2   =   2y^2+2[^2-^2-1sin^22x]ρ^2-2 (^2-1)sin(4 x)
 sinh(2 ρ) ρz
 

       + 2 ^2x^2+2[^2cosh(4ρ)-(^2-1)cos^2(2x)sinh^22ρ]z^2-4 ^2cosh(2 ρ) xz ,
which reduces to eq. metric2 when =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 penalty2 for most of the following. 

As in the unpenalized case, this metric enjoys the four Killing vectors (k̂_I)^i given in eq. (<ref>).[Again, these arise from the right-invariance of the expression in the first line of eq. penalty1. The fifth “accidental” Killing vector ∂_x in eq. Killed5 no longer gives rise to a symmetry for the penalized metric, as is clear from eq. penalty2.] For the metric (<ref>), the associated conserved quantities ĉ_I=(k̂_I)^ig_ijẋ^i are

ĉ_1≡    2ẏ ,


ĉ_2≡     -ẋ[2 ^2 sinh(2 ρ)  cos(2 z)]

   +ż [cos(2 z) sinh(4 ρ)  (2^2-(^2-1) cos^2 (2 x)) -(^2-1)sin(2 z) sinh(2 ρ) sin(4 x) ]

   + ρ̇[2sin(2 z) (^2-(^2-1) sin^2 (2 x))-(^2-1) cos(2 z)cosh(2 ρ) sin(4 x) ] ,


ĉ_3≡    ẋ[2 ^2 sinh(2 ρ)  sin(2 z)]

   - ż [ sin(2 z) sinh(4 ρ)  (2^2-(^2-1) cos^2 (2 x))+(^2-1) cos(2 z) sinh(2 ρ) sin(4 x) ]

   + ρ̇[2cos(2 z) (^2-(^2-1) sin^2 (2 x))+(^2-1) sin(2 z) cosh(2 ρ) sin(4 x) ] ,


ĉ_4≡     -2 ^2 cosh(2 ρ)  ẋ
+2 (^2 cosh(4 ρ)-(^2-1) cos^2 (2 x)sinh^2(2 ρ) ) ż

    - (^2-1)  sinh(2 ρ) sin(4 x) ρ̇ .


One can check that these quantities indeed reduce to those given in eq. (<ref>) when =1. Solving the first equation for y is trivial, and we simply recover eq. (<ref>), y=y_1 s. The next three equations may be solved for ρ̇, ẋ, and ż:
ρ̇=   
1/4 ^2{(^2-1) cosh(2 ρ) sin(4 x) (ĉ_2 cos(2 z)-ĉ_3 sin(2 z))

   +2(^2-^2-1cos^2 (2 x)) (ĉ_2 sin(2 z)+ĉ_3 cos(2 z))-(^2-1) ĉ_4 sinh(2 ρ) sin(4 x)}


ẋ=   1/4 ^2{1/sinh(2 ρ)[2cosh(4ρ)+^2-1cos^2(2x)cosh^22ρĉ_2cos(2z)-ĉ_3sin(2 z)

   +^2-1sin(4x)cosh2ρĉ_2sin(2z)+ĉ_3cos(2z)]
-2(1+(^2-1) cos^2 (2 x))cosh2ρĉ_4}


ż=   1/4 ^2{1/sinh(2ρ)[2(1+(^2-1) cos^2 (2 x))cosh2ρ(ĉ_2 cos(2 z)-ĉ_3 sin(2 z))

   +(^2-1) sin(4 x) (ĉ_2 sin(2 z)+ĉ_3 cos(2 z))]
-2(1+(^2-1) cos^2 (2 x))ĉ_4}



Now recall that in the unpenalized case, the expression for θ̇ in eq. trun diverged at the origin ρ(s=0)=0 unless the conserved quantities were properly tuned. This divergence is simply the usual angular momentum barrier at the origin, and the tuning amounts to setting the angular momentum to zero. The same issue arises here, as reflected in the fact that both ẋ and ż have the same pole structure as ρ→0, in this limit, θ̇=ẋ+ż diverges but τ̇=ẋ-ż does not. Taking the limit ρ→0, and setting x_0=z_0 (since τ(s=0)=0) in eq. revel1, one finds that this divergence is avoided by choosing
ĉ_2=ĉ_3/^2tan2z_0 .

Substituting eq. lov2 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 penalty2, the normalization of the tangent vector k^2=g_ij ẋ^iẋ^j becomes
k^2=    2 y_1^2+2[^2-^2-1sin^22x]ρ̇^2-2 (^2-1)sin(4 x)
 sinh(2 ρ) ρ̇ ż

   + 2 ^2ẋ^2 +2[^2cosh(4ρ)-(^2-1)cos^2(2x)sinh^22ρ]ż^2-4 ^2cosh(2 ρ) ẋ ż .

In principle, one can substitute the expressions in eq. revel1 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 solver5, and as we now verify, this trajectory remains a geodesic in the penalized geometry penalty2. Recall that the first constraint in eq. revel0 yielded the desired behaviour for y, y(s)=y_1 s, as in eq. yrun. The straight-line 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. lov2 then yields ĉ_2=0. Substituting these values of x and z into the last two expressions in eq. revel1, we obtain
ẋ=-ĉ_4/2 cosh(2ρ)     and    ż=-ĉ_4/2 .

Hence consistency with the condition ẋ=ż=0 demands that we set
ĉ_4=0. Finally, the ρ̇ equation yields
ρ̇=-ĉ_3/2^2    ĉ_3=-2^2 ρ_1 ,

and we arrive at the desired solution: ρ(s)=ρ_1 s. 

Having shown that this simple trajectory remains a geodesic in the penalized geometry penalty2, we substitute this geodesic into eq. (<ref>) to obtain
k^2=2y_1^2+^2ρ_1^2≡k_0^2 ,
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. (<ref>) is the natural generalization of eq. (<ref>) to the case with >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 penalty2 more closely, we see that 
g_ρρ=2[^2-^2-1sin^22x]

depends on the x coordinate. This contrasts with the unpenalized metric metric1 (or the lattice metric metric6a) where we found that the geometry of the normal-mode subspace (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. rats6, it is clear that there should be short-cuts for the motion along ρ if we move away from the normal-mode subspace, away from x=π/2.[As an amusing observation, let us add that it is also clear that there are no such short-cuts if <1. That is, if we weight the scaling gates more heavily than the entangling gates, then the straight-line circuit solver5 will remain the optimal circuit.] For example, we might consider the following simple path consisting of two segments:
a)  0≤s̅≤1 :      y=0 , ρ=ρ_1 s̅ , x=π/4 , z=π/4 ;


b)  1≤s̅≤2 :      y=y_1(s̅-1) , ρ=ρ_1 , x=π/4 s̅, z=π/4 ,
where s̅ 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
k_s=√(2)ρ_1+√(2π/4^2^2+2y_1^2) ,
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 . 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, ≫ρ_1,y_1 and ρ_1,y_1≫1. Then we may approximate the two lengths with
k_s   ≃   π/2√(2) +√(2)ρ_1+ 2√(2)/π y_1^2/^2 +⋯ ,


k_0   ≃   √(2)ρ_1+ y_1^2/√(2)ρ_1 +⋯ ,      k_0/ k_s≃4/π ρ_1≫1 .
Thus in the penalized geometry penalty2, 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. eq:kPenFull, obtaining the general solution seems out of reach. However we can make progress with a simplifying assumption: if we examine the penalized metric penalty2, we see that as the radius ρ increases, the most rapidly growing component of the metric is g_zz∼^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 origin. Therefore we simplify our problem by considering trajectories confined to a constant-z submanifold, for which the relevant metric is given by
s^2=2y^2+2[^2-^2-1sin^22x]ρ^2+2 ^2x^2 .
We will return to justify our assumption of no (or little) motion in the z direction below. Obviously, eq. (<ref>) 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 <ref> 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 , but we must add that we have not provided an irrefutable proof of this result. Furthermore, using the results of appendix <ref>, we also explicitly show below that the segmented path segments provides a good approximation of this optimal geodesic in the regime where the penalty factor is large.

Our analysis in appendix <ref> suggests we use the following quantity to more conveniently compare the lengths of the various paths:
k̅^2≡k^2/2-y_1^2 .

For the straight-line geodesic, we have simply
k̅_0= ρ_1 ,

while for the segmented path, we have
k̅_s   =   π/4[1+8ρ_1/π2ρ_1/π+√(1+4y_1/π^2)]^1/2


   ≈   π/4+ρ_1+8/π^2 y_1^2/^2 ρ_1+… ,

where as above, the expansion in the second line assumes ≫ρ_1,y_1. Now for the optimal geodesic, we have from eq. eq:kbarCoords
k̅   ≈   /2tan^-1√(^2-1)+ρ_1


   ≈   π/4+ρ_1-1/2-1/12^2+1/^4 .

Comparing these results, we see that in this large  regime, k̅ for the optimal geodesic is much smaller than k̅_0 for the straight-line geodesic, and extremely close to k̅ for the segmented path. Furthermore, we note that both k̅_0 and k̅ are completely independent of y_1, while it only appears in k̅_s at order  y_1^2/^2.

In figure <ref>, we plot k̅, k̅_0, and k̅_s as functions of  for fixed values of the variable  (see eq. eq:kEps), as well as for various values of y_1 in the case of k̅_s. Recall the definition newer 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 .
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. kbars and kbaropt. In fact, if we take the difference of these two equations in the large  limit, we find
k̅_s-k̅≃1/2 .
-2ex


Given that in this limit, both k̅ and k̅_s are diverging, it is impressive to find the simple (1) difference shown above. Figure <ref> 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 segments is not itself a geodesic, it must have a longer length. However, the impressive agreement in eq. eq:overhead 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 <ref>, 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 . We might note that the growth in y is uniform throughout the entire span 0≤ s≤ 1. However, for large , 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 segments for large . 



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 for ẋ and ρ̇ cross in figure <ref>, ẋ(s)/ẋ|_max =ρ̇(s)/ρ̇|_max at s=s_trans. In the limit ≪1, this point can be well-approximated by[This expression is obtained by equating the normalized quantities ẋ/ẋ|_max and ρ̇/ρ̇|_max (from eqs. (<ref>) and (<ref>), respectively) to solve for the critical point x_crit at which the curves in the right plot in figure <ref> cross. Upon substituting this into (<ref>) for s(x), and using (<ref>) for k̅ and (<ref>) for c̅_2, one obtains an expression that depends only on  and , s,, in which we then take ≪1.]
s_trans≈1-Π1-^2, h(), 1-/Π1-^2, 1- ,
where
h()≡^-1[^2-1+√(5^2-4)/2^3--√(5^2-4)]^1/2 .
We plot s_trans() in figure <ref>. In conjunction with eq. (<ref>), 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 →∞. We can verify this by first approximating h() for large  as
h()≈^-1[2/√(5)-1]^1/2 ,
and then expanding eq. (<ref>) in the limit →∞, →0:
s_trans()≈1/[1/πlog1/+^0]+1/^2 ,
where ^0≈1.11502. Comparing this expression with eq. (<ref>) for large , we see that the leading-order term in s_trans() above may be written as
s()≈4 ρ_1/π  .
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  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 eq:metricSimple 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. pen99, the complexity of the ground state becomes
=√(2)[π/4 +ρ_1-1/2+2 y_1^2/π  +1/^2] .

in the regime ≫ρ_1,y_1.



 §.§ New optimal circuit


Given the optimal geodesic for the penalized geometry, we would now like to examine the properties of the corresponding circuit. For simplicity, we rely on the fact that for ≫1 the optimal geodesic is well approximated by the segmented path described in eq. segments and explicitly build the circuit for the latter path. First however, let us rewrite eq. segments in terms of θ=x+z and τ=x-z:
a)  0≤s̅≤1 :      y=0 , ρ=ρ_1 s̅ , θ=π/2 , τ=0 ;


b)  1≤s̅≤2 :      y=y_1(s̅-1) , ρ=ρ_1 , θ=π/4(s̅+1), τ=π/4(s̅-1) .
Recall that s̅ 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 <ref>. This is essentially a comparison of the optimal geodesic in the original geometry metric1 to that in the new penalized geometry penalty1. 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  in eq. (<ref>), we write the segmented circuit U_s(s̅) as
U_s(s̅)=
U_a(s̅)=[ e^-ρ_1 s̅         0;         0   e^ρ_1s̅ ]     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 ,

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
U_b(s̅)=e^y_1(s̅-1)R̅(s̅) U_a(s̅=1) , 
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) ] .
The interpretation of eq. (<ref>) 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 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, s̅=2, the rotation matrix reduces to
R̅(s̅=2)=1/√(2)[  1 -1;  1  1 ] ,

which is closely related but distinct from the rotation matrix R defined previously in eq. (<ref>). 

At its end-point, this new circuit becomes
U_s(s̅=2)=e^y_1/√(2)[ e^-ρ_1 -e^ρ_1; e^-ρ_1  e^ρ_1 ] ,

which we might compare to the end-point of the straight-line circuit solver5,
U_0(s=1)=e^y_1[  coshρ_1 -sinhρ_1; -sinhρ_1  coshρ_1 ] .

These are clearly different, in accordance with our comment about the geodesics ending at different points in  figure <ref>. However, it is straightforward to show that the transformations implemented in the segmented and straight-line circuits are related by the rotation in eq. dead, U_s(s̅=2)=U_0(s=1) R̅(s̅=2) .

Both of these transformations act on the reference state to produce the target state (see eq. eq:stateMatrix) as A_T =U A_R U^T. Since the reference state is proportional to the identity, the additional rotation in eq. eq:U2phys leaves this state invariant, A_R=R̅(s̅=2)  A_R R̅^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 eq:segU in the normal-mode space using eq. rotateU, Ũ(s)=R U(s) R^T, which yields
Ũ_s(s̅)=
Ũ_a(s̅)=[  coshρ_1 s̅ -sinhρ_1 s̅; -sinhρ_1 s̅  coshρ_1 s̅ ]     for 0≤s̅≤1 ,


Ũ_b(s̅)=e^y_1(s̅-1)[  cosπ/4(s̅-1)  sinπ/4(s̅-1); -sinπ/4(s̅-1)  cosπ/4(s̅-1) ] Ũ_a(s̅=1)     for 1≤s̅≤2 ,

where we have expressed Ũ_b(s̅) in a form analogous to eq. eq:seg2. 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, gates which entangle (or disentangle) the normal-mode coordinates. Since we know that the reference state Atilde and the target state Atilde2 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 (for 0≤s̅≤1), we find
Ã(s̅)=
ω_0 Ũ_a(s̅) Ũ^T_a(s̅)=ω_0 [  cosh2ρ_1 s̅ -sinh2ρ_1 s̅; -sinh2ρ_1 s̅  cosh2ρ_1 s̅ ]  .

Here we see the entanglement (the off-diagonal terms) begins at zero and steadily grows to a maximum at  s̅=1 at the end of the first segment. Subsequently, along the second segment (for 1≤s̅≤2), we find
Ã(s̅)   =   ω_0 Ũ_b(s̅) Ũ^T_b(s̅)


   =   ω_0 e^2y_1(s̅-1)[ cosh2ρ_1 -sinh2ρ_1  sinπ/2(s̅-1)          -sinh2ρ_1  cosπ/2(s̅-1);          -sinh2ρ_1  cosπ/2(s̅-1) cosh2ρ_1 +sinh2ρ_1  sinπ/2(s̅-1) ]  .
Here, the entanglement shrinks steadily back to zero as s̅ runs over this second interval. Recall that y_1 and ρ_1 are given in terms of the normal-mode frequencies in eq. fin2.

To reiterate our key observation, with eq. pen99 penalty factors are introduced to increase the cost of the entangling gates in the position space. As a result (for large ), the optimal geodesic is deformed to be close to the segmented path described in eqs. segments and segmentsX. 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 pen99 to the normal-mode basis using the orthogonal matrix  defined in eq. bark2. The new metric then becomes
G_ĨJ̃ = _Ĩ I  G_I J ^T_JJ̃=1/2 
[ 1+^2    0    0 1-^2;    0 1+^2 1-^2    0;    0 1-^2 1+^2    0; 1-^2    0    0 1+^2 ]
 .

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 <cit.>, 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 eq:gates, in particular, scaling and entangling gates. We also chose our reference state to be a factorized Gaussian state refX, 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 targetX 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 metric1 on . 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 .



 §.§ 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 <ref> 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 Cholo in both the CV and CA proposals are extensive, they are proportional to the volume of the time slice on which the boundary state is evaluated, as shown in <cit.>. 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. Dalpha, all of which yield an extensive complexity for the scalar field theory.[Again, we remind the reader that one can continue to work with the F_2 measure if the cost of the individual gates is set proportional to V/δ^d-1^1/2—see footnote <ref>.] 

With the new cost functions Dalpha, 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 refX. However, this precisely matched an ambiguity in the holographic complexity <cit.> found for the CA construction ccaction. In the latter case, the logarithmic factor came from joint terms <cit.> 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 volver. However, as explained in section <ref>, the additional logarithmic factor in the leading term calpha 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. calpha3 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 <cit.>.

As further noted in section <ref>, our scalar field complexity emulates the CA proposal ccaction most closely if we choose the F_1 cost function, κ=1 in eq. Dalpha. 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 state-surface conjecture in <cit.>. 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. scotch. Clearly, even though it is easy to set out interesting questions for complexity (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.[Of course, the tolerance does not play a role in the geometric framework adopted here because the gates are no longer discrete—see discussion around eq. success.] The ambiguity in our reference state, 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, <cit.>. 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. croot and calpha), 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, 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 <cit.> 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. eq:Cd1, where we are examining the case d=2 and κ=1. In this instance, we may isolate this constant in[Of course, we are inspired to formulate this quantity by the constructions using entanglement entropy to examine RG flows of three-dimensional QFTs <cit.>.]
(L ∂_L-1)=-a_0 .

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 <ref>.

1.5ex

Let us also remark here that the reference state refX 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 <cit.>—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⟩:
⟨ψ_R| H |ψ_R⟩-⟨ψ_T| H |ψ_T⟩   = 1/4 N^d-1 ω_0+1/4ω_0 ∑ω_k⃗^2-1/2 ∑ω_k⃗

   ≈V/δ^d  [ω_0δ/4+1/ω_0δ-1] .
We therefore see that that generically, if ω_0δ≪ 1 or ωδ≫1 (which were advocated to be the natural choices in section <ref>), the “renormalized” energy density of |ψ_R⟩ diverges as δ→0.[Note that with some fine-tuning of ω_0, we could arrange the difference of energy densities to be finite.] Hence this would not be a state that would be considered to be part of the standard 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. 

1.5ex

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.[We thank Edward Witten for making this observation and raising the following question.] 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 <ref>, 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 pen99. 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 (), as shown in eq. loft. 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 <ref>.

However, in our simple experiment in section <ref>, 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 <ref> 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. (<ref>) 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
Q_13=Q_12^-1Q_23^-1Q_12Q_23^1/ .
In other words, the next-to-nearest neighbour entangling gate Q_13 is equivalent to 4/ nearest-neighbour entangling gates. A simple generalization of the above result is
Q_14=Q_13^-1 Q_34^-1 Q_13 Q_34^1/ ,
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
c(n)=2/1+c(n-1)
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 <cit.>, 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. tsoc suggests increasing the penalty factors as a power law to match the growth of the nonlocality, the directions corresponding to I=(a,a+1+n) and (a+1+n,a) would be assigned a penalty factor ^2n. Note that this does not penalize the nearest-neighbour gates at all, in contrast to our experiment in section <ref>. Of course, for a periodic chain of oscillators, the maximum penalty factor would be ^N-2 and ^N-3 for even and odd N, respectively.

We can gain further insight by translating eq. eq:Q13 into a macroscopic circuit described by a path-ordered exponential pathA. In particular, consider the following path:
Y^23(s)   =    α [ 1-Θ(s-1) -Θ(s-2) +Θ(s-3)] ,


Y^12(s)   =    α [ Θ(s-1)-Θ(s-2) +Θ(s-3) -Θ(s-4) ] ,

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 
U_1(s)=𝒫 exp∫_0^xs̃ Y^I(s̃) M_I

then yields U_1(s=4)=exp[α^2 M_13], following the same calculation that yields eq. eq:Q13. Hence we could accomplish the same transformation with
U_2(s)=𝒫 exp∫_0^xs̃ Y^13(s̃) M_13
    where    Y^13(s)=α^2/4   for  
0≤s≤4 .

Now let us compare the costs of these two circuits using the F_q measure eq:Fmetrics, where the nearest neighbour gates are assigned cost 1 while the next-to-nearest neighbour gates are assigned cost . The cost functions are then easily evaluated to be
D(U_1)=∫_0^4s √(δ_IJY^I(s)Y^J(s))=4α ,    
D(U_2) = 4|Y^13|= α^2 .
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 <cit.>. This proposal suggests that the MERA (Multiscale Entanglement Renormalization Ansatz)  tensor network <cit.> provides a discrete representation of a time slice of (three-dimensional) AdS space. As illustrated in figure <ref>, 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 <cit.>. Furthermore, in these discussions, it was implicitly considered the optimal circuit for the preparation of the CFT ground state.



There has been some progress towards developing a continuum version of MERA, however, these constructions are limited to describing very simple QFTs <cit.>—see also <cit.>. 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 <ref>. The key difference between the two circuits is that the cMERA circuit is organized to systematically introduce entanglement scale-by-scale, to order the amplification of the normal modes according to their wavelength <cit.>. However, by almost all of the measures considered in section <ref>, including the F_2 cost function and the κ cost functions in eq. Dalpha, 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,[Let us add here that the F_1 measure also exhibits some unusual properties under a change of basis, as discussed in appendix <ref>.] 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, 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=-∑λ_ilogλ_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 <cit.>. The latter applies familiar tools (path integrals) in a novel setting (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 <cit.>—see also <cit.>.



§ ACKNOWLEDGMENTS

It is a pleasure to thank Micha Berkooz, Eugenio Bianchi, Shira Chapman, Adrián Franco-Rubio, Lucas Hackl, Markus Hauru, Michal Heller, Qi Hu, Steve Jordan, Hugo Marrochio, John Preskill, Djordje Radicevic, Grant Salton, Joan Simon, Guillaume Verdon-Akzam, Guifré Vidal, Edward Witten, and Beni Yoshida for helpful conversations.  Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development and by the Province of Ontario through the Ministry of Research & Innovation. RCM is also supported by an NSERC Discovery grant, as well as research funding from the Canadian Institute for Advanced Research and from the Simons Foundation through the “It from Qubit" Collaboration. RJ is supported by the Δ-ITP consortium and the Foundation for Fundamental Research on Matter (FOM), both of which are parts of the Netherlands Organization for Scientific Research (NWO) funded by the Dutch Ministry of Education, Culture, and Science (OCW). RJ also gratefully acknowledges the support of the Perimeter Institute Visiting Graduate Fellows program. Finally, RCM would also like to thank the organizers of the It-From-Qubit “Complexity and Black Holes” workshop at Stanford University, the “Tensor Networks for Quantum Field Theories II” workshop at Perimeter Institute and the “Strings 2017” conference in Tel Aviv for the opportunity to present this work.

1.5ex
While this paper was in preparation, we were informed of <cit.>, which seems to have significant overlap with the present work.




§ EXAMPLE CIRCUITS

In this section, we analyze the circuit depth of a few discrete circuits, as  introduced in section <ref>. In particular, let us consider the example given in eq. (<ref>),
ψ_T=U_1 ψ_R= Q_22^α_3 Q_21^α_2 Q_11^α_1 ψ_R ,

where the target wave function is given in eq. eq:targetPhys and the reference wave function, in eq. eq:refPhys. The question then is to determine the exponents α_i in this equation, 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 gate applied is Q_11, which rescales the coefficient of x_1^2 in the exponent of the Gaussian wave function. Next, by applying Q_21, the two oscillators become entangled. Finally, the application of Q_22 rescales x_2 to ensure that x_2^2 appears with the correct coefficient. Hence let us begin the quantitative analysis by considering:
Q_11^α_1 ψ_0(x_1,x_2)
=e^α_1/2ψe^α_1x_1,x_2=√(ω_0/π)e^α_1/2exp[-ω_1/2x_1^2-ω_0/2x_2^2] ,
where
ω_1≡e^2α_1ω_0α_1=1/2logω_1/ω_0 .

Next, applying the Q_21 gates yields
Q_21^α_2 Q_11^α_1 ψ_0(x_1,x_2)
   =√(ω_0/π)e^α_1/2e^ix_2p_1exp[-ω_1/2x_1^2-ω_0/2x_2^2]

   =√(ω_0/π)e^α_1/2exp[-ω_1/2x_1+α_2 x_2^2-ω_0/2x_2^2]

   =√(ω_0/π)e^α_1/2exp[-ω_1/2x_1^2-1/2ω_0+^2α_2^2ω_1x_2^2-α_2ω_1 x_1x_2] .
Note that the x_1x_2 cross-term will be rescaled in the next step, so we cannot fix any of the coefficients quite yet. Finally, we rescale x_2 with the Q_22 gates:
Q_22^α_3 Q_21^α_2 Q_11^α_1 ψ_0(x_1,x_2)
=√(ω_0/π)e^(α_1+α_3)/2exp[-ω_1/2x_1^2-ω_2/2x_2^2-βx_1x_2] ,
where α_2 and α_3 are determined by matching the second and third coefficients in the exponent, ω_2=ω_0+^2α_2^2ω_1e^2α_3 ,
   
 β≡α_2ω_1  e^α_3 .
Solving the above constraints then yields
α_2=1/√(ω_0/ω_1)β/√(ω_1ω_2-β^2) ,    α_3=1/2logω_1ω_2-β^2/ω_0 ω_1 .

As a consistency check, note that with these identifications, the normalization factor of the final wave function becomes
√(ω_0/π)e^(α_1+α_3)/2=ω_1ω_2-β^2^1/4/√(π) ,
which correctly preserves the unit norm. 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 Gaussian wave functions.

Hence the total number of gates in the circuit U_1 in eq. exam1 is given by 
𝒟(U_1)=|α_1|+|α_2|+|α_3|
=1/2logω_1ω_2-β^2/ω_0^2+1/√(ω_0/ω_1) |β|/√(ω_1ω_2-β^2) ,
where we have assumed here that ω_1>ω_0 and ω_1ω_2-β^2>ω_0^2.[See further comments at the end of this appendix.] As in the main text, we refer to 𝒟(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. eq:Deg1Phys in order to give an equal complexity cost for the inverse gates Q_ij^-1 as for the original gates Q_ij, we count the appearance of Q_ij^-1 as one gate in a circuit. At a pragmatic level, this is required because α_2 is negative in our example, β=(ω_+-ω_-)/2<0. 

Note that in evaluating the exponents α_i in eqs. exp1 and exp2, 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
|∫d^2x ψ_T^† Q_22^α_3 Q_21^α_2 Q_11^α_1 ψ_R|^2 =1-χ ,

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 eq:pathPsi. This discussion is related to the choice of a tolerance  in eq. scotch, 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 applying the entangling gate before either of the scaling gates:
ψ_T=U_2ψ_R= Q_22^α̃_3 Q_11^α̃_1 
Q_21^α̃_2 ψ_R .
Note that for comparison purposes, our numbering of the exponents is such that they are associated with the same gates as appear in eq. exam1. The calculation proceeds essentially as above; in the end, we must match the coefficients 
ω_1=ω_0e^2α̃_1 ,    ω_2=1+^2α̃_2^2ω_0 e^2α̃_3 ,    β=α̃_2ω_0 e^(α̃_1+α̃_3) .
Solving for the exponents α_i then yields
α̃_1=1/2logω_1/ω_0 ,    α̃_2=1/β/√(ω_1ω_2-β^2) ,    α̃_3=1/2logω_1ω_2-β^2/ω_0ω_1 ,

and hence the circuit depth becomes
𝒟(U_2)=∑|α̃_i|=
1/2logω_1ω_2-β^2/ω_0^2+1/|β|/√(ω_1ω_2-β^2) .


Comparing the results in eqs. exp1 and exp2 with those in eq. exp3, we see that the exponents for the scaling gates are identical, α_1=α̃_1 and α_3=α̃_3, and only the exponent for the entangling gate has changed. Hence the circuit depth is almost identical to (<ref>), except that the second term lacks the factor √(ω_0/ω_1). If we assume ω_1>ω_0 as before, this implies that the present circuit will be slightly longer, 𝒟(U_2)>𝒟(U_1).

As a third simple example, let us consider instead applying the entangling gate after both of the scaling gates:
ψ_T=U_3ψ_R= Q_21^α̂_2 Q_22^α̂_3 
Q_11^α̂_1 
ψ_R .

Again we skip over the details of the calculation; we find that we must match the coefficients 
ω_1=ω_0e^2α̂_1 ,    ω_2=e^2α̂_3+^2α̂_2^2 e^2α̂_1ω_0  ,    
β=α̂_2ω_0 e^2α̂_1 .
Solving for the exponents α_i then yields
α̂_1=1/2logω_1/ω_0 ,    α̂_2=1/β/ω_1 ,
    α̂_3=1/2logω_1ω_2-β^2/ω_0ω_1 ,

and hence the circuit depth becomes
𝒟(U_3)=∑|α̃_i|=
1/2logω_1ω_2-β^2/ω_0^2+1/|β|/ω_1 .


Again, comparing with the exponents in eqs. exp1 and exp2 or in eq. exp3, we see that only the exponent for the entangling gate has changed. Hence the circuit depth here is similar to those for the previous two circuits, and whether the present circuit is longer or shorter depends on the values of the parameters ω_0, ω_1, ω_2 and β

Let us consider one more general example. Another interesting circuit would be
ψ_T=U_4ψ_R= Q_22^α̅_3 Q_21^α̅_2 
(Q_21^-1 Q_11)^α̅_1 
ψ_R .

Note that
(Q_21^-1 Q_11)^nψ(x_1,x_2)
=e^n /2ψe^nx_1-e^1-e^n/1-e^ x_2,x_2 ,
the derivation of which is as follows: first, consider
Q_11ψ=e^/2ψe^x_1,x_2 Q_21^-1Q_11ψ=e^/2ψe^x_1-e^x_2,x_2 .
Then acting with this combination twice yields
Q_21^-1Q_11^2ψ   =e^2/2Q_21^-1ψe^2 x_1-e^ x_2,x_2=e^2/2ψe^2x_1-x_2-e^ x_2,x_2

   =e^2/2ψe^2 x_1-e^e^+1x_2,x_2 .
And a third time:
Q_21^-1Q_11^3ψ=e^3/2Q_21^-1ψe^3 x_1-e^e^+1x_2,x_2=e^3/2ψe^3 x_1-e^e^2+e^+1x_2,x_2 .
Now the pattern is clear, and we deduce
Q_21^-1Q_11^nψ=e^n/2ψe^n x_1-e^∑_k=0^n-1 e^k x_2,x_2 .
Since
∑_k=0^n-1e^k=1-e^n/1-e^ ,
this becomes (<ref>), as claimed.

Now, acting with the circuit U_4 and matching coefficients as before, we find
α̅_1   =1/2logω_1/ω_0 ,

α̅_2   =1/√(ω_0/ω_1)β/√(ω_1ω_2-β^2)+
1/e^--11-√(ω_0/ω_1) ,

α̅_3   =1/2logω_1ω_2-β^2/ω_0ω_1 .
Expanding α̅_2 near ≈0, we have
α̅_2   =1/√(ω_0/ω_1)β/√(ω_1ω_2-β^2)-1/+1/2+1-√(ω_0/ω_1)+

   =1/[√(ω_0/ω_1)1+β/√(ω_1ω_2-β^2)-1]-1/21-√(ω_0/ω_1)+() .
We therefore find that the circuit depth for U_4 is
𝒟U_4=

1/2logω_1ω_2-β^2/ω_0^2+|1/[√(ω_0/ω_1)1+β/√(ω_1ω_2-β^2)-1]
-1/21-√(ω_0/ω_1)+()| .
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 of the geometric approach in the main text, it is useful to express the present results in terms of the normal-mode frequencies via eq. eq:omega12pm. If we focus our attention on the first circuit U_1 in eq. exam1, the exponents given in eqs. exp1 and exp2 become
α_1=1/2log_++_-/2ω_0 ,
  α_2=-1/√(ω_0/_++_-)_–_+/√(2_+_-) ,  α_3=1/2log2 _+_-/ω_0 (_++_-) .

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 <ref>, there are two natural hierarchies to consider: (i)  _+<_-<ω_0 or (ii) ω_0<_+<_-.[Implicitly, we chose the second hierarchy above in presenting our results in eqs. eq:Deg1Phys, horse0, and horse.] 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 a 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. Evaluating 𝒟(U_1)=|α_1|+|α_2|+|α_3| in these two cases yields:
𝒟(U_1)=1/√(ω_0/2 _++ ω_0/2 _-) _--_+/_++_-+1/2 ×{
   
logω_0^2/_+_-  for  ω_0>_+,_- ,

   
log_+_-/ω_0^2   for  ω_0<_+,_- .

 
.

Recall that _->_+ from eq. qm3. We may now compare this result with those derived using the geometric approach of section <ref>. In particular, if we recall the F_1 measure given in eq. (<ref>), the complexity would be given by
C=1/2|log_+/ω_0|+1/2|log_-/ω_0|={
   
1/2logω_0^2/_+_-  for  ω_0>_+,_- ,

   
1/2log_+_-/ω_0^2   for  ω_0<_+,_- .

 
.

Furthermore, recall that we should compare this with the coefficient of the 1/ϵ factor in the discrete calculations (see footnote <ref>). Hence we see the second contribution in eq. jjj precisely matches the complexity above. However, there is an additional positive term in 𝒟(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 eq:pathPsi. In this case, U_1 as given in eq. exam1 becomes
0≤s≤|α_1|/𝒟(U_1)    :        Y^11= 𝒟(U_1) , Y^22=Y^12=Y^21=0 ,


|α_1|/𝒟(U_1)≤s≤|α_1|+|α_2|/𝒟(U_1) :        Y^11= Y^22= Y^12=0 , Y^21=𝒟(U_1) ,


|α_1|+|α_2|/𝒟(U_1)≤s≤1    :      
  Y^11= 0 , Y^22=𝒟(U_1) , Y^12=Y^21=0 .
This form makes clear that the circuit consists of three separate “straight” segments, and so U_1 does not correspond to a geodesic path or an optimal circuit.



§ KILLING VECTORS AND MORE GEOMETRY


Inspecting the metric in eq. metric1, we can see three obvious Killing coordinates: y, τ, θ. When the penalty factors were introduced in section <ref>, we found that this is reduced to two Killing coordinates, y and z=(θ-τ)/2, in the geometry described by eqs. penalty1 or penalty2. However, by construction, all of these metrics are right-invariant, and hence the corresponding geometries must have one Killing vector for each generator Msimple, namely, four.[We thank Lucas Hackl for discussions on this point.] Furthermore, 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 it is a constant matrix), but rather is determined by the structure of eq. eq:metricRightInv.

One way to think of a Killing vector k^i is as providing a coordinate transformation
x^i →x^i+ k^i

which leaves the geometry or line element invariant. (Note  is just an infinitesimal parameter.) For example, eq. metric1 is certainly invariant under δτ = and so we write the corresponding Killing vector as
k^i∂_i = ∂_τ or k^i =δ^i_τ. 

So let us identify the coordinate transformations which generally leave eq. eq:metricRightInv invariant. For a general coordinate shift in eq. Umatrix, we have
δU = ∂_i U δx^i  .

As long as G_IJ is a constant matrix, all of the coordinate dependence is hidden in the one-forms U(s) U^-1(s)M^T_I, eq. penalty1. However, it is clear that these expressions are invariant if we right-multiply U by a global  transformation. Hence let us make the infinitesimal transformation: U→ U exp[^I M_I], where the ^I are (infinitesimal) constants. To leading order in these parameters, this reduces to
δU = U   M_I ^I .

Equating eqs. shift1 and shift2, we have
U M_I^I=_iUδx^i 
^I=U^-1_i U M_I^T δx^i ,

where we have assumed that we are working with the orthogonal basis of generators satisfying M_I M_J^T=δ_IJ, eq. Msimple. 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
δx^i=[U^-1_i U M_I^T]^-1^I
=k_I^i^I .
Thus we obtain four independent Killing vectors k_I= k_I^i∂_i.

Given our basis of generators in eq. Msimple and our parametrization of the circuit space in eq. (<ref>), we can easily compute k_I^i.  We can then identify the Killing vectors by simply reading off this matrix row-by-row:
k_1=   1/2_y-1/2sin(2z)_ρ-cos(2 z)/2sinh(2 ρ)  _x-cosh(2ρ)/2sinh(2 ρ)cos(2 z) _z ,

k_2=   1/2 cos(2 z)_ρ-sin(2 z)/2sinh(2 ρ) _x+1/2 (1-cosh(2ρ)/sinh(2ρ)  sin(2z) )_z ,

k_3=   1/2 cos(2 z)_ρ-sin(2 z)/2sinh(2 ρ)_x-1/2 (1+cosh(2ρ)/sinh(2ρ)  sin(2z) )_z ,

k_4=   1/2_y+1/2sin(2z)_ρ+cos(2 z)/2sinh(2 ρ)  _x+cosh(2ρ)/2sinh(2 ρ)cos(2 z) _z ,

where we are using the pseudo-lightcone coordinates of eq. (<ref>), with θ=x+z, τ=x-z.
One can explicitly verify that these indeed satisfy the Killing equations, 
0=∇_i k_I_j+∇_jk_I_i
=g_jℓ∇_i+g_iℓ∇_jk_I^ℓ ,
for either of the metrics in eqs. metric2 or penalty2. However, it is clear that eq. Killed does not organize the Killing vectors in the simplest way, so we define:
k̂_1≡      k_1+k_4 = _y ,

k̂_2≡   -k_1+k_4=sin(2z)_ρ+cos(2 z)/sinh(2 ρ)  _x+cosh(2ρ)/sinh(2 ρ)cos(2 z) _z ,

k̂_3≡      k_2+k_3 = cos(2 z)_ρ-sin(2 z)/sinh(2 ρ)_x-cosh(2ρ)/sinh(2ρ)  sin(2z) _z ,

k̂_4≡      k_2-k_3 =_z .


However, a simple inspection of the first metric (<ref>) reveals that _x is also an independent Killing vector, hence:
k̂_5≡_x .
This is an accidental symmetry that emerges with the choice G_IJ=δ_IJ.[We thank Lucas Hackl for discussions on the Killing symmetries.] However, as noted above, the four Killing vectors in eq. Killed2 apply for any (constant) choice of 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, c_I≡ (k̂_I)^i g_ij ẋ^j. We make use of these momenta in solving for the optimal circuits in sections <ref> and <ref>.



 §.§ AdS_3 geometry:
 

In section <ref>, we noted the appearance of a three-dimensional anti-de Sitter geometry in discussing the parametrization of U∈= ℝ×SL(2,ℝ).  Of course, the appearance of AdS_3 is natural since it is the universal cover of the SL(2,ℝ) subgroup. Here we would like to show how the AdS_3 geometry can be realized using the formalism introduced in section <ref>. In particular, we consider the geometry that results from the choice 
G_IJ=[ 1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0 ]
=η_IJ .

We have designated G_IJ=η_IJ because the I=3={12} and I=4={21} directions are null and hence the metric has a Minkowski signature. With this choice, eq. metric1 is replaced with
s^2   =η_IJ   U(s) U^-1(s)M^T_I U(s) U^-1(s)M^T_J

   =2y^2+2ρ^2-2cosh^2ρ τ^2+2sinh^2ρ θ^2 .

Hence we have produced precisely the AdS_3×ℝ geometry anticipated in section <ref>. Eq. adds describes the natural group invariant metric for , the left- and right-invariant metric, whereas the Euclidean metric metric1 is a less symmetric metric with only right-invariance. 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, zero cost. This would allow the construction of arbitrarily low-complexity circuits simply by deforming the circuit along the null directions.[Of course, moving in a timelike direction also yields a negative cost.] 




 §.§ Alternate basis of generators:
 

The basis of matrix generators in eq. Msimple is natural in the sense that it straightforwardly extends from the problem of complexity in the case of two coupled harmonic oscillators to the case of N coupled oscillators. However, this is not the most convenient basis for certain calculations in section <ref>. Hence for the interested reader, we describe here an alternate basis of generators which simplifies some of the calculations. In particular, consider the Pauli-like basis:
M_1   =1/√(2)[ 1 0; 0 1 ]=1/√(2) 1 ,      M_2=1/√(2)[  1  0;  0 -1 ]=1/√(2) σ_3 , 


M_3   =1/√(2)[ 0 1; 1 0 ]=1/√(2) σ_1 ,    M_4=1/√(2)[  0  1; -1  0 ]=-i/√(2) σ_2
 .

The normalization of the generators is still given by (M_I M_J^T)=δ_IJ. In fact, the new generators are easily related to the original generators in eq. Msimple by an orthogonal transformation: M_I = R_I^J  M_J with R_I^J∈ O(2)× O(2)∈ O(4). In this new basis, the M_2,3,4 generators naturally form the sl(2,ℝ) subalgebra, with
[_2,_3]=√(2) _4 ,   
[_2,_4]=√(2) _3 ,   
[_3,_4]=-√(2) _2 ,
while M_1 describes the remaining fibre over ℝ in the  group.

With this new basis, the Killing vectors which emerge from the right-invariance of the metric naturally appear in the form given in eq. Killed2. 
One can easily show that working with these new generators, the metrics appearing in eqs. metric1 and penalty1 are unchanged, the corresponding G_IJ are left unchanged by the rotation R_I^J introduced above. Additionally, the AdS_3 geometry in eq. adds now results from the choice G_IJ=η_IJ=diag(1,1,1,-1).



§ NORMAL-MODE FREQUENCIES _K
 

The derivation of the normal-mode frequencies in eq. (<ref>) – or eq. freakd for a lattice of coupled oscillators – is straightforward, and can be found in a number of different sources, any elementary condensed matter textbook. For completeness, we briefly review the result (<ref>) for the periodic one-dimensional lattice discussed in section <ref>. Essentially, we need only apply the inverse Fourier transform 
x_a≡1/√(N)∑_k=0^N-1exp2πi k/N ax̃_k ,

to re-express the Hamiltonian qm88 in terms of the normal modes, qm288. In particular, we focus on the potential
V=1/2∑_a=0^N-1[ω^2x_a^2+Ω^2x_a-x_a+1^2] .


Considering the second term involving the coupling between the oscillators, we find

    Ω^2∑_a=0^N-1 x_a-x_a+1^2
       =Ω^2∑_a=0^N-1[1/√(N)∑_k=0^N-1exp2π ik a/Nx̃_k1-exp2π ik/N]^2
    
       =Ω^2/N∑_a,k,k'exp2π i(k+k')a/Nx̃_kx̃_k'1-exp2π ik/N1-exp2π ik'/N
    
       =Ω^2∑_kx̃_kx̃_-k1-exp2π ik/N1-exp-2π ik/N
       =2Ω^2∑_kx̃_kx̃_-k1-cos2π k/N
    =4Ω^2∑_k|x̃_k|^2 sin^2π k/N ,

where in going to the third line we applied the normalization condition (<ref>), and in the last step we used x̃_kx̃_-k=x̃_kx̃^†_k. Here all sums run from 0 to N-1. The Fourier transform of the first term in the potential (<ref>) is trivial, and thus we find
V=1/2∑_k=0^N-1[ ω^2+4Ω^2 sin^2πk/N]  |x̃_k|^2
=1/2 ∑_k=0^N-1_k^2 |x̃_k|^2 .
Hence we have identified the desired normal-mode frequencies,
_k^2=ω^2+4Ω^2 sin^2πk/N ,

(<ref>). If instead we were examining a d-dimensional free scalar field, the lattice would be extended to d-1 (spatial) dimensions, whereupon the corresponding normal-mode frequencies become
_k⃗^2=ω^2+4Ω^2 ∑_i=1^d-1sin^2πk_i/N ,

where k_i are the components of the momentum vector 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.



§ CHANGE OF BASIS
 

In this appendix, we would like to extend the discussion around eq. bark2 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 <ref>, where we presented a tentative argument that the metric on the normal-mode subspace is flat, s^2_n-m=|ỹ_0|^2+|ỹ_1|^2+⋯+|ỹ_N-1|^2 ,

(<ref>). However, we also noted that, at the time, this conclusion was somewhat premature, since implicitly we applied eq. metric5, which defines the metric in the position basis, to a calculation with the diagonal circuit raffle4 written in the normal-mode basis. That is, in eq. metric5, 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. raffle4, we would write 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 coordinates x̃_k. In other words, in eq. raffle4, where we are restricting our attention to the normal-mode subspace, we are considering the diagonal generators Ỹ^kℓ=δ^kℓ ỹ_k. 

Hence to show that the result metric6xx is correct, we must take care to translate between the two bases of generators discussed above. As in eqs. nmU and straight9, we can transform from generators acting in the normal-mode basis to those in the position basis via[Note that this transformation removes the tilde from M but not from the index. For example, the new generator M_kk acts on Gaussian wave functions written in terms of the oscillator position coordinates x_a, but still has the effect of scaling the kth normal mode x̃_k.] 
M_Ĩ = ^†  M̃_Ĩ   .

Implicitly, the normal-mode generators M̃_Ĩ have the same form as that given in eq. matrix, namely
[M̃_kℓ]_pq=δ_kpδ_ℓq ,

where we have denoted Ĩ=(kℓ) with momentum labels k, ℓ. 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[Note that the complex conjugation appears on the first factor in =^*⊗ because our convention is that written in terms of the normal modes, the Gaussian wave functions involve both x̃_k and x̃_k^†, the appearance of |x̃_k|^2 in eq. targetk.]
[M_kℓ]_ab    =    [^†]_ap [M̃_kℓ]_pq []_qb
=[^†]_a k   []_ℓb


   =    [^†]_c k   []_ℓd [M_cd]_ab


       M_Ĩ = []_Ĩ  J  M_J     with  =^*⊗ .

In going from the second to third line, we used eq. matrix and identified Ĩ=(kℓ) and J=(cd). This equation generalizes eq. bark2 from  in the previous section to the case of  studied here. Furthermore, given the properties of , one can easily see that the matrix  is a unitary matrix. Hence we can invert the transformation in eq. lot2 to write M_I =[^†]_IJ̃  M_J̃. 

Similarly, we can invert the transformation in eq. lot0, transform from generators acting in the position basis to the normal-mode basis via
M̃_ I= M_I ^† ,

and combine this expression with eq. matrix, [ M_ab]_cd =δ_acδ_bd, to write
[M̃_ab]_kℓ    =    []_kc [M_ab]_cd [^†]_dℓ


=[]_ ka   [^†]_bℓ


   =    []_pa   [^†]_bq  [M̃_pq]_kℓ


       M̃_I = [^†]_I  J̃  M̃_J̃     with  ^†=^T⊗^† .

In going from the second to third line, we used eq. lot1 and identified I=(ab) and J̃=(pq). As before, we can easily invert the transformation in eq. lot2a to write M̃_Ĩ =[]_Ĩ J M̃_J. As our notation indicates,  is precisely the unitary matrix appearing in eq. lot2, and hence it also plays a role in transforming the generators acting in the normal-mode space. 

Hence by using the special structure of the generators in eqs. matrix and lot1, we have re-organized the  transformation acting on the matrix indices in eqs. lot0 and lot0a to a transformation acting on the generator labels in eqs. lot2 and lot2a, respectively.

With these tools in hand, let us consider re-expressing the cost function cost5 or the metric metric5 in terms of the normal-mode basis, using eq. lot0. Here, we show the calculation for the metric; the transformation of the cost function follows in a similar manner. Beginning with the differential Y^I= U U^-1  M^†_I defined in eq. metric5, we transform the circuit to the normal-mode space via U=^† Ũ, which yields 
Y^I=Ũ  Ũ^-1   M^†_I  ^†=Ũ  Ũ^-1 M̃^†_I =[]_Ĩ  I Ỹ^Ĩ

where Ỹ^Ĩ=Ũ Ũ^-1 M̃^†_Ĩ,
and we have employed eqs. lot0a and lot2a in the second and third equalities. Hence the metric metric5 transforms as
s^2=[]_Ĩ  I δ_IJ  [^†]_J J̃  Ỹ^Ĩ 
(Ỹ^J̃)^* =δ_Ĩ J̃  Ỹ^Ĩ 
(Ỹ^J̃)^*  .

Note that here we are using the fact that  is a unitary matrix. Thus 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.[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.] Of course, the same is true of the cost function cost5, it can also be written as
𝒟(U)=∫_0^1s√(δ_ĨJ̃ Y^Ĩ(s) ( Y^J̃(s))^*) ,     where  Y^Ĩ(s)=_sŨ(s) Ũ^-1(s)M̃^†_Ĩ  .



Note that this transformation is slightly different than that expressed in eq. metric1a for the metric for two coupled oscillators. In the latter case, we are considering the metric to still be in the position basis but evaluated with a different basis of generators. The same invariance holds here for a lattice of oscillators, as can be seen by applying eq. lot2 directly to the metric metric5 to produce
s^2=δ_Ĩ J̃ Y^Ĩ 
(Y^J̃)^* ,     where    Y^Ĩ=U U^-1  M^†_Ĩ .

Of course, the same change of basis could also be performed with eq. lot2a when working in the normal-mode space.



 §.§ General cost functions
 

In eq. Dalpha, 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 that arise from using the original oscillator position basis. That is, we could equally well define cost functions with
𝒟_κ=∫_0^1s∑|Y^I(s)|^κ .
In the discussion of the F_2 cost function in the previous section, 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. solver5, for which the velocity components in the position basis become
Y^11=Y^22=y_1 ,    Y^12=Y^21=-ρ_1 .

These two factors are written in terms of the normal-mode frequencies in eq. fin2, but we can re-express these results as
y_1=1/4(log_-/ω_0+log_+/ω_0) ,    ρ_1=1/4(log_-/ω_0-log_+/ω_0) .

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. DalphaP for a few values of κ,[In the case that ω_0>_±, one should replace
_±/ω_0→ω_0/_∓ in these formulae. Note this substitution only really changes the results for odd κ.]
𝒟_κ(U_0)   =2y_1^κ+2ρ_1^κ=
log_-/ω_0   for   κ=1 ,

1/4(log^2_-/ω_0+log^2_+/ω_0)   for   κ=2 ,

1/16(log^3_-/ω_0+3log_-/ω_0 log^2_+/ω_0)   for   κ=3 ,

1/64(log^4_-/ω_0+log^4_+/ω_0+6log^2_-/ω_0 log^2_+/ω_0)   for   κ=4 .


Hence we see that it is only for κ=2 that we reproduce the cost found using eq. Dalpha in the normal-mode basis,
𝒟_κ(U_0)≃log^κ(_-/ω_0)+ log^κ(_+/ω_0).

These differences in the cost can be understood using the approach developed to implement a change of basis for a lattice of oscillators in section <ref>. In particular, transforming from the position basis to the normal-mode basis can be described in terms of the unitary matrix  defined in eq. lot2. Given the definition of the velocity components in eq. cost5, we then have
Y^Ĩ=[^*]_Ĩ J Y^J ,

or, inverting this expression, Y^I=[^T]_IJ̃ Y^J̃. Furthermore, recall that the quadratic construction δ_IJ Y^I (Y^J)^*= δ_ĨJ̃ Y^Ĩ (Y^J̃)^* 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 general metric G_IJ, which transforms nontrivially under the change of basis, G_ĨJ̃= []_J̃ J G_IJ [^†]_IĨ ,

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. DalphaP to allow for general penalty factors by writing
𝒟_κ=∫_0^1s∑_I_1,I_2,⋯,I_κ
G_I_1I_2⋯I_κ 
|Y^I_1(s)|  |Y^I_2(s)|⋯|Y^I_κ(s)| ,
where G_I_1I_2⋯ I_κ is a symmetric tensor with κ indices.  In eq. Dalpha, we are implicitly considering simple “penalty” tensors of the form
G_I_1I_2⋯I_κ   =   δ_I_1I_2 δ_I_2I_3⋯δ_I_κ-1 I_κ     for  κ≥2 ,


G_I   =   1                    for  κ=1 .

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. In addition to the simplicity of this case, recall that this was also the cost function favoured in the comparison to holographic complexity in section <ref>.

Let us begin with the case of N=2, in which case the transformation matrix =_2 takes the simple form given in eq. bark2. For κ=1, the penalty tensor tens2 becomes the four component vector
G_I=(1,1,1,1) ,

which is actually an eigenvector of . Hence if we transform as in eq. chan1, we find the rather surprising result that
G_Ĩ = _Ĩ J  G_J= (2,0,0,0) .

That is, expressing our κ=1 cost function DalphaP in terms of the normal-mode basis, we are only penalizing the scaling gate  associated with x̃_+! The other (normal-mode) gates can be inserted in the circuit at zero cost. However, we must add that the transformation in eq. tens3 is slightly naive since it assumes that the absolute values in the cost function DalphaP play no role, 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 Ĩ=+-,-+,[This is the case in eq. veloc for the optimal circuit with _±>ω_0. Hence in eq. trials, the κ=1 cost function only depends on _-.] one finds G_Ĩ = (0,0,0,2), only the scaling gate associated with x̃_- is penalized. 

Similar results arise if we begin with the κ cost functions Dalpha 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 tensors become 
G_Ĩ=(1,1,1,1)  ⟶  G_I = ^T_IJ̃  G_J̃= (2,0,0,0) .

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. tens5 assumes that the absolute values in the cost function Dalpha 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, all of the scaling components have a definite sign and all components in the entangling directions vanish. 

Furthermore, using eq. Dalpha, we might note that the cost of our straight-line circuit is simply
𝒟_κ=1(U_0) =ỹ_++ỹ_-=1/2log_+/ω_0 +
1/2log_-/ω_0 ,

again assuming _±>ω_0. Here we emphasize that since only two of the velocity components were non-vanishing, namely Y^++ and Y^–, we would arrive at the same cost for a family of penalty tensors of the form
G_Ĩ=(1,_1^2,_2^2,1) .

In this case, transforming to the position basis as in eq. tens5 yields
G_I = ^T_IJ̃  G_J̃=1/2 (2+_1^2+_2^2, _2^2-_1^2,   _1^2-_2^2, 2-_1^2-_2^2) .

At first sight, this result seems to yield a more reasonable penalty tensor relative to eq. tens5. 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, _1=_2=, whereupon eq. tens7 becomes
G_I = (1+^2, 0,   0, 1-^2) ,

which still requires that ^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 in section <ref>; however, we should say that some of these issues arise because we focused on the simple case of κ=1. For example, if we consider instead the κ=2 cost function DalphaP with our penalized metric pen99, then transforming to the normal-mode basis yields eq. pen99a. 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, and since these have not been changed by the transformation, 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 <ref> using the κ cost functions Dalpha constructed in the normal-mode basis remain unchanged when working with eq. DalphaP in the position basis. For a (d-1)-dimensional spatial lattice of oscillators, the κ=1 penalty tensor in eq. tens2 becomes
G_Ĩ=(N^d-1,0,0,⋯,0)

in terms of the normal modes. Hence as in eq. tens4, 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.[Note that in eq. tens4, the penalty associated to the M_++ was increased to 2, the number of oscillators.] The cost for the straight-line circuit then becomes
D_κ=1(U_0)=1/2 N^d-1 |logm/ω_0|=V/2 δ^d-1 |logm/ω_0| .

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, ω_0=e^-σ/δ, but it does not appear if the reference frequency is chosen as an IR frequency. We observe that this is the opposite of the situation discussed in section <ref>.



§ APPROXIMATING THE COMPLEXITY


In section <ref>, we compared our result complexityNd for the complexity of the ground state of 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, mδ≪1. In particular, this led us to consider the generalized family of κ cost functions given in eq. Dalpha, which yields
=1/2^κ∑_{k_i}=0^N-1|logω̃_k⃗/ω_0|^κ .

To identify the leading contribution to either eq. complexityNd or nothing, 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 nothing in more detail. Our result (<ref>) is still an approximation, but it allows us to consider the subleading contributions to eq. calpha. In particular, we will determine the leading corrections involving the mass. 

First we substitute the normal-mode frequencies freakd into the above expression for general κ to find
=1/4^κ∑_{k_i}=0^N-1|logm^2/ω_0^2+2/ω_0δ^2 ∑_j=1^d-1sin^2πk_j/N|^κ
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. atall then becomes
_IR=V/4^κ∫_0^Λ_IR
^d-1p/(2π)^d-1 [logm^2+p^2/ω_0^2]^κ .

Note that we have dropped the absolute value symbol here since 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. atall. 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 <ref>. Choosing this cut-off to be Λ_IR∝ m, the IR contribution takes the general form
_IR=V m^d-1 ∑_a=0^κc_a [logm/ω_0]^a ,

where the numerical coefficients c_a are independent of m and ω_0, but will depend on the spacetime dimension d. The leading contribution then takes the form
_IR≃ c_κ V m^d-1[log m/ω_0]^κ.

Having isolated the IR contribution, we return to the UV contributions to eq. atall. In these remaining terms, we can consider m^2/ω_0^2 to be a small correction to the argument of the logarithm, and so we perform a Taylor series expansion and keep only the first correction in mδ:
_UV   ≃   1/4^κ∑_{k_i}>IR
log[2/ω_0δ^2 ∑_j=1^d-1sin^2πk_j/N]^κ

                × {1+κ (mδ)^2/4(∑_j=1^d-1sin^2πk_j/N 
log[2/ω_0δ^2 ∑_j=1^d-1sin^2πk_j/N]
 )^-1} ,
where the notation ∑_{k_i}>IR indicates that the summation begins at the IR cut-off, |k⃗|≥Λ_IR Nδ/(2π). The leading term will produce the expected leading contribution in eq. calpha 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 proceed further, we focus on the case κ=1, which simplifies the calculations slightly and was also the case that we found best emulated the holographic complexity. The following analysis is essentially unchanged for larger values of κ. Substituting κ=1 into eq. (<ref>) yields
_UV≃1/4∑_{k_i}>IR{log[2/ω_0δ^2 ∑_j=1^d-1sin^2πk_j/N]
+(mδ)^2/4[∑_j=1^d-1sin^2πk_j/N]^-1} .

We now examine the two sums separately. First, we break the leading sum into two:
log4/ω_0^2δ^2∑_j=1^d-1sin^2πk_j/N=2log2/ω_0δ

+log∑_j=1^d-1sin^2πk_j/N .
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 Λ_IR^d-1). 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∑_{k_i}>IRlog∑_j=1^d-1sin^2πk_j/N=a_d-1N^d-1+a_d-3N^d-3+⋯+a_0 ,
where the a_i are fixed numerical coefficients. Note that the constant term a_0 appears for both odd and even d, and the former may also have a logarithmic correction (log N). Turning now to the second sum in eq. (<ref>), we again carried out numerical fits to find
1/4∑_{k_i}>IR[∑_j=1^d-1sin^2πk_j/N]^-1≃b_d-1N^d-1+b_d-3N^d-3+⋯+b_0 .


Collecting results, our approximation to the total complexity for κ=1 (and assuming ω_0<m) is therefore
   ≃   V m^d-1 (c_IR 1 logm/ω_0+ c_IR 0 )+N^d-1/2log2/ω_0δ

+a_d-1N^d-1+a_d-3N^d-3+⋯+a_0

              +mδ^2b_d-1N^d-1+b_d-3N^d-3+⋯+b_0 .

To make a comparison with holographic complexity, as in section <ref>, we substitute N^d-1=V/δ^d-1, and introduce L=V^1/(d-1) as the linear size of the lattice. The complexity exit then becomes
   ≃   
1/2V/δ^d-1log2/ω_0δ
+V/δ^d-1a_d-1+a_d-3δ^2/L^2+⋯+m^2V/δ^d-3b_d-1+b_d-3δ^2/L^2+⋯

            +V m^d-1 (c_1 logm/ω_0+ c_0 ) .

Hence we find the expected leading term, which corresponds to the result in eq. calpha with κ=1. We also find a subleading term proportional to V/δ^d-1, as is found in holographic complexity in the CA proposal <cit.>. Additionally, we would highlight the correction proportional to m^2V/δ^d-3, for which analogous results can again be found in holographic calculations—see section <ref> for further comments. It is interesting that we also see corrections, 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, 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. exit2 reduces to
=1/2L/δlog2/ω_0δ
+a_1L/δ+a_0+L m (c_1 logm/ω_0+ c_0 ) ,

where we have replaced V=L to emphasize that the volume is only a linear length here.



§ OPTIMAL GEODESIC FOR PENALIZED GEOMETRY


We would like to find the optimal geodesic in the penalized geometry penalty2, but as commented below eq. eq:kPenFull, 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 solver5 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. segments was shorter than this geodesic when the penalty factor was large, ≫ρ_1,y_1; eq. compar2.
 
To make progress towards finding the optimal geodesic in the new geometry, we make a simplifying assumption. To begin, we examine the penalized metric penalty2 and observe that as the radius ρ increases, the fastest growing component of the metric is g_zz∼^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 motion on a constant-z submanifold:
s^2=2y^2+2[^2-^2-1sin^22x]ρ^2+2 ^2x^2 .
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 this appendix. 

Working with this simpler geometry (<ref>), 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
2 c̅_1≡2ẏ ,    2a  c̅_2≡2[ ^2-^2-1sin^22 x]ρ̇ ,
where the factors of 2 and  were chosen to simplify expressions below, and we have used the notation c̅_i to avoid confusion with the ĉ_i in eq. revel0. As before, the first constraint gives the usual solution yrun for y, y(s)=y_1 s with c̅_1=y_1. The second constraint yields
ρ̇= c̅_2/^2-^2-1sin^22 x .

The normalization of the tangent vector then becomes
k^2   =2ẏ^2+2[ ^2-^2-1sin^22 x]ρ̇^2+2^2ẋ^2

   =2y_1^2+2^2 c̅_2^2/^2-^2-1sin^22 x+2^2ẋ^2 .
It is possible to integrate this equation to find s(x). To simplify the subsequent equations, we shall define k̅ via
2k̅^2≡k^2-2y_1^2 .

Isolating ẋ in eq. eq:ktemp then yields
s/x=√(2)a[ 2k̅^2-2^2  c̅_2^2/^2-^2-1sin^22 x]^-1/2 .

Upon integrating and choosing the constant of integration such that s(x_1=π/2)=1, the result can be simplified to
s(x)=1-^2/2√(k̅^2-c̅_2^2) Π-(^2-1);-f(x) | ^2-1c̅_2^2/k̅^2-c̅_2^2 ,
where
f(x)≡1/sin^-1(2x)√(^2+tan^2(2x)) ,
and Π is the incomplete elliptic integral of the third kind, which we write here as
Π(-n;-z|m)=-∫_0^zt/1+nsin^2t√(1-msin^2t) .
We now combine the expressions for ρ̇ and 1/ẋ in eqs. revel2 and revel3 to find 
ρ/x=ρ̇/ẋ=√(2)^2 c̅_2/^2-^2-1sin^22 x[k̅^2-2^2  c̅_2^2/^2-^2-1sin^22 x]^-1/2 .

This expression can likewise be integrated to obtain
ρ(x)=ρ_1-ic̅_2/2√(k̅^2-c̅_2^2) Fi g(x)| 1-^2-1c̅_2^2/k̅^2-c̅_2^2 ,
where
g(x)≡1/sinh^-1(2x) ,

and F is the incomplete elliptic integral of the first kind,
F(z=ix|m)=i ∫_0^xτ/√(1+msinh^2τ) ,    x∈ℝ .
Note that in eq. eq:rX, we have fixed the integration constant via the boundary condition ρ(x_1=π/2)=ρ_1. 

Now, unfortunately (<ref>) cannot be inverted to find an analytical expression for x(s), which we could then use to obtain ρ(s) via eq. (<ref>). However, we can still study the behaviour of these geodesics numerically. Before doing so, it remains to relate the parameters c̅_2 and k (or k̅) 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. (<ref>) is always complex; hence in order for ρ(x) to be real, the coefficient must also be imaginary, which requires
k̅^2>c̅_2^2  k^2>2y_1^2+c̅_2^2 .
Turning now to s(x) in eq. (<ref>), the elliptic integral Π in this case is always real, and the coefficient will also be real in precisely the same regime (<ref>). 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, namely s(x=π/2)=1 and ρ(x=π/2)=ρ_1. For the optimal geodesic, we further choose the boundary condition x(s=0)=π/4, which minimizes the cost of motion in the ρ direction, eq. rats6.[Alternatively, we could choose x(s=0)=3π/4, but the resulting trajectory is simply of a copy of the present geodesic rotated 180^o around the (θ,τ)=(π,0) axis—see figure <ref>.] However, we must be careful in evaluating eqs. eq:sX and eq:rX at this value of x; in particular, we must consider the limits
lim_x→π/4^+s(x)   =1-^2/2√(k̅^2-c̅_2^2) Π-(^2-1)| ^2-1c̅_2^2/k̅^2-c̅_2^2 ,

lim_x→π/4^+ρ(x)   =
ρ_1- c̅_2/2√(k̅^2-c̅_2^2) K^2-1c̅_2^2/k̅^2-c̅_2^2 ,
where F and K are the complete elliptic integrals of the first and third kind, respectively, defined via
K(z)=Fπ/2 | z ,    Π(n | m)=Πn;π/2 | m .
The parameters c̅_2 and k̅ 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 is not required: we seek only a valid case to compare with the length  (<ref>) of the simple straight-line geodesic.

To that end, observe that the elliptic integral K is of order 1 almost everywhere, except when the argument approaches 1 in where it diverges, lim_w→1K(w)=∞. Since we want ρ to be large, let us choose
^2-1c̅_2^2/k̅^2-c̅_2^2=1-    k̅^2=c̅_2^2 ^2-/1- ,
where 0<≪1. Note that this is within the reality domain (<ref>) since >1. The boundary condition that the limits (<ref>) should vanish then allows us to solve for ρ_1 and c̅_2; one finds:
ρ_1=/2√(1-/^2-1) K1- ,    c̅_2=^2/2√(1-/^2-1) Π1-^2 | 1- .
One can then make ρ_1 arbitrarily large by taking →0; note that c̅_2 becomes arbitrarily large in the same limit. In fact, the divergence in both cases is logarithmic:
ρ_1   =/√(^2-1)1/4log1/+log2+() ,

c̅_2   =1/4√(^2-1)log1/+1/2tan^-1√(^2-1)+log2/√(^2-1)+() ,

where higher-order terms vanish in the limit →0. 

We may now numerically compare the length of this geodesic to the proposed minimum (<ref>) associated with the straight-line circuit. Substituting c̅_2 and ρ_1 from eq. (<ref>) into k̅ given eq. (<ref>) and the analogous quantity k̅_0 from eq. (<ref>), we find
k̅=^2/2√(^2-/^2-1) Π1-^2 | 1-    and    k̅_0=^2/2√(1-/^2-1) K1- ,
where 2k̅_0^2≡ k_0^2-2y_1^2. 

Of course, while these expressions are well-suited to numerics, we would also like to express k̅ in terms of the coordinates ρ_1, y_1, so as to compare with (<ref>) on more physical footing. We can obtain an approximation of this form by first replacing c̅_2 in (<ref>) by its expression in (<ref>), and then expanding for →0:
k̅=/√(^2-1)1/4log1/+log2+/2tan^-1√(^2-1)+() , 
where as above the () terms vanish as →0, and we shall drop them henceforth. Comparing this expression to eq. (<ref>), we observe that we can equivalently write this as
k̅≃/2tan^-1√(^2-1)+ρ_1
 ≃ π/4+ρ_1-1/2-1/12^2+1/^4 ,
where in the second approximation we have performed an expansion in the limit →∞.

Additionally, it will be interesting to compare these geodesics against the segmented trajectory described in eq. segments. The length of this path is given by eq. eq:ks, and so as in eq. newerX, we define 2k̅_s^2 =k_s^2 -2 y_1^2, 
k̅_s   =1/4[π^2^2+8ρ_12ρ_1+√(π^2^2+ (4y_1)^2)]^1/2

   =/4[π^2+41-/^2-1^1/2K1- 1-/^2-1^1/2 K1-+π1+16y_1^2/π^2^2 ^1/2]^1/2 ,

where in the second line we have replaced ρ_1 using eq. (<ref>). Note that unlike k̅ and k̅_0 in eq. (<ref>), the parameter y_1 still appears in this expression—although this contribution is suppressed for ≫ y_1. Again however, the second line above is more suited to numerics than physical inspection; to compare with (<ref>), we shall expand with ≫ρ_1,y_1 (as well as ρ_1,y_1≫1, and assuming ρ_1 and y_1 are roughly the same order of magnitude). Hence:
k̅_s=π/4[1+8ρ_1/π2ρ_1/π+√(1+4y_1/π^2)]^1/2
≃π/4+ρ_1+8/π^2y_1^2/^2 ρ_1+… .
We mentioned above that the segmented path constitutes a remarkably good approximation to the geodesic. Comparing k̅ in (<ref>) and k̅_s in (<ref>), 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 . We discuss this point further in the main text—see eq. eq:overhead. We also explicitly confirm that the two paths are very close to one another in the large  regime by examining x(s) and ρ(s) numerically, as shown in figure <ref>.

In closing this appendix, we remind the reader that in order to make progress, we confined our attention to motion in the constant-z subspace given by the simpler metric (<ref>). Hence for completeness, we should go back and examine whether or not this was a reasonable assumption.  In particular, we wish to argue that, at least in the limit ≫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 (<ref>). 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. (<ref>). Recall that τ_0=0, and hence z_0=x_0=π/4. Then the finiteness condition (<ref>) requires that we set ĉ_3=0.[Note that ĉ_2 is still free, since we can rewrite eq. (<ref>) as ĉ_3=^2ĉ_2(2z_0)|_z_0=π/4=0.] Along the initial segment, where x=π/4, the derivatives (<ref>) then reduce to
ẋ|_x=π/4=-ĉ_4cosh2ρ/^2 ,    ż|_x=π/4=-ĉ_4/2^2 ,    ρ̇|_x=π/4=ĉ_2/2 .
Therefore, in the large  limit under consideration, 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 (1) in , but ż is still exponentially suppressed in ρ_1 relative to ẋ. (Meanwhile ρ̇ decreases sharply to 0 on this segment in the limit ≫1.) Thus geodesics in the full spacetime (<ref>) can indeed be approximated by those in the constant-z subspace (<ref>), at least in the limits that we are considering. 



ytphys

