Circuit Complexity for Coherent States

We examine the circuit complexity of coherent states in a free scalar field theory, applying Nielsen's geometric approach as in [1]. The complexity of the coherent states have the same UV divergences as the vacuum state complexity and so we consider the finite increase of the complexity of these states over the vacuum state. One observation is that generally, the optimal circuits introduce entanglement between the normal modes at intermediate stages even though our reference state and target states are not entangled in this basis. We also compare our results from Nielsen's approach with those found using the Fubini-Study method of [2]. For general coherent states, we find that the complexities, as well as the optimal circuits, derived from these two approaches, are different.


Introduction
In recent years, a new bridge has begun to develop connecting quantum information theory to quantum gravity and quantum field theory. In particular, understanding the relation between quantum entanglement and the emergence of semi-classical spacetime geometry [3][4][5] has become an active field of research. Gauge/gravity duality [6][7][8] has been the central arena for the exploration of these connections and much of the understanding of the connection between entanglement and geometry has come from investigations of holographic entanglement entropy [9][10][11]. However, it has become clear that holographic entanglement entropy is not able to probe the bulk spacetime far behind the event horizon of black holes [12,13]. Inspired by this problem, Susskind [13][14][15] proposed the study of new bulk observables, which he conjectured should be the gravitational dual of the circuit complexity in the boundary theory. In particular, there are two proposals for 'holographic complexity': complexity=volume (CV) [15,16] and complexity=action (CA) [17,18]. The CV conjecture states that the complexity of the boundary state is proportional to the volume of an extremal codimension-one surface extending the boundary time slice into the bulk. The CA conjecture identifies the complexity of the boundary state with the gravitational action evaluated at special bulk region called the Wheeler-DeWitt patch, i.e., the causal development of the bulk surface identified in the previous approach. Both conjectures bring to our attention new gravitational observables which contain information about the spacetime region deep behind the black hole horizon and they have been vigorously investigated in the recent literature, e.g., .
An obstacle in this program is that precise comparisons with the boundary theory are not yet possible because we still do not have a precise definition of circuit complexity for states in quantum field theory. However, some preliminary steps towards developing such a definition have been taken in the past year, e.g., [1,2,[40][41][42][43][44][45][46][47][48][49][50]. In particular, refs. [1] adapted the approach of Nielsen and his collaborators [51][52][53] to translate the task of finding the complexity of the ground state of a free scalar field theory into a geometric problem of finding optimal geodesics in an associated geometry. As similar geometric approach was developed for this question in [2] based on the Fubini-Study metric. 1 In fact, both approaches produced the same simple circuit to prepare the vacuum state for a simple (unentangled) reference state and assigned the same complexity to the vacuum. In these calculations the field theory must be regulated since the complexity is dominated by contributions from the high energy modes and the result is UV divergent. However, an interesting agreement was found in comparing the structure of these divergences with those appearing in holographic complexity. In particular, the leading divergence found for holographic complexity [21] takes the form where V is the spatial volume, δ is the short-distance cutoff, and d is the spacetime dimension of the boundary theory. The scale is undetermined and arises because of ambiguities in defining the gravitational action on regions with null boundary segments [23,38]. An analogous ambiguity appears in evaluating the complexity for the scalar field theory because an undetermined scale must be introduced to define the reference state, and the leading divergence of the vacuum has precisely the same form as shown above for C A [1,2]. Of course, in either calculation, the interesting logarithmic factor can be removed by choosing ∼ δ and so this does not rule out the CV conjecture.
In this paper, we are extending the investigations of complexity in refs. [1,2] by examining the complexity of excited states in the free scalar field theory. In particular, we develop the additional techniques needed to evaluate the complexity of coherent states in which the scalar field acquires a nonvanishing expectation value. An exploratory investigation of the complexity of coherent states already appears in [47], however, the analysis there differs in many essential ways from our approach and there is no substantive overlap between the previous work and the present paper, as we will describe in more detail in section 6. Here, we might add that the complexity of excited fermionic states was considered in [41]. But this was a special case where the excited states were still Gaussian states and so no new ingredients were needed beyond those needed to evaluate the complexity of the vacuum. Further, refs. [44,45] examined the complexity of excitations for the free scalar produced by a quench of the mass term. However, again these excited states could be assessed using the same techniques used to evaluated the vacuum complexity. To prepare the coherent states, we must introduce a new class of gates in our circuits and in particular, this requires introducing a new (undetermined) scale into our model for the complexity. We develop the extended geometry associated with this larger gate set for both the Nielsen and Fubini-Study approaches and the resulting optimal circuits and complexities exhibit a number of interesting features. For example, we find that the optimal circuits introduce entanglement between the normal modes at intermediate stages even though our reference state and target states are not entangled in this basis. Further for general coherent states, we show that the complexities, as well as the optimal circuits, derived by Nielsen and Fubini-Study approaches are different.

Nielsen, geometry and complexity
In this section, we briefly review the salient ideas required to apply Nielsen's geometric approach to circuit complexity [51][52][53] to evaluate the complexity of state in a quantum field theory, as developed in [1]. In this setting, complexity is a measure of the difficulty or cost to prepare the particular target state |ψ T starting with a certain simple reference state |ψ R . We are using a quantum circuit model where the preparation is accomplished by applying a series of elementary unitaries, chosen from a particular set of gates {g 1 , · · · , g N }. That is, 2 |ψ T = U T |ψ R = g in · · · g i 2 g i 1 |ψ R , (1.2) Now in general, we must expect that there are a large (e.g., infinite) number of circuits or sequences of elementary gates which will accomplish the above transformation. The complexity of the target state |ψ T is then defined as the minimum number of gates needed to construct a unitary U T satisfying eq. (1.2). We stress that this optimal number will depend on the choices for the reference state |ψ R and for the gate set {g 1 , · · · , g N }, however, one can still obtain interesting physical insights by comparing the complexities for families of target states. Nonetheless, given a particular set of choices, the main challenge is to identify the optimal circuit from amongst the infinite range of possibilities to prepare a certain target state. To overcome this challenge, Nielsen and collaborators [51][52][53] developed a geometric method. Adapting this approach to evaluate the complexity of QFT states [1], one begins with a continuum construction of the unitary transformations acting on the states U (σ) = P exp −i When working with discrete gates as discussed here, we will typically only prepare |ψ T within some tolerance ε, e.g., |ψ T − U T |ψ R 2 ≤ ε. However, with the continous construction of unitaries introduced in eq. (1.3), we are always able to exactly prepare the target states with a finite cost, and so we will not need to introduce a tolerance.
would be an infinitesimal parameter). The control functions Y I (s) then specify which gates (and how many times they) are applied at a particular point s in the circuit. Note that eq. (1.3) specifies not only the full transformation U T in eq. (1.2) but also a trajectory U (σ) through the space of unitaries, or through the space of states with |ψ(σ) = U (σ)|ψ R , where 0 ≤ σ ≤ 1. The circuits of interest are then the trajectories satisfying the boundary conditions U (σ = 0) = 1 , U (σ = 1) = U T , (1.4) i.e., we start from the identity and end with the desired unitary U T producing the desired transformation in eq. (1.2). From this perspective, the Y I (s) can be understood as the tangent vectors to the trajectories with which will play a important role later. Then Nielsen's approach for identifying the optimal circuit is to minimize the cost defined as where F is a local cost function depending on the position U (s) and the tangent vector Y I (s). This question is then similar to the physical problem of identifying a particle trajectory by minimizing the action with Lagrangian F (U (s), Y I (s)). While the precise form of the cost function F is not fixed, there are a number of desirable features for reasonable cost functions [53]: 1) Smoothness, 2) Positivity, 3) Triangle inequality and 4) Positive homogeneity -see also [1,59]. Some simple examples of cost functions that satisfy the above constraints are Given the role of the Y I (s) as control functions, the F 1 measure comes the closest to the original concept of counting the number of gates. The F 2 measure can be recognized as the proper distance in a Riemannian geometry, and this choice reduces the problem of identifying the optimal circuit to finding the shortest geodesic connecting the reference and target states in this geometry. Another class of cost functions introduced in [1] take the form (1.8) These κ cost functions can be thought of as a generalization of the F 1 cost function. The corresponding vacuum complexity compares well with the results from holographic complexity but these cost functions do not satisfy the 'homogeneity' property above, i.e., the cost (1.6) is not invariant under reparametrizations of s. We also note that the κ = 2 cost function will yield exactly the same extremal trajectories or optimal circuits as the F 2 cost function. An interesting suggestion in [41] was to construct a family of new cost functions using the Schatten norm (e.g., see [60][61][62]) where V = Y I (s)O I is the tangent vector defined as an operator which transforms the states (see details in section 3.2). These cost functions satisfy all of the desired properties and further are independent of the particular choice of basis for the O I -a issue for the F 1 measure and the general κ cost functions (for κ = 2) [1]. Before closing this short review, we must mention the group theoretic structure that naturally appears in applying this approach to evaluate the complexity of QFT states. For the problem to be tractable, one only considers a limited basis of operators O I to constructing the unitaries (1.3). A practical restriction is that this basis should then form a closed algebra, and hence in many examples, the O I provide a representation of a Lie algebra g, i.e., For example, in examining the complexity of fermionic Gaussian states, an O(2N ) group structure emerges [41]. In [1], a GL(N, R) algebra appeared in evaluating the complexity of the ground state of a free scalar field, and the latter was extended to an Sp(2N, R) algebra examining the corresponding thermofield double state in [48] -see also [41]. In the following, we will find that an R N GL(N, R) algebra plays a central role in evaluating the complexity of coherent states. The utility of this group theoretic perspective is that the physical details of the basis operators O I can be pushed to the background. Instead, the generators in eq. (1.3) are simply elements of the Lie algebra g, and we can choose the most convenient representation for the particular calculations of interest.
The rest of the paper is organized as follows: In section 2, we study the complexity for coherent Gaussian states in a system of two coupled oscillators. For general states, we must resort to numerical methods to evaluate the complexity in section 2.4, however, we also produce some analytic results for simple cases in which only one mode is excited in section 2.3 or in which the excitations have small amplitudes in section 2.5. This initial analysis is based on the F 2 and κ = 2 cost functions, and so we discuss analogous results for the F 1 cost function in section 3.1 and the p = 1 Schatten norm in section 3.2. We extend our results to a free scalar field theory by a lattice regularization in section 4.
In section 5, we apply the geometric approach based on the Fubini-Study metric [2] to reinvestigate the complexity of our coherent states for two coupled harmonic oscillators. The results for this simple system are also compared with our results in section 2 using Nielsen's approach with the F 2 cost function. We conclude with a discussion of our results and possible future directions in section 6.

Complexity of coupled harmonic oscillators
Our goal is to evaluate the complexity of coherent states in a free scalar field theory, applying the techniques of [1]. However, as a warm-up exercise, we begin here by considering coherent states in the simpler system of two coupled harmonic oscillators. In this section, our focus will be on the F 2 cost function (1.7), and also on the κ = 2 cost function (1.8) which are extremized by the same trajectories. We will turn to consider other cost functions in section 3. Our approach here closely follows that in [1] and we refer the reader there for a more detailed discussion.

Gate set and group structure
Let us consider two coupled harmonic oscillators with the Hamiltonian where in the second line, the two oscillators were decoupled by introducing the normal modes, Given the decoupled Hamiltonian, the ground state wave function is easily written as While the normal modes are completely unentangled here, from the perspective of the physical coordinates {x 1 , x 2 }, the ground state is an entangled state. Ref. [1] developed the techniques needed to evaluate the complexity of this state relative to a factorized Gaussian state as the reference state, where the reference frequency ω R , which characterizes the width of this state, is left as a free parameter. 3 We note that this reference state is unentangled in both the physical and the normal-mode bases. Now we would like to extend the calculations in [1] to evaluate the complexity of coherent states of the form These coherent states are characterized by the expectation values The coherent states in eq. (2.5) are written in terms of the normal modes since this simplifies the calculations below, as shown in [1], and this will be our working basis throughout the rest of this paper. 4 The next step is to identify the set of elementary unitary gates with which we will construct the desired unitary U , which implements (2.8) With the new shift parameters a ± , we need additional gates than those described by the GL(2, R) algebra in [1]. However, the full complement of gates required to construct an arbitrary Gaussian state were discussed in [1] and for the coherent states of the form (2.5), we only need three types of elementary gates: scaling gates : shift gates : Q 0i = e iεx 0 p i , 3 Note that our notation is slightly different from that in [1]. In particular, the latter had mω 0 in place of ω 2 R in eq. (2.4). 4 With the states chosen in eq. (2.5), we are focusing on real wavefunctions with ψ T |x i |ψ T = 0 but ψ T |p i |ψ T = 0. In principle, by considering complex wavefunctions, one could examine more general states which also have ψ T |p i |ψ T = 0, as would naturally arise from the time evolution of the wavefunctions in eq. (2.5). We note that this would require the extending the GL(2, R) algebra appearing below to Sp(4, R), e.g., see [41,48]. We thank Lucas Hackl for a discussion on this point.
where the i, j can be either {1, 2} or {+, −}, but as mentioned above, we will work in the normal mode basis, i.e., i, j ∈ {+, −}. Further we recall that ε is a small (dimensionless) parameter which ensures that these gates only make small changes to the states on which they act. The dimensionful parameter x 0 appearing in the shift gates is another free parameter (a c-number) which characterizes our complexity model. As we discuss below, we might simplify the model by setting x 0 ∼ 1/ω R (or x 0 ∼ δ in the QFT calculations). The action of these gates is illustrated with the following examples: (2.10) Note that our set of elementary gates (2.9) can be summarized by where a ∈ {+, −, 0} and i ∈ {+, −}. We have also introduced the notation O ai to denote the Hermitian generators of these elementary gates. Now following [1] to make further progress, next, we construct a matrix representation of these gates. In general, we are interested in coherent states of the form where again the sums over a, b run over {+, −, 0}, and A is a symmetric 3×3 matrix with A 00 = c. We introduced the term cx 2 0 above to eliminate this c-number contribution from the exponent and hence N is the normalization constant. It will be convenient to keep A 00 in the following calculations, but we stress that its value will be unimportant since the wavefunction (2.12) is independent of this coefficient. 5 Of course, the matrix A completely determines the wave function, and so instead of working with these wavefunctions directly, we focus our attention on the fivedimensional space of A's. Again, the full space of symmetric 3×3 matrices would be six-dimensional but since as explained above, the wavefunctions are independent of A 00 , we have a five-dimensional space of distinct wavefunctions. With this matrix form, we can represent the reference state (2.4) as 13) and the target state (2.5) is represented by where a ± ≡ a ± /x 0 . We emphasize once more that the values of c R and c T are completely undetermined since they do not affect the wavefunctions. By considering the action of the elementary gates (2.10) on the general wavefunctions (2.12), we produce a 3×3 matrix representation of the gates which transforms the A as follows An easy way to verify this result is to consider the action of the matrices Q ai on the vectorx T = (x + , x − , x 0 ) and confirm that the result agrees with the transformation by the original gates (2.9), e.g., we can comparẽ (2.18) The convenience of this representation is that we can define a simple inner product of these matrix generators (2.16), where I, J ∈ {++, −−, −+, +−, 0+, 0−}. We will use this inner product in a moment in constructing a metric on the six-dimensional space of unitary transformations generated by our elementary gates (2.9). Now we recall from [1] that the four generators M ij for the scaling and entangling gates (appearing in the first line of eq. (2.18)) form a GL(2, R) algebra. More generally if we consider the action of a string of the elementary gates on x + and x − , we find that they are transformed as . That is, our gates produce affine transformations of the coordinates. Hence the full group generated the six gates Q ai has a structure similar to that of the Poincaré group. The GL(2, R) of the scaling and entangling gatess plays the role of the Lorentz group O(1, 3) and the R 2 of translations generated by the Q 0i is analogous to the translations in Minkowski space. Hence, the structure of our algebra here 6 is the semidirect product of R 2 by general linear transformations GL(2, R), i.e., (2.20)

Six-dimensional geometry and its geodesics
The group structure (2.20) is manifest by any transformation U generated by the M ai taking the form where u T = (u + , u − ) ∈ R 2 and U 2 ∈ GL(2, R). It will be convenient to parametrize the latter with the following polar decomposition where R denotes a rotation matrix and S is a 'squeezing' matrix. This construction then introduces the coordinates y T = (y, ρ, x, z, u + , u − ) on the group of affine transformations (2.20). There is also a surjective, but not injective, map that associates a wavefunction of the form (2.12) to every group element, given by with the reference state given in eq. (2.4). The corresponding transformation using the matrix representation (2.21) becomes 24) where A R is given in eq. (2.13) and In fact, one can see that Λ = U 2 ·u where we have assembled the vector Λ T = (Λ + , Λ − ). This observation is useful because it allows us to see that the final matrix A(y) is independent of z in the following sense: First, it is obvious that the upper-left 2×2 block in eq. (2.24) is invariant under arbitrary shifts z → z = z + δz, i.e., this block is completely independent of z. Now given the form of U 2 in eq. (2.22), it is also evident that Λ is invariant as long as we accompany the shift of z with a rotation u → u = R(−δz) · u. Finally, this rotation also leaves invariant the component [A(y)] 00 , as can be seen by writing this term as u T · u. This result reflects the fact that the map from the space of unitary transformations to Gaussian states is surjective but not injective, i.e., the space of unitaries which we are considering is six-dimensional while our space of Gaussian states is only five-dimensional. 7 7 Further, the fact that this mismatch appears as A(y) being independent of z is a reflection of the rotation invariance of the reference state (2.13). This symmetry can be made more explicit by reparametrizing the group elements as We then find where x T = (x + , x − ) and A 2 = A 2 (y, ρ, x) is the 2×2 upper-left matrix in eq. (2.24). The wavefunction is then clearly independent of z. We chose not to use this parametrization because the metric in these coordinates is much more complicated. The metric in these coordinates is Further, for the parametrization of the group elements in R 2 GL(2, R) in eq. (2.21), we can define a metric on the space of unitary transformations as 8 , ds 2 = dv 2 + + dv 2 − + 2v + dv + (dy + cos 2x dρ − sin 2x sinh 2ρ dz) + 2v + dv − (−dx + sin 2x dρ + (cosh 2ρ + cos 2x sinh 2ρ) dz) + 2v − dv − (dy − cos 2x dρ + sin 2x sinh 2ρ dz) + 2v − dv + (dx + sin 2x dρ − (cosh 2ρ − cos 2x sinh 2ρ) dz) . 8 More generally, one could replace δ IJ → G IJ in constructing this metric. However, the present choice assigns the same cost to all of elementary gates and further it corresponds to the F 2 cost function introduced in eq. (1.7). Following a construction analogous to that in [63] (see also [41,48]), we can also construct the metric by defining In this case, the metric takes the form Of course, this metric agrees with eq. (2.28) when we choose c R = ω 2 R , i.e., with A R ∝ 1. Recall that up to this point, c R was a spurious parameter but with the above construction, it plays an essential role in defining the geometry. In particular, if we were to adopt this approach, we would have to restrict our attention to c R ≥ 0 to produce a sensible geometry.
An intuitive cost function in this context is the F 2 norm (1.7), which becomes D 2 (U ) = 1 0 ds g abẋ aẋb , (2.29) i.e., this simply corresponds to evaluating the geodesic distance in geometry defined by eq. (2.28). However, as was argued in [1] (see also [2]), this cost function does not reproduce the expected UV divergences found in holographic complexity [21]. However, this situation can be remedied by using the κ = 2 cost function (1.8), which corresponds to D κ=2 (U ) = 1 0 ds g abẋ aẋb . (2.30) Of course, from a physicist's perspective, this can be seen as the action of a test particle moving in the same geometry and so it yields the same extremal trajectories. We will also consider two alternative cost functions in section 3, the F 1 and Schatten measures, but in the following we will focus on finding the circuits that minimize the distance (1.6) using the cost functions (2.29) and (2.30). Now the complexity is given by the shortest unitary connecting the reference and target state, i.e., C (A T ) = min U D (U ) with and With the cost functions in eqs. (2.29) and (2.30), this corresponds to finding a geodesic from the origin in the geometry (2.28) to the point corresponding to the desired transformation U (σ = 1). However, as we described for the transformation in eq. (2.24), the target state is independent of one of the coordinates in U or alternatively, the reference state is invariant under a family of transformations (known as the stabilizer group, e.g., see [41,48]). Hence for any target state A T , there exists a one-parameter family of transformations satisfying the boundary conditions in eq. (2.31). Thus, there is a one-parameter family of geodesics connecting the reference state to the target state and the complexity will be determined by the length of the shortest geodesic in this family. For simplicity, we describe the determination of the geodesics in terms of extremizing eq. (2.30), which takes the form of a particle action with Lagrangian We solve the resulting 'equations of motion' analytically for simpler target states in section 2.3 and provide numerical solutions for general target states of the form (2.14) in section 2.4.

Solving for simple geodesics
While we were not able to find analytic solutions for the geodesics to arbitrary target states (2.14), for some simpler set of target states, the optimal path between the reference and target states remains in a H 2 × R slice of the full geometry (2.28). We begin by describing these simple geodesics which have an analytic solution. In section 2.4, we confirm numerically that these are indeed globally the shortest geodesics for the particular target states of interest. First, we can use the freedom to reparametrize s in the cost function (2.29) to fix k = g abẋ aẋb (2.33) where k is a (positive) constant. This means that when evaluated for the optimal trajectory, the complexity with this cost function is simply given by C 2 = k. Similarly, with the κ = 2 cost function (2.30), the complexity is given by C κ=2 = k 2 . Now to identify simple geodesics, we begin by looking at the equations of motions for x(s) and z(s): Now ifu ± = 0 (e.g., if we were simply preparing the vacuum state as in [1]), these equations would be solved by x(s) =x and z(s) =z, i.e., setting both of these coordinates to be constant along the trajectory. These are the geodesics that we will focus on below.
To pick appropriate values forx andz, we must examine the boundary conditions.
The initial boundary condition U (s = 0) = 1 and comparing to eqs. (2.21) and (2.22) then gives where the subscript notation indicates y a0 = y a (s = 0). 9 The first restriction implies that we must choosex =z for our simple geodesics. Similarly for the final boundary conditions, comparing (2.14) and (2.24), we see that sin(2x 1 ) = 0 is required to produce a target state which is unentangled in normal mode basis. Hence this end point condition requires sin(2x) = 0, i.e.,x = nπ/2. Combining these conditions for the final state from eq. (2.24) (with c R = 0) gives at s = 1, A(s = 1) = U (s = 1) A R U T (s = 1) (2.37) 9 We will use a similar notation for the final boundary conditions at s = 1.
With these choices, the z equation (2.35) reduces to Together with the initial boundary conditions (2.36), this is satisfied for ρ = 0 or u + = 0 or u − = 0. The first of these possibilities is inconsistent with the final boundary condition if the normal mode frequencies given in eq. (2.2) are not equal, i.e., Ω = 0.
Hence we must choose one of the latter two possibilities. That is, the consistency of our simple geodesics (with constant x and z) demands that we only shift one of the normal modes to produce either a nonvanishing expectation value x + or x − , but not both! 10 We continue our discussion here with the choice u − (s) = 0, i.e., we consider states with x + = 0 and x − = 0. To ensure that the choicex =z, sin(2z) = 0 and u − = 0 is still a geodesic of the full geometry (2.28), we verify that the equation of motion for u − is satisfied, i.e., which is indeed the case. The induced geometry on the slice given by these choices becomes ds 2 = dy 2 + + e −2y + du 2 where we have introduced y ± = y ± ρ. Our analysis has guaranteed that finding geodesics (y + (s), y − (s), u + (s)) in the induced metric (2.40) will give us geodesics (y + (s), y − (s), x = nπ, z = nπ, u + = 0 = u − ) in the full six-dimensional geometry described by eq. (2.28). It is straightforward to see that the three-dimensional geometry (2.40) is the direct product of two-dimensional hyperbolic space (covered by y + and u + ) with the real line (covered by y − ). Since two components of this geometry are maximally symmetric, it would be straightforward to evaluate the distance between any two points in this geometry using standard formulae. However, it is useful for us to have explicit expressions describing the geodesics and so we proceed by evaluating the equations of motion in the H 2 × R geometry. Of course, from eq. (2.36), the initial conditions for the geodesics are simply: y +0 = 0 = y −0 = u +0 . To determine the final boundary conditions, we begin with eq. (2.37) for the final state, which simplifies to Requiring that this state matches the target state A T with a − = 0 in eq. (2.14) yields the following boundary conditions at s = 1: where for convenience, we are using the following dimensionless ratios: 11 Now to find the path which these geodesics follow, we extremize the cost function (either eq. (2.29) or (2.30)) subject to the restriction that the motion is only in the three-dimensional subspace described by eq. (2.40). Each of the three equations of motion can be integrated to yield the following first order equationṡ where the three integration constants correspond to the velocities at s = 0, e.g., y + | s=0 = A. These equations can be integrated and after imposing the initial conditions, the solution becomes . Furthermore, the final conditions (2.42) fixes the integration constants as

46)
11 Of course, a ± are the same quantities which already appeared in eq. (2.14).
where the sign of B is chosen to match the sign of u +1 (i.e., the opposite sign as a + ). 12 Further, it follows that For these simple geodesics, the constraint (2.33) reduces to As will be shown in section 2.4, these geodesics are indeed the shortest ones connecting the reference state (2.13) to the target states (2.14) with a − = 0 in the full geometry. Therefore eq. (2.48) yields the complexity of the coherent states with the F 2 and κ = 2 cost functions, i.e., C 2 = k sim and C κ=2 = k 2 sim . As a check, one can easily verify that in the limit a + → 0, this result (2.48) yields the complexity of the ground state found in [1], i.e., As expected, the difference in complexity between the coherent states and the ground state comes from the normal mode which has been translated (x + in this case). It is interesting to examine the difference in various limiting cases. 13 That is, let us consider the asymptotic behavior of Expanding for small |a + |, eq. (2.50) yields Here, we assume the definition of 'arccosh' is such that it always yields a positive result. 13 In the following, we focus on the complexity for the κ = 2 cost function rather that the F 2 measure. There are two motivations for doing so: First, the κ = 2 complexity reproduces the expected UV divergences of holographic complexity as was found in [1]. Second, as we will see in section 4.2, the change in F 2 complexity ∆C 2 vanishes when generalizing our results to free scalar field theory. In contrast, the change in κ = 2 complexity ∆C κ=2 remains finite when generalizing to field theory.

Numerical results in full geometry
In this section, we describe our results for numerical solutions of the geodesic equations.
Our approach was to derive the second order differential equations from the variation of eq. (2.32) and then use the pseudo-spectral method where Chebyshev polynomials were used in the s direction. Combining the Dirichlet boundary conditions at s = 0 and s = 1, the solutions can be uniquely determined. One subtlety is that in the initial conditions (2.36), the value of x 0 = z 0 is not fixed. However, with our method, this parameter is easily fixed by scanning through a range of values for x 0 and demanding that the solution is smooth in the vicinity of s = 0. Our first application was to verify our numerical results by comparing them with the analytic solutions for the simple geodesics found in the previous section. An example is shown in figure 2. As expected, if u + (u − ) ends at zero, it remains zero along the entire trajectory, and the scale coordinate y − (y + ) follows a straight line. The other two coordinates follow curved paths, as expected from eqs. (2.45) and (2.46). In every case, we found excellent agreement between the numerical and the analytic solutions.
Next we considered the family of geodesics connecting the reference state to a specific target state with a −1 = 0 (or a +1 = 0), as shown in figure 3. Recall that as described below eqs. (2.24) and (2.25), the final state was independent of z 1 (as long as the final values u ±1 were rotated appropriately). In the figure, we see that the shortest geodesic is that for which z 1 = 0. That is, for all the examples that we examined, our numerics confirm that the optimal geodesics correspond to the simple geodesics derived in the previous section. Hence these numerical studies provide strong evidence for the claim that the simple geodesics are indeed the shortest ones for the target states in which only one of normal modes is shifted. Figure 4 shows an example of the optimal geodesic to a target state with both a + and a − non-zero. In this situation, we do not have an analytic solution and we can see in our numerical solution that these geodesics do not take a simple form, e.g., none of the coordinates follow a straight path. Similar to the previous discussion, to determine the optimal geodesic, we vary z(s = 1) while keeping the final state (2.24) fixed, evaluate the lengths of the corresponding geodesics and then choose the trajectory with the minimal length. Let us note in passing that generally these optimal geodesics pass through regions where x(s) and ρ(s) are nonvanishing. Therefore even though both the initial (2.13) and final (2.14) states are unentangled, the intermediate states (all) along the optimal trajectory are entangled when preparing a target state with both a + and a − are nonvanishing -see further discussion in section 6 Recall that for the simple geodesics with only a single excitation, we found x(s) = 0 = z(s). In contrast, when both normal modes are excited in the target state, the optimal geodesic has nonvanishing profiles for both x(s) and z(s), as illustrated in figure 4. While at the final point, x 1 = 0 in order to ensure that the target state is unentangled, as can be inferred from eq. (2.24), in general we have z 1 = 0 for the optimal geodesic. For comparison purposes, we can also consider the geodesic with z 1 = 0, which we will denote as the "simpler" geodesic, which is also shown in figure  4. There we can see that the biggest difference between these two geodesics is in the profiles of x(s) and z(s). In fact, the profiles for y ± (s) are indistinguishable in the figure. It is also interesting to compare the length of these geodesics, which we do in  . Example of geodesics preparing a target state with both a ± nonvanishing. We compare the optimal geodesic and the "simpler" geodesic with z 1 = 0. In this example, both geodesics have the final boundary conditions: y +1 = 0.1, y −1 = 1.1, Λ + = 1.105, Λ − = 6.008. For the optimal geodesic, we also have x 0 = z 0 = −0.0976π, x 1 = 0, z 1 = −0.0167π, while for the "simpler" one, x 0 = z 0 = −0.0749π, x 1 = 0 = z 1 . Note that y ± essentially coincide in both geodesics, as shown in the far left panel.
figure 5. There we introduce the new parameter ∆k = k sim − k opt , 14 In the figure, we show the results for ∆k for geodesics where the boundary conditions are fixed as in figure 4 except that Λ − varies from 0 to 12.017. With Λ − = 0, ∆k = 0 because the two geodesics coincide with the simple geodesics found analytically in the previous section. However, we see in the figure that as a − increases, the difference in lengths increases 14 Of course, ∆k > 0 because the optimal geodesic is the shortest geodesic. Figure 5. A comparison of the lengths of the optimal geodesic and the "simpler" geodesic with z 1 = 0. ∆k = k sim −k opt and k opt are shown as functions of Λ − (s = 1) = − mω + x 0 ω 2 R a − . This example is characterized by the boundary conditions y +1 = 0.1, y −1 = 1.1 and Λ +1 = 1.105, while Λ − (s = 1) varies from 0 to 12.017. We note that while ∆k grows as |a − | increases, it represents at most an increase of 0.2% over k opt for the geodesics shown here. monotonically.

Small excitations
We began in section 2.3 by considering simple geodesics for states where only one normal mode is excited. Then in section 2.4, we applied numerical techniques to examine the geodesics for target states where both normal modes are excited. In particular, we noted that the resulting geodesics are driven away from the space of states with no entanglement between the two normal modes. While we cannot find the geodesics for these general target states analytically, we can at least find the leading order contributions to the length of the geodesics for small shifts, i.e., when both a ± are nonvanishing but |a ± | 1. To examine this situation, we consider small perturbations from the optimal geodesics connecting the reference state to ground state state. It was already shown in [1] that the optimal circuit connecting the reference state (2.13) to the ground state, is the 'straight line' circuit: In terms of the six-dimensional geometry, the corresponding geodesic is given by where y +1 = 1 2 log w + and y −1 = 1 2 log w − , with y ± = y ± ρ as before. The κ = 2 complexity (2.30) is then given by the expression in eq. (2.49).
Now we want to evaluate the leading order change to the above circuit depth (2.49) evaluated with eq. (2.30) when we introduce small shifts for both normal modes, i.e., a + , a − ∼ O(ε). In particular, the final boundary conditions are then modified for u ± but it will be true that u ± (s) ,u ± (s) ∼ O(ε) all along the new geodesic. This follows because the second line in (2.32) is positive definite, so havingu ± = O(1) would increase the length by an O(1) factor. The x and z equations of motion take the form . Now if we expand the cost function, i.e., determine the leading corrections to eq. (2.32), we find Effectively, the modified geodesics move on a four-dimensional submanifold of the full geometry (2.28) which takes the form H 2 × H 2 . Hence to leading order, the complexity becomes where ∆ + is the expression in eq. (2.47) and ∆ − is the same expression after substituting w + → w − and a + → a − . The leading order change in the complexity then becomes where we are dropping terms of the form a 4 + , a 4 − and a 2 + a 2 − . The key result here is that the leading order corrections to complexity factorize into contributions from the individual normal modes, i.e., there are no second order contributions involving the cross-term a + a − .
We would like to go further and so that in fact, these geodesics on the effective H 2 × H 2 geometry are optimal, i.e., that we are correctly evaluating the leading corrections to the complexity in eq. (2.59). We can argue for a proof by contradiction of this result as follows: Imagine that we find a geodesic where the deviations of x(s) and z(s) from the straight-line geodesic (2.55) are the same order as a ± , i.e., x, z ∼ O(ε). Examining the cost function, we see that the motion in x and z will introduce a strictly positive contribution of order ε 2 , from the terms in the first line of eq. (2.32). Similarly the Figure 6. Comparison between the quadratic approximation for small excitations, i.e., the second line of eq. (2.59), and the complexity found numerically for various target states. In the figures, we choose w + = 1.221 and w − = 9.025 and a − = γa + . We let a + range from 0 to 0.01 and then compare the two results with γ = 1, 1.5 and 2. The three figures show excellent agreement between the quadratic approximation and the true complexity for a ± 1.
second and third terms in the second line will make contributions of O(ε 4 ) and O(ε 3 ), respectively. There is no definite sign of these contributions but being higher order, they will never be able to make up for the O(ε 2 ) increase generated by moving in the x and z directions. One might consider even stronger deviations, e.g., x, z ∼ O(1), but then theẋ,ż terms in eq. (2.32) will only increase the length of the geodesic by order one whileu ± terms will still only contribute at order ε 2 . We can also use the numerical approach developed in the previous section to find evidence that eq. (2.59) correctly provides the leading corrections to the complexity. In particular, we looked at families of states where a − = γa + with γ being some fixed numerical constant. As a + increased from zero, we found that the numerical results matched the approximation provided in the second line of eq. (2.59) when a ± 1 in all of the cases that we examined. Figure 6 provides an example of our numerical analysis.

Complexity with alternate cost functions
In [1], the UV divergences in complexity of the ground state of the free scalar field were compared to those in holographic complexity (see also [2]). In particular, it was found that the F 1 cost function in eq. (1.7) gave the most promising comparison. In particular, the leading divergence for the F 1 cost function took the form V /δ d−1 log( /δ) where is some undetermined length scale. This is precisely the same form as the leading divergence in holographic complexity evaluated with the CA proposal [21]. However, an apparent shortcoming of the F 1 cost function is that the complexity depends on the basis used for the gates, e.g., the results will change upon rotating between the physical basis and the normal mode basis. However, in [41], it was suggested that we would recover the same essential results of the F 1 measure using the Schatten norm (with p = 1 -e.g., see [60,61] and further details below). The advantage of the Schatten norm is that the results are basis independent. Hence in the following, we will re-examine the complexity of the coherent states (2.5) for the system of two coupled harmonic oscillators introduced in section 2.1 using these two alternatives for the cost function.

F 1 cost function
First, we turn to the task of studying the F 1 cost function introduced in eq. (1.7) that is, we want to study the circuits U (s) that optimize the cost function However, we re-iterate that this measure is not invariant under changes of the basis, and therefore the results depend on the choice of basis I which we choose in its definition. As we saw in the previous section, the normal modes provide a natural basis to work in for the circuit optimization problem and so in the following, we simply define the F 1 metric in this basis. Hence for the problem of two coupled harmonic oscillators, which we focus on in the following, the index I in eq. (3.1) runs over {++, −−, −+, +−, 0+, 0−}.
Using the results of section 2, we find the components Y I appearing in the cost function (1.6) to be We will not attempt to find the extremal trajectories in complete generality here. Rather we will focus on the analog of the simple geodesics found in section 2.3, which prepare coherent states where only one of a ± is nonvanishing. We will also consider the case of small excitations, i.e., |a ± | 1, to parallel the analysis in section 2.5. To begin let us consider constraining the trajectories with x(s) = 0 = z(s), in which case the F 1 cost function (3.2) takes a simple form where once again we used y ± = y ± ρ. A key feature here is that the motions for the + and − coordinates have decoupled, which is reminiscent of the trajectories studied in sections 2.3 and 2.5. We will proceed with examining the geodesics in the x = 0 = z subspace which extremize eq. (3.4) in a moment. However, first imagine that we have these solutions and then we wish to show that they also extremize the full cost function (3.2) by considering perturbations away from this subspace. Let us construct a perturbative expansion with x(s), z(s) ∼ O(ε). Then keeping the leading perturbations in eq. (3.2) yields First, let us consider the case where the excitations are small, i.e., |a ± | = O(ε). In this scenario, we expect u ± = O(ε), and therefore the x, z dependent terms in the first two lines are O(ε 2 ). Therefore, we drop the latter and the only O(ε) terms involving these variables are in the third line and the action is minimized when they vanish, which yieldsẋ These expressions in turn are solved by x tanh(ρ) = constant. The initial condition ρ 0 = 0 implies that the constant is zero, which in turn leads to x(s) = 0 = z(s) as necessary conditions to minimize the action. The remainder of the cost function then takes the simple form in eq. (3.4), in which the ± coordinates decouple, but recall that here we assumed that both |a ± | 1. Instead let us assume that we have a coherent state where a + is large but that a − = O(ε) or zero. In this situation, we again expect u − = O(ε), which means that the x, z dependent terms in the first line are O(ε 2 ) and can be ignored again. The O(ε) terms involving x and z are the second term in the second line and the two terms in the third line of eq. (3.5). The action will be minimized if we can find a solution where they vanish. The terms in the third line vanish with x = 0 = z, with the analysis following eq. (3.6). Hence the action again reduces to the form given in eq. (3.4), although we must keep in mind that the term involvingu − is O(ε). 15 We now turn to the problem of finding the geodesics in the simple "geometry" appearing above in eq. (3.4). That is, we consider (3.7) With the "flat measure" D = ds (|ẏ| + |u|) the minimal trajectories are simply those which traverse between the initial and final endpoints without backtracking in u or y. However, with the addition of the scaling factor e −y in theu term in eq. (3.7), there is a balance between backtracking in y and attempting to reduce the scaling factor by going to a larger y before turning back to the final value. This leads to two possible classes of paths that can minimize the distance (3.7), illustrated in figure 7. We call these the L and J paths. We will assume y 1 > y 0 and u 1 > u 0 to simplify our discussion, but the other cases are very similar to these. The length of the L-shaped path that is a straight line from (y 0 , u 0 ) to (y 1 , u 0 ) and then a straight line from (y 1 , u 0 ) to (y 1 , u 1 ) is terms in the third line vanish, and the two modes would decouple in determining the optimal path. That is, the optimal trajectories would effectively be determined by the product geometry H 2 × H 2 even when |a ± | 1.
where ∆y = y 1 −y 0 and ∆u = u 1 −u 0 . For the J-shaped paths, there are three straightline sections: (y 0 , u 0 ) to (y 2 , u 0 ), (y 2 , u 0 ) to (y 2 , u 1 ) and then (y 2 , u 1 ) to (y 1 , u 1 ). The length of this path is D J (δy) = ∆y + 2δy + e −(y 1 +δy) ∆u , (3.9) where δy = y 2 − y 1 > 0. This cost is minimized with δy min = log ∆u 2 − y 1 . The optimal J path therefore goes up to y 2,min = log ∆u 2 , then over to u 1 and then back to y 1 , and has length D J (δy min ) = ∆y + 2 log ∆u 2 − 2 y 1 + 2 . (3.10) The L paths are shorter for e −y 1 ∆u 2, while the J paths are shorter for e −y 1 ∆u 2. Finally, we can express the boundary conditions in terms of parameters of the target state (and the reference state). For the reference states (2.4) and the target states (2.5), we find 16 Substituting these expressions into eqs. (3.8) and (3.10), we find Using eq. (3.11), we also find e −y 1 ∆u = a and so we see that D L is smaller for a 2 while D J is smaller for a 2.
Note that since we assumed that y 1 > y 0 above, we are implicitly considering w > 1. Carefully going through the argument above for the case y 0 > y 1 , we find that the following lengths for the two types of paths for w < 1. Here, the J paths are shorter for √ wa 2 and otherwise the L paths are shorter. Note that D J has precisely the same form in eqs. (3.12) and (3.13). So in general, we can write the cost of these two families of extremal paths as where the J geodesics are defined (and shorter) for min(1, √ w) |a| 2. Hence in general, the F 1 cost function (3.7) yields the following complexities To conclude here, let us note that if we are considering small excitations of the ground state, i.e., |a| 1, then the optimal circuit will be described by a L-shaped geodesic. Further, in this case, the change in the complexity will be linear in |a|, i.e., ∆C 1 = C 1 − C 1,vac ∝ |a|. The latter behaviour contrasts with our previous results for the κ = 2 cost function, where we found ∆C κ=2 ∝ a 2 in eq. (2.51). Further, for large excitations with |a| 1, the optimal circuit is descriped by a J-shaped geodesic and we find ∆C 1 log a 2 . Recall that we found in eq. (2.52) that ∆C κ=2 (log a 2 ) 2 and so there is again a difference in the power of the leading contribution.
For target states where both modes are excited but with a ± 1, the resulting "geometry" is a product space of two copies of the above geometry. The geodesics will therefore correspond to the L-shaped geodesics and following eq. (3.15), the total complexity is then where D L (w, a) is defined in (3.14). With these small excitations, the increase in complexity above the vacuum complexity is given by which is linear in |a ± |. We can also consider the complexity of states where one excitation is large, e.g., |a + | 1 (but the other is still small). This contribution dominates and C 1,tot D J (w + , a + ). Hence the change in the complexity becomes 17 Let us note that unlike the case of the κ = 2 complexity of section 2 and of the Schatten complexity of the next section, ∆C 1 is independent of the excited frequency w, unless w < 1 and √ w|a| 2, in which case it is proportional to √ w.

Schatten cost function
A suggestion put forward in [41] is that we might use the p = 1 Schatten norm (e.g., see [60][61][62]) rather than the F 1 cost function. The observation was that with this new cost function that we would recover the same leading divergence as with the F 1 measure for the complexity of the vacuum, 18 however, the results are now basis independent 17 As described in footnote 15, if both of excitation parameters a ± are large, we expect that two modes still decouple in the optimal preparation of most such states. In this situation, the change in complexity would scale as ∆C 1 log a 2 + + log a 2 − . 18 In fact, the vacuum complexity was identical for both measures, but we will see below that this does not carry over for the coherent states studied here.
when described in terms of the Schatten norm. This norm actually provides a family of measures based on computing the singular value decomposition of the desired transformation. Given a transformation A, this norm takes the form Note that with p = 2, this reduces to the standard Frobenius-Hilbert-Schmidt norm, i.e., we recover the F 2 measure which we were studying in the previous section. As with the F 2 cost function, the results are basis independent for the Schatten measure for any value of p. Another property worth noting is that the Schatten p-norms are non-increasing in p, which means we have In the present case, the transformation of interest is the velocity tangent to the path of unitaries, namely In particular then, V 1 = k s k . Given eq. (3.20), we can explicitly write out the self-adjoint operator where implicitly we are summing over a ∈ {+, −, 0} in each component. Hence we can immediately see that in the special case of interest, the third singular value is zero and we simply need to find the eigenvalues of the upper 2×2 block. The latter is a simple exercise, which yields In general, we can write V = R 1 DR 2 where R 1 and R 2 are two independent rotation matrices while D = diag(s 1 , s 2 , s 3 ) with s k ≥ 0. The s k are the singular values of V , which only agree with the eigenvalues of V in special cases. For example, the two agree when V is symmetric and non-negative. We note that the singular values and eigenvalues do not agree for the case of interest in eq. (3.20).
Substituting these expressions into eq. (3.21) for p = 2, we recover 24) in agreement with the F 2 cost function in eq. (1.7), as expected.
Turning instead to the Schatten cost function with p = 1, we find It is useful to consider some simple examples, i.e., some simple trajectories. First imagine that we are only scaling the two normal modes, as we would in preparing the ground state. Then from eq. (3.3), we have Y ++ =ẏ +ρ =ẏ + and Y −− =ẏ −ρ =ẏ − , and eq. (3.25) reduces to Of course, this expression has the same form as the F 1 cost function for these trajectories, and so both measures would yield the same complexities in situations where these simple scaling circuits were the optimal ones. But now let us consider trajectories where there is also a displacement for, say, the + mode, i.e., whereu + = 0. Then another component of the tangent vector is also nonvanishing, namely, Y 0+ = e −y +u + . The cost function (3.25) then becomes Hence the 'Schatten' cost of this simple trajectory is already different from the F 1 cost. 20 Interestingly, because the motions associated with the ± modes are decoupled in the above cost function (3.27), we can easily find the optimal trajectory is a geodesic in the product geometry H 2 × R. The optimal trajectory which extremizes eq. (3.27) is precisely the 'simple geodesic' discussed in section 2.3. However, we have restricted the motion of the trajectories in constructing the expression in eq. (3.27) and so next we would like to show that our 'simple geodesics' also extremize the full Schatten norm (3.25). Towards this goal, we consider a new Lagrangian (or cost function) which is the square of Schatten cost function, and if we find trajectories which extremize L 0 (and yield V 1 = 0), then they will also extremize V 1 , the desired cost function. Now we have divided the result in 20 Further, we can anticipate that for small displacements of u + , i.e., small excitations of a + , the total cost will have a contribution proportional to ∆u 2 + ∼ a 2 + . eq. (3.28) into the sum of two pieces: L 0 = γ 1 + γ 2 , which corresponds to the κ = 2 cost function, 21 and L 1 = γ 1 γ 2 . Now we wish to consider the simple geodesics given by eqs. (2.45) and (2.46), as well as Now the analysis in section 2.3 showed that these trajectories extremized the κ = 2 cost function (2.32). Hence we already know that the simple geodesics will extremize the first part of eq. (3.28), and we need only examine the variations of the L 1 term. These equations of motion are generally very complicated but they simplify enormously when we substitute eq. (3.29). The simplied variations are (3.30) However, one can easily show that the three remaining variations will vanish upon substituting the corresponding equations of motion derived from the κ = 2 cost function: Therefore we arrive at the desired conclusion that the simple geodesics in the H 2 × R geometry also describe the optimal circuits for the (full) Schatten p = 1 cost function (3.25).
Hence the coherent states in which a single normal mode is excited are prepared in precisely the same way as in section 2.3. Recall that the boundary conditions for these trajectories are given by eq. (2.42). Further, using the subsequent analysis in section 2.3, it is then straightforward to show that the complexity measured by the Schatten cost function is then given by where C and ∆ are given in eqs. (2.46) and (2.47), respectively. The increase in the complexity above that of the vacuum state is given by (3.33) 21 As we saw above in evaluating the Schatten norm with p = 2, Expanding for small |a + |, eq. (2.50) yields while for large |a + |, we find In analogy to section 2.4, one might attempt to study numerically the full equations of motion resulting from eq. (3.28) to investigate the complexity of states where both of the normal modes are excited. However, we do not pursue this direction here. Instead, we turn to an analysis for such states in the regime where the excitations are small, i.e., a ± 1, in analogy to section 2.5. We will assume that in the excited state that a ± ∼ O(ε) where ε is a small expansion parameter in the following perturbative construction. Assuming the variation of the geodesics is smooth and starting from the geodesic with u ± = 0 = x = z, i.e., a ± = 0, the perturbed geodesic line for to leading order in our ε expansion. Therefore we define the perturbative solution with u ± (s) = εu (1) ± (s) + O(ε 3 ), and similarly for x(s) and z(s). Substituting these expansions into the expressions in the cost function (3.28), we find (3.37) So let us consider the solutions extremizing L 0 first: It is consistent to solve with x (1) (s) = 0 = z (1) (s). 22 With this choice, the ± modes decouple at this order and the solutions correspond to geodesics on H 2 × H 2 . Further note that for either mode on these geodesics,ẏ 2 + e −2y ε 2 (u (1) ) 2 = ∆ is a constant of the motion. That is, this is like a conserved energy, which corresponds to ∆ 2 in eq. (2.47) for the simple geodesics. Now we move to consider variations of L 1 . Again it is straightforward to show that x (1) (s) = 0 = z (1) (s) is a consistent solution. To facilitate the discussion, we can then write where we should only really pay attention to the terms to O(ε 2 ). But given this form, the variations with respect to y ± and u ± are all proportional either to equations of motion from L 0 or to ∂ s (ẏ 2 + e −2y ε 2u2 ), both of which vanish for the perturbative solutions of the equations of motion from L 0 . Therefore to leading order, the two modes decouple and we can just consider geodesics on H 2 × H 2 . From eq. (3.34), the resulting change in complexity is then just

Complexity of coherent states in QFT
In the previous section, we examined the complexity of coherent states in a system of two coupled harmonic oscillators. In this section, we extend the results to the quantum field theory describing a free scalar. In particular, we consider a free scalar theory in d spacetime dimensions with the Hamiltonian Following [1], we regulate the theory by putting it on a lattice with lattice spacing δ, in which case the regulated theory becomes a family of coupled harmonic oscillators. The lattice Hamiltonian can be written as 23 where in the second line, we have defined X( n) = δ d/2 φ( n), P ( n) = δ (d−2)/2 π( n), m = 1/δ, ω = µ and Ω = 1/δ. Hence as noted above, the lattice Hamiltonian describes a system of the coupled harmonic oscillators on an (d-1)-dimensional lattice. For 23 We approximate the spatial derivatives as ∂ i φ(x) 1 δ (φ(x) − φ(x − δx i )) and we designate the lattice sites with n = n ix i , wherex i are unit normals along the spatial axes.
simplicity in the following, let us consider the case of d = 2. That is, we will consider N oscillators on a one-dimensional lattice with the Hamiltonian and periodic boundary conditionsx a+N =x a . The Hamiltonian is straightforwardly rewritten in terms of normal modes, where the normal modes and the corresponding frequencies are given by with k ∈ {1, ..., N }. The conjugates are x † k = x −k = x N −k and similarly for p k . In the normal mode basis, the ground state wave function becomes The complexity associated with preparing the ground state from an unentangled reference state 24 was analyzed in [1] and, e.g., with the κ = 2 cost function (1.8), the complexity is given by where in the latter expression, we substituted the notation introduced in eq. (2.43). We now consider the complexity of coherent states in the (regulated) scalar field theory of the form Of course, this question is a simple extension to N modes of that examined in the previous section for two coupled harmonic oscillators. Hence the construction of the circuits preparing ψ T (x k ) from the reference state ψ R (x k ) also only requires a straightforward extension of the previous discussion. For example, to define the gates, we need only extend the range of the indices in eq. (2.11), i.e., i ∈ {1, . . . , N } and a ∈ {0, . . . , N }. 25 With this set of gates, the group structure in eq. (2.20) is generalized to R N GL(N, R), and eq. (2.21) is replaced by a representation of (N + 1)×(N + 1) matrices taking the form where u T = (u 1 , . . . , u N ) ∈ R N and U N ∈ GL(N, R). In principle, we can then construct a metric on the corresponding N (N + 1)-dimensional space of unitaries, analogous to eq. (2.28), and the geodesics would be solutions extremizing the following particle action, analogous to eq. (2.32), However, parametrizing the transformations in eq. (4.10) (and in particular, the GL(N, R) transformation U N ) is far more complicated. In any event, given our experience in the previous section, we do not expect that we will be able to find analytical solutions for geodesics preparing general states of the form in eq. (4.9). Instead then, let us focus on the special case where only a single mode in the lattice model is shifted. In particular, we focus on target states of the form (4.12) where only the i'th mode is excited by shifting to x i = a i . Motivated by the results in last section and also in [1], we are led to conjecture that the optimal circuit preparing this state corresponds to a geodesic in the geometry H 2 ×R N −1 . That is, for the optimal transformation preparing the above state (4.12), eq. (4.10) reduces to The y j (s) with j = i would simply grow linearly with s, while y i (s) and u i (s) satisfy the geodesic equations on the hyperbolic space H 2 . This suggestion generalizes the geodesics on H 2 × R found for two coupled harmonic oscillators in the previous section, and if we set a i = 0 (and hence u i (s) = 0), the motion would be restricted to the R N parametrized by y i , which was dubbed the normal mode subspace in [1]. In fact, we can prove that eq. (4.13) do indeed yield a family of simple geodesics in the full N (N + 1)-dimensional manifold described by eq. (4.10). We save the proof for the next subsection where we consider the more general geodesics necessary to prepare coherent states where more than one of the a i are nonvanishing. Given these simple geodesics describing coherent states (4.12) with a single excited mode, we can easily find their complexity as in section 2.3.

Perturbations of simple geodesics
Here we would like to examine the effect of exciting some subset of the normal modes with a shift producing x i = a i . To motivate the conjecture which we will prove in the following, let us review: First, with no such excitations at all, it was found in [1] that the optimal geodesics preparing the ground state were confined to an R N submanifold of the full GL(N, R) geometry. In the previous section, it was found that for two coupled oscillators that the geodesics preparing a coherent state in which a single normal mode was excited were confined to a H 2 ×R submanifold of the full R 2 GL(2, R) geometry. That is, the motion of the geodesic was still confined to the normal mode subspace for the second unexcited oscillator. However, when both normal modes were excited, we had to consider the geodesic motion in the full six-dimensional geometry, as described in section 2.4. Now the geodesics describing optimal circuits to prepare coherent states in our (regulated) scalar field theory are governed by the N (N + 1)dimensional geometry R N GL(N, R). However, given the previous observations, it is natural to conjecture that if we are considering coherent states (4.9) where only K of the normal modes are excited, then the motion is confined to the normal mode subspace R N −K for the unexcited modes, while the geodesics explore the full R K GL(K, R) subspace describing all of the gates acting on the remaining normal modes. That is, for these states, the optimal geodesics are confined to a (K 2 + N )-dimensional submanifold of the full geometry, described by where U K ∈ GL(K, R), d T is a K-dimensional vector (with entries u i (s)), and D is an (N − K)×(N − K) diagonal matrix (with entries e ya(s) ). For convenience, we have arranged the basis of normal modes so that the first K modes are being excited with a i = 0.
Given the ansatz (4.14), we can use eq. (4.11) to describe geodesics restricted to move on this (K 2 + N )-dimensional submanifold. However, we would like to show that geodesics lying within this subpace are in fact geodesics of the full R N GL(N, R) geometry. Hence we consider perturbing the above trajectories as followŝ where V , X, Y and Z represent small first-order excursions away from the submanifold described by eq. (4.14). Here, V , X and Y fill out the three 'zero' blocks on the left-hand side of U and Z comprises the off-diagonal components of the central (N − K)×(N − k) block. We have also introduced a small expansion parameter ε here and so that if we substituteÛ into eq. (4.11), the particle action can be expanded as If we set δU = 0, the variation of L 0 (U ) yields the geodesic equations on the submanifold of interest, i.e., R K (GL(K, R) × R N −K ). The order ε and higher order terms will contribute to determine the geodesics in the full geometry as they move away from the submanifold. However, the terms of order ε 2 and higher will vanish in the equations of motion if we simply set the components of δU to zero. The dangerous terms are those linear in ε since they may yield nonvanishing terms which do not vanish in the equations derived from variations of the components of δU , i.e., these terms may produce source terms which drive the geodesics away from the submanifold. Therefore our goal is to verify that in fact L 0 vanishes. Towards the latter goal, let us begin by writing the inverse ofÛ to first order in perturbations: We then examine the expansion of the tangent vector 18) and let us explicitly write out the zero'th order term Given this last expression and the form of the action (4.11), we can conclude that a nonvanishing order ε term will arise in eq. (4. 16) if and only if ∂ s δU U −1 − ∂ s U U −1 δU U −1 has contributions proportional to the same matrix generators as appear in ∂ s U U −1 , i.e., the O(ε) term in eq. (4.18) has nonvanishing components in the same entries as eq. (4.19). However, given our explict expressions above, it is straightforward to show that all of these entries vanish. For example, where the only potential overlap with eq. (4.19) is in the central block. However, since D is a diagonal matrix while Z has only off-diagonal components, these contributions are orthogonal in the sense of the inner product (2.19) on the matrix generators. As we argued above, since we were able to show that L 0 vanishes in eq. (4.16) above, we can conclude that the geodesics determined with L 0 (U ) on the R K (GL(K, R) × R N −K ) submanifold are in fact geodesics in the full geometry R N GL(N, R). In particular, notice that if we choose K = 1, i.e., our target state is only excited in one normal mode, this proof shows that there is a simple geodesic in an (N +1)-dimensional slice of the full geometry which takes the form R (GL(1, R) × R N −1 ) = H 2 × R N −1 , discussed in the previous section. In this case, the present argument is a generalization of that presented in section 2.3, in which we showed this geometry plays a role in determining simple geodesics for N = 2.
We will not examine here the geodesics in the more general case where K ≥ 2, as it seems that this will demand rather intensive numerical work. For example, the numerical results in section 2.4 are easily extended to the case of K = 2 for the present discussion with larger values of N . However, we would remark that if we excite K normal modes but all with small amplitudes, it is straightforward to show that to leading order the optimal geodesics can be evaluated using an (N + K)-dimensional submanifold of the form (H 2 ) K × R N −K -see discussion in the next section. Hence, for example, eq. (2.59) would be easily extended here to give the change in the κ = 2 complexity for a QFT state of this form, as we discuss in the next section.
Recall that it seems that the simple geodesics found in section 2.3 (i.e., K = 1 and N = 2) actually seem to provide the optimal geodesics for the corresponding family of target states. Strong evidence for this claim came from our numerical studies in section 2. 4. An interesting open question is whether the generalization of these simple geodesics found here for larger values of N and K will actually provide the optimal geodesics.

Complexity for simple target states
Now we would like to evaluate the complexity of coherent states in the free scalar field theory using various cost functions. We will focus on two situations: a) where a single mode is excited and b) where many modes are excited but all with small amplitudes. Recall that the complexity of the ground state (4.6) is divergent because the complexity is dominated by contributions of the UV modes [1,2]. In particular, with the F 2 cost function (1.7), this leading divergence takes the form where V is the spatial volume, δ is the short-distance cutoff (i.e., the lattice spacing) and d is the spacetime dimension of the scalar field theory. Further we have introduced m = 1/δ (as in eq. (4.2)) and ω k ∼ 1/δ for a typical UV mode. The form of this divergence did not match the leading divergence (1.1) found for holographic complexity [21] and hence the κ measures (1.8) were introduced in [1] to ameliorate this problem. With these cost functions, the leading divergence becomes Let us add that following the reasoning presented in [41], one can show that the Schatten p = 1 cost function yields the same leading divergence in the vacuum complexity as for the κ = 1 complexity (or the F 1 cost function).
a) Single excitation: In the previous section, we have argued that the simple geodesics found in section 2.3 also describe the optimal circuit preparing QFT coherent states (4.12) with a single excited mode, for the F 2 and κ = 2 cost functions. Hence we can apply our earlier results to evaluate the complexity of these states. For example with the κ = 2 cost function (1.8), we would have where, in analogy to eqs. (2.46) and (2.47), we have with a i = a/x 0 and w j = mω j /ω 2 R . In particular, we can evaluate the difference between the complexity of this coherent state and the complexity (4.8) of the ground state, which yields precisely the same result as for two coupled harmonic oscillators in eq. (2.50) with the substitution a + , w + → a i , w i , (4.25) Further, we can consider various limits of this result in analogy to those presented at the end of section 2.3. For example, to leading order for a i 1, we have ∆C κ=2 ∝ a 2 i as in eq. (2.51). The arguments in the previous section can also be extended to the F 1 metric (3.1) and p = 1 Schatten cost function (3.19). This would combine the reasoning given in section 4.1 and extending the perturbative arguments given below for the case of small amplitudes, but we do not present the details here. Hence the results for the simple geodesics can be extended to give the Schatten complexity for QFT coherent states with a single excited mode, with where again ∆ i and C k are given in eq. (4.24). Similarly, for this class of states, the F 1 cost function would be extremized by the L-or J-shaped paths described in section 3.1. Hence from eq. (3.15), the F 1 complexity becomes where D L and D J are the costs given in eq. (3.14). Of course, as described for eq. (3.15), the L cost applies for w > 1, |a| ≤ 2 or w < 1, √ w|a| ≤ 2, and the J cost applies otherwise.
One may also be tempted to extend the previous analysis to the higher κ cost functions (1.8). In this case, we would have where C i are given in eq. (4.24). This expression can be evaluated numerically for the simple geodesics given the expressions in eqs. (2.45) and (2.46). In the limit a i → 0, the simple geodesic reduces to a straight-line geodesic with u i (s) = 0 and of course, ∆C κ → 0. It was shown in [1] that these straight-line geodesics were still optimal geodesics preparing the vacuum state for general κ measures, not just for κ = 2. Hence the vacuum complexity is correctly given by C κ,vac = N k=1 (C k ) κ , but for the coherent states with a i = 0, eq. (4.28) only provides an upper bound on the complexity. Hence the above expression provides an upper bound on the increase of the complexity for these special states.
An interesting feature of the result in eq. (4.25) is that ∆C κ=2 is finite, i.e., there is a single contribution from the excited mode. In particular, the sum over the contributions from the UV modes causes C κ=2 to diverge for both the vacuum and the coherent state (in the limit δ → 0), but these UV divergences cancel in the difference. In contrast, we might carry out the analogous calculations with the F 2 cost function but in this case, the leading contribution takes the form where C 2,vac and ∆C κ=2 are given by eqs. (4.21) and (4.25), respectively. Hence combining these expressions, we find that ∆C 2 vanishes as δ b) Small amplitudes: Another interesting case which is relatively easy to analyze is that of excited states where a number of modes are excited but all with small amplitudes, i.e., with a k 1 for every mode. As alluded to in the previous section, the final result will be that, to leading order, the modes decouple and ∆C is simply given by the sum of the leading results found when simply exciting a single mode. In the following, we will demonstrate that this result applies for the F 1 , κ = 2 and p = 1 Schatten measures, using the techniques developed in section 4.1. Let us summarize the (leading) result for the increase in the complexity for each of these cost functions here, where the sums run over the excited modes. Let us begin with the κ = 2 cost function, where we are generalizing the arguments made for two coupled harmonic oscillators in section 2.5. Imagine that we have a number of modes excited but that a k 1. We wish to construct a perturbation expansion in which we designate the excitations as first order, i.e., a + ∼ O(ε). Then using the formalism of section 4.1, we assume that u T = εu (1)T + O(ε 2 ) and Expanding the tangent then yields 32) and the particle action (4.11) in the κ = 2 cost function becomes The O(1) part of the Lagrangian is the precisely same as for considered in [1], and the solution which prepares Gaussian states (with a k = 0) is a diagonal matrix and u T = 0. Therefore in our perturbative expansion for small excitations, we may assume where D is a diagonal matrix and Z (i) are completely off-diagonal. Furthermore, since the diagonal is a local minimum of the zeroth order Lagrangian, substituting this into eq. (4.33) gives an expression of the form where F (Z (1) ,Ż (1) ; D,Ḋ) is quadratic in Z (1) , positive semidefinite and vanishes if and only if Z (1) =Ż (1) = 0. Because there is no term linear in ε, the optimal solution for Z (1) is just zero since any nonvanishing Z (1) would only increase the cost at second order of ε.
With this choice, the separate modes simply decouple and the resulting cost function indeed describes motion on (H 2 ) N . That is, the motion in the full R N GL(N, R) geometry is restricted to a (H 2 ) N submanifold to leading order when the excitations are small. According to the previous result (2.51) for the simple geodesics in the hyperbolic geometry, we find the leading order change of complexity ∆C κ=2 is given by the expression in eq. (4.30) above. Of course, extremizing the κ = 2 cost function (1.8) also extremizes the F 2 cost function (1.7). Hence using eq. (4.29), we can evaluate ∆C 2 for only small excitations as Now we turn to the p = 1 Schatten norm, where at the end of section 3.2, we already argued that for two coupled harmonic oscillators, the leading order result for ∆C is simply the sum of those for the individual modes. It is straightforward to extend the above perturbative argument to the Schatten norm for the free scalar field theory, i.e., for N coupled modes. In order to proceed, we need to consider the eigenvalues of the square matrix as whereÛ is defined as in eq. (4.10) and then the matrix M is given by The eigenvalues of M are labeled as γ i with i = 1, · · · N . The general p = 1 Schatten cost function is then defined as Our perturbative construction again begins with the small excitations where a k ∼ O(ε). It is straightforward to show that the zeroth order solution is then given by u(s) = 0 and which means that here the straight-line geodesics also provide the optimal circuit for the Schatten cost function. Perturbing around these solutions as above, we have where D is a diagonal matrix with D ii = e 2y i and Z is a completely off-diagonal perturbation. The leading perturbations of the eigenvalues are now given by where substituting eq. (4.41) into eq. (4.38) has produced an expansion M = M (0) + εM (1) +ε 2 M (2) +O(ε 3 ) and we combine all of the higher order terms as ∆M = M −M (0) . Further, since U N is diagonal to leading order, M (0) is a diagonal matrix as well and so the eigenvectors take the special form: (v i ) a = δ ia . As a result, the perturbations δγ i all come from the diagonal components of ∆M , i.e., eq. (4.42) yields δγ i = (∆M ) ii . However, it is straightforward to show that M (1) has only off-diagonal components and so there is no O(ε) contribution to δγ i . Hence we may focus on M (2) to find the leading perturbations of the eigenvalues. First of all, we can easily see the second term in (4.38) provides us with the term like ε 2 e −2y iu (1)

iu
(1) i in δγ i , which implies the leading corrections contain the hyperbolic geometry H 2 for every mode. Secondly, the corrections on δγ i from the first term in (4.38) provide terms which are quadratic in Z (1) andŻ (1) , and hence we solve the corresponding equations of motion with Z (1) = 0. Finally, as eq. (3.28), we consider the square of Schatten norm Now the first term is precisely the κ = 2 cost function, which in our perturbative expansion describes motion in the restricted subspace (H 2 ) N of the full geometry, as above. Combining these observations with the arguments at the end of section 3.2, it is straightforward to show that to leading order, the simple geodesics for each of the individual modes extremize the above squared cost function and then the p = 1 Schatten cost function. Hence we may simply sum the leading order results for the change in complexity given in eq. (3.34) for each of the decoupled modes to find the expression for ∆C Schat given in eq. (4.30). 26 Lastly, to close this section, we consider the case of small excitations with F 1 cost function. Firstly, we re-iterate that the F 1 cost function depends on the choice of the basis of generators M I . Here we work with the normal mode basis where the M take the simple form given in eq. (2.16), i.e., [M ai ] cd = δ ac δ id . Again we construct a perturbative expansion with the a i ∼ O(ε) and at zeroth order, we begin with the simple straight-line solution (without any excitations). We then consider perturbations of the F 1 cost function, 26 Let us note that this discussion can be easily adapted to show that for states in which a single mode is excited, the full result of the simple geodesics can be applied. That is, the increase in the complexity is given by eq. (3.33) with the substitution w + , a + → w i , a i , where the subscript i indicates the mode which is excited. The discussion is almost the same. We only need to replace the original solutions by u(s) = (0, · · · u i · · · , 0) and notice the eigenvector for this mode is also (0, · · · 1 · · · , 0). While the simple geodesic is clearly a solution of the restricted cost function analogous to eq. (3.27), we can perturb around this trajectory to find that it is also extremizes the full cost function.
where the index I ∈ {ij, 0i} with i, j = 1, · · · N , and given the simple form of the generators, the components Y I are read off from the entries of V = ∂ sÛÛ −1 . As above, we assume u T = εu (1)T + O(ε 2 ) and expand U N as in eq. (4.34), which yields Now the leading perturbation of eq. (4.44) comes from the second term above which produces O(ε) contributions with |Y ij | with i = j and |Y 0i |. Here we are using the fact that the original simple solution, i.e., the first term, only contains Y ii components. Now because of the absolute value for all of the terms in (4.44) and the boundary conditions Z (1) (s = 0) = 0 = Z (1) (s = 1), we minimize the |Y ij | (with i = j) contribution by setting Z (1) (s) = 0. Finally the measure of the optimal path should have N copies of the analogous structure in eq. (3.7), which are extremized by the L-shaped paths (for small a i ). Hence to leading order in our expansion, the F 1 complexity becomes the sum of the D L costs in eq. (3.14) for the individual modes and then ∆C 1 is given by the expression in eq. (4.30).

Fubini-Study approach for circuit complexity
In this section, we apply the Fubini-Study approach proposed in [2] to examine the complexity of coherent states (2.5) for a pair of coupled harmonic oscillators. In contrast to the Nielsen approach, which defines a geometry on the space of unitaries (1.3), this method makes use of the Fubini-Study metric to define a geometry on the space of states. First, to introduce the basic definitions, let us imagine that the space of states of interest is covered by some convenient set of coordinates λ µ -we will be explicit about the coordinates in our calculations but for the time being one might think of the coordinates in eq. (2.24). In the following, we focus on a family of pure states |ψ(λ) and then we can consider the quantum fidelity as the inner product between two such states, e.g., [64,65], The quantum information metric then measures the distance between nearby states as The quantum information metric is also known as the fidelity susceptibility since it encodes the response of the fidelity to small changes in one of the states. 27 In the present case of pure states, eq. (5.3) also corresponds to the desired Fubini-Study metric. This metric may also be evaluated with the following expression Then following [2], we consider curves λ µ (σ) on the space of states parameterized by σ ∈ [0, 1] which take us from the reference state to the desired target state, i.e., We then assign a cost to each of these trajectories as the distance as measured by the Fubini-Study metric (5.3), whereλ µ (s) = dλ µ (s) ds specifies the tangent vector to the trajectory. The complexity assigned to the target state is then the minimal distance according to this measure, i.e., the complexity is the length of the geodesic in the state space equipped with the Fubini-Study metric.
Before proceeding with our calculation of the Fubini-Study complexity for coherent states, it is interesting to express this approach in a way that is closer to the circuit construction introduced in eq. (1.3). In particular, given a trajectory described by a particular choice of λ µ (σ), we may express the corresponding states as where O µ (λ) is the set of Hermitian operators which generate the evolution of state |ψ(λ) in the λ µ direction, i.e., i∂ µ |ψ(λ) = O µ (λ) |ψ(λ) . (5.8) 27 We might add that in the context of the AdS/CFT correspondence, the information metric or fidelity susceptibility for boundary states deformed by a marginal operator was proposed to be described by the volume of maximal time slice in AdS spacetime in [66]. Of course, the latter is also the conjectured dual of complexity according to the CV proposal [15,16]. Different proposals for the holographic dual of information metric are also discussed in [67][68][69][70].
Note that we may think of the operators O µ (λ) as being linear combinations of the O I appearing in eq. (1.3). We show a λ dependence to indicate that these linear combinations vary as we move through the space of states. However, this leaves the definition of the O µ (λ) ambiguous since, at any particular point, there will be degenerate operations which leave the state unchanged, i.e., O 0 (λ)|ψ(λ) = 0. Therefore, in general, one finds that the space of states has a smaller dimension than the space of unitaries, as will be illustrated by the example discussed below. Given eq. (5.8), we can also rewrite the Fubini-Study metric as connected correlation functions of the operators O µ , Let us also add that ref. [2] also proposed an alternative formulation where only one gate acts at any given point in the circuit. In preparing the vacuum state of the scalar field theory, this formulation gave a result similar to that of the F 1 measure. However, this formulation was developed in [2] with a limited gate set and so it would be interesting to extend it to the more general setting discussed here. 28

FS complexity of two harmonic oscillators
We would like to describe the coherent states discussed in section 2 as where i, j ∈ {+, −}. The 2×2 coefficient matrix A 2 is given by and U 2 is the GL(2, R) matrix given in eq. (2.22). The explicit form of A 2 is given by the upper left 2×2 block found in eq. (2.24), and as noted there, A 2 is independent of z. 29 Let us parametrize the displacements of the coherent states in terms of the dimensionless coordinates v ± with a ± ≡x 0 v ± , where we have introduced a convenient length scalex 0 in this definition. 30 Then our family (5.10) of coherent states is described by five dimensionless coordinates λ µ = {y, ρ, x, v ± }, and by construction, the origin of this coordinate system corresponds to the reference state (2.4). Of course, this is one less coordinate than described the unitary transformations in section 2. Now by the methods introduced above, we can define the Fubini-Study metric for the space of states |ψ(y, ρ, x, v ± ) . The metric can be constructed with eq. (5.3) by evaluating the integrals where the wave function ψ(x + , x − ; y, ρ, x, v ± ) is defined in eq. (5.10). Alternatively, we can calculate the fidelity (5.1) 13) and then evaluate the metric with eq. (5.4).
Using either method, we find the Fubini-Study metric is given by If we begin by focusing on Gaussian states with a ± = 0, we expect that the optimal trajectories will not involve motion in the v ± directions and hence we focus on the first three terms in eq. (5.14). This three-dimensional subspace has the geometry R × H 2 . As noted above, the reference state corresponds to the origin, i.e., y = 0 = ρ (while the angle x is unspecified). Hence the geodesics are simply lines moving along the R and radially outward in the hyperbolic space, i.e., y = y 1 s, ρ = ρ 1 s and x = x 1 where (y 1 , ρ 1 , x 1 ) is the position specifying the target state [2]. However, the complexity or the length of the geodesic is precisely the same as found using the Nielsen approach [1], except for the overall constant factor. 31 The full Fubini-Study metric (5.14) has a form similar to the Nielsen metric in eq. (2.28) defined on the space of unitaries, although the dimension of the geometry differs by one as we already noted. In order to define complexity with this metric (5.14), we would need to solve the corresponding geodesic equations, but generally the only tractable approach is to find numerical solutions, as we did in section 2.4 for the 31 Note that our conventions were such that the metric (2.28) for the Nielsen geometry had an extra overall factor of 2 compared to the Fubini-Study metric (5.14), i.e., ds 2 Nielsen = 2dy 2 + 2dρ 2 + · · · while ds 2 FS = dy 2 + dρ 2 + · · · .
Nielsen geometry. However, as in section 2.3, we can find a simple analytic solution here for states with a single excitation, e.g., a + = 0 and a − = 0. Examining the full geodesic equations, we find it is consistent to set x = 0 and v − = 0 in this case. Hence we are simply solving for the geodesic equations in the reduced Fubini-Study metric where as before y ± = y ± ρ and for convenience, we have setκ = 1. We note that this geometry again has the familiar form H 2 × R but comparing to the corresponding geometry in section 2.3, we see that to identify this metric with eq. (2.40) and the corresponding geodesic equations, we must set The initial boundary conditions are simply y +0 = 0 = y −0 = v +0 and to match the final target state (2.5) with a − = 0, the final boundary conditions are where w ± are the same dimensionless ratios as in eq. (2.43), whileã ± ≡ ω R a ± . 32 We note that, of course, the boundary conditions for y ± are the same here as in eq. (2.42), but u + and v + are different coordinates and so their boundary conditions do not match.
Using the above observations, we can use the solution found in section 2.3 given by eqs. (2.45) and (2.46) to produce the simple geodesics for the Fubini-Study metric (5.15), which take the form where the sign of B FS is chosen to match that ofã + . The coefficients above can also be derived from eq. (2.46) by replacing a 2 + w + →ã 2 + and w + → 1/w + at the same time. Further, one can verify the above solutions satisfẏ 20) and as expected, this combination of the velocities is constant along the geodesic. Certainly the new trajectories in the Fubini-Study geometry (5.15) should be different from those in the Nielsen geometry (2.40) because of the differences in A FS , B FS compared to A, B in eq. (2.46). Of course, the y − part of the trajectory is identical in both cases. However, to make clear that the simple geodesics describe distinct circuits in the Nielsen and Fubini-Study geometries, we compare the evolution of the states as described by the 3×3 coefficient matrix in eq. (2.24), which for the simple geodesics reduces to

The final boundary conditions (2.42) fixes the integration constants as
(5.21) The scale x 0 appears in A FS because by definition this 3×3 matrix is contracted with x a = (x + , x − , x 0 ) to construct the wave function. In general, the comparison depends on the combination ω R x 0 appearing in [A FS ] 0+ , 33 however, to simplify the comparison we might simply set ω R x 0 = 1. With this choice, figure 8 illustrates an example comparing the components of A Niel and A FS for a fixed target state (and reference state). As expected, the evolution of [A] −− ∝ ω − is identical in both approaches because this component is controlled entirely by y − (s), which we already noted is the same in the two cases. The evolution of [A] 0+ ∝ Λ + distinguishes the two trajectories, but in both cases, this component is monotonically increasing from zero to the final value in both cases. The difference between the two circuits is shown most dramatically in [A] ++ ∝ ω + . In the Nielsen approach, the evolution of this component is concave down, i.e., it begins by increasing but it overshoots the final value and so it must decrease again towards the end of the trajectory. In contrast, for the Fubini-Study approach the evolution is concave up, i.e., this component begins by decreasing but this direction is reversed in the latter part of the geodesic so that it can reach the final positive value. This reversal of the concavity might be expected from the fact that the metrics look identical under the identification (5.16), i.e., where the sign of y + is reversed. We should also examine the length of the simple geodesics in the Fubini-Study geometry. Towards this end, we first evaluate ∆ FS = A 2 FS + B 2 FS using the expression in eq. (5.19) to find where we are again usingx 0 = 1/ω R in expressing the result in terms of a + . Comparing this result to ∆ in eq. (2.47) for the Nielsen construction, we see that the two expressions generally differ because of the factor of ω 2 R x 2 0 in the final expression above. However, quite remarkably if we set ω R x 0 = 1 (as above), we find that ∆ FS = ∆. Now the Fubini-Study complexity is given by This result is naturally compared with the F 2 complexity in the Nielsen approach, and while the two complexities differ in general, there is a remarkable agreement between the two complexities if we simplify the analysis by choosing ω R x 0 = 1. We emphasize that with this choice, the Fubini-Study and Nielsen complexities agree even though we have shown above that the two approaches are constructing different optimal circuits. We can highlight the difference between the Nielsen and Fubini-Study geometries by introducing the same coordinate systems for both. In fact, the coordinates (y, ρ, x, v ± ) match those introduced in footnote 7. If we focus on the subspace x = 0 = v − , we can compare the coefficient matrices in eq. (5.21) to find u + = −e y + v + (where we also set ω R x 0 = 1 as before). Then the Nielsen geometry (2.40) becomes Comparing this expression to the Fubini-Study metric (5.15) we see that even though both describe an H 2 × R geometry, the physical states are assigned to the geometries in different ways. That is, if we choose a particular state described by particular values of the coordinates (y + , y − , v + ) in eq. (5.21), then the distances to nearby states (y + + δy + , y − + δy − , v + + δv + ) are very different in the two metrics in eqs. (5.15) and (5.24). Of course, if we fix our attention on the plane v + = 0, e.g., to evaluate the ground state complexity, the metrics are the same (except for an overall factor of 1/2). However, when we move away from this 'normal mode subspace' [1], we should not expect that the optimal circuits between two states (or the corresponding complexities) to be the same in the two geometries. ) that Λ ± = [A] 0± /ω 2 R , and we defined mω ± = [A] ±± .) From this figure, we see that the optimal circuits from the two approaches follow different trajectories even though they begin and end at the same states. This difference appears most dramatically for mω + (s)/ω 2 R and [A] +− (s)/ω 2 R in the upper two plots.
The comparison above highlights that when considering coherent states, the Nielsen and Fubini-Study approaches to complexity are really different systems. That is, for a given reference state and target state, the optimal preparation chosen by these two approaches moves through different families of intermediate states. This stands in contrast to the preparation of simple Gaussian states, where the optimal trajectories were the same for both approaches. We have seen that these differences already arise for the simple geodesics preparing states where only one of a ± is nonvanishing. We provide another example in figure 9 for a general target state in which both a ± are nonzero (where we again set ω R x 0 = 1). We solved for the optimal trajectory for the Fubini-Study and the Nielsen approaches numerically and then translated to trajectories into the corresponding (physical) components of the A matrix in eq. (2.12). As shown in the plots, the two approaches prepare the same target state through different families of intermediate states. Further for general states like this in which both a ± are nonvanishing, the complexity derived by the two methods also differs. In the example shown in the figure, C FS = 1.518 while C 1 = 1.511. 34 For comparison purposes, let us note that the ground state complexity is C FS,vac = C 1,vac = 1.221 for this example.

Discussion
Refs. [1,2] provided the first calculations of complexity in quantum field theory. In the present paper, we extended this analysis, which examined the ground state of a free scalar field theory, to evaluate the complexity of excited states in the same theory.
In particular, we considered coherent states with a nonvanishing expectation value of the scalar field (but for which the expectation value of the conjugate momentum was vanishing). Following the analysis of [1], we began by examining in detail the complexity of the analogous coherent states in a pair of coupled harmonic oscillators in sections 2 and 3, and then extending the results to the free scalar in section 4. The generators of the gates preparing our coherent state naturally gave rise to a group structure R N GL(N, R), which is a simple extension of the GL(N, R) structure found in [1]. 35 While this analysis focused on Nielsen's geometric approach [51][52][53] for evaluating circuit complexity, we also considered the Fubini-Study approach proposed by [2] in section 5. Before proceeding, let us remind the reader that a brief discussion of the complexity of coherent states appeared in [47]. This recent work was one of the first investigations 34 The difference is small but significant, i.e., we are confident that the accuracy of our numerical calculations goes well beyond the fourth significant digit here. 35 We reiterate that for the most general coherent states (e.g., for which the expectation values of the momenta are also nonvanishing), this group structure is enlarged to R 2N Sp (2N, R), which extends the Sp(2N, R) for general bosonic Gaussian states (with vanishing expectation values) [41,48]. of the application of Nielsen's geometric approach to evaluate state complexity in a quantum field theory, and as an example in a free scalar field theory, they consider coherent states where both x and p can be nonvanishing for a single mode. However, their analysis differs from ours in a number of essential ways: First of all, rather than considering an unentangled reference state, [47] considers preparing their coherent states beginning with the vacuum state of the field theory. Further, the gate scale introduced for the shift gates in eq. (2.9) is implicitly set by the frequency of the excited mode in [47]. In particular, x 2 0 = 2/(mω k ) is chosen there. Finally, we would add that the complexity is evaluated there by optimizing a somewhat unconventional cost function and the circuits considered are generally not unitary. Hence there is no substantive overlap between our work and the discussion in [47].

Optimal Trajectories/Circuits:
When applying the Nielsen or the Fubini-Study approach to coherent states for the system of two coupled harmonic oscillators, we could only find the desired geodesics numerically for states in which both normal modes were excited, i.e., both a ± = 0. One of the interesting features of these geodesics was that generally they pass through nonvanishing values of x. The physical significance of this feature appears in eq. (2.24), where we see that [A] +− = 0 ( = [A] −+ ) with x = 0 (and also ρ = 0). Therefore, even though the two normal modes are unentangled in both the reference state (2.13) and target state (2.14), they become entangled in the intermediate states that appear in the optimal circuit joining these states. This behaviour is illustrated schematically in figure 10. It is also exhibited by the explicit examples shown in figures 4 and 9. We emphasize again that this behaviour is common to both the Nielsen and Fubini-Study approaches.
In section 4, we showed that the 'complexity' of determining the optimal trajectories grew with a larger number of excitations. In particular, determining the optimal circuit for states with K normal modes excited, required studying the geodesic equations on a R K GL(K, R) manifold. The remaining unexcited modes decouple and they are simply prepared with the linear application of the corresponding scaling gates. It maybe interesting to use numerical methods to investigate the general properties of optimal circuits and corresponding complexity for states where K ≥ 3.
However, a particularly simple case is K = 1, i.e., only one normal mode was excited. In this case, we found analytic solutions for a class of simple geodesics for the Nielsen approach in section 2.3. These geodesics moved in a H 2 slice of the full geometry, involving the coordinates corresponding to scaling and shift gates for the excited mode, e.g., y + and u + in eq. (2.40). However, we still had to rely on numerical tests to support the claim that these simple geodesics were the optimal geodesics connecting Figure 10. Schematic diagram of the paths followed by the optimal circuits connecting the unentangled product state A R to three different target states A Ti . The states A T1 and A T2 have a + = 0 and a − = 0, respectively. The optimal circuits preparing such states in which only one of the normal modes is shifted remain in the x = 0 plane, i.e., [A] +− = 0. Therefore the normal modes are unentangled for all of the states along these trajectories. In contrast, the trajectory preparing A T3 begins and ends with [A] +− = 0 but this component is nonvanishing everywhere away from these endpoints. That is, both the reference state and target state are unentangled but the optimal circuit introduces entanglement in the intermediate states when both a ± = 0 in the final state. the reference state (2.13) to a target state (2.14) with only one normal mode excited. These simple geodesics played a role not just for the F 2 and κ = 2 cost functions as discussed in section 2.3, but also for the Schatten p = 1 cost function as described in section 3.2 and also with the Fubini-Study approach in section 5. The analysis of these states showed a similar behaviour for the F 1 cost function, in that the optimal circuit only involved the scaling and shift gates for the excited mode, while the other modes decoupled.

Schatten measures:
In section 3.2, we investigated the complexity of our coherent states for a cost function constructed from the p = 1 Schatten norm (3.19). This cost function was first suggested in [41] as a replacement for the F 1 cost function. There it was observed that in preparing the ground state, this Schatten cost function would produce the same optimal circuit and complexity as the F 1 measure constructed in the normal mode basis. However, these results are basis independent when described in terms of the Schatten norm. In comparing the results in sections 3.1 and 3.2, one of the most striking results is that the F 1 and Schatten measures no longer give the same circuits or complexities when considering coherent states. For example, the increase above the complexity produced by a small amplitude excitation produced by the F 1 measure gave ∆C 1 ∝ |a| (see eq. (3.17)) while the Schatten norm gave ∆C Schat ∝ a 2 (see eq. (2.51)). We return to these different behaviours for these two cost functions below.
The examination of the Schatten cost function in section 3.2 focused on the complexity of coherent states for two coupled harmonic oscillators. However, this is easily extended to the (regulated) scalar field theory where the circuits act in the group R N GL(N, R). For example, with p = 1, eq. (3.25) is replaced by where the γ i are the eigenvalues of V T V . Note that the range of i implicitly indicates that γ N +1 = 0, i.e., V T V is represented by a square (N +1)×(N +1) matrix but one of the eigenvalues automatically vanishes (because the last column of V is filled by zeros, as in eq. (3.20)). We note that the number of eigenvalues matches the number of types of gates that are applied to prepare the ground state, i.e., the optimal circuit only uses the scaling gates for each of the N normal modes. This match is why the p = 1 Schatten complexity agrees with the F 1 complexity for the ground state. However, the N eigenvalues encode information about the shift gates, as well as the scaling gates, when preparing the coherent states, and so as noted above, this agreement does not extend to these states. Generally the p = 1 Schatten cost function also involves a complicated coupling between the different modes, e.g., as is implicit in the singular values given in eq. (3.23). However, the modes seem to decouple when evaluating the complexity of coherent states where a single mode is excited, and the optimal circuit follow the same simple geodesics described above for the F 2 or κ = 2 cost function. We were able to prove these geodesics extremized the full Schatten norm (3.25) by considering a new cost function L 0 = V 2 1 , see eq. (3.28). This could be decomposed into two parts: L 0 = γ 1 + γ 2 and L 1 = γ 1 γ 2 . The first coincides with the κ = 2 cost function and so the simple geodesics extremized this term. It was then straightforward to show that they also extremized L 1 .
If we recall that there is a family of Schatten norms (3.19) labeled by a positive integer p, it is interesting that the previous reasoning can be extended to the higher p norms. That is, we can argue that the simple geodesics extremize the Schatten cost functions for general p as follows: First it is straightforward to show the recursion relation Now we have shown that the simple geodesics extremize V 1 and L 1 and therefore if they also extremize V p and V p−1 , then the same geodesics will extremize V p+1 . 36 Since the p = 2 norm corresponds to the F 2 norm, it is also extremized by these simple geodesics. Hence beginning with p = 1, 2, we can work iteratively to show that our simple geodesics are in fact also geodesics for the general p Schatten cost functions. Given the previous result, we can apply the interesting property of Schatten norms that A p ≥ A q for 1 ≤ p ≤ q ≤ ∞ [62]. This leads us to conclude that given a particular simple geodesic describing the optimal circuit for a particular state, the complexity of the same circuit increases if we increase the index of the Schatten norm with which the complexity is evaluated, i.e., We stress that our discussion above focused on the simple geodesics describing the optimal circuits for states with a single excitation (and this discussion easily generalizes to the case of N normal modes but only a single excitation). General geodesics of the F 2 or κ = 2 measures, i.e., for states with multiple excitations, will not extremize the auxillary functional L 1 and so they will not be optimal trajectories for any of the Schatten norms except p = 2. However, we did argue at the end of section 4.2 that it is possible to consider multiple excitations as long as the amplitudes are small, i.e., a i 1. In this case, the different normal modes can be decoupled at least to first order in a perturbative expansion.
To close here, we would like to point out that we can use a modified Schatten cost function of the form, i.e., we eliminate the overall p'th root in eq. (3.19). These cost functions are rather analogous to the κ cost functions (1.8) with κ = p, i.e., optimizing these new cost functions would yield the same optimal circuit and complexity as the κ = p cost function constructed in the normal mode basis when considering Gaussian states. Therefore, the divergence structure of the ground state complexity would match for these two sets of cost functions. Of course, the advantage of using eq. (6.4) would be that the results are basis independent. However, as with the case of κ = p = 1, this agreement would not extend to the coherent states considered here. We should also note that like the κ measures, these modified Schatten cost functions are not homogeneous, i.e., the total cost associated with a path is generally not invariant under reparametrizations of s.

Fubini-Study approach:
In Section 5, we examined the Fubini-Study approach developed in [2] in some detail.
In particular, we applied this approach to examine the complexity of coherent states for a pair of coupled harmonic oscillators, the same problem that we studied using the Nielsen approach in section 2. Both the Nielsen and the Fubini-Study approaches identify the complexity of a state as the distance from a simple reference state in some geometry. Nielsen's method [51][52][53] is motivated by the definition of complexity as the number of elementary gates in the optimal circuit, and so in this case, a metric is defined on the space of quantum circuits or unitary transformations, e.g., as in eq. (2.28). Optimizing the trajectory in this space then has a direct interpretation as minimizing the number gates used in the circuit preparing the desired target state (or at least, optimizing this number according to some cost function). The Fubini-Study approach instead accounts for the complexity by keeping track of the changes of the state throughout the preparation of the target state. As its title indicates, this method makes use of the Fubini-Study metric, which defines a geometry directly on the space of states. An important difference is then that the latter geometry assigns a variable cost to specific gates, i.e., the cost depends on the details of the state on which they act, whereas the gates are assigned fixed costs in the Nielsen approach. Further, at any point in the space of states, there will be degenerate operations which leave the state unchanged, i.e., |ψ = U 0 |ψ . Therefore, in general, one finds that the space of unitaries has a larger dimension than the space of states, as illustrated by comparing the geometries in sections 2 and 5. 37 For a more detailed discussion comparing these two approaches, the interested reader is referred to [59]. However, we want to stress that the definition of Fubini-Study metric only depends on the physical parameters which characterize the states. This is clear from the definitions in terms of the fidelity in eqs. (5.2) and (5.4). For example, even though the coordinates λ µ may be dimensionful, producing a dimensionful metric, the cost is dimensionless due to the appearance of the compensating factors ofλ in eq. (5.6). Hence, the parameterx 0 , which was introduced to define the dimensionless coordinates v ± = a ± /x 0 and which appears in the metric (5.14), will never appear in the complexity or in the distance along any trajectories. Instead it will be absorbed by the boundary conditions which would be defined in terms of the dimensionful displacements a ± . In contrast, the parameter x 0 is an essential ingredient in the definition of the shift gates (2.9), which must have dimensionless generators. 38 This parameter reflects a true freedom in the choice of the fundamental gates and it will affect the final complexity evaluated using the Nielsen approach. For example, it implicitly appears in eqs. (2.51) and (2.52) through the definition of a ± = a ± /x 0 Hence we see that the Fubini-Study and Nielsen approaches must define different complexities for the optimal circuit with the same target and reference state. However, we remind the reader that the ground state complexities, and in fact the optimal circuits, were found to agree with these two different approaches [1,2]. In this case, the optimal circuits only involved GL(N, R) gates and so no additional scale was needed to define the corresponding generators. In fact, in this case, the Fubini-Study geometry can be embedded in the corresponding Nielsen geometry. However, in the case of coherent states, we saw in section 5.1 that the Nielsen and Fubini-Study approaches produced different optimal circuits for a fixed pair of reference and target states. We were able to show this analytically for the simple geodesics where only one of a ± is nonvanishing. However, even though the optimal circuits are clearly different (see figure  8), a somewhat surprising result was that the Fubini-Study complexity still matched the Nielsen complexity (measured with the F 2 cost function) if we make the choice x 0 = 1/ω R . It would be interesting to better understand this agreement. Nevertheless, when we explored the geodesics for coherent states with both a ± nonvanishing, we found that the optimal circuits produced by the Nielsen and Fubini-Study approaches were again different (see figure 9) and that the corresponding complexities were also distinct.
Complexity for free scalar field: As we described in section 4.2, the complexity of coherent states (or any state) in the free scalar field theory is UV divergent. However, considering the difference ∆C = C coh − C vac yields an interesting UV finite quantity. Hence in the following, we focus on discussing this difference, i.e., the increase of the complexity of the coherent state over the complexity of the vacuum state. However, we must add that as explained with eq. (4.29), this difference vanishes for the F 2 complexity. This same reasoning would apply for the complexity evaluated with the Schatten cost functions (3.19) with p ≥ 2. Further, this difference would also vanish for the Fubini-Study complexities if we were to extend the result of section 5 to the quantum field theory. However, we can still consider this difference when evaluating the complexity with F 1 cost function, κ = 2 cost function and the p = 1 Schatten norm, and as we will discuss below the QFT complexities produced with these cost functions are most closely aligned with the result of holographic complexity.
If we only excite a single mode of the field theory, we can use the analytic results for the simple geodesics found for the κ = 2 cost function or the p = 1 Schatten norm. That is, eqs. (2.50) and (3.33) would produce ∆C κ=2 and ∆C Schat for the full field theory with w + , a + corresponding to the frequency and shift of the excited mode. Similarly, eq. (3.15) could be used to evaluate ∆C 1 for a field theory state with a single excitation. In principle, one could use numerical methods, e.g., as in section 2.4, to study the increase in complexity for coherent states in which more than one mode is excited.
However, a simpler and more interesting situation is one where many modes are excited in the coherent state but with small shifts, i.e., a k 1 for all of the modes. As we argued in section 4.2 for these three cost functions, to leading order, the shift in the complexities for each of the individual modes can be added together to produce where the sums run over the excited modes. We would like to stress that verifying these results required a nontrivial analysis and relied on the special form of the simple trajectories for the individual modes. Here we might recall the definitions of the dimensionless ratios from eq. (2.43) where we have also substituted m = 1/δ (i.e., the inverse of the lattice spacing) and a k = δ d/2 φ k from the discussion of the lattice regularization of the scalar field theory at the beginning of section 4. While the full dispersion relation for arbitrary modes is given in eq. (4.5), we would typically only be interested in exciting low energy modes, i.e., with ω k 1/δ, and so the dispersion relation would be well approximated by ω 2 k = | k| 2 + µ 2 (where µ is the mass of the scalar in eq. (4.1)). One interesting difference here is that the leading contribution for the F 1 complexity scales as ∆C 1 ∝ |a k | while in the other two cases, we have ∆C ∝ a 2 k . We return to this point below.
fermionic states has two disconnected components, i.e., states with odd and even particle number, where the component with even particle number contains the vacuum. It is the excited states in this component whose complexity was evaluated in [41]. The precise increase in the complexity depends on the details of the excited state, but generally ∆C is finite and larger for lower energy modes. For example, considering the class of states with n particle excitations and n antiparticle excitations, but where the momenta of all of these excitations are different, ∆C κ=2 = n π 2 − i tan −1 | k i |/µ 2 , (6.9) where µ is the fermion mass. For these states, we see that with the κ = 2 cost function, ∆C κ=2 n π 2 if all of the excitations have low energy (i.e., | k i | µ) whereas ∆C κ=2 1 2 n π 2 with all high energy excitations (i.e., | k i | µ). Even more dramatically, p = 1 Schatten cost function yields ∆C Schat n π if all | k i | µ and ∆C Schat 0 if all | k i | µ. Hence the behaviour of the fermionic states (with even particle number) contrasts with the bosonic coherent states above since for the latter, excitations in the higher momentum modes generally produces a larger ∆C.
Let us conclude with a few comments on possible future extensions. One obvious extension would be to consider more general coherent states with expectation values for both the field modes and their conjugate momenta. As we commented before, this would require extending the R N GL(N, R) group structure found here to R 2N Sp(2N, R). In particular, this would allow us to follow the time evolution of the coherent states. An obvious question would be to then to examine if the complexity increases, decreases or remains constant as a coherent state evolves. Coherent states also provide an interesting forum to compare to the QFT complexity with holographic complexity. Recall that the leading divergences appearing in the QFT calculations of complexity compared well with those appearing in holographic complexity (1.1) with an appropriate choice for the cost function [1,2]. The holographic analog of our coherent states would be a bulk configuration where a bulk scalar has excited in the vacuum AdS spacetime. Here we observe that to leading order, modification of the bulk geometry will be proportional to the square of the scalar amplitude since the bulk scalar backreacts on the geometry through the stress tensor in Einstein's equations, which is quadratic in scalar field. Hence we expect that the change in the holographic complexity must also be quadratic in the scalar amplitude, which is in agreement with our results in eq. (6.5) for the κ = 2 and p = 1 Schatten cost functions. However, the F 1 cost function does not exhibit this behaviour. We plan to return to this topic and make a detailed comparison between our results for the complexity of QFT coherent states and holographic complexity in [71].