Complexity and entanglement for thermofield double states

Motivated by holographic complexity proposals as novel probes of black hole spacetimes, we explore circuit complexity for thermofield double (TFD) states in free scalar quantum field theories using the Nielsen approach. For TFD states at t = 0, we show that the complexity of formation is proportional to the thermodynamic entropy, in qualitative agreement with holographic complexity proposals. For TFD states at t>0, we demonstrate that the complexity evolves in time and saturates after a time of the order of the inverse temperature. The latter feature, which is in contrast with the results of holographic proposals, is due to the Gaussian nature of the TFD state of the free bosonic QFT. A novel technical aspect of our work is framing complexity calculations in the language of covariance matrices and the associated symplectic transformations, which provide a natural language for dealing with Gaussian states. Furthermore, for free QFTs in 1+1 dimension, we compare the dynamics of circuit complexity with the time dependence of the entanglement entropy for simple bipartitions of TFDs. We relate our results for the entanglement entropy to previous studies on non-equilibrium entanglement evolution following quenches. We also present a new analytic derivation of a logarithmic contribution due to the zero momentum mode in the limit of vanishing mass for a subsystem containing a single degree of freedom on each side of the TFD and argue why a similar logarithmic growth should be present for larger subsystems.


Introduction
In the context of holography [1], thermofield double (TFD) states play an especially important role. From the perspective of the boundary theory, such states are formed by entangling two copies of a conformal field theory (CFT) in such a manner that tracing out either copy produces the thermal density matrix at inverse temperature β > 0 for the other, i.e., where |E n L,R are the energy eigenstate of the left/right CFT, and Z β is the canonical partition function at inverse temperature β. The TFD state (1) plays a special role in holography due to the fact that it is dual to an eternal black hole in AdS [2]. Hence, it provides a particularly well-controlled setup for studying various aspects of entanglement, black holes, and quantum information, e.g., the time-evolution of entanglement entropy [3], scrambling and quantum chaos [4][5][6], firewalls [7][8][9], ER=EPR [8,10], and emergent spacetime [11]. The left and right sides of the geometry associated with the black hole dual to the TFD state are connected by a wormhole, or Einstein-Rosen bridge (ERB), whose volume increases for a time which is exponential in the number of degrees of freedom of the boundary theory [12,13]. This time is much larger than other characteristic times in holography, e.g., the time for the mutual information to saturate or the scrambling time t * ∼ β 2π log S. This implies that "entanglement (entropy) is not enough" [14] to capture the dynamics behind the horizon, and in particular that there must be some other quantity in the dual field theory which encodes this late-time evolution of the wormhole interior. It has been suggested that this quantity is the quantum computational complexity of the boundary state [10].
The concept of circuit complexity is rooted in theoretical computer science [15,16]. In classical computing, one is interested in implementing a given algorithm, i.e., a function that maps an arbitrary set of input bits to specific output bits, using the minimum number of elementary operations or gates. In quantum computing, this function becomes a unitary operation U which maps an input quantum state for some number of qubits to an output quantum state on an equal number of qubits [17][18][19]. Here, one may adopt a circuit model in which U is constructed from elementary gates -which typically act on some limited number of qubits -chosen from some fixed set of universal gates { V I }. The circuit complexity of  [23]. the unitary U is then given by the minimal number D of elementary gates V I k required to construct the desired unitary, i.e., U = D k=1 V I k . We note that with discrete elementary gates { V I }, most unitaries U can only be obtained up to some accuracy ε > 0 (in operator norm), and therefore the circuit complexity is defined up to some tolerance ε. In the context of the preceding holographic conjectures, these notions should be adapted to define the circuit complexity for the states in the boundary theory [10,[12][13][14]. That is, we seek the minimum number of elementary gates required to prepare a desired target state |ψ T from a simple (unentangled) reference state |ψ R , i.e., (2) Hence the complexity of a family of target states will be defined with respect to a particular reference state |ψ R and a choice of the gate set { V I }, as well as a tolerance ε.
In terms of the gravitational description of the holographic theory, Susskind and collaborators have put forward two proposals for quantifying the size of the ERB dual to holographic complexity. The "complexity=volume" (CV) [20] proposal states that the complexity is given by the codimension-one volume of the maximal spacelike slice that connects the left and right CFTs through the bulk. The "complexity=action" (CA) [21,22] conjecture instead suggests that a more appropriate bulk dual is given by the gravitational action of the bulk region known as the Wheeler-deWitt (WdW) patch bounded by light sheets. These two proposals are illustrated in fig. 1.
Over the past few years, a wide variety of aspects of these proposals have been investigated on the gravitational side of the duality, e.g., see refs. , but this research program is still in its nascent stages. A particular obstacle to further progress is that a precise definition of complexity is still lacking for the boundary CFT. That is, it remains to construct a concrete formulation of eq. (2) for strongly coupled CFTs, or more broadly for general quantum field theories. Clearly, such a definition will be necessary in order to firmly establish any holographic complexity proposal as a new entry in the holographic dictionary. To that end, some initial steps towards precisely defining complexity in quantum field theory have recently been made in refs. [23,[51][52][53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68][69]. In particular, ref. [23] has considered the complexity of the ground state of a free scalar field theory. As alluded above, one must first identify both an appropriate reference state, and a suitable set of gates with which to construct the target state (i.e., the ground state). However, evaluating the complexity still requires identifying the optimal circuit out of the infinite number of possible circuits which prepare the ground state of the theory. To overcome this challenge, the authors of ref. [23] has adapted a geometric approach developed by Nielsen and collaborators [70][71][72]. The essential idea in the latter is that the elementary gates acting on a chain of n qubits form a representation of the Lie algebra su(2 n ), and hence one can define a natural geometry on the associated Lie manifold. This allows one to translate the question of finding the optimal circuit to that of finding the minimal geodesic in the space of unitaries U equipped with a suitable metric. This idea has been applied to free scalar field theories in ref. [23], and subsequently extended to free fermionic theories in refs. [55,56] (see also ref. [57]). Another geometric definition of complexity based on the Fubini-Study metric has been explored for the case of a free scalar theory in ref. [51]. Despite the restriction to free scalar field theories, the results of refs. [23,51] exhibit some surprising similarities with holography in the structure of the UV divergences. While this is not sufficient to either confirm or rule out either the CV or CA proposals, it may suggest some degree of universality.
The purpose of the present work is to extend the construction of ref. [23] to evaluate the complexity of the TFD state in a free scalar field theory. That is, the state (1) will serve as our target state |ψ T , while the reference state |ψ R will consist of two copies of the reference state used in refs. [23,51], namely the product state in which all of the positions are disentangled from each other. It will be sufficient to consider circuits U generated by quadratic operators of the canonical variables to effect the necessary changes in the entanglement structure to produce the desired state (1) via eq. (2). The key feature which permits this construction is that both the reference and target states are Gaussian; furthermore, since all of the gates are generated by quadratic operators, all of the intermediate states will be Gaussian as well. We will use an approach based on the covariance matrix to manipulate these Gaussian states. This will allow us to study the time-dependence of complexity, as well as the complexity of formation of the thermal state, and to compare our results with those of the holographic complexity conjectures. We will find that the complexity of formation of the TFD state is proportional to the thermodynamic entropy, in qualitative agreement with holographic complexity proposals, but that it does not exhibit the late-time growth characteristic of holographic theories/fast scramblers [73] due to the Gaussian nature of the TFD state for free scalar theories. We note that this problem has been studied previously in refs. [52,53]. Our methods differ however, in particular in the choice of the cost function, and hence certain key features of our respective results differ as well, e.g., the time-dependence of the complexity. We will compare the different approaches and results in detail in section 7.
In addition, for QFTs in 1+1 dimensions, we study the time evolution of entanglement entropy for subregions containing an equal number of sites on each side of the TFD state. We adapt the analytic formula of the time evolution of entanglement entropy following a global quench which was constructed based on a quasi-particle picture in [74,75] to the case of the TFD state, and find a good agreement up to the effect of a zero mode which leads to a logarithmic growth for the case of a scalar with vanishing mass. We present an analytic derivation of the effect of this zero mode for a subsystem containing one degree of freedom on each side of the TFD and argue why a similar logarithmic growth is obtained for larger subsystems. We also comment on the comparison between the time evolution of complexity and that of entanglement entropy.
The remainder of this work is organized as follows: In section 2, we introduce the required preliminaries. We briefly review the approach of ref. [23] to circuit complexity, and then demonstrate how the TFD state of two harmonic oscillators can be generated by quadratic operators. In section 3, we introduce a formalism for manipulating Gaussian states based on their covariance matrices, which serves as an alternate -and indeed, more natural -characterization of Gaussian states. We then proceed to evaluate the complexity for TFD states comprised of a pair of harmonic oscillators in section 4, by considering circuits that form a representation of Sp(4, R). We then explain how to extend this construction to evaluate the complexity of TFD states of a free scalar field theory in section 5. In section 6, we study the dynamics of entanglement entropy of simple subsystems involving an equal number of sites on each side of the TFD in 1 + 1-dimensional QFTs and compare it to the dynamics of complexity. We close in section 7 with a brief discussion and outlook. Various technical details have been relegated to a series of appendices. In appendix A, we collect information about our conventions and notation. In appendix B we review how to construct the TFD state for a single harmonic oscillator in a unitary way. In appendix C, we collect the matrix generators of Sp(4, R) which are relevant for the discussion in section 4. Appendix D comments on a few details with respect to the basis of operators use to describe the circuits and the complexity in the scalar field theory. In appendix E, we provide further details on the complexity of formation of the TFD state in the diagonal basis. In appendix F, we present the proof that the shortest unpenalized circuit in the full Sp(2N, R) geometry (for the case λ R = 1 using the L 2 norm) amounts to an independent squeezing of the normal modes. We conclude with appendix G where we derive compact matrix functions for the time-dependent covariance matrix in position space, which we use for the efficient numerical evaluation of the entanglement entropy.

Preliminaries: circuit complexity and thermofield double states
In this section, we provide the relevant preliminaries for the construction of both the TFD state and the relevant quantum circuits. We start with a brief overview of Nielsen's geometric approach to evaluating quantum circuit complexity. In the second part of this section, we will present the thermofield double state for a simple harmonic oscillator in preparation for our studies of its complexity in the rest of the paper.

Circuit complexity from Nielsen geometry
As mentioned in the introduction, in applying circuit complexity to quantum field theory, we are essentially quantifying the effort needed to prepare a target state |ψ T beginning from a certain reference state |ψ R . That is, we wish to construct the shortest circuit which performs the transformation where the unitary is constructed as a sequence of elementary unitary gates, U T = D k=1 V I k , and the complexity corresponds to the number of gates D in the optimal construction. A generic gate V I can be written as an exponential V I = e −iε K I , where K I is some Hermitian operator which only acts on a few degrees of freedom, and ε is a small parameter which ensures that each gate only produces a small change on the state. In general however, there may be arbitrarily many different circuits, i.e., different sequences of elementary gates, which yield the same target state, and so the primary challenge in defining complexity for the latter is identifying the optimal circuit from amongst the infinite family of possible constructions of U T .
To overcome this challenge, Nielsen and collaborators [70][71][72] introduced a geometric approach to identify the optimal circuit, which ref. [23] adapted to define the complexity of Gaussian states in free scalar field theory. While the above discussion phrased the construction of U T as a string of discrete gates, Nielsen's approach begins by introducing a continuous parametrization of the unitary operators, namely where the Hermitian operators K I , which appeared in the unitary gates above, now form a basis for the interactions appearing in the "time"-dependant Hamiltonian H(s). The pathordering symbol P indicates that the circuit is constructed from right to left as s increases.
Hence we can think of the control functions Y I (s) as indicating which gates are inserted at a given point in the circuit, i.e., at a given "time" s. One can then extend this construction such that the circuit (4) represents a trajectory through the space of unitaries, From this perspective, Y I (σ) becomes the tangent vector to the trajectory U (σ), with Our problem is then to find the "shortest" path which starts at U (σ=0) = 1, and ends at U (σ = 1) = U T , where the latter effects the desired transformation in eq. (3). 1 Of course, there are still (infinitely) many paths which will produce the desired unitary U T , and so the question remains how to identify the optimal trajectory. Nielsen's approach [70][71][72] is to define the circuit depth where the cost function F (U (s), Y I (s)) is a local functional of the position U (s) and the tangent vector Y I (s). This functional must satisfy a number of conditions: Note that this last condition is the requirement for the cost function to be an asymmetric norm. 2 These conditions still leave us with an enormous freedom in choosing the cost function.
In the following, we will focus on some simple choices, e.g., which represent the L 1 -and L 2 -norms with respect to the basis K I , and satisfy the conditions (1)-(4) above. Another simple set of cost functions introduced in ref. [23] are where R κ > 1, which however do not satisfy the condition (4). In any event, Nielsen's approach reduces the technical problem of identifying the optimal trajectory to the familiar physics problem of studying the motion of a particle along geodesics, albeit perhaps governed by an unusual Lagrangian.
To make this problem tractable, one typically focuses on a limited basis of operators K I to construct the unitary circuit (4). Of course, this limits the family of target states that can be constructed from a given reference state |ψ R in eq. (3), but it admits a powerful group-theoretic structure if these operators form a closed algebra, i.e., a Lie algebra g with [ K I , K J ] = if I,J K K K . Recent studies of the complexity in free quantum field theories (see, e.g., refs. [23,55,56,58,59]) focused on circuits which remain within the space of Gaussian states, and hence the group of Bogoliubov transformations played an important role. For example, a GL(N, R) algebra appeared in the construction of the free scalar ground state using a lattice of N bosonic degrees of freedom in [23], while the analogous group structure for free (non-interacting) fermions was O(2N ) in refs. [55,56]. The group-theoretic perspective proves to be quite powerful in evaluating the circuit complexity as seen in these examples.
One key benefit of this perspective is that the specific physical details of the operators K I become unimportant once the underlying group structure is identified. Rather, we can think of these generators as the elements of the corresponding Lie algebra g, and of the circuits (5) as trajectories on the corresponding group manifold G. In particular, we are free to choose whichever representation is most convenient for our calculations. In the present paper, we again consider bosonic Gaussian states, but it will be necessary to expand the approach of [23] to the full group of Bogoliubov transformations, namely the symplectic group Sp(2N, R). Additionally, it will turn out to be very useful to characterize the Gaussian states by their covariance matrices G a,b , as explained in section 3, whereupon the unitary operators U are represented as symplectic matrices U a b acting on these covariance matrices. Similarly, in this representation, we can construct matrix generators K a b associated to the operators K. 3

TFD state for the simple harmonic oscillator
As motivated in the introduction, we wish to study the complexity of the TFD state (1) by applying the notion of circuit complexity for Gaussian states introduced in ref. [23]. The complexity of states in QFTs is in general divergent, due to the need to introduce correlations up to arbitrarily short length scales when building the states. In order to study the complexity of the ground state of a free bosonic QFT, ref. [23] regulated the theory by discretizing it on a spatial lattice. The theory then takes the form of a set of coupled simple harmonic oscillators. We will follow the same approach here, and begin by considering the TFD state (1) constructed from two simple harmonic oscillators. We will later explain how this allows us to evaluate the complexity of the TFD formed from two copies of a free bosonic field theory in section 5.
Let us introduce some basic notation before we proceed. The Hamiltonian for a single oscillator is given by 4 where M is the mass of the oscillator, ω is its frequency, and Q and P are the position and momentum operators satisfying the standard commutation relations. In terms of the canonical annihilation and creation operators, which satisfy [a, a † ] = 1, the Hamiltonian (10) becomes where N ≡ a † a is the number operator. In the following, it will be useful to invert (11) via The (normalized) energy eigenstates are given by acting on the vacuum |0 with creation operators, where a † |n = √ n + 1|n + 1 , a|n = √ n|n − 1 , N |n = n|n .
Of course, one then has H |n = ω n + 1 2 |n , as usual.
generators K ∈ sp(2N, R). When we want to refer to them explicitly as quantum operators acting on the Hilbert space H, we will use hats, e.g., U and K. 4 Note that in ref. [23], the mass of the oscillators was set to unity. We avoid doing so here, for reasons that will become clear below.
In the following, we consider two identical copies of a simple harmonic oscillator, which we denote by subscripts L and R in analogy to the left-and right-copies of the CFT in the Penrose diagram in fig. 1. We start by considering the case of the TFD state (1) at t L +t R = 0 and explain how to express it as a quadratic operator acting on the tensor product of the vacua of two harmonic oscillators. We then consider the generalization to the time-dependent TFD.
2.2.1 TFD state at t = 0 The TFD state (1) at t L = 0 = t R can be constructed from two copies of the vacuum state by acting with creation operators in the following manner, e.g., see ref. [76], where we have taken E n = ω(n + 1 2 ). However, since we wish to construct our circuit -and hence the state -from unitary operators, the form (17) is not quite the desired result. Rather, we wish to express the TFD state by acting with a unitary operator on the vacuum |0 L |0 R . Such a rearrangement of eq. (17) can be achieved, as explained in appendix B, using the decomposition formula (252), which yields where we have defined tanh α := exp(−βω/2) .
For later purposes, it is convenient to express α in the form where of course α > 0. Note that the normalization Z −1/2 β = 1 − e −βω 1/2 in eq. (18) has been absorbed into the exponential. Now, observe that by using eq. (11), we can re-express the generator in eq. (18) as Thus we see that, in the language of ref. [23], the TFD state (18) can be interpreted as the result of acting with an entangling operator on the product state of the vacua of the two oscillators.
As we will explain below, a key feature of our approach is that the TFD factorizes in a particular basis, which we refer to as the diagonal basis. 5 One can see this by expressing the generator (21) in terms of the diagonal coordinates for the combined system In other words, the generator can be re-expressed in terms of the scaling operators of the individual diagonal modes. The advantage now is that these two components commute, and hence eq. (18) factorizes as where |0 ± denotes the vacuum of the Hamiltonian of each diagonal mode. Of course, these Hamiltonians look the same as in eq. (10), cf. (25) below. Note that we have rearranged the expression in eq. (24) in such a way that the operator in each factor is unitary. At this point, the set-up for the time-independent (t = 0) TFD is essentially complete. To facilitate our circuit computations below however, let us examine the corresponding wavefunctions that serve as our reference and target states. We began with two independent harmonic oscillators, and hence the total Hamiltonian is given by (cf. (10)) where in the second equality, we have used the diagonal basis (22). The ground-state wavefunction for this Hamiltonian takes the simple form We might characterize the analysis in refs. [23,51] as defining the complexity of this ground state given the reference state 6 We can think of this reference state as being the ground state of a Hamiltonian as eq. (25), but where the frequency is fixed to be µ, the reference frequency of our complexity model -compare with eq. (182) in our analysis of the scalar field theory. In the present work, we extend these results to define the complexity of the TFD state relative to the same reference state |ψ R . Examination of eq. (24) indicates that the wavefunction for the TFD is obtained by acting on the vacuum wavefunction (26) with the appropriate scaling operators in the diagonal basis. For example, using the results of ref. [23], we have Therefore the desired wavefunction is simply given by Thus, the target state is indeed a simple Gaussian wavefunction of the same form studied in refs. [23,51].

Time-dependent TFD state
It is relatively straightforward to extend the manipulations of the previous subsection to the full time-dependent case. For simplicity, we shall follow the common convention in holography (see, e.g., refs. [20,28]) and set t L = t R = t/2 in eq. (1). Then the generalization of eq. (17) becomes Note that in what follows we will simply drop the global time-dependent phase given by the e − i 2 ωt factor, since this does not change the physical state. 7 As above, we wish to express this state in terms of a unitary operator acting on the vacuum state. As explained in appendix B, using eq. (257), one finds where we have defined z = α e −iωt (32) and α is given by eq. (19). Note that this reduces to eq. (18) upon setting t = 0. As before, we use eq. (13) to re-express the generator in eq. (31) as At this point, we make the important observation that producing the time-dependent TFD requires moving beyond the scaling and entangling operators of the gl(2, R) algebra considered 7 One could consider adding a contribution to the complexity to account for this phase using a phase gate as in [23]. This contribution to the complexity would just be some weight or penalty factor times a simple function of the phase, i.e., one would take the absolute value of the portion of the phase between −π and π. However, when considering the vacuum state, one would acquire an identical phase contribution to the complexity as it evolves in time (recall we are evolving with tL = tR = t/2). Note that no such contribution arises in holographic complexity. Furthermore, if we compare the complexity of the TFD to that of the vacuum at general times, this contribution would cancel out. in ref. [23]. Rather, we must build our unitaries using generators from the full algebra formed by all possible bilinears of Q L , Q R , P L , P R , which form the algebra sp(4, R) [77]. As we demonstrate in the following section, such unitaries describe the set of transformations of two-mode Gaussian states with vanishing one-point functions among themselves. In fact, we will be mostly using the diagonal basis and the corresponding sp(2, R) subalgebras for each diagonal mode formed by bilinears of P + , Q + or P − , Q − respectively.
Rotating to the diagonal basis (22) as above, the generator (33) becomes where we see that as in the time-independent case, the generator can be decomposed into two (anti-Hermitian) operators, which act separately in the '+' or '-' Hilbert spaces. Therefore the time-dependent TFD continues to factorizes as This state can also be expressed as a Gaussian wavefunction of the general form (29), and serves as our target state for the time-dependent case. In this case the explicit expression for the wavefunction is somewhat complicated, so we will not write it out here. Rather, below we will introduce a more elegant formalism for characterizing these Gaussian states by their covariance matrices, which then allows us to represent the symplectic transformations of sp(4, R) (or of the sp(2, R) subalgebras) acting on them as matrix operations.

Gate scale
Looking ahead, we note that in constructing the quantum circuits (4), our generators, i.e., the K I , will consist of all quadratic combinations of the Q's and P 's -see section 3.2. However, since as usual these operators are dimensionful, we will need to include an additional scale in our complexity model. The simplest approach, which we follow here, is to introduce the dimensionless position and momenta where we refer to ω g as the gate scale. In particular, we can then construct dimensionless generators K I as all quadratic combinations of these dimensionless q's and p's. For example, we might think of some typical elementary gates represented as where ε is a dimensionless infinitesimal parameter. 8 Hence the gate scale implicitly ensures that all of the components of the control functions Y I (s) are dimensionless and are readily combined in cost functions, such as those given in eqs. (8) and (9). In terms of these dimensionless variables, the reference (27), vacuum (26), and timeindependent TFD (29) states become, respectively, where we have defined the dimensionless ratios 3 Covariance matrix approach So far we have described Gaussian states of bosonic systems in terms of their associated wavefunctions. It turns out, however, that this representation is more complicated than necessary; in particular, as alluded above, the wavefunction for the time-dependent TFD is rather unwieldy. Fortunately, an equivalent and simpler representation is available at the level of covariance matrices, which greatly facilitates our analysis. In this section, we review the relevant aspects of this approach.

From quantum states to covariance matrices
A bosonic system with N degrees of freedom can be described by 2N linear observables ξ = (q 1 , · · · , q N , p 1 , · · · , p N ) consisting of canonical coordinates. Given an arbitrary quantum state |ψ , its two-point function may be expressed as where G a,b = G (a,b) is symmetric and Ω a,b = Ω [a,b] is antisymmetric. Such a decomposition can always be performed for an arbitrary matrix, but we included an i in front of Ω a,b to ensure that both G a,b and Ω a,b are real linear forms. This follows directly from the fact that ξ a ξ b + ξ b ξ a is a Hermitian operator with real eigenvalues, while ξ a ξ b − ξ b ξ a is anti-Hermitian with purely imaginary eigenvalues. In fact, the latter is completely fixed by the canonical commutation relations to with respect to the coordinates ξ introduced above. If |ψ is a pure Gaussian state with ψ|ξ a |ψ = 0, then it is completely characterized by its symmetric two-point function, often referred to as its (symmetric) covariance matrix G with entries For a mixed state ρ with vanishing first moments, one can in the same way define [78,79] In particular, we can use Wick's theorem to compute any higher order n-point functions from G a,b . Hence we can without ambiguity use G as label for these Gaussian states. For example, let us consider a general pure Gaussian state in the Hilbert space of a single degree of freedom. Such a state can be parameterized by its wavefunction as where a, b ∈ R, and a > 0 for normalizability. The two-point function (43) may be explicitly evaluated, and the results encapsulated in the covariance matrix One sees that it is straightforward to extract the parameters a and b of the wavefunction from the covariance matrix G: a = 1/G 1,1 and b = −G 2,1 /G 1,1 . 9 This illustrates our claim above that the covariance matrix provides a complete characterization of Gaussian states, and can therefore be used as an alternative description thereof. In particular, the wavefunctions (39) each decouple into a product of wavefunctions for the ± modes, and hence the associated covariance matrices for the total states are block-diagonal. The covariance matrices for the + mode read and similarly for the − mode, with α → −α.

Trajectories between states and their generators
The power of the covariance matrix formalism lies in the fact that we can study trajectories in state space purely in terms of G a,b , provided that we do not leave the sector of Gaussian states. Restricting to this class of states can be easily achieved by focusing on a natural subgroup of unitary transformations which evolves Gaussian states into Gaussian states, i.e., the Bogoliubov transformations. This subgroup is generated by Hermitian operators that are quadratic in the canonical coordinates ξ introduced above. The most general quadratic operator can be written as where k = k (a,b) is chosen to be a real symmetric form. This is because any antisymmetric part in k does not affect K due to ξ a ξ b being symmetric. 10 Such a general operator K generates unitaries that map Gaussian states into Gaussian states [80,81]. To find the covariance matrix associated with the new state |G σ , we start by computing the operation of U (σ) on ξ a given by where we have used the well-known Baker-Campbell-Hausdorff formula. Using the commuta- where we have defined the matrix generator associated with K, One can check that K ∈ sp(2N, R), and that U (σ) = e σK belongs to the symplectic group Sp(2N, R), namely 11 Satisfaction of the last condition is most transparent when expressed in terms of the algebra, i.e., which can be verified by the use of eq. (52) and the fact that Ω T = −Ω.
The symplectic group Sp(2N, R) is the group of elements satisfying the relations (53), which amounts to linear transformations on the canonical variables ξ a that preserve the canonical commutation relations; sp(2N, R) is the associated algebra of generators. We can express the operation of U (σ) on ξ a as 10 Note that K is a positive operator if and only if k is a positive definite bilinear form. In this case it can be viewed as a physical Hamiltonian which is bounded from below. 11 Note that we have dropped the hats in moving from the operators to the matrix representation, and that despite the suggestive notation, U (σ) is not a unitary matrix.
U reflects the so-called metaplectic representation of U [77]. With the above in hand, the covariance matrix of |G σ may be computed as Note that the transformation of the covariance matrix is linear, i.e., each of the components of G σ is a linear combination of the entries in G 0 . This allows for the very compact notation Of course, we can express any Gaussian state |G also as a Gaussian wavefunction of the form with q := (q 1 , · · · , q N ), where A and B are real bilinear forms that can be computed from G as explained in ref. [81]. However, the action of U (σ) on these bilinear forms -namely A(σ) and B(σ) for the sequence of states |G σ -is much more cumbersome than the simple expression (57) above. Only if one enforces B = 0 with respect to a specific splitting of the classical phase space into positions q i and their conjugate momenta p i does one find that A(σ) transforms in a simple way under the subgroup GL(2N, R). This was the case studied in ref. [23], but if we want to extend our analysis to the full symplectic group (as required for the TFD state), then it is much more convenient to label Gaussian states by their covariance matrices rather than by the parameters A and B in the wavefunction. Let us consider the special case N = 1 in order to gain some intuition for the generators of the corresponding group Sp(2, R). In this case, eq. (48) yields only three independent generators, which we denote These close to form the algebra sp(2, R), with commutation relations Note that these are also the generators that serve as building blocks for the operators in eq. (35). The associated matrices k (a,b) in eq. (48) are given by We can then use eq. (52) to obtain the relevant matrix generators which satisfy Finally, exponentiating these generators yields the group elements that will serve as the elementary gates used in the construction of our quantum circuits below: where is a real parameter with | | 1. To see how these Sp(2, R) gates modify the state, we evaluate the change in the covariance matrix effected by these gates via eq. (57). Suppose we start with the generic pure Gaussian state |ψ given in eq. (45). Then, denoting the state after acting with the U Z gate byψ Z , i.e., |ψ Z = e −i Z |ψ , one finds 12 The parametersã,b of the transformed wavefunctionψ Z are therefore Similarly, the changes effected by the other two gates are These last two expressions are relatively simple; and indeed, in U W we recognize the action of a scaling gate as in ref. [23]. 13 In contrast, one sees from (66) that, while the state remains Gaussian under the action of U Z , this gate produces a nonlinear change in the parameters of the wavefunction, 14 in contrast to the simple transformation of the covariance matrix in eq. (56). At a practical level, the fact that this action can be deduced straightforwardly from the change in the two-point function is the main advantage of working with the covariance matrix. Now, returning to the quantum circuits introduced in eq. (5), we can use the covariance matrix language to replace the circuits U (σ) by their matrix representation U (σ). In particular, we can ask how a state changes under the evolution U (σ) = Pe −i σ 0 K(s)ds of a varying quadratic operator 12 Alternatively, one may use the Baker-Campbell-Hausdorff formula directly to obtain the same relations. 13 This was also referred to as a squeezing gate in ref. [51], since it shrinks the variance of the position operator q at the expense of increasing the variance of the momentum operator p, while keeping the expectation value of the cross-product fixed.
14 In fact, the action in eq. (66) is somewhat akin to a special conformal transformation in conformallyinvariant systems, insofar as it can be obtained from an inversion a = a/ a 2 + b 2 and b = b/ a 2 + b 2 , followed by a translation in the b-direction (i.e., applying UV ), a = a and b = b − 2 , and finally by another inversion a = a / a 2 + b 2 and b = b / a 2 + b 2 .
We can use the same arguments as before to find where the matrix representation of the path U (σ) satisfies the equation whose solution is given by the path-ordered exponential that acts on the covariance matrix (or the state) as in eq. (57).

Covariance matrix for the time-dependent TFD state
In this section, we demonstrate how the above machinery may be employed to evaluate the covariance matrix for the time-dependent thermofield double state. We saw in eqs. (35)- (36) that this state can be written as a tensor product decomposition in the diagonal basis. Hence we focus our attention on only one of these Gaussian factors, e.g., the state formed by acting with O + (t) in eq. (36). The Gaussian state formed by acting with O − (t) is obtained simply by replacing α → −α.
The most straightforward way to obtain the covariance matrix of the time-dependent TFD state is by evolving the covariance matrix of the TFD state at the t = 0 in eq. (47) forward in time. More concretely, we would like to apply to the latter state the unitary where we have used eqs. (1), (25), (37), and (40), such that the state whose covariance matrix we wish to obtain is given by This problem falls precisely within the formalism of the last subsection, where the evolution is given by the operator Translating this to the level of matrix operators, (52) and (53) become whereupon the action of U (t) on the covariance matrix (57) is found to be G + TFD (t) = U (t) G + TFD U (t) = 1 λ (cosh 2α + sinh 2α cos ωt) − sinh 2α sin ωt − sinh 2α sin ωt λ (cosh 2α − sinh 2α cos ωt) = cosh 2 α 1 λ 1 + 2 tanh α cos ωt + tanh 2 α −2 tanh α sin ωt −2 tanh α sin ωt λ 1 − 2 tanh α cos ωt + tanh 2 α , where G + TFD was given in (47). In the second line, we have presented the result in a way that is easily related to the physical variables via (19), i.e., tanh α = exp(−βω/2) and cosh α = (1 − e −βω ) −1/2 . The time-dependence of this expression is of course periodic. Note also that we recover G TFD from eq. (47) upon setting t = 0, as expected.
Alternatively, one may obtain this covariance matrix by using the operation of O + (t) on the vacuum state according to eqs. (35)- (36), i.e., |ψ(t) : and where we have rewritten the generator O + (t) in eq. (35) in terms of the rescaled variables q, p defined in eq. (37) and the parameter λ which was defined in eq. (40). To simplify the notation, we drop the + subscript in the following. The relevant matrix operator that obtains the TFD at time t from the vacuum state, which we denote U TFD (t), is again obtained according to eqs. (52), (53), and (57): We can then obtain the covariance matrix of the time-dependent TFD state by acting on the vacuum covariance matrix G 0 in eq. (47) with U TFD (t) according to eq. (57), i.e., which of course reproduces eq. (76) above.
To close this section, let us also give the transformation U vac that brings the reference state G R to the vacuum state G 0 (also given in eq. (47)), since we will need this in the following sections. This transformation was studied in ref. [23], and is related to the following quantum operator (again focusing on the + mode and dropping the subscripts): Following the same steps as above, one finds and one can readily verify that

Relative covariance matrices and stabilizer group
The stabilizer subgroup Sta G ⊂ Sp(2N, R) associated to a Gaussian state G is defined as The importance of the stabilizer group lies in the fact that it allows one to relate different unitaries that map a given reference state to the same target state. Explicitly, if we have a matrix U satisfying then for any U R ∈ Sta G R the operator U U R will also obtain the same state G T . Thus when minimizing over circuits to compute the complexity, we must also minimize over this family of transformations. As a group, we have Sta G U(N ), but different choices of G will lead to different embeddings of this subgroup within Sp(2N, R). 15 For the case N = 1, we consider the stabilizer group which leaves invariant the reference state G R in eq. (47). This is an SO(2) U(1) subgroup of Sp(2, R) of the form (85) In this case, if U T achieves the desired transformation between the reference and target states, then we must minimize over the family of circuits given by the transformation U T U φ (subject to the same boundary conditions) for all values of the rotation angle φ.
If we are interested in the relation of two Gaussian states G and G, then we can express the relation between them in a basis independent way in terms of the relative covariance matrix where g is the inverse of G, such that G a,c g c,b = δ a b . Any quantity which is invariant under the Sp(2N, R) group is necessarily a pure function of the spectrum of ∆. 16 To show this, we may first use the full group Sp(2N, R) to choose a basis, such that the matrix representation of G becomes the identity, i.e., G = 1. We can then diagonalize the covariance matrix G using only transformations within the stabilizer subgroup Sta G , which provides the freedom to change basis without affecting G. This is equivalent to finding the spectrum of the relative covariance matrix. An example of such an invariant function is the inner product | G| G |, which can be computed as [81] | G| G | 2 = det √ In the case of complexity, we can make a choice for the cost function that is defined in terms of the reference state G R . 17 As this implies that the complexity only depends on G R and G T , we will find a simple formula for the F 2 complexity in terms of ∆, namely 15 In fact, up to a deformation, the elements of the group are nothing but the elements of the passive subgroup Sp(2N, R) ∩ O(2N ). This, in turn, is isomorphic to U (N ), as this subgroup reflects unitary transformations of vectors of bosonic operators. 16 The eigenvalues of ∆ come in pairs where each eigenvalue is accompanied by its inverse. We refer to the set containing the first of each pair as the spectrum. 17 This is not actually the choice that we make in our complexity calculations, see discussion in section 4.1.
Generally, we distinguish the reference state (47) from the state G = 1 appearing in the definition of the cost functions, e.g., see eq. (108). Therefore, the complexity expression in eq. (88) only applies for the choice λR = 1, as we make clear in appendix F. or for the κ = 2 complexity, we have C κ=2 (G R , G T ) = [C 2 (G R , G T )] 2 = 1 8 Tr[(log ∆) 2 ]. Both of these expressions are derived in appendix F -see eqs. (322) and (323). However, we will also consider other cost functions that explicitly depend on a choice of basis, in which case knowledge of ∆ does not suffice to compute the complexity.
To make the above more concrete, let us write the relative covariance matrix between the reference state G R and the time-independent TFD state G + TFD given in eq. (47): Hence in this case, any cost function (i.e., definition of complexity) which is invariant under the Sp(2N, R) group must only depend on the combination e 2(α R +α) = λ R λ e 2α . For example, the inner product is which indeed depends only on the spectrum of ∆.

Generators of Sp(2N, R)
In examining the complexity, we replace the circuits (5) with their matrix-valued counterparts (71), but we will still need to decompose the exponential in terms of some basis in order to evaluate the appropriate cost function, e.g., in eq. (8) or (9). Hence we will find it useful to have explicit expressions for the generators of the group Sp(2N, R), in particular for N = 1 and N = 2. Accordingly, here we give a list of the generators of Sp(2N, R) for general values of N . We may start with the quantum generators, i, j, k ∈ (1, . . . , N ) We note that the number of W i,j , V i,j , and Z i,j operators is N 2 , 1 2 (N 2 + N ), and 1 2 (N 2 + N ), respectively, for a total of N (2N + 1) generators. In these expressions, ω g is the gate scale introduced in eq. (37) in order to render the generators V i,j and Z i,j dimensionless. This scale does not enter W i,j , since it is invariant under a rescaling of Q i and P j . The generators W i,j span the subalgebra gl(N, R), which was analyzed in ref. [23], and hence the gate scale did not enter into the complexity calculations considered therein. However, as shown above, the preparation of the time-evolve TFD states will also involve the V i,j and Z i,j generators, and so a priori the complexity of these states will depend on the choice of the gate scale.
The above generators have been chosen to be orthonormal according to the Frobenius inner product where g a,b is the inverse of G a,b as above. We will choose G = 1 in the basis spanned by our dimensionless variables (q i , p j ). In particular, this normalization is responsible for the extra factors of 1/ √ 2 appearing in eq. (91) for the diagonal generators. This inner product will play an important role in what follows.
We can translate these generators into the relevant matrix representation by using eq. (48) and eq. (52), which leads to Following eq. (52), the associated matrix generators are obtained by multiplying with the symplectic form which yields For N = 1, these expressions reproduce the generators W, V, Z in eq. (62). For the purposes of this paper, we will mainly use the generators of Sp(4, R) which we list explicitly in appendix C. When computing the complexity of a particular circuit described by eq. (71), we may need to expand a given generator K(s) with respect to different bases of generators, say K I and K I . Provided that these bases are both orthonormal with respect to the Frobenius inner product (92), i.e., we can accomplish this by computing Y I (s) = K I |K(s) and Y I (s) = K I |K(s) .
In the following, we will work with two bases that have already appeared in section 2.2. In particular, we have the L, R basis referring to the two copies of the physical degrees of freedom entangled in the TFD state, and the diagonal or ± basis in which the TFD state factorizes.
As indicated in eq. (22), these two bases are related by a simple rotation 18 but it will be useful to systematize the transformation for general expressions. For example, eq. (98) extends to an analogous equation for the momenta, and so the transformation on the full phase space reads The two-point function (43) then transforms as It is now straightforward to see that if we have a circuit acting in the diagonal basis as G T = U G R U , then the same circuit in the physical basis is Furthermore, the matrix generators K I in eq. (95) are also simply transformed to 19 For example we may use these transformation rules to transform the covariance matrix from the (q + , q − , p + , p − ) basis to the (q L , q R , p L , p R ) basis. We start from the covariance matrix in the (q + , q − , p + , p − ) basis where the direct sum inputs the + and minus components in a 4 by 4 combined matrix and where the G + TFD (t) was defined in eq. (76) and G − TFD (t) is a similar matrix with α → −α. The time dependent covariance matrix with respect to the (q L , q R , p L , p R ) basis is given according to the rotation (100) or explicitly This expression will come handy later on in section 6. 18 We use the subscript 2 here to indicate that this is a 2×2 matrix and distinguish this rotation matrix from R4 below. Furthermore, while as a numerical matrix R2 is symmetric, we nonetheless distinguish R2 and R 2 in the following to emphasize the fact that R2 provides a mapping from the physical to the diagonal coordinates, while R −1 2 = R 2 provides the inverse mapping. In other words, the columns of R2 are labeled L, R while the rows are labeled +, −, and vice versa for R . 19 The transformation (99) is special since it is orthogonal. Generally, such coordinate transformations on the phase space have a similar effect, except that eqs. (101) and (102) are replaced with U = R −1 U R and K = R −1 KR, respectively.

Complexity of TFD states
In this section, we apply the tools developed above to study the complexity of the TFD state comprised of two harmonic oscillators. We will later use the results of this section to evaluate the complexity of the TFD state of a free scalar field theory in section 5.

Circuit geometry for the TFD state
In subsection 2.1, we introduced the definition of complexity for general groups. In the present subsection, we specialize to the group G = Sp(2N, R), and the associated algebra g = sp(2N, R) of matrices acting on the covariance matrix. We will explain how the general definitions presented in section 2.1 can be applied in this specific case. In general, our matrixvalued circuits will take the form given in eq. (71), i.e., where K I are the generators of the algebra sp(2N, R) in a given basis. We will start with the covariance matrix associated with our reference state G R , and follow a path in the space of Gaussian states represented by the covariance matrices We will be interested in trajectories that end on a given Gaussian target state G T , i.e., our circuits satisfy the boundary conditions The length of a given circuit will be given by integrating certain cost functions along the path, as we have discussed in subsection 2.1. We introduced some examples with the F 1 , F 2 , and D κ cost functions in eqs. (8) and (9) above. Note however that we still have an enormous amount of freedom, since different choices of basis vectors K I will in general lead to different results for the total cost [23]. The F 2 cost function, as well as D κ=2 , is invariant when the two bases are related by an orthogonal transformation. 20 However, even these cost functions are implicitly defined in terms of a particular reference state [82]. 21 We can also view the F 2 cost function as arising from a natural construction using a positive-definite matrix G a,b (for instance, by taking it from the covariance matrix G of a Gaussian state), namely where again g a,b denotes the inverse of G a,b , i.e., G a,c g c,b = δ a b . This cost function coincides with the norm induced by the inner product (92). Extending this inner product to the full 20 For instance, the transformation between the qL,R basis and the diagonal basis q± in eq. (22) is such an orthogonal transformation, and hence the F2 and Dκ=2 cost functions will be invariant under this change. 21 This state-dependence does not occur for fermions [56].
group turns Sp(2N, R) into a Riemannian manifold, whose metric at the point U ∈ Sp(2N, R) can be computed as If we choose G = 1, this metric reproduces that given in eq. (8). This is the choice that we will make in the following. Recall from eq. (97) that we can use the inner product to extract the coefficients Y I in a given generator K = Y I K I with respect to an orthonormal basis K I by simply evaluating If we are interested in the circuit complexity defined with respect to a given reference state G R , then a great simplification occurs when the matrix G used to define the geometry above and the covariance matrix G R of the reference state coincide. This is the case when we choose G = 1, and the gate scale ω g is chosen to be equal to the characteristic scale of the reference state √ M µ, equivalently when setting λ R = 1, cf. eq. (40). In this case the covariance matrix of the reference state is simply the identity, and is hence equal to the matrix G used in defining the geometry above. For the cost function F 1 , we also have considerable freedom in choosing the basis of generators K I . We will impose that our generators be orthonormal under the inner product inducing the F 2 cost function. Even then we retain quite a bit of freedom, e.g., the rotation between the LR basis and the ± basis in the group Sp(4, R). This change of basis does not affect the F 2 cost function, but it does affect F 1 .
Let us now consider the case of a single degree of freedom. In this case, we have the group Sp(2, R) = SL(2, R), whose algebra is given by the traceless matrices K I ∈ {W, V, Z} given in eq. (62), where we have chosen the coordinates ξ = (q, p) (which will later represent the conjugate pairs (q ± , p ± ) that mix the left and right sides of the TFD). With respect to these coordinates, our spatially unentangled reference state G R is given in eq. (47). To obtain the complexity of a given target state G T , we will then study geodesics in the geometry (109) which satisfy the boundary conditions (106). More precisely, we will have to minimize over the family of geodesics that end at U (σ = 1) = U T U φ , cf. eq. (85), all of which transform the reference state to the same target state.
The relevant geodesics for this metric were already worked out in ref. [23]. In particular, if we denote the geometric coordinates x(σ) = (ρ(σ), θ(σ) , τ (σ)), then the initial condition U (σ = 0) = 1 in eq. (107) fixes cf. (111), where θ 0 := θ(σ = 0). The freedom in specifying the initial angle θ 0 is simply the freedom to specify the initial direction in which the geodesic moves away from the origin. Now, denote the coordinates at the endpoint of the geodesic by These will be fixed in terms of the physical quantities using eq. (107) when we consider specific cases below. Of course, there will still be some residual freedom from the stabilizer group of the reference state as discussed in subsection 3.4, in particular around eq. (85). The geodesics with these general boundary conditions are Eq. (114) then allows us to extract the values of (θ 0 , ∆θ, c) in terms of the boundary condition (ρ 1 , θ 1 , τ 1 ); in general, this inversion has to be done numerically. The geodesics in eq. (115) are affinely parameterized such that the line element is constant and equal to Given a reference state G R and a target state G T , we wish to find the minimal geodesic in the group that takes us from 1 to an element of the family In order to perform this minimization, we need to understand this family F R→T for different target states, so that we can use them as boundary condition for our solutions of the geodesic equation. In particular, we will be interested in the family of minimal geodesics starting from the reference state G R in eq. (47). Recall from subsection 3.4 that the stabilizer subgroup Sta G R associated with the reference state G R is given in terms of U φ in eq. (85). Given a group element U T that prepares the target state G T = U T G R U T , we find In order to better understand the geometry of F R→T for different U T , we can plot a selection of them in fig. 3 for various values of λ R . Minimizing the geodesic distance over all end points U T (φ) := U T U φ ∈ F R→T is a difficult task for general λ R , and we will discuss it further in subsection 4.5 below.
As we have already mentioned, a significant simplification occurs for the special case of λ R = 1, in which the reference and gate scales are equivalent, √ M µ = ω g . The reference-state covariance matrix becomes G R = 1, and the transformation U (ρ, θ, τ ) takes us to the target state with covariance matrix We immediately observe that this expression only depends on the two coordinates θ and τ through the combination θ + τ , leading to a one-parameter family of solutions. This is precisely the U(1) invariance of the stabilizer group of the reference state (note that in this case the matrix U φ in eq. (85) is a simple rotation matrix). Let us define χ := θ + τ , so that we can label the covariance matrix G T (ρ, χ). Then the equivalence class of group elements that prepare the same state G T (ρ, χ) is a spiral, which can be parameterized by τ with We illustrate this equivalence class in fig. 3. As alluded above, obtaining the explicit parameters (c, ∆θ, θ 0 ) for a general boundary condition (ρ 1 , θ 1 , τ 1 ) at σ = 1 is difficult, but for ∆θ = 0 the expressions simplify to This means that the geodesics in the plane with τ = 0 are just given by straight lines. The line element associated with the trajectory is given by eq. (116), which for this specific case reads ds = ρ 1 dσ, Ref. [23] showed that in certain cases the optimal trajectory is indeed the one associated with ∆θ = 0. The analysis is based on a series of inequalities which appear in eqs. (3.39)-(3.44) therein. For completeness, we give a sketch of the derivation here: first, recall that evaluating the norm of the velocity along the geodesics using the metric (112) yields a constant, i.e., |∂ σ x| 2 = k 2 . Using the geodesic solution (115) to evaluate ∂ σ θ and ∂ σ τ , and averaging the resulting expression for k 2 over the geodesic, we can then show that We may now use the two inequalities to prove the geodesic inequality It is then obvious that in cases where ρ 1 is a constant independent of ∆θ (as in ref. [23]), the geodesic with the minimal value of k 2 indeed has ∆θ = 0. In the present case, with the F 2 or D κ=2 cost functions, we have simply F 2 = k or D κ=2 = k 2 , and so in both cases minimizing k 2 corresponds to minimizing the corresponding circuit depth. Therefore we conclude that the optimal geodesic also corresponds to that with ∆θ = 0. It is straightforward to extract the value of ρ 1 (φ) associated to the end point of the trajectory U U φ -see eq. (168) below. In the particular case when λ R = 1, the matrix U φ in eq. (85) is simply a rotation matrix, and we will find that the value of ρ 1 (φ) does not depend on φ, i.e., ρ 1 becomes a constant with λ R = 1. We are therefore able to apply the above result and conclude that for λ R = 1, the straight-line geodesics with ∆θ = 0 describe the optimal circuits (where again we have chosen either the F 2 or D κ=2 cost functions).
In appendix F, we prove a similar result for Sp(2N, R) for general N , again for the special case that λ R = 1. In particular, we demonstrate that the optimal circuit that prepares the target state G T from the reference state G R whose covariance matrix is the identity is given by the straight-line geodesic, generated by a single generator determined by the relative covariance matrix introduced in eq. (86), i.e., The proof of this general result requires some Lie group techniques together with a well-known decomposition of group elements of Sp(2N, R). In particular, this Lie group can be represented as a U(N ) fiber bundle over the symmetric space Sp(2N, R)/U(N ). Here the fiber is nothing but the stabilizer group (83), and the base manifold can be interpreted as the space of Gaussian quantum states. We build on this decomposition to produce a generalized cylindrical foliation of the group manifold of the form e A u with u ∈ U(N ), where A plays the role of the radius. One can show that all geodesics that prepare the desired state will end on the cylinder of radius K . The final step is then to show that the minimal geodesic connecting 1 with this cylinder is the one which moves in a purely radial direction, i.e., it is the geodesic given in eq. (126), as we saw in the previous discussion of the special case N = 1.
It now remains to find the geodesics that produce the particular target states of interest, namely the TFD states introduced above. We shall first consider the special case of the timeindependent TFD with t = 0 given by eq. (47). This case is relatively straightforward, and also enables us to make contact with the holographic results on the complexity of formation [27]. We shall then move on to the full time-dependent TFD (76) in the subsequent subsection.

Complexity of the TFD at
In this subsection, we focus on the complexity of one of the diagonal modes comprising the TFD at t = 0. Specifically, the covariance matrix G + TFD in eq. (47) for the mode associated with the x + coordinate will serve as our target state. Our reference state is given by G R in eq. (47). As mentioned above, setting λ R = 1 provides a significant simplification, so we will focus on this case first. For convenience, we restate the relevant covariance matrices here: Using eq. (107) together with eq. (119), we obtain the boundary conditions for our circuit 1 λ e 2α = cosh 2ρ 1 − sinh 2ρ 1 sin(θ 1 + τ 1 ) , λ e −2α = cosh 2ρ 1 + sinh 2ρ 1 sin(θ 1 + τ 1 ) , The determinant of G T in eq. (119) is one, and so the above equations only represent two independent relations for θ 1 + τ 1 and ρ 1 in terms of the physical parameters of the problem (i.e., e −2α λ). Of course, this leaves the linear combination θ 1 − τ 1 unfixed. This remaining freedom is due to the equivalence class of circuits which produce the same target state, as discussed in the previous subsection. Additionally, as explained above, when ρ 1 is a fixed constant, the optimal geodesic is the straight line moving at a fixed angle in the τ = 0 plane, as given in eq. (121). Hence we have τ 1 = 0, whereupon solving eq. (165) yields Substituting into eq. (111) and using eq. (121), we thus obtain the optimal circuit where W + is the generator for the scaling gate acting on (x + , p + ) in eq. (62). Repeating the analysis or the x − mode, 22 we simply get U − (σ) = exp − 1 2 log λ + α σ W − , where W − is the scaling generator of the − modes. The simple form of these circuits, each containing only a single generator, allows us to easily compute the complexity C according to the cost functions in eq. (8). Since the two circuits commute, we can think of 1 2 log λ ± α as the number of times each scaling gate was applied, and simply combine these numbers using the chosen norm. The F 2 cost function yields the complexity 23 or, in terms of the physical variables, using eqs. (20) and (40), Of course, the complexity for the κ = 2 cost function is simply related to the above result, i.e., We can also evaluate the length of this circuit with the F 1 cost function, which yields where we have included the superscript ± to indicate that this complexity was evaluated with respect to the diagonal basis of generators, see eq. (22). Note that we did not prove that this trajectory was optimal for the F 1 cost function. It may be that a proof can be formulated, 22 Recall that the state associated with the x− coordinate is obtained by replacing α → −α in the + state. 23 Alternatively, we can use eq. (109) to read the line element directly, which for a circuit generated by a constant generator reduces to ds 2 = 1 2 tr(K · K )dσ 2 .
but it seems likely that with the F 1 cost function, many different circuits will be assigned the same, minimal circuit depth-see discussion in ref. [23]. However, we have not pursued this possibility in any detail at this point, and hence we may only state that our results for the F 1 cost function are an upper bound on the complexity of the target state. Additionally, we noted above that the F 1 cost function depends on the basis chosen for the generators-again, see discussion in ref. [23]. It is therefore interesting to explore how the result (134) changes when we consider the basis of generators which naturally act on the physical (left (L) and right (R)) modes, rather than the diagonal (±) modes-see eq. (22). For this purpose, we must first combine the two circuits U ± (σ) into a single 4 × 4 matrix (rather than two independent 2 × 2 matrices). The relevant transformation is then given by the combined U (σ) acting on the covariance matrix describing the two ± oscillators. Organizing the diagonal modes asξ = (q + , q − , p + , p − ), the transformation is block-diagonal, and takes the form Organizing the physical modes as ξ = (q L , q R , p L , p R ), the transformation of the circuit to the left-right basis via eq. (102) yields To compute the complexity in this basis, we need to decompose this generator K in terms of the generators of Sp(4, R) defined in eq. (95), see also appendix C. It is straightforward to extract the coefficients of these basis generators by taking the inner product (110) of eq. (136) with each of the generators. In this way, we find that and so only four components of the tangent vector Y I are nonvanishing. Evaluating the complexity with the F 1 cost function (8) then yields which clearly differs from the result in eq. (134) in the diagonal basis. It is also interesting to compare the complexity of the entangled TFD state of the two oscillators with that of the unentangled vacuum state, i.e., the β → ∞ or α → 0 limit of the TFD. The difference between these complexities for the present case of two oscillators serves as a precursor to the complexity of formation for the free scalar field in the next section. This quantity was originally defined in the context of holographic complexity, see, e.g., ref. [27]. For the F 2 and κ = 2 complexities above, one finds Hence we see that there is a more complete cancellation in the complexity with the κ = 2 cost function. In particular, the difference only depends on the combination βω through α from eq. (20), and is independent of the reference state scale µ or the mass M . For the F 1 complexity, we can compare the difference for the diagonal and physical bases, respectively: We note that a similar cancellation arises for ∆C as in ∆C κ=2 , where the result only depends on βω (but not µ or M ).
In the following section, we will apply the above results to evaluate the complexity of formation for a free quantum field theory. These calculations will be based on the fact that the field theory can be represented as a collection of momentum modes, where each momentum mode is essentially entangled with its counterpart in the TFD to form a product of twooscillator TFD states of the form studied in this section. We will find that the there are some interesting similarities between the F 1 result in the physical basis and previous results obtained for holographic complexity [27].
Of course, as in eq. (165), the right-hand side only depends on ρ 1 and θ 1 + τ 1 , while the combination θ 1 − τ 1 is left undetermined. We can solve eq. (141) explicitly to obtain Now, for the minimal straight-line geodesic (121) in the τ = 0 plane, the circuit (111) is given by where in the last expression we have used the matrix generators in eq. (62). Now, for the F 2 cost function (8), we can use eq. (122) to evaluate the complexity, i.e., the length of this optimal circuit, which yields which is independent of θ 1 , as expected from our analysis above. Of course, for the κ = 2 cost function (9), we find C where ρ 1 was given in eq. (142). Alternatively, the F 1 cost function yields 24 As we commented above, we are not assured that the straight-line trajectory (121) is the shortest geodesic for this cost function, but it does at least provide an upper bound on the complexity. Substituting in the boundary conditions ρ 1 and θ 1 given in eq. (142) (with τ 1 = 0) into the above expressions then yields the complexity of |TFD(t) in terms of the physical parameters of the state, i.e., ω and β. The contribution to complexity from the − mode are obtained from the above by simply replacing α → −α in eqs. (141) and (142). Next, we would like to look at the F 1 complexity in the physical LR basis. By construction, Now we apply the transformation (102) to obtain the relevant generator M = R 4 M R 4 in the LR basis. We can then use eq. (97) to extract the coefficients of the SP(4, R) generators in this basis (see appendix C for the full list of generators). We finally obtain the following decomposition: If we measure the complexity with the F 2 cost function (8), combining the two modes, we find 24 We might note that the results here and in the following are simplified somewhat if we replace {W, V, Z} → {W, 1 √ 2 (V ± Z)}. With this new basis, eq. (146) becomes C1 = ρ1(| cos θ1| + | sin θ1|), but the results are qualitatively unchanged.
while the κ = 2 cost function (9) instead yields Both these results are invariant under the orthogonal basis transformation (102), i.e., these complexities are the same in the diagonal and the physical bases, as well as being independent of θ 1± . In contrast, using the F 1 cost function in the physical basis, we arrive at a different result C Translating eq. (142) for the ± boundary conditions then allows us to express the coefficients above in terms of the physical parameters, namely cosh 2ρ 1+ = cosh 2α cosh 2α − sinh 2α sinh 2α cos ωt , cosh 2ρ 1− = cosh 2α cosh 2α + sinh 2α sinh 2α cos ωt , where to simplify these expressions, we have introduced In fig. 2, we plot the results for C κ=2 and C is initially linear, while that for C 2 and C κ=2 is quadratic. The examples depicted in fig. 2 exhibit this behaviour. We also see that the complexity oscillates in time, but that the amplitude of the oscillations decreases as βω increases. From eq. (20), we see that for large βω, α exp[−βω/2], and expanding the expressions in eq. (152) shows that the oscillations are indeed exponentially suppressed in this regime. This fact will allow our results for the field theory to be integrated with respect to the frequency in the next section.

A simple limit
Recall that above we have set λ R = 1, which amounts to setting the gate scale and referencestate scales equal, i.e., ω 2 g = M µ. Hence from eq. (40) we have Now let us consider the limit in which µ is much bigger than any other scale, i.e., λ → 0 or α → −∞. In this limit, we have cosh 2α 1 2 in which case eq. (152) simplifies to cosh 2ρ 1± 1 2 µ ω (cosh 2α ± sinh 2α cos ωt) , , and 2 (green). We note that the complexity oscillates in time with periodicity δt = π/ω as expected from the explicit expressions and with an amplitude that decreases (approximately exponentially for large βω) for increasing βω.
We can also substitute the simplified expressions from eq. (157) into the F 2 complexity in eq. (149) to find where we have again subtracted the zero-temperature contribution, which in this case is We have also dropped terms which are suppressed by inverse powers of log µ ω . Finally, substituting the simplified result (157) into the F 1 complexity (151), we find where the subtracted zero-temperature contribution is C is maximal when | cos ωt | = 1, i.e., t = nπ/ω. In particular, we have a maximum at t = 0.
We can also consider the opposite limit where µ is much smaller than any other scale, i.e., λ → ∞ orα → ∞ from eq. (154). Again, we have set ω 2 g = M µ. Then in this limit, we have cosh 2α 1 2 in which case eq. (152) simplifies to These expressions then produce the simple solution which is quite similar to that found above in eq. (157). We can substitute these simplified expressions into all of the various expressions for the complexity from the different cost functions above, but let us focus on the F 1 complexity in eq. (151). In this case, we find an identical result to eq. (160), namely ∆C (LR) 1 = log (cosh 2α + sinh 2α |cos ωt|) .
Thus, even though we are considering the opposite limit here (i.e., λ → ∞ rather than λ → 0), ∆C is unchanged. In particular, we still find that the complexity has a maximum at t = 0. The results for ∆C 2 and ∆C κ=2 are very similar to those obtained using the previous limit.

Complexity of the TFD at general t with λ R = 1
For the general case λ R = 1, the boundary conditions (107) for the circuits leading to the time-dependent TFD state (76) read Here we have used the parametrization in eq. (111) in terms of the values (ρ 1 , θ 1 , τ 1 ) at the end point of the trajectory, and as usual G R is given in eq. (47). We have also focused on the + mode in the TFD. The TFD state at t = 0 with general λ R is obtained by simply setting t = 0 in the results of this section. Of course, since the determinant of the matrices on either side of eq. (107) is one, there are again only two independent relations in eq. (165), and as a result there will be a one-parameter family geodesics (i.e., circuits) that produce the desired target state. As before, this family of end points maps out a spiral in the space of unitaries spanned by (ρ, θ, τ ). However, an important feature is that for λ R = 1, ρ 1 is no longer a fixed constant, but rather varies as we move along the spiral. Hence our previous arguments that the optimal circuit corresponds to ∆θ = 0 = τ in section 4.1 no longer apply, and we must undertake a more extensive analysis of all possible geodesics ending on the spiral. Given an end point, we may use the geodesic solution (115) evaluated at σ = 1 to solve for c, ∆θ, and θ 1 for a given boundary condition (in general, this solution is found numerically below). The optimal circuit will be given by eq. (115) for general values of σ, and its length can be evaluated according to eq. (116). We must then minimize this length over all possible solutions of eq. (165) in order to find the optimal circuit. One particular solution to eq. (165) arises naturally from our construction of the timedependent TFD state, namely, that which corresponds to the unitary where U vac→TFD(t) is given by eq. (78), and U R→vac by eq. (81). Other end points can then be obtained using the stabilizer group of the reference state by multiplying this transformation on the right by U φ in eq. (85) (see also eq. (118)), More explicitly, one extracts from the unitary U T (φ) the corresponding end point (ρ 1 (φ), θ 1 (φ), τ 1 (φ)) using eq. (111) via the relations Alternatively, we may derive differential equations for the end points by varying the two sides of eq. (165) with respect to φ, whereupon we find These differential equations provide an efficient way of generating end points for geodesics in our numerical solutions below. Once we have the family of end points (ρ 1 (φ), θ 1 (φ), τ 1 (φ)), we can compute the parameters (θ 0 (φ), c(φ), ∆θ(φ)) of the corresponding geodesic by inverting the boundary condition numerically and then evaluate the corresponding length (i.e., complexity) of the geodesic ending at U T U φ as We can gain some intuition for this set-up by plotting the curves of equivalent end points in the three-dimensional space spanned by (ρ cos(θ), ρ sin(θ), τ ). These curves have a spiral like shape, and their projection on the τ = 0 plane is a closed curve-see the plots in fig. 3. We may use the point where the curves cross the τ = 0 plane as a representative point for each spiral. In the figure, we have chosen two such representative end points, one on the ρ sin(θ) axis in the upper half-plane, and a generic point in the upper-right quadrant.
Let us now focus on the example of the TFD state at t = 0. Using the explicit expression for the symplectic transformations which bring us to this state, i.e., eq. (166) together with eqs. (78), (81), and (85) evaluated at t = 0, as well as the inversions (168), we find a family of end points (ρ 1 (φ), θ 1 (φ), τ 1 (φ)) associated with this state. We do not write them explicitly here since the expressions are rather lengthy and uninformative. One of these end points, corresponding to τ 1 = 0, turns out to be associated to the value φ = 0 in the unitary U φ . This means that the natural symplectic transformation coming from our construction of the TFD state was associated with τ 1 = 0. Now, the corresponding value of θ 1 can be read from the inversion formula (168) for the case of φ = 0 in eq. (166), to get This means that θ 1 = π/2 for e −2α λ/λ R > 1, while for e −2α λ/λ R < 1 we have θ 1 = −π/2. Now let us consider the expression for ρ 1 (φ). This determines the projection of the spiral of equivalent end points on the τ = 0 plane. We obtain the following expression: where we have defined e −2y := λ R and e −2x := λe −2α . The above curve is a closed shape with a long and a short axis where and the short axis is either aligned with φ = 0 or π/2 depending on the signs of x and y. In particular, when they both have the same sign (i.e., xy > 0), the short axis is at φ = 0, and hence given the previous observations, it is aligned with the point at which the spiral crosses the τ = 0 plane. Recall that the geodesic ending at τ 1 = 0 is still the straight-line geodesic with ∆θ = 0, which remains in the τ = 0 plane for the entire trajectory. 25 Since for other values of φ we will have ρ 1 (φ) ≥ ρ 1 (0) and ∆θ 2 ≥ 0, we can use the inequality (125) to argue that the shortest geodesic is in fact still the simple straight-line geodesic (121), where the previous analysis indicates that θ 1 = ±π/2. Alternatively, if x and y have opposite signs, then the short axis is aligned with φ = π/2, and we should expect that the optimal geodesic will no longer be the simple one which remains in the τ = 0 plane. Numerical testing reveals that this is indeed the case. This means that even for the TFD at t = 0, there exist values of the parameters for which the straight line geodesic does not represent the optimal trajectory. This is illustrated in fig. 4. This conclusion continues to hold even in the limit of low temperatures where the two sides of the thermofield double decouple from one another and our target state becomes two copies of the vacuum state. This is simply achieved by setting α = 0 in the definition of x. Exploring the relevant ranges of x and y reveals that we deviate from the straight line trajectories when λ > 1 and λ R < 1 or when λ < 1 and λ R > 1. This is different from the conclusion of [23] where only W i,j gates were used, see eq. (91), and indeed we see that when exploring the full symplectic group, shorter trajectories exist using the full set of gates in eq. (91).
Numerical testing reveals that the optimal geodesics deviate from the τ = 0 plane in all cases where the spirals do not intersect τ = 0 at θ = ± π 2 . Of course, the alignment with θ = ± π 2 happens not only at t = 0 but more generally for ω k t = πn where n ∈ N (since our solutions are periodic in time). When n is even, the conditions are identical to the previous case. When n is odd, one has to substitute α → −α to obtain the relevant conditions for the trajectories to move in the τ = 0 plane. Of course, if we consider the TFD for oscillators with different values of ω, 26 the time for this alignment will differ for the different oscillators and therefore, in the QFT calculations, the time in which all trajectories move in the τ = 0 plane for all values of ω k can only be achieved at t = 0.
We have plotted the time dependence of complexity (using the κ = 2 cost function, i.e., eq. (150)) for the TFD state (76) in fig. 5. We have chosen to plot it as a function of ωt to account for the periodicity of the result. One striking feature of all of these plots is that the complexity (or rather the difference C κ=2 (t) − C κ=2 (0)) decreases as λ R is varied away from one, i.e., , the complexity decreases as both λ R is increased or decreased. We see that as before the complexity decreases exponentially as we increase ωβ. We observe that as we increase µ/ω, the λ R < 1 result becomes very close to the one for λ R = 1. We do not understand what is the reason for this behaviour at present and leave this issue for future study. The low temperature limit is obtained by taking βω → 0 while keeping µ/ω fixed. Even in this limit, 25 The analysis of the geometry and the geodesics is not changed by choosing λR = 1, only the positions of the end points. 26 For example, the regulated scalar field theory in the next section. In this plot we have fixed ω 2 R β/M = 10 and βω = 1 and this means that λ R /λ = 10 and focused on the + mode. In this case we expect that for 1 < λ R < 40.84 we will have deviations from the straight line trajectories. Outside this range, we recover straight line trajectories and since the ratio of λ/λ R is fixed ((x − y) is fixed), eq. (172) predicts no dependence on λ R outside this range and this is indeed what we see in the figure when comparing to the λ R = 1 value. We remark that this plot does not have the resolution to show exponentially small deviations around ∆C κ=2 = 0.
we see that significant deviations are obtained from the previous λ R = 1 results. Further, we observe that for certain values of λ R , the complexity is decreased compared to that at t = 0 for all other times.
In the next section we explain how to use the results of this section for the complexity of the TFD state of two simple harmonic oscillators for the purpose of computing the complexity of the TFD state of two copies of a free scalar field theory. We will discretize our field theory on the spatial lattice in such a way that the TFD state becomes a product of these harmonic oscillator TFD states for the normal modes on the lattice.

Complexity of TFD states in quantum field theory
In this section, we combine the results of section 4 with the lattice discretization and normal mode decomposition of free scalar QFTs explained below in subsection 5.1 in order to obtain the complexity of TFD states of free scalar QFTs. The results of subsection 5.1 also serve as background for the entanglement calculations presented in section 6. In subsection 5.2, we present results for complexity obtained using lattice regularization, which will turn out to be particularly relevant for the comparison with the physics of entanglement entropy. In subsection 5.3, we focus on the complexity of formation of TFD states at t = 0 with λ R = 1. In subsection 5.4, we use a similar approach to study the complexity of the TFD of a free scalar QFT at t = 0 with λ R = 1. Finally, in subsection 5.5 we comment on the case λ R = 1. The plot is always periodic as a function of ωt with periodicity π as previously, and so we have only plotted one period. The result for the complexity is suppressed as we increase βω. When increasing the reference state scale encoded in the parameter µβ, the curves for λ R < 1 approach the one of λ R = 1. We do not fully understand what is the reason for this behavior and leave it for future study. In some instances the (regularized) complexity is negative for all times.

Normal mode decomposition for a free QFT
Let us start by focusing on a 1+1-dimensional free QFT living on a (Lorentzian) cylinder of circumference L, defined by the Hamiltonian We regulate this theory in the UV by putting it on a lattice with lattice spacing δ > 0. If we assume that the lattice has N sites arranged on the spatial circle, we then have The Hamiltonian then takes the form of N coupled harmonic oscillators, with the redefined canonical variables 27 where Q N +1 := Q 1 and P N +1 := P 1 . Passing to Fourier space gives where we have the Fourier transformed variables If N is even, we also haveQ † N/2 =Q N/2 andP † N/2 =P N/2 , which are therefore also real. The canonical commutation relations are given by [Q k ,P † ] = i δ k . If we define the frequency our Hamiltonian (177) takes the form which is a sum of N independent harmonic oscillators 28 of equal mass M = δ −1 (not to be confused with the physical mass m), but k-dependent frequencies ω k , cf. eq. (10). It is 27 The scaling here ensures that Qa and Pa have the usual dimensions of positions and momenta, respectively. 28 Note that our variablesQ k andP k are complex, except for k = 0 and k = N/2 (for even N ). The two complex degrees of freedom labeled by (Q k ,P k ) and (Q N −k ,P N −k ) only contain two real degrees of freedom, because they are subject to the constraintsQ † k =Q N −k andP † k =P N −k . When expressing the Hamiltonian in terms of the real and imaginary parts of these variables, we find two independent real harmonic oscillators with common frequency ω k . This is discussed in more detail in appendix D.
important to point out that when m = 0, the frequency ω 0 vanishes and, as a result, the zero mode (k = 0) Hamiltonian does not have a normalizable ground state. Of course, the case of m = 0 represents the conformal limit, and as such is the case of primary interest in our paper, since we want to compare with holographic CFTs. The most straightforward means of dealing with the potential problems stemming from the presence of this zero mode is to regulate its behaviour by introducing a small but non-vanishing mass (e.g., m 1/L), which we shall do below. We can then obtain results with decreasing values of this IR regulator.
We take the reference state to be the ground state of the ultralocal Hamiltonian where the spatial derivative term is absent, cf. eq. (173) and where µ > 0 takes the role of the mass of this fictitious Hamiltonian. After performing the above discretization and Fourier transform, we arrive at Hence the momentum modes remain decoupled and the reference state is a simple product over the momentum modes of Gaussian wavefunctions, all with a fixed width set by µ. Of course, the ground state of the physical Hamiltonian (180) has the same form where the variance of each mode is set by ω k . Since each mode k is decoupled from the other modes, the respective TFD state will be the product of TFD states for each of the oscillators. This brings us to the setup of subsection 2.2.2, after we make the following identifications for each mode: The dimensionless ratios λ and λ R in eq. (40) for each mode now take the form It may seem curious that these coefficients in eq. (184) seem to depend on the lattice spacing δ. However, it turns out that the natural gate scale must be modified when working with the field theory as follows: for a general spacetime dimension, the relation (176) becomes Q a = δ d/2 φ(x a ) and P a = δ d/2−1 π(x a ). Following the structure in eq. (38), we express the quantum circuits of interest as Note that in going from the lattice to the continuum expressions, we have absorbed a factor of 1/δ d−1 into the control functions, and with this choice, the three functions y,ŷ, andỹ all have the same dimensions, i.e., length −(d−1) for the field theory in (d − 1) spatial dimensions. 29 In the continuum, we have also defined the gate scale as which naturally appears in the φ 2 and π 2 gates. This also absorbs the lattice spacing δ in eq. (184). That is, the dimensionless ratios in eq. (184) reduce to Additionally, note that our most heavily analyzed case above, namely λ R = 1, now equates the new gate scale with the reference scale, i.e., µ g = µ, and gives λ = ω k /µ. We can also take the limit of a large chain, L δ (i.e., N 1 while keeping δ fixed), in which case the system becomes infinite and we obtain a lattice-regularized quantum field theory living on an infinite line, where In these expressions, we have introduced the continuous label p ∈ [− π δ , π δ ], defined as forQ p ,P p . The range of p corresponds to the Brillouin zone familiar from the physics of crystals, see e.g., ref. [83]. What we have done is introduce a UV-regularization and mode decomposition in which the Hamiltonian of a continuous quantum-many body system becomes a sum over independent harmonic oscillators (bosonic modes). In subsequent parts of this section, we will calculate the complexity for a quantum field theory by simply adding up (or integrating over) contributions for each bosonic mode. Before we conclude this subsection, it is worth commenting on a different regularization scheme introduced for the purposes of continuous multiscale entanglement renormalization ansatz (cMERA) in ref. [84] (see also refs. [85,86]), which was used in the context of QFT complexity in ref. [51]. In this approach, the UV regulator is introduced by modifying the Hamiltonian of the free quantum field theory in momentum space, such that the ground state behaves at large momenta |p| > Λ as a product state, i.e., the ground state of the ultralocal Hamiltonian (181). In particular, taking the reference state to be the ground state of an ultralocal Hamiltonian (182) and working in momentum space, one can truncate the momentum integrals at |p| = Λ when computing the complexity, since for higher momenta the state already has the form of the target state.
We extend this approach to the complexity of the thermofield double by requiring that we only reproduce the TFD state up to a cutoff scale Λ in momentum. This means that we can again cut our momentum integrals at p = |Λ|. The result with this regulator can be easily obtained by starting with the previous lattice regularization and placing the new regulator Λ far below the scale of the lattice spacing Λ π/δ, such that upon sending δ → 0, the result remains finite and regulated by the new scale Λ. In light of eq. (188), and using the fact that Λ π δ , we may linearize the sine function in eq. (189) for ω p as We may also demand that the frequency of the oscillator is continuous at the transition point, i.e., ω p=Λ = µ where ω p was defined in eq. (189).
Of course, as alluded in the discussion around eq. (185), the analysis discussed here generalizes in a straightforward manner to higher dimensions. In the Hamiltonian (188), one then replaces the integral over the range [− π δ , π δ ] by an integral over the corresponding (d−1)dimensional hypercube. When calculating complexity in the cMERA-inspired approach, one replaces the one-dimensional integral over momenta by an integral in momentum space confined to the ball of radius Λ centered at the origin.
Lastly, we would like to comment that most of the observables in which we will be interested are regularized quantities, and will exhibit exponential suppression for momenta larger than the temperature scale. Hence in this case the details of the regularization scheme will not matter as long as β δ (or β 1/Λ).

5.2
Warm-up: complexity as a function of time with λ R = 1 on the lattice As a warm-up for both the complexity calculations in field theory and the analysis of entanglement in section 6, we present here representative lattice results in (1 + 1)-dimensions. For concreteness, we again choose the gate scale such that λ R = 1 and focus on the κ = 2 cost function, see eqs. (9) and (150). We will consider three cases: keeping the total size L of the system fixed, i.e., working on a circle, see figs. 6 and 7; keeping the lattice spacing fixed and increasing the total number of lattice sites; and working with infinite system at fixed lattice spacing-see fig. 8 for these last two cases. The motivation to study these three cases is, first, to isolate finite size effects and, second, to see how well the cMERA-type regularization in eq. (191) agrees with the answer on an infinite lattice. As mentioned above, when λ R = 1 the parameter λ for each mode becomes ω k /µ, see eq. (187). We now have to add contributions from all modes with each mode contributing as in eq. (150), see also eq. (142). For a finite lattice this gives where for compactness we have used the identity ( Here, α k is given by eq. (183), ω k by eq. (179), and µ is the QFT gate scale defined in terms of eqs. (181) and (183), see also the discussion below eq. (187). In order to make contact with holography, we are primarily interested in CFTs (i.e., massless models). However, the massless limit of eq. (192) is ill-defined, since the zero mode gives a divergent contribution. As explained above, we regulate this divergence by instead working with a small but non-zero mass. In most of the analyzed cases, we chose it to be m = 10 −6 /L, where L is the circumference of the circle. Finally, the continuum limit on a circle can be approached by keeping the size of the circle and other parameters fixed while increasing the number of points. In order to UV regulate our result for the complexity, we subtract its value at the initial time t = 0.
In fig. 6, we look at complexity as a function of time for theories close to the conformal limit for different temperatures. We measure time in units of the inverse temperature β, since ultimately we want this to be the dominant scale. We observe that at low temperatures the intermediate late time behaviour is dominated by a logarithmic growth, i.e., by a term proportional to log 2 t/β (which we associate with the zero mode below), whereas for higher temperatures we see initially propagating wave packets on a circle superimposed with the logarithmic growth and ultimately, at large temperatures, saturation. Subsequent results in subsection 5.4 applicable to a line are consistent with the finding that the zero mode becomes subdominant in the large temperature limit.
In fig. 7, we investigate the effect of the IR regulator mass on the time-dependence at low and intermediate temperatures. What we see is that decreasing the mass allows us to recover more and more of the logarithmic growth, and the transition to the oscillatory regime occurs at times inversely proportional to the mass. Therefore, we interpret the logarithmic growth as a feature associated with the zero mode. Let us analyze how the logarithmic growth arises primarily from the zero mode contribution in eq. (192) (i.e., only counting the contribution from k = 0). In the limit where both m/µ and βm are small parameters, the coefficients f ± 0 simplify to Now, if we expand for large times (for instance, compared to the temperature or the ultralocal mass µ), but still such that mt is small, we have Notice that as long as the condition mt 1 is satisfied, the effective dimensionless quantity associated with the time dependence becomes to leading order µt 2 2β . Therefore, in this particular regime, eq. (192) for the zero mode reads which corresponds to the behaviour shown in figs. 6 and 7.
The last thing that we investigate in this subsection is the comparison between the κ = 2 complexity on a circle with fixed lattice spacing δ, the analogous lattice calculation on the line, and a calculation using the cMERA inspired regularization, see subsection 5.4. As one sees in fig. 8, all of the approaches beautifully agree in their overlapping domains of applicability, i.e., for times which are not so large that one becomes sensitive to finite size effects.
Let us briefly comment on why the dependence of complexity on the zero mode is not present in the decompactification limit, where the circumference of the circle becomes infinite. 30 The maximum complexity is reached at times π/(2m), which from eq. (194) implies 30 We can also think of this as a high temperature limit, because in the continuum with an arbitrarily small IR regulator mass, the physics is controlled by the combination L/β = L T .  Figure 6: Grey curves represent time dependence of κ = 2 complexity with the initial value subtracted for the TFD on a circle with circumference L with µ = 1/L, m = 10 −6 /L and increasing temperatures T = 1/β. We use 1601 lattice sites on each side. For smaller temperatures (top two plots), we observe a late-time behaviour proportional to log 2 (t/L). The green dotted curve demonstrates an excellent fit of a 2 log 2 (t/L) + a 1 log(t/L) + a 0 to the full function (upper left plot) and its maxima (upper right plot) for later values of the time. For higher temperatures (bottom two plots), we observe saturation. Transitioning between these two regimes are oscillations, which occur with a period of half of the circle's circumference, as if two wave packets were propagating on a circle with the speed of light in opposite directions. One should think of the saturation as resulting from the presence of many modes non-trivially contributing to the sum (192) at high temperatures, where the zero mode contribution becomes subdominant when dividing by the thermal entropy which scales with the temperature, see the discussion around eq. (197).  that the contributions from f ± 0 should be essentially the same. The complexity difference at times 0 and π/(2m) then diverges logarithmically with the small masses m, Naively, this equation suggests that the zero mode should dominate the growth of complexity indefinitely. However, for a decompactified system, the contribution from the other modes in the sum given by eq. (192) will result in a scaling with the volume of the system, which in the decompactification limit is taken to infinity. Therefore, if we evaluate a ratio of quantities such as the difference in complexities divided by the thermal entropy of the system, we expect that the zero mode contribution will be subleading with respect to the other modes. Of course, the thermal entropy of the system increases with the temperature and so we expect the same suppression of the zero mode in the limit of high temperatures, as we see in fig. 6. We will explicitly confirm in section 5.4 that for decompactified theories in the continuum limit, there is not an unbounded growth.

Complexity of formation with λ R = 1
In this subsection, we evaluate the complexity of formation [27]. That is, we evaluate the extra complexity required to prepare the two copies of the scalar field theory in the TFD state compared to simply preparing each of the copies in the vacuum state. We consider the scalar in d spacetime dimensions and with mass m, as well as the m → 0 limit. For simplicity, we focus on the case λ R = 1. We will start with the L 1 -norm in the left-right (LR) basis in eq. (140), for which the complexity of formation in the continuum limit reads Here we have used the cMERA inspired regularization discussed at the end of subsection 5.1. One may worry that this integral will produce UV divergences. However, these potential divergences are eliminated because we have (positive) powers of Λ competing against powers of exp[−βΛ] in these contributions, and so they actually vanish in the limit that Λ → ∞.
In fact, we can therefore remove the UV regulator altogether and integrating all the way to infinite momenta still leaves a finite result. As a result, one obtains with where Ω d−2 is the volume of a (d − 2)-sphere, i.e., Ω d−2 = 2π The complexity of formation for holographic CFTs in a flat decompactified space is directly proportional to the entropy, for d ≥ 3 in both CA and CV proposals [27]. Therefore, it is natural to normalize the complexity of formation by the entropy. Let us then consider a gas of free bosons. Its partition function in d dimensions reads and hence the thermodynamic entropy is given by We can rewrite this expression as where s th (x) is defined as The ratio of the complexity of formation and entropy is then simply the ratio of the two functions of βm in eqs. (200) and (204), For the massless theory, the ratio has a simple analytic expression, As argued previously, we find a similar agreement with holography, where the complexity of formation scales like the entropy, with a dimensionless coefficient that increases with the dimension of spacetime. For the CA and CV proposals in holography, the coefficients of proportionality of complexity and entropy were [27] ∆C The CA coefficient increases essentially linearly with the dimension, while the CV coefficient increases at first but very quickly saturates to a constant. The result of the massless free scalar in eq. (206) also increases with dimension, but it grows exponentially fast. Generally, comparing the QFT expression (206) with the holographic results (207), we see that in both frameworks, the complexity of formation is UV finite (i.e., independent of the cutoff δ), positive, and independent of the reference state scale µ (or of the normalization constant α or of the counterterm scale ct , in the case of holography). Of course, in both cases, we also have ∆C ∝ S th to leading order. However, one difference is that for d = 2, the holographic complexity of formation is a constant (i.e., independent of the temperature), while eq. (206) is still proportional to the entropy for d = 2. Next, we evaluate how the ratio of complexity of formation to entropy behaves for a massive scalar, given by eq. (205). We show in fig. 9 the numerical evaluation of how the coefficient of the complexity of formation in eq. (205) increases once the theory is massive, for various dimensions. Both the complexity of formation and the entropy go to zero as the parameter βm increases, but the ratio increases exponentially as a function of βm. We can expand eq. (205) for large masses, resulting in

Comments on the diagonal basis
We now turn our attention to evaluating the complexity of formation with the L 1 -norm in the diagonal basis with λ R = 1, where each mode contributes according to eq. (140). The complexity of formation reads In order to evaluate the above expression, one has to study carefully how the sign changes inside each argument of the logarithms, and this behaviour depends strongly on the reference scale µ. In appendix E, we break down carefully how one computes the complexity of formation in the diagonal basis for µ smaller than, equal to, and larger than the cutoff Λ in the massless scalar theory (where ω k = k). For now, let us focus on the case where the reference scale is higher than the cutoff (µ > Λ) for a massless theory. In this regime, and for integration variable that ranges from 0 to Λ, 1 2 log k µ + α k changes sign at a value k f given by There are two important limits to the above transcendental equation: when k f is very small and when k f is close to the cutoff scale Λ. Solving for the temperature in these two regimes, we find For T < T c1 , the arguments of all the absolute values in eq. (209) for k in the range [0, Λ] are negative. As a consequence, the complexity of formation is identically zero,  211) and (212)). This is the result obtained when holding the temperature fixed while sending the cutoff to infinity. For temperatures of the order of the UV cutoff, the complexity of formation in the diagonal basis develops a dependence on the temperature and the cutoff scale Λ. Of course, the results in this unphysical regime are sensitive to the UV regulator (and are only presented for illustrative purposes). The fact that the complexity of formation is either zero or not proportional to the entropy contrasts with the holographic results of ref. [27].
Next, if the temperature is within the range T c1 < T < T c2 , there is a single solution k f to the transcendental equation (210) in the range [0, Λ]. We find that for k < k f the argument of the first absolute value is negative, and for k > k f the argument is positive. The complexity of formation in this situation is found by integrating only over modes larger than k f , Finally, if the temperature is bigger than T c2 , we find that 1 2 log k µ + α k is always positive in the range of momenta [0, Λ]. Therefore, we have We show in fig. 10 the integrated complexity of formation divided by the entropy for several dimensions. For small temperatures with respect to the cutoff scale, we find that the complexity of formation is exactly zero. For higher temperatures, there are some nontrivial cancellations between the circuits that introduce some dependence on T and Λ that contrasts with the LR basis and the holographic results of ref. [27]. Of course, the high temperature regime with T ∼ Λ should be regarded as unphysical and we are only presenting the results for illustrative purposes. Furthermore, we note that in this regime, the results are sensitive to the details of the UV regulator (which we simply chose as a hard cutoff on the momentum integral for the calculations here). In appendix E, we also analyze the other cases µ = Λ and µ < Λ. Generally we still find that the complexity of formation in the diagonal basis vanishes for small temperatures and it is only a nontrivial function of the temperature with T ∼ Λ.

Comments on different cost functions
Let us briefly comment on the integrated complexity of formation for different cost functions with λ R = 1. For the κ = 2 cost function, 31 the one-mode complexity is simply given by eq. (139). The integrated complexity of formation then reads We present the dependence on the dimension of the ratio of the complexity of formation to the entropy for different cost functions in fig. 11. We focus on the massless theory. The ratio decreases when the dimension increases for the κ = 2 cost function, in contrast to the exponential increase in the L 1 norm using the left-right basis, see eq. (206). In addition, the complexity of formation for the F 2 cost function in eq. (139) resembles the structure of the complexity of formation for the L 1 norm in the diagonal basis. We find again that the integrated complexity of formation in general depends on the cutoff scale Λ and in the physically relevant regime where the temperature is much smaller than the cutoff scale, the ratio of complexity to entropy vanishes which does not match with holographic expectations.

Time dependence with λ R = 1
Next, we investigate the time dependence of the TFD state in the continuum limit. The full time dependence in the L 1 norm in the LR basis with λ R = 1 is given by integrating eqs. (151) and (152) over momenta. 31 Recall that the complexity is basis independent with this cost function, i.e., the same complexity results using either the left-right or diagonal basis in this case.

"Simple" limit
Let us start by studying the time dependence in the "simple" limit considered in eqs. (160) and (160). These simple results correspond taking either the limit λ → 0 or λ → ∞ (as well as setting λ R = 1). For the first limit, it seems that in the corresponding field theory calculation, we would need to considering the case where the reference scale µ is much larger than the frequency ω k of all possible modes. However, as we will see in a moment, very high frequency contributions are exponentially suppressed and so in fact, we need only consider µ T = 1/β. For the second limit, we would be considering µ ω 0 = m, i.e., the reference scale is much lower than all frequencies and so is much lower than the minimum frequency, which corresponds to the scalar field mass. In either case, we find Recall that in this result, we are subtracting off the zero-temperature complexity (i.e., subtracting off the contributions coming in the limit α k → 0). Examining the above expression, we note that as in the complexity of formation calculation, the integrand in eq. (216) is exponentially suppressed for high energy modes. Effectively, this means we can consider the UV cutoff to be much higher than any other scale, and simply integrate up to infinity. Furthermore, this also means that the UV divergences in the complexity of the time-dependent TFD state still match those of (two copies of) the vacuum -see the discussion around eq. (198). We present the time dependence of the complexity in this simple limit in figs. 12 and 13 for a massless theory in d = 2 and d = 3, respectively. In these plots, we present the difference between the complexity of the TFD at a given time t and that at t = 0 by subtracting from eq. (216) its value at t = 0, i.e., the complexity of formation, see eq. (209). We find that this difference of complexities actually decreases as a function of time, which can be understood from eq. (160). The maximum is at t = 0 since the contributions of all of the individual modes will take their maximum value at this time. However, the oscillating factors in eq. (160) will all become misaligned after this initial time. Hence the complexity begins by decreasing and it never recovers the maximum value. Instead, the complexity saturates at a new value that is still of the order of the entropy (for m = 0) on a time scale which is of the order of the inverse temperature. When summing over all modes, the exponential suppression (see discussion after eq. (152)) implies that the modes with frequencies of order β or less are dominant in the summation. This means that the oscillations of different modes (with periodicity 2π/ω k ) will be maximally dephased after times of order t ∼ β, which explains the saturation that we observe.
Of course, this is very different than what we see for holographic complexity, i.e., where we see a linear growth at late times. However, this difference is not unexpected. The intuitive argument for the linear increase of the holographic complexity is that the effect of the chaotic/fast-scrambling Hamiltonian is like throwing random gates at the state, and in the large-N limit, it is very unusual for a new gate to reduce the complexity. Since we have (almost) the simplest of possible Hamiltonians with the free theory in our QFT calculations, we should not expect to find analogous behaviour here -see further discussion in section 7.
Next, we investigate the influence of turning on the mass of the scalar field on the time dependence in the "simple" limit given by eq. (216). We show the time dependence of complexity for different masses in figs. 14  temperature. These oscillations arise because of the simple dependence on time cos ωt in eq. (216), from which follows an oscillatory behaviour with period ∆t ≈ π/m.

General reference scale
We now turn our attention to the effect of the reference scale on the time evolution of complexity. That is, we will consider the time dependence for general values of λ. In order to take into account the reference scale, one has to integrate the complexity contributions in eqs. (151) and (152) over all momenta. For concreteness, let us focus on the massless case. We define the dimensionless variables k := βk ,t := t β ,γ := 1 βµ , whereγk is simply the parameter λ as in eq. (154). The two simple limits above are theñ γ → 0 andγ → ∞. In figs. 16 and 17, we investigate the time evolution of complexity for different values of the dimensionless parameter related to the reference scale. For both small and large values ofγ, we recover a similar behaviour to the one discussed previously in figs. 12 and 13, with the complexity decreasing at first, then saturating to a constant. For intermediate values of γ, the complexity increases with respect to its value at t = 0, and then it again saturates to a constant value. Of course, despite increasing at first, the complexity does not continue increasing linearly with the energy for long times, as found for holographic complexity.
Let us comment further on another possible lesson from holography that is manifest in figs. 16 and 17. The time dependence of complexity for the TFD state dual to an eternal black hole exhibits non-universal behaviour at early times due to the normalization of the null normals to the boundary of the Wheeler-DeWitt patch in the complexity=action proposal. As discussed in ref. [28], the transient behaviour of the time derivative is controlled by a dimensionless parameter α, as even under affine parametrization there is freedom in the overall scale of the null normals. Even if we fix reparametrization invariance by a boundary counterterm, as introduced in ref. [26] and argued to be necessary in order to reproduce   desired properties of complexity in refs. [38,39], the transient behaviour is also controlled by an arbitrary dimensionless parameter ct /L. 32 In this sense, figs. 16 and 17 reproduce the intuition of holography, that at early times the evolution of complexity is dominated by non-universal effects, such as the scale of the reference state. The late time growth rate of complexity in holography was independent of these ambiguities in the null boundaries, which also seems to be a property in figs. 16 and 17. However, unlike holography, the growth rate of the complexity for the massless Gaussian TFD vanishes after times of the order of the inverse temperature. We will return to this point in the discussion in section 7.

κ = 2 cost function
We now investigate the time evolution for the κ = 2 cost function given by eq. (150). The integrated complexity has the simple expression where ρ 1 + and ρ 1 − where defined in eq. (152). We show the time evolution for different reference scalesγ in fig. 18. Although the complexity still saturates to a constant at times of the order of a few inverse temperatures, the reference scale changes significantly the rate of change during the transient regimes at early times, and the time derivative grows as the reference scaleγ becomes very large or very small. This dependence on the reference scale in the time derivative is something that holographic complexity also exhibits when a counterterm is added to the null boundaries of the Wheeler-DeWitt patch (e.g., see ref. [39]). Nonetheless, we should add that the UV structure for this cost function has an extra logarithmic factor compared to the UV divergences observed in holographic complexity, i.e., C κ=2 ∼ V δ d−1 log 2 ( /δ) compared to C holo ∼ V δ d−1 log( /δ).   Figure 19: Integrated regularized κ = 2 complexity as a function of time for general λ R in 1 + 1 and 2 + 1 dimensions with µβ = 10 and m = 0. Note that the λ R < 1 curves are rather close to the λ R = 1 curve and that the λ R = 5 curve in 2 + 1 dimensions is completely negative. The saturation time in 1 + 1 dimensions is longer, similarly to what happened in fig. 18 for λ R = 1 andγ = 0.1, however we stress that this is not a logarithmic growth.

Comments on λ R = 1
Finally, we investigate the influence of varying the gate scale (λ R = 1) on the D κ=2 cost function for a massless theory in d = 2 + 1 and d = 1 + 1 dimensions with reference scale µβ = 10 in fig. 19. Varying the gate scale does not change the vanishing rate of change after a time of the order of the inverse temperature, however, it does change the growth rate at times t ∼ β, as well as the final value of the complexity. The maximum value for the complexity is found for λ R = 1, which suggests that the results we have obtained for λ R = 1 can be understood as an upper limit on the complexity for other possible values of the gate scale. It would still be interesting to understand how the gate scale influences other properties of the complexity, such as the structure of the UV divergences, and whether there is a simple intuition as to why the circuit complexity has a maximum at λ R = 1. We can also observe in the figure that for some values of λ R (in this case λ R = 5 in 2 + 1 dimensions) we obtain a complexity that is always lower than its value at t = 0 since the relevant curve is everywhere negative. When the reference scale is large, the curves are similar to those of λ R = 1, for reasons that we do not fully understand at present. In 1+1 dimensions, the transient regime is longer, similar to what happened for λ R = 1 andγ = 0.1, see fig. 18. However, we stress that the plot does not exhibit a logarithmic growth and the complexity finally saturates, which we tested by plotting the exponent of the regularized complexity against time.

Entanglement production in TFD states
Among other motivations, circuit complexity was proposed as a quantity that can probe aspects of field theory states that cannot be captured by entanglement entropy. 33 In particular, when studying the thermofield double state of an interacting field theory with a holographic dual, it is expected that circuit complexity can probe the degrees of freedom which are encoded in the dual geometry deep in the interior of black holes. In holography two proposals have been made for the holographic dual of complexity -the CV and CA proposals and they both predict that the complexity keeps increasing for a very long time as long as the geometry is still well-approximated classically, e.g., [13,14]. On the other hand, the entanglement entropy will not keep increasing for such a long time, and instead simply saturates at times of the order of the subregion size [3]. In this section, we will be primarily concerned with the time evolution of the entanglement entropy in the (1+1)-dimensional free field TFD state on a circle. We will consider the entanglement entropy of a subregion of the entire quantum system which contains parts in both the left and right QFTs. However, our free field theory does not have a holographic dual and so we do not expect to see the same contrast in the behaviours of the circuit complexity and the entanglement entropy, as seen in holography. Still we believe that it is illuminating to analyze the entanglement dynamics where it can be done semi-analytically, and to compare our results with those for complexity.

Entanglement entropy for Gaussian states from covariance matrices
For bosonic Gaussian states, the knowledge of the covariance matrix G A for canonical coordinates supported in a subregion A is equivalent to the knowledge of the reduced density matrix.
Of course, G A itself admits a decomposition into normal modes, 34 and the von Neumann entropy and higher Rényi entropies of the reduced density matrix are sums of independent contributions from each of these modes. The entropies can be then computed by viewing each mode as an auxiliary thermal system.
The remaining question is then how to efficiently implement the above logic to evaluate the entanglement entropy. One can indeed easily find the expression for the entanglement entropy [87] S(ρ A (t)) of the reduced density matrix ρ A (t) formed from the time-dependent global state ρ(t) = |ψ(t) ψ(t)|, which is now a mixed quantum state. Before we perform the computation, let us start by re-emphasizing that the total subsystem consists of two spatially disconnected regions: an interval on the left (I L ) and the corresponding (identical) interval on the right I R , together constituting At the level of the covariance matrix of A, we can clearly see this decomposition in the associated block structure. The relevant covariance matrix G A can be viewed as consisting of four blocks with respect to the LR decomposition in real space. The time independent blocks G L,L A and G R,R A originate from the thermal density matrix (describing each side individually), reduced with respect to I L and I R (respectively), and are therefore time-independent and equal to each other because we assumed symmetric intervals. The time-dependence of the TFD state then manifests itself in the mixed L(eft)-R(ight) blocks G R,L A (t) and G L,R A (t). In this way, one has [78,79,87] Here, Ω A is the symplectic matrix of the degrees of freedom belonging to A, | · | is the matrix absolute value, and the function s : [1, ∞) → [0, ∞) is defined as G A is the principal submatrix of the covariance matrix G that corresponds to A. In the same way, the Rényi entropies S q (ρ A (t)) for q > 0 can be computed by replacing the above function by s q : Of course, one recovers the familiar von Neumann entanglement entropy in the limit S(ρ A (t)) = lim q→1 S q (ρ A (t)). Furthermore, we find S 2 (ρ A (t)) = 1 2 log det G A (t), which follows from (222) with s 2 (λ) = log λ and tr log = log det. Appendix G provides some useful closed form expressions for the relevant covariance matrices of TFD states. 34 Note that in general, these are different from the normal modes of the full system. 35 This means taking the absolute values of the eigenvalues λi of G

Bounds on the entanglement entropy
Since our approach to computing entanglement entropy ultimately uses numerics to find the relevant symplectic eigenvalues, and we can therefore scan only a finite number of values of the underlying parameters including time, it is important to have additional guiding principles that constrain the problem at hand. A relevant question one can ask in this context is if it is possible to bound the maximal amount of entropy produced as the TFD is evolved in time.
We can use the subadditivity of the entanglement entropy applied to the decomposition (219) to upper bound S(ρ A (t)) as The discussion in the previous subsection clarifies that S(ρ I L (t)) = S(ρ I R (t)), and that it is given by calculating the von Neumann entropy of a single interval in a finite temperature thermal state. Therefore, the upper bound is time-independent and depends only on the values of m L, β/L, the total number of sites and the size of the subsystem. It should be also clear that when the number of sites N is sufficiently large (such as 1001 or 2001 sites, as we use in our numerics) and for large enough subsystems, subtracting from both sides of eq. (223) the initial (t = 0) entanglement entropy gives rise to a bound that is to a good degree independent of the total number of sites and the size of the circle L. The main advantage of the bound (223) in the context of our studies of the TFD state is that even if our numerics show a growing behaviour over a large range of times, we know that this growth will have to ultimately terminate since otherwise it would violate the inequality (223). For completeness, we would also like to mention two other useful inequalities for entanglement entropy involving the Rényi entropy S 2 for q = 2. The first inequality, similarly to eq. (223), holds for any quantum state and is a lower bound The second inequality is specific to Gaussian states and takes the form where N A is the number of bosonic degrees of freedom in the subregion A on each side of the TFD [88]. It can be derived using s(λ) ≤ log(λ) + (1 − log 2) for all λ ≥ 1. Note that, as opposed to inequality (223), the two bounds above are in general time-dependent. Note also that they can be viewed as providing a band in which the entanglement entropy necessarily resides, which is easier to calculate than the entanglement entropy itself, because we find a simple determinant in eq. (222) rather an expression involving eigenvalues. This will also be relevant for the asymptotic analysis in section 6.5, which is easier with a determinant than with the eigenvalues, for which we lack simple analytical formulas.

Quasiparticle picture of entanglement production
In our numerical studies in the next subsection we focus exclusively on the difference between the entanglement entropy at a given instant of time and the initial entanglement entropy normalized with respect to the thermodynamic entropy S th defined in eq. (202), Of course, the upper bound in eq. (223) now becomes ∆S(ρ A (t)) ≤ S(ρ I L (t)) + S(ρ I R (t)) − S(ρ A (0)) .
Let us also note that ∆S(ρ A (t)) is a UV-finite quantity. Our setup, following ref. [3], can be regarded as an unusual quantum quench involving two decoupled, albeit entangled via their initial conditions, subsystems. An important regime in quenches is the linear growth regime that occurs when the lattice spacing δ > 0 is much smaller than the correlation length set in our case by β > 0, which itself is much smaller than the subsystem size , i.e., δ β .
In the context of TFD states for holographic CFTs in (1+1)-dimensional Minkowski space, 36 ref. [3] observed in this regime a linear growth of ∆S(ρ A (t)) with the slope set by the associated thermodynamic entropy density. The growth terminates after a time equal to half the size of the interval, t = /2, with ∆S(ρ A (t)) remaining constant afterwards. These features follow from the fact that the holographic entanglement entropy [89,90] (see, e.g., ref. [91] for a review) is given by surfaces of minimal area (in the relevant case of AdS 3 these reduce to geodesics), and as t approaches /2 there is an exchange of dominance from geometric objects penetrating black hole interior to ones that remain outside. The latter case leads necessarily to a saturation, since the black hole exterior in the case of interest is static and equivalent geometric configurations contribute at every instance of time t ≥ /2. Of course, in this regime, the holographic entanglement entropy has precisely saturated the subadditivity bound (223).
In standard quenches involving a rapid global change in a local Hamiltonian and a single interval of length , the corresponding initial linear growth in ∆S(ρ A (t)) can be accurately captured using the quasiparticle picture [92][93][94][95][96][97]. In this picture, independent entangled pairs propagate in a ballistic fashion following an effective group velocity, and which is bounded from above by a Lieb-Robinson bound (which also holds in the regime of harmonic infinite dimensional constituents discussed in [95]). The linear increase in the entanglement entropy arises from quasiparticles that move out of the interval while their partner is still inside (or alternatively, from quasiparticles created outside the interval entering while their partner remains outside).
Indeed, Alba and Calabrese in ref. [74] (see also ref. [75]) go as far as providing a formula for the change in the entanglement entropy as a function of time -which is largely independent of the underlying model as long as it is integrable -by counting the quasiparticles with a given weight contributing to the entanglement entropy for a given interval. Within the time scale of /2, all the pairs created within the interval are distributed with one partner outside and the other inside. After this time, for some of these pairs, the second partner also leaves reducing the entanglement, in principle, but this effect is precisely counteracted by the increase produced by new quasiparticles coming into the interval. The net effect is a saturation of the entanglement entropy. At a quantitative level, this formula for standard quenches involving only one copy of a quantum many body system and homogeneous quenches is derived as follows: Suppose that we have a set of quasiparticles with labels n (for the case of free QFTs, this is simply the momentum label) with a function s n characterizing their density times their contribution to the entanglement entropy. The non-trivial insight of refs. [74]- [75] was that the function s n can be simply related to the thermal entropy density s th characterizing the global equilibrium state with all conserved charges equal to the ones characterizing the postquench set-up. Furthermore, it is assumed that pairs of quasiparticles are created, moving in opposite directions, each with velocity v n (as dictated by the conservation of momentum), without interacting with one another and that they are pairwise entangled. As a result, s n are treated as additive contributions to the entanglement entropy of a given interval with the rest of the system when one of the quasiparticles constituting a pair leaves the interval. Of course, the quasiparticle production is uniform through the interval because of homogeneity. With this set of assumptions, the problem becomes now a counting problem in which one needs to keep track of which quasiparticles are located within the interval at any given time.
Focusing on a single kind of quasiparticles of type n, the flux of quasiparticles across either end of the interval is proportional to v n t, and because the interval has two ends, the relevant total number is proportional to 2 v n t. 37 Given the above logic, the contribution to the entanglement takes the form t (2 v n s n ). This process lasts until the time t = /(2 v n ), when the pairs created at the center the interval reach the ends and the entanglement saturates. Afterwards, with the entanglement constant, the contribution of this quasiparticle species to the total is just s n . Since by assumption quasiparticles are independent, the formula proposed is Of course, the contributions linear in time in the above equations are the ones driving the linear growth of entanglement. This simple formula was successfully tested in refs. [74,75] in a variety of integrable models including free boson quenches. In the latter case, the authors considered theories where the boson mass is very large, i.e., comparable to the inverse lattice spacing (which was set to one there). As a result, eq. (229) does not account for gapless or nearly gapless zero modes. In addition, the equation does not account for finite size effects when the system is confined to a circle, rather than to an infinite line. It is possible to adapt the formula (229) to the case of the entanglement dynamics in the TFD state of free bosons living on a circle. In this case, the quasiparticle pairs are excitations of eigenmodes of a Hamiltonian with energies given by eq. (179) which move with group velocities given by 38 where n is an integer running between 0 and N (and of course, we assume N 1). In the previous sections we have regarded the TFD state as consisting of two copies of the field theory living on two geometrically distinct spaces. Alternatively one can regard the TFD as a state living on a single space, but in which two copies of the fields reside. The relevant quasiparticle excitations would be combinations of excitations of the left and right degrees of freedom, 39 however, the simple count presented above still holds for these (more complicated) combined excitations. The only change needed is to recall that the time we used in the full 37 The additional factor of proportionality needed in order to obtain the quasiparticle number is in fact the quasiparticle density which is encoded as a multiplicative factor in the function sn. 38 More precisely the group velocity is defined by differentiating ωp in equation eq. (189) with respect to p in eq. (190) and this introduces the extra factor of L 2π in eq. (230). Note that in the massless limit vn = 1 for all species of quasiparticles. 39 Naturally, those would be excitations of the L + R and L − R combinations of the fields, see eq. (22).
Hamiltonian evolution H L + H R for TFD states in the previous sections is not the parameter t but rather t/2, see eq. (30). Evaluating eq. (229) with this change yields where the relevant function s TFD n can be again deduced from the thermodynamic entropy density of this (double field) configuration In the above equation, the portion within the large parentheses represents a contribution to the thermodynamic entropy of the free boson system at the inverse temperature β from the mode n, cf. eq. (202), and the factor of 1/L renders this an entropy density, and the overall factor of 2 stands for the two copies (Left and Right) of the fields. Of course, this formula is valid only for times t < L − before finite size effects appear. We can very easily include these by carefully tracking quasiparticles on a circle leaving and re-entering the interval, which leads to the final formula which is periodic in time with period L/v n for any given species of quasiparticles 40 where frac denotes the fractional part, e.g., frac(3/2) = 1/2. Of course in the almost-massless case v n = 1 for all the different species and the entanglement production is approximately periodic, while for non vanishing value of the mass the different modes dephase and we obtain a picture that is not perfectly periodic (see fig. 21). Since eq. (233) is a phenomenological relation, it is interesting to see to what extent it quantitatively agrees with our set-up. Among other things, we will test it in the next section, where we study the dynamics of entanglement entropy numerically.

Numerical results
In the following, we evaluate ∆S(ρ A (t)), defined in eq. (226), numerically by using the techniques of subsections 6.1 for several values of the mass and temperature. 41 We work with each side being a circle of circumference L discretized into 1001 or 2001 lattice sites. We consider subsystems consisting of intervals on each side with lengths varying from = 0.1 L to = 0.5 L in increments of 0.05 L. Since ∆S(ρ A (t)) is a UV-finite quantity, comparing the results for different numbers of lattice sites, but for the same relative subsystem size, provides us with a test of the numerical convergence of our calculations. In fig. 20, we explore the linear regime and saturation of entanglement evolution on a circle for N = 2001 and β = 0.01L, for masses ranging from mL = 0.001 to 100, where we increase the mass by a factor of 10 in each step. We plot times smaller than L/2 where the system is 40 We would like to thank Pasquale Calabrese for pointing this out to us. 41 In fact we fix the values of the dimensionless parameters m L and β/L. not yet very sensitive to finite size effects (which will be explored in fig. 21). As expected in eq. (231), we observe a regime of linear growth that lasts until approximately t ∼ , followed by a saturation. For m β 1, the slope is equal with a good accuracy to the total thermal entropy density of two copies of the QFT (at inverse temperature β). Of course, this matches with the prediction of eq. (233), since for small m the group velocity is v n = 1 for all species of quasiparticles and the slope becomes the sum of entropy densities over all species, which is then simply the total entropy density. That is, with v n 1, the prefactor in eq. (231) is given by eq. (202) divided by the spatial volume, i.e., L. For the larger masses, where v n varies significantly in the mode sum, the prefactor should not have this simple form and in fact, fig. 20 shows the slope is significantly less than 2S th /L for the largest mass m = 100/L. For times larger than L − , the entanglement entropy becomes very sensitive to finite size effects, which manifest themselves as oscillations of periodicity L rather than a saturation. 42 This is depicted in fig. 21. For small masses, we observe a clear structure consisting of linear rise for a time t ∼ , followed by a plateau of width approximately L − 2 and then a linear decrease for another t ∼ . For large masses, the features of ∆S(ρ A (t)) change as time progresses, i.e., the oscillations are rather irregular. We observe that for large masses and not too large intervals, the maxima of the oscillations lie close to the bound (227), see fig. 22.
Another interesting feature of the entanglement growth is hidden in the plateaus for small values of the masses. If we restrict our attention to these flat regimes, we observe a logarithmic growth with, to an excellent degree, unit coefficient, see fig. 22, where we show the relevant fit to the quasi-plateau regions after the first period, i.e., for t > L. Fits to the data shown in figs. 20 and 21 at m = 0.001/L for < t < L − point towards logarithmic growth with a smaller coefficient. In the next subsection we will develop an analytic understanding of this logarithmic growth and demonstrate that for < t < L − it has a prefactor of 1/2. Let us also recall that a similar logarithmic growth was observed for the complexity on the circle in section 5.2.
A similar logarithmic growth of the entanglement entropy was reported earlier in ref. [98] for a global quench to a free massless bosonic quantum field theory on a circle, where it was attributed to the presence of an approximately gapless zero mode, namely the momentum mode with k = 0 in the massless limit. 43 The authors of ref. [98] claimed that an almost gapless zero mode does not lead to a ballistic propagation, as in the quasiparticle picture, but is instead of a diffusive nature. Since the times studied there were too short for the finite size effects to take effect, the logarithmic growth was observed with prefactor 1/2. Indeed, in the next subsection, we confirm this statement and demonstrate how one is able to predict the logarithmic growth (with prefactors 1/2 and 1 for the different regimes) by studying analytically simple subsystems with a single degree of freedom on each side of the TFD. We also observed that this logarithmically growing contribution is absent if we break translational invariance by considering the system on an interval with, for example, Dirichlet instead of periodic boundary conditions.
In fig. 23, we compare our numerical predictions with the Alba-Calabrese formula [74,75]  Larger intervals are related to those already in the plot since the TFD is a pure state. The dashed black lines correspond to linear fits passing through the origin. We see that for small masses the slope is approximately 2. The linear regime terminates (albeit not in a sharp way for larger values of the mass) at times approximately equal to the size of an interval. The smooth transition for larger masses is due to dephasing of the different kinds of quasiparticels with different group velocities. See also fig. 21 where a longer period of time is plotted in order to study finite size effects.  fig. 20 plotted over a larger range of times. One sees finite size effects, which induce a periodic behavior with periodicity L. The plateaus for small values of the mass follow a logarithmic growth due to the presence of a zero mode, which we study further analytically in the next subsection. For m = 10/L and with < 0.4 L, we observe that the entanglement entropy reaches rather close to the upper bound in eq. (227) after approximately t ∼ . The entanglement does not, however, saturate at this value for the times we considered, but rather keeps oscillating. In the right plot, for l < 0.4 L, we also plot dotted black curves presenting very good fits to the quasi-plateau regimes for t > L that follow the behaviour (constant + log m t)/S th with unit coefficient in front of logarithm and with different constants for the different subsystem sizes. We interpret the latter behaviour as a manifestation of an approximately gapless zero mode on the circle as we explore further in the next subsection.  [74,75] adapted to our setup, see eq. (233). In the most massive case we considered (left), the formula works remarkably well, especially for the smaller subsystems. In particular, it predicts the slope of the linear regime to be approximately 1.30, a fraction of a percent from the observed value. For the larger subsystems we suspect that the mismatch is partially due to numerical errors. However, for lower values of masses all the way to m = 0.001/L depicted here, the formula predicts the slope only to about 5% accuracy (already for m = 10/L), and does not lead to a logarithmic growth, due to not properly accounting for the zero mode which does not ballistically propagate. It would be interesting to further explore how to incorporate the zero mode effect into the Alba-Calabrese formula. adapted to our setup, i.e., with eq. (233). We observe a remarkably good fit for the highest value of the mass we consider, namely m = 100/L. For smaller values of the mass, even for m = 10/L, we observe that the Alba-Calabrese formula, as it stands, misses the slope of the linear regime by about 5% and of course, due to the built-in saturation, it does not fit the logarithmic growth for t , since it does not account for a nearly gapless zero mode behaving non-ballistically. It would be very interesting to be able to account for the zero mode behaviour in a generalization of eq. (229). Note that just adding ∼ log t would not work, because it behaves badly near t = 0. We leave these interesting issues for future work.
We can also use our derivation on the circle to obtain entanglement dynamics when the QFT lives on a line. This can be achieved by numerically studying a system on a circle with fixed lattice spacing δ and increasing the number of lattice sites N , while keeping the mass m, inverse temperature β and subsystem size fixed ratios of the lattice spacing as was done for the complexity in subsection 5.2. Numerics with increasing N reproduce more and more parts of a single curve which corresponds to the limit L/δ → ∞. Note that this pushes the first oscillation period on the circle to larger and larger times. Fig. 24 shows these results for small masses. Note that here we observe a logarithmic growth of the entanglement entropy even when the theory lives on the line. In contrast for the complexity in section 5, we saw that the logarithmic growth initially found for the TFD on a circle did not survive in the decompactification limit. This obviously happens also when N = 8001, but at later times (this is the dip in the grey curve). Taking larger N pushes finite size effects to later and later times, recovering more and more of the curve for a QFT on a line. As in fig. 22, we also observe a logarithmic growth, but now with coefficient equal to approximately 1/2 (dashed red curve), as reported earlier for standard quenches in ref. [98]. We interpret this as the effect of an almost gapless zero mode. This growth was harder to see before in fig. 21 since 1 2 log (t/β) becomes a good description only right before finite size effects kick in. Finally, we refer the reader to subsection 6.5 for further studies of the logarithmic growth with coefficients 1/2 and 1.

Logarithmic contributions to the entanglement from the zero mode
Both refs. [74] and [98] argue that the logarithmic growth of the entanglement entropy for intermediate times could arise from logarithmic corrections due to the presence of a zero mode, i.e., a mode of frequency m L −1 . In this subsection, we present an analytical argument which explains the logarithmic growth arising from the zero mode. The argument is based on overwhelming numerical evidence that the logarithmic growth with unit coefficient, as described in the previous subsection, persists for smaller subsystems, in particular for a subsystem A consisting of a single local degree of freedom (single site) on each side. We comment on more general situations later. In this case of a single site A, we derive below the full asymptotic behaviour of the entanglement entropy in the limit β, m −1 L, which will turn out to be given by where the second case offers a simplification of the third case under the condition t m −1 . The large time asymptotics given by oscillations with frequency m is not surprising, because we already knew that the entanglement entropy cannot grow arbitrarily due to the upper bound provided by the thermal state, cf. eq. (223). Only for times t m −1 will we see the logarithmic asymptotics, while for longer times oscillations with frequency m become visible. Another important regime is t < L, which shows the characteristic behaviour for a field theory on a line (L → ∞) rather than circle (L fixed). We can see that these predictions match the different regimes in our numerical results in fig. 25. We begin our argument by recalling from inequality (225) that for highly entangled states the von Neumann entropy S(ρ A ) rapidly approaches S 2 (ρ A ) + 2N A (1 − log 2), where S 2 is the second Rényi entropy. In fact, we know that our error scales as exp(−2S 2 (ρ A )), which means it becomes exponentially small for highly entangled states, as explained in refs. [88,102]. The Rényi entropy of order 2 can be written as the determinant which we already discussed in the context of equation (220). We will apply this formula to a single degree of freedom on each side of the TFD. To do so, let us derive the covariance matrix for a subregion consisting of a single site at position x on both sides of the TFD with respect to the dimensionless variables ξ a k = (q L k ,q R k ,p L k ,p R k ), which are related to the discretized field variables in momentum space according to eqs. (37) and (176). The covariance matrix is derived by first decomposing the complex modes into real modes, as explained in appendix D. Then for each the real degree of freedom, we write the covariance matrix (104) with the substitutions ω → ω k in eq. (179) and λ → ω k L -see eq. (268), 44 It may be surprising that this covariance matrix is real, even though the underlying operator ξ a kξ †b k +ξ †b kξ a k is not Hermitian. However, if we rememberξ †a k =ξ a −k , it follows from the translational invariance of the state that the expectation value should be invariant under the swap k ↔ −k, which implies thatG a,b k is real. Moreover, translational invariance also ensures that there cannot be any other cross correlations, so that the total covariance matrix is block diagonal. We can go to the position basis ξ x = (q L x , q R x , p L x , p R x ) by applying the inverse Fourier transformation toξ a k given by Conjugating this equation gives rise to an extra minus sign in the complex exponent. With this in hand, we can write the covariance matrix in position basis as where we have exploited the fact that the covariance matrix in momentum space is block diagonal. When setting x = y, this 4 × 4 matrix describes the quantum correlations and entanglement in a subsystem consisting of a single local degree of freedom on each side (left and right) of the theory. We will use this explicit representation to study the asymptotics of the entanglement entropy analytically, which complements our complexity studies of the TFD state. In particular, we will be able to identify the characteristic contribution of the zero mode (m → 0 limit) to the production of the entanglement entropy.
Restricting to x = y yields the following 4 × 4 covariance matrix (239) 44 Here we use the indices a, b ∈ {L, R}. Further, let us note that, in contrast to the complexity, the entanglement entropy has no connection to the auxiliary gate scale. This means that we can choose any convenient scale to make the covariance matrix dimensionless. Effectively, our substitution λ → ω k L replaces ωg → 1/ √ L δ or equivalently µg → L -see eqs. (184) and (187).
where we have defined λ k = ω k L. Schematically, the structure of this matrix is whose determinant can be explicitly computed as Leading order One can observe that only the first two terms are leading order for t δ = L/N in the limit under consideration, namely N → ∞, mL 1, and β/L 1. The three ingredients of the leading order piece are then given by We may use eq. (20) to show that sinh(2α k ) = 1/ sinh(βω k /2) and for k N and β L we may effectively approximate sinh(2α k ) ∼ 2/(βω k ). Using these equalities we are able to approximate the function F 1 (t) for t > 0 as follows where we have used the fact that ω k = ω N −k ∼ 2πk/L for k N and m δ −1 , see eq. (179), and ignored differences introduced for larger values of k where both the terms in the above sum as well as sinh(2α k ) are highly suppressed (the latter exponentially) in the limit N → ∞. Furthermore, the sum over cos(2πkt/L) can be recognized as the Fourier series expansion of t(t − L) + L 2 /6 over the interval [0, L]. Due to the periodicity of the Fourier series, we define the function P L (x) as the periodic function that repeats its values from the interval [0, L] on all other intervals [nL, (n + 1)L]. Putting things together and using we find where we have neglected the square of P L (t(L − t)), and we approximated cos (mt) ∼ 1 in its product with P L (t(L − t)). These approximations are consistent as long as we are not too close to points where mt is an integer multiple of π, since we are working in a regime where mL 1. Consequently, the entanglement entropy is given by We can simplify this expression in three time regimes, namely δ t < L, L < t m −1 , and L t. In these regimes, we find the asymptotic behaviour of the entanglement entropy to be given by This expression is well-matched by the numerical results in fig. 25. Note that the last case becomes inaccurate around its singularities, where mt is a multiple of π. For large N , we can use the asymptotics c 2 ∼ 2L/β to obtain further simplifications. We should emphasize that this expression will fail to describe the entanglement accurately in the regime t < δ, in which the precise oscillations of the various modes will determine the time evolution. This discussion gave a precise account of the logarithmic entanglement growth as a function of time for a single mode in I R and I L , respectively, constituting the subsystem A. To argue about more general subsystems composed of many modes, let us first again state that due to the aforementioned upper and lower bounds to the von-Neumann entanglement entropy, and since we are in what follows not that much interested in pre-factors, we can basically express the entanglement entropy freely in terms of the first and second Rényi entropies, respectively. Let us denote with σ(t) the reduced state of a single pair of sites in R and L, so that eq. (248) applies to S(σ(t)). Hence, for a general region A = I L ∪ I R , invoking the subadditivity of the von-Neumann entropy and making use of the above bound, we have as a rigorous upper bound Hence, we find that S(ρ A (t)) for the entire multi-mode region A is again upper bounded by a constant plus a function that grows logarithmically in time, until saturation. At the same time, since the initial state contains short-ranged correlations, one expects to be able to decompose the system A in largely uncorrelated regions, giving rise to each logarithmically growing contributions, with little error. That again has the consequence that one expects the behaviour S(ρ A (t)) ∼ c(constant + log m t) for some constant c > 0. Note the different normalization of the time axis on each plot. Note also that for such small subsystems finite size effects manifest themselves differently than in the quasiparticle regime covered by fig. 21, i.e., there are no pronounced oscillations. Note that we do not divide by the thermodynamic entropy S th here, since the entanglement entropy of two sites does not have a well-defined continuum limit.
To close this discussion, let us emphasize that contributions similar to eq. (248) have already appeared earlier in section 5.2 in our discussion of the κ = 2 complexity on a circle, see eq. (196). There, the contribution of the zero mode had a fixed numerical prefactor and we argued that in the limit where we decompactify the circle, it is subdominant with respect to the other modes contributing on the order of the thermodynamic entropy. As a result, unlike for the entanglement entropy, the zero mode played no role in the dynamics of complexity for quantum field theories in Minkowski space.

Discussion
In this work, we have studied the entanglement dynamics as well as the definition of circuit complexity for time-dependent TFD states (1) in free, non-interacting scalar QFTs. In our investigations, we have been directly motivated by the results of the proposed holographic duals of complexity and entanglement evaluated in eternal black hole backgrounds in asymptotically anti-de Sitter spacetimes, which are the dual manifestation of TFD states in a class of strongly coupled quantum field theories with large number of degrees of freedom.
In defining complexity, one needs to specify the elementary operations or gates whose number one counts according to some cost function which assigns a length to the circuit. The crucial insight in this respect is that for free systems, the TFD state introduced in section 2.2, the vacuum state, and simple, spatially disentangled reference states considered previously in refs. [23,51] -are Gaussian. As a result, the natural choice of gates is the set of all Gaussian transformations (i.e., the full set of Bogoliubov transformations), which for N bosonic modes 45 is the symplectic group Sp (2N, R). This is a larger set of elementary operations than considered previously in refs. [23,51], but at the same time it is a small subset of all bosonic unitary transformations. A technical novelty of our work with respect to earlier works is phrasing the complexity calculation in the language of symplectic transformations acting on covariance matrices, which is the most efficient way of dealing with the full group of symplectic transformations on the Gaussian states, see section 3. An additional new feature with respect to ref. [23] is the need to make some of the new generators of the symplectic transformations dimensionless, see eqs. (37) and (91), which leads to the appearance of an additional scale in the problem -the so-called gate scale ω g .
Following the ideas of Nielsen and collaborators [70][71][72], we have represented the circuits which prepare the desired state as paths in a Riemannian geometry of the underlying Lie group, i.e., Sp(2N, R) with appropriate boundary conditions. The optimal circuits then become geodesics in this geometry. An important point to mention is that we are not just interested in representing or approximating a particular unitary operator as a circuit, but rather in considering a whole class of circuits, all of which will give rise to the desired target state. There is no unique choice of a circuit geometry. 46 In the bulk of the text, we focused on the distance measure induced by the Hilbert-Schmidt inner product between two infinitesimally separated group elements, see eq. (109) with G = g = 1. As we demonstrate in appendix F, when λ R = 1 we can take advantage of a special property of this geometry to show that the optimal circuit does not mix different normal modes, which significantly simplifies the complexity calculation. Furthermore, in this situation where the scale controlling the on-site correlations in the reference state (i.e., the reference state scale) is associated to the gate scale, the expression for complexity is a particularly simple function of the relative covariance matrix between the target and the reference state, see eqs. (86) and (88). In this case, in the basis in which both the reference and target state covariance matrices are diagonal, the optimal circuit amounts to squeezing each of the individual normal modes at a constant rate. This is reminiscent of earlier results obtained in refs. [23,51]. However, we stress that while both of those references only optimized their respective circuits within some restricted gate set (i.e., GL(N, R) in ref. [23] and SU(1, 1) N in ref. [51]), we have shown that our shortest geodesics really optimize the corresponding circuits within the full family of Sp(2N, R) generators, as proven in appendix F. Unfortunately, the simplicity of the optimal circuits is lost when the reference state scale and the gate scale are not simply related, i.e., λ R = 1. In this case, for each normal mode, one needs to find the optimal path numerically, see section 4.5. We must add here that to make the numerical problem tractable, we have still assumed that the optimal circuits did not mix the normal modes (although this was only rigorously proven for the case of λ R = 1). Working in the normal mode basis implies that our circuits are constructed using spatially non-local gates. This is not unlike the action of MERA tensor networks [103], where entanglement at different scales is injected by gates in different layers of the circuit. 47 We might also add that a recent holographic work [36] suggests that if the holographic complexity proposals indeed capture some variant of a circuit complexity, then the optimal circuits must be utilizing spatially non-local gates.
Before we summarize the main physics lessons, let us remark that the complexity defined by the Riemannian geometry can be thought of as using the L 2 norm, i.e., the F 2 cost function, when counting the number of gates, see eq. (8) (right). The results of refs. [23,51] show UV behaviour closely aligned with the holographic results [29] when counting the gates using the L 1 norm, i.e., the F 1 cost function, see eq. (8) (left) and the discussion below. Such counting leads to a much more complicated minimization problem, and in the present article we adopted the proviso from refs. [23,51], i.e., we used the L 1 norm to express the length of a circuit optimized with respect to the L 2 norm. In fact, for the preparation of the ground state, ref. [23] showed that this circuit also optimized the L 1 norm. But their proof does not extend to the present situation and so our results should be seen as an upper bound on the F 1 complexity. Setting aside the problem of optimality in the L 1 norm, there is also the issue of basis dependence, i.e., the length of the path will depend on the basis chosen for the generators. While this may be seen as a shortcoming of this norm, we actually took advantage of this basis dependence in evaluating the complexity of formation in section 5.3. There we saw that using the physical left-right (LR) basis (but not the diagonal basis) produced a complexity of formation that compared well with the holographic results [27]-see further discussion below.
Having sketched the mathematical underpinnings of our work, let us now discuss the physics of complexity in the Gaussian TFD states that we have uncovered. We begin by reiterating some of the key lessons stemming from studies of the holographic complexity proposals. First, the complexity of the vacuum state diverges as the spatial volume occupied by a boundary quantum field theory measured in the units of the UV cut-off with a non-universal prefactor [29]. Furthermore, for the CA proposal, there is a possibility that this leading divergence is also multiplied by a | log( ct /L)| factor arising from the detailed structure of null boundary terms in the gravitational action. 48 Second, calculating holographic complexity for the TFD state (1) at t L = 0 = t R , and then subtracting twice the vacuum complexity, one finds that the remainder is proportional to the thermodynamic entropy, i.e., the entanglement entropy between the left and right CFTs, see eq. (207) here and ref. [27] for details. 49 This remainder term is called the complexity of formation, i.e., the extra complexity needed to prepare the left and right CFTs into the entangled TFD state, instead of a product state of two vacuum states (i.e., the β → ∞ limit of the TFD state). Finally, for the holographic TFD states both proposals predict a linear growth of the holographic complexity at late times, which is perhaps their most well-known feature.
Following refs. [23,51], we always take as our reference state a spatially disentangled state. Furthermore, we focus on the L 1 norm, i.e., F 1 cost function, since evaluating the complexity of the vacuum state in the normal mode basis then yields a result with the same leading UV divergence as found in holographic complexity. However, it is also interesting to consider the complexity evaluated with the κ = 2 cost function since in this case the geodesics correctly describe optimal circuits. 50

Complexity of formation and time evolution
Let us begin by discussing the complexity of formation in the case when the gate scale is equal to the reference state scale, i.e., λ R = 1, for which the results have been presented in section 5.3. As we noted above, using the L 1 norm and physical LR basis, our result for the complexity of formation compared well with the holographic results [27]. Three general features shared with the holographic result are, ∆C (LR) 1 is UV finite, it is independent of the reference scale µ, and it is positive. 51 Furthermore, if we set m = 0 to emulate a CFT, then ∆C (LR) 1 is proportional to the entanglement entropy between the left and right CFTs as found with holographic complexity when d ≥ 3. However, there is a difference in the dependence of the proportionality constant on the number of spacetime dimensions compared to holography, e.g., see eqs. (206) and (207). For two dimensions, our QFT complexity of formation is also proportional to the entropy and so grows with increasing temperature, i.e., ∆C (LR) 1 ∝ V T for d = 2. Of course, this deviates from the holographic results where, as we explained above, ∆C holo ∝ c for d = 2, and so we might also classify this as part of the previous discrepancy in the constant of proportionality (between the complexity and the entropy), i.e., this constant vanishes for d = 2 in holography while it remains finite in our QFT calculations. We might note that this coefficient had a very different dependence when the κ = 2 measure was applied, as shown in fig. 11, but that it still remained nonvanishing for d = 2 in this case as well. As a final comment, let us reiterate that when evaluating the complexity of formation to produce the results which compared well with holography, we have taken advantage of the basis dependence of the L 1 norm and in particular, evaluated it using the LR basis. In section 5.3 -with further details in appendix E -we have shown that using the diagonal basis (along with reasonable choices for the reference and gate scales) leads to the complexity of formation being suppressed with respect to the entropy at low temperatures (with respect to the cut-off Λ), peaking and decreasing as the temperature approaches the cut-off, 52 see figs. 10, 26, 27 and 28. Therefore, the good agreement with holography produced with the LR basis might actually be a clue as to the microscopic rules that implicitly define the holographic proposals for complexity.
Turning to the time-dependence of the TFD state, the complexity in our quantum field theory calculations is a combination of contributions from individual normal modes, each of which exhibits oscillatory behavior with a frequency set by the normal mode frequency, e.g., see eqs. (149), (151) and (152). While the oscillations are all aligned at t = 0, they quickly dephase afterwards and as a result, summing over the modes yields a result in which contributions from different normal modes average out at times of the order of the inverse temperature β. This leads to a quick saturation in the TFD complexity of free QFTs, a feature which is very different from the late-time linear growth found for holographic complexity. Depending on the choice of parameters, the complexity can either grow or decrease after t = 0 before saturating with main features of the process independent of the dimensionality of the QFT. The dependence of the transient early time behaviour on µ is comparable to the the CA proposal. However, this issue is evaded if we choose the reference frequency µ to be proportional to the cutoff scale, e.g., see the discussion in refs. [23,29]. 51 A priori, ∆C > 0 is not an obvious result since while in preparing the TFD state, one introduces extra entanglement between the left and right copies of the QFT, one does not introduce as much long-range entanglement in either copy as one would in the vacuum state. For further discussion, see appendix E of ref. [27]. 52 With an exception of d = 2, where the ratio is the largest for temperatures approaching the cut-off; one should note that this is an unphysical regime to consider. sensitivity of the initial transients to the ambiguous scale arising in the null boundary terms in the CA proposal, similar to the aforementioned logarithmic divergence for the vacuum complexity [38,39].
Since the late-time growth of holographic complexity was one of the hallmark features of the CA and CV proposals, it may seem somewhat disappointing that this feature is not recovered in our present complexity calculations. However, we must consider that the QFT for which we are calculating the complexity is a free theory, while the boundary CFTs which are described by holographic complexity are strongly coupled theories. In fact, the time evolution of the complexity of the TFD in these theories is predicted to initially show linear growth, i.e., ∂ t C ∼ M , and to saturate only on a time scale of the order e S BH [13,105,106]. This expectation is based on the fact that these CFTs are fast scramblers and the evolution of the TFD state explores a larger and larger region within the Hilbert space of the boundary theory. Of course, this is clearly not the situation in the present case of a free scalar field theory. In fact, our complexity calculations have relied on the fact that the initial TFD state and time-evolved TFD are Gaussian states. Hence, the time evolution of the TFD state in our free scalar theory is confined to the submanifold in the full Hilbert space describing Gaussian states. From this perspective, it is perhaps not so surprising that we found that the complexity saturates on the thermal time scale. One could construct far more "complex" states by replacing the regularly spaced phases, (n + 1/2)ωt in eq. (30), with random phases θ n . Furthermore upon integrating out one of the copies in these new entangled states, we would still recover the desired (time-independent) thermal density matrix. Unfortunately, while we expect the random phases increase the complexity of these states, we do not yet have the technology required to actually perform concrete calculations of their complexity.
One is tempted to draw here an analogy to the early studies of quark-gluon plasmas using holography, e.g., see ref. [107] for an overview of these efforts. Such studies have revealed that the thermodynamics of N = 4 super Yang-Mills theory at vanishing 't Hooft coupling constant differ only by 25% from the thermodynamics of the same model in the holographic regime of infinite coupling [108]. In our case, the complexity of formation for free QFTs behaved in a qualitatively similar way to the results of the holographic complexity proposals. However, the time-dependent properties of complexity in free QFTs turn out to be very different from the predictions of holographic complexity proposals. To close the analogy, it is well known that the time-dependence of holographic QFTs is very different from their weakly-coupled cousins, for example the value of the famous shear viscosity to the entropy density ratio is parametrically different between the two [109]. 53 Altering the ratio of the reference state scale to the gate scale away from unity, i.e., choosing λ R = 1, does not change our main conclusions here in a significant way (e.g., as expected on general grounds, the complexity still saturates at times of the order of β), but it nevertheless leads to certain interesting effects. As shown in fig. 19, when λ R is either increased or decreased away from one, the complexity generally appears to decrease. This decrease seems to be a universal behaviour, which was already observed in the contributions of the individual modes in fig. 5. Hence, it appears that our results for the time dependence 53 One might wonder what features are expected in holographic complexity once more general features are added to the spacetime, such as higher curvature corrections. The investigation in this direction so far suggests that for holographic CFTs in Minkowski space as considered here the CA late time growth rate does not receive direct corrections for the Lovelock class of theories [47] (see refs. [25,44] for discussion focused on Gauss-Bonnet gravity). It is an interesting question as to whether for more generic higher curvature theories this feature persists, or whether that is due to the special properties of the Lovelock action. of the complexity with λ R = 1 set an upper bound on the complexity with these more general parameters. Furthermore, we should keep in mind that fig. 19 shows the complexity at time t relative to the complexity at the initial time. Given the results in fig. 4 for a single mode, we should also expect that varying λ R will reduce the initial complexity of the QFT by reducing the contributions for modes in a particular range of frequencies.

Comparison with entanglement dynamics
One of the most interesting aspects of the holographic proposals is the difference in the time evolution of the complexity and the entanglement entropy, as we have described in the introduction. This has motivated us to compare the time-dependent complexity with the behaviour of the entanglement entropy in time-dependent TFD states in free QFTs, where we have focused on 1+1 dimensional systems. We have studied the entanglement entropy associated to a subsystem consisting of two equal intervals on each side of the timedependent TFD state. This setup corresponds to the free-field, non-interacting analogue of the holographic studies reported in ref. [3]. For convenience, we have worked with a theory on a spatial circle and considered two cases: first, keeping the size of the circle and physical parameters fixed, while making the lattice denser, and second, keeping the lattice spacing and physical parameters fixed, while increasing the number of lattice sites in order to approach the field theory on a line.
Our investigations have revealed an expected linear growth and saturation pattern, usually present in entanglement of quenched systems, as well as oscillations for systems on the circle with periodicity given by the circle size. In addition, for small masses we have observed a logarithmic growth associated with the presence of a zero mode, both on the circle and on the line. This is in contrast with our findings for complexity on the line where a logarithmic growth was absent, see section 5. The complexity of the TFD on the circle does exhibit a logarithmic growth. This growth, however, does not survive the decompactification limit. We note however that the logarithmic growth of the entanglement cannot last forever given the upper bounds (223) (universal) and (225) (specific to Gaussian states) and the entanglement entropy has to eventually saturate.
In order to understand the linear phase of entanglement production, we related our results to previous studies of the entanglement entropy after global quenches in a wide class of bosonic systems. Building on earlier work [92][93][94], it has become clear that the linear growth of the entanglement entropy in time can be largely explained using a quasiparticle picture where pairs of entangled quasiparticle excitations are created and propagate ballistically at their group velocity in opposite directions [74,[95][96][97]. Initially, the flux of quasiparticles coming in or going out of the interval produces to an increase in the entanglement entropy and a careful accounting shows that the entanglement rises linearly until times t ∼ where it saturates at a value which is proportional to the thermal entropy of the system. When studying the entanglement entropy on a circle, we find that it exhibits recurrences with periodically equal to the circle circumference. This can be related to quasiparticles going around the circle and reducing the correlations between the subsystem and its complement again once they meet on the opposite side of the circle. These features becomes particularly sharp in the massless limit where the quasiparticles effectively move with the speed of light [95], while for larger masses, the different quasiparticles have different group velocities and their contributions to the entanglement entropy quickly dephase. As a result, the variation of the entanglement entropy over time periods of ∆t ∼ L are greatly reduced.
Overall, we find that the effective model that has been proposed in ref. [74] provides a surprisingly accurate description of the linear slope and quasi-periodic behavior for large masses when adapted to our setting, see eq. (233). For smaller masses, we identified a logarithmic correction to the periodic behavior. Such a behavior was previously observed in refs. [74,98] and attributed to the zero momentum mode in the massless limit. In fact, this contribution was purposefully removed through appropriate boundary conditions in ref. [74]. However, in the present work, we are able to present a largely analytical formula of this logarithmic growth, based on the determinant of submatrices of the covariance matrix, a contribution that we are convinced is interesting in its own right.
Specifically, to make an analytical analysis of the zero mode contribution feasible, we first focused on a subsystem consisting of a single lattice site on each side of the TFD state. We then argued that the observed asymptotics contribute additively to larger subsystems, due to the fact that the zero mode is completely non-local and therefore affects local subsystems in a similar way, regardless of their size. Our analytical study is based on a well-known relation between the von Neumann entropy and the Rényi entropy of order 2 that can be efficiently computed as the determinant of the covariance matrix for Gaussian states. In the case of a single site on each side of the TFD, this covariance matrix is a 4-by-4 matrix whose time dependence can be analyzed analytically. For a circle of circumference L and a small mass m, we find three distinct regimes, namely an initial logarithmic regime given by S(ρ A (t)) ∼ 1 2 log t/L + const for t < L, a second logarithmic regime given by S(ρ A (t)) ∼ log t/L + const for L < t m −1 and finally an oscillating regime given by S(ρ A (t)) ∼ log sin mt + const when t is of the same order as m −1 , see eq. (248). Note that we intentionally suppress the linear regime by choosing a subsystem size that vanishes in the continuum limit, so we can fully focus on understanding the logarithmic contribution analytically. In the non-compact limit, L → ∞, we expect to only see the initial contribution of S(ρ A (t)) ∼ 1 2 log t/L + const. Let us emphasize that to the best of our knowledge, this asymptotic analysis provides one of the first analytical insights into the zero mode contributions. However, understanding its additive contribution for larger subsystems, i.e., those that do not vanish in the continuum limit, will require further work if one wants to go beyond the rough analysis and heuristic arguments presented at the end of subsection 6.5.

Relation to other works
We note that the question of the complexity of the TFD state within a free scalar field theory has been studied previously in refs. [52,53]. However, our methods differ in a variety of ways. For example, these works have not considered unitary circuits to prepare the target state; they used a different reference state, i.e., their reference state was two unentangled copies of the vacuum state; they used a very restricted gate set, i.e., their generators formed an [Sp(2, R)] N algebra compared to Sp(2N, R) here; 54 the gate scale was implicitly set by the frequency of the individual modes, i.e., ω 2 g = M ω k ; and finally, they chose an unconventional cost function which was different than any of the cost functions considered in our paper. Given these extensive differences, it may seem remarkable if any of our results were to match with those in these earlier references. However, one finds an interesting parallel between our complexity of formation for a massless field (see eqs. (200) and (206)) and the complexity of the TFD relative to an unentangled product of two vacuum states evaluated in [52], in that both results are proportional to V T d−1 . We might add that for the special case of t L = 0 = t R , i.e., the TFD without any complex phases, the measure in ref. [52] coincides with the F 1 cost function, but their results do not match ours because of the other differences in our two approaches described above, e.g., they would not produce the same complexities for a massive scalar field. Rather, the simple behaviour with ∆C ∝ V T d−1 seems to be a result of dimensional analysis. If we assume that the complexity is extensive, then certainly in a situation where the temperature is the only other dimensionful parameter, we must find ∆C ∝ V T d−1 (recall that we found similar results for both the F 1 and κ = 2 cost functions).
Our results for the complexity of the time-dependent TFD state differ quite dramatically from those found in ref. [53]. The latter reported a linear growth of complexity at late times (when applying the Nielsen approach), a behaviour which stands in stark contrast with our findings, i.e., the saturation of complexity after roughly the thermal scale. The key difference in our approaches leading to this discrepancy is that ref. [53] adopts the unusual cost function introduced in ref. [52]: compare eq. (B.9) in ref. [53] with our eq. (8). With this choice for a unitary involving exponentiation of a phase factor, the cost function becomes arbitrarily large as the phase increases beyond 2π, which directly leads to the aforementioned growth. We regard it as physically incorrect to assign complexity growth to a periodic behaviour, and indeed in the approach pursued in the present paper the complexity for each normal mode is simply periodic, as are the relevant unitaries.
We might add that it is straightforward to adapt our analysis to consider the complexity of the (time-dependent) TFD state relative to the reference state being two unentangled copies of the vacuum state, as in refs. [52,53]. For each momentum mode, the reference state analogous to eq. (27) would be replaced by the vacuum, as in eq. (26). As a result, we would simply substitute ω 2 R → M ω k for each mode. The single-mode calculations would proceed as in section (4.5) with the substitutions (λ, λ R ) → (1, λ = M ω k /ω 2 g ). It would be interesting to investigate this relative complexity further, both for the individual modes and for the full field theory. However, as noted above, refs. [52,53] go one step further and set the gate scale for each individual mode to be ω 2 g = M ω k . 55 If we added this choice in our approach, we would also set λ R = 1 and our analysis would reduce to that in section 4.4 with λ = 1. In this case, eq. (142) yields ρ 1 = α and θ 1 = π 2 − ω k t (with τ 1 = 0), and hence the complexity for the individual modes is simply a constant, e.g., C κ=2 = ρ 2 1 = α 2 . Furthermore, integrating over all the modes would yield a total complexity which precisely matches the expressions which we gave to the complexity of formation, e.g., C κ=2 would be given by eq. (215). Hence the total complexity would be time-independent, and of course, our analysis still does not yield the linear growth found in ref. [53].
Ref. [53] also uses the Fubini-Study approach of ref. [51] to assign a complexity to the TFD state of a free scalar, again relative to the product of two vacua. With this approach, they report that the complexity initially grows but saturates in roughly the thermal scale. However, we still disagree with their calculations in this case. This result should correspond to considering eq. (144) together with eq. (142) evaluated for λ = 1 (and λ R = 1), which again leads to a time-independent answer. While circuit complexity is a priori a different notion from the complexity of states, in the case of free QFTs (or more generally, for Gaussian states which are prepared with only squeezing gates 56 ) they agree, see refs. [23,51,111], which motivated us to look into ref. [53] in more detail. And indeed, upon a closer examination one notices that eq. (4.32) in ref. [53], and following from it eq. (4.34), should have a phase factor pulled out from them. When this correction is included, the complexity remains constant throughout the time evolution, as suggested by the argument above.
The time-dependence of the TFD state discussed in the present manuscript is similar to a quantum quench in a free quantum field theory (in this context, evolution of a ground state under a Hamiltonian whose mass term becomes time-dependent for a certain period). Such set-ups have been discussed in the context of cMERA in refs. [112,113] and it was reported there, in the language of ref. [51], that the Fubini-Study distance along the RG scale direction of quantum circuit modelling a state after a quench exhibits a linear growth, in contrast with our claims. However, there is no contradiction since the results of refs. [112,113] can be regarded as measuring the circuit depth of a particular unitary rather than measuring the circuit depth of the optimal unitary and, as we demonstrated here, it does not exhibit any long-time growth. 57 When it comes to the complexity of quench setups in free QFTs, studies reported in recent ref. [58] define complexity using translationally-invariant Gaussian circuits and yield results that exhibit qualitative similarity to our results, i.e., a short period of transient effects following a quench and a subsequent saturation or mild oscillations of complexity. Finally, ref. [59] demonstrated that complexity in free quantum field theory defined using some of the machinery introduced in the present manuscript can obey scaling relations as a function of the quench rate in quenches passing through a critical point, similar to an earlier analysis of one-point functions [114,115] and entanglement entropy [116].

Open questions
Perhaps the most interesting open question that our work raises lies in defining circuit complexity in interacting QFTs in a calculable manner. This is, of course, related to a rather quick saturation of circuit complexity in free QFTs that we observed here, which is in stark contrast with the results of the holographic complexity proposals. It would be very interesting to understand if recent attempts of refs. [117,118] to include interactions in the continuous version of MERA (cMERA) [84] can help in addressing this problem. See also the recent work [65].
Another important problem is to understand if there are relations between circuit complexity as considered here and complexities defined using the Fubini-Study metric [51] and the path-integral optimization [60,61,119]. Free QFTs provide what we believe is a fruitful testing ground in which all three approaches become computable.
In the realm of Gaussian circuit complexity, it would be fascinating to explore more systematically the properties of optimal circuits in the presence of non-trivial penalty factors. An interesting question is whether insisting on (quasi-)locality of gates (even if it might not be the case in holography, see again ref. [36]) can significantly alter the results of existing complexity calculations, such as reported here or in refs. [23,51,55,56,58]. Such results can also have an interesting application in the field of tensor networks, namely if MERA [103] or cMERA [84,117,118], as tensor networks are in some sense an optimal way of preparing critical ground states from product states.
Taking a more quantum information perspective, one can quantitatively relate the notions of complexity with those of the entangling power of quantum gates. The complexity counts the 57 In particular, the time dependence of such unitaries is heavily dependent on the time dependence of a phase θ k (t) that does not change the state, see, for instance, eq. (3.31) in ref. [113] number of gates from a specific gate set needed to prepare a quantum state. The entanglement over any cut can be upper bounded in terms of the number of gates that are supported nontrivially on both sides of the cut, discounted by the entangling power of these specific quantum gates. In a continuum limit, the latter can be captured in terms of entanglement rates [120]. The quantitative connection between upper bounds to entanglement over cuts quantified in terms of entanglement rates and circuit complexity will be explored in detail elsewhere.
Finally, our work raises many interesting open problems regarding nontrivial gate scales, i.e., λ R = 1. For instance, since every mode has a more intricate minimization procedure, it is possible that the vacuum UV divergence will get modified as well. Understanding whether the case λ R = 1 still provides an upper bound, and exactly how the vacuum divergent terms get modified, could shed light on the nature of the gate scale, as well as whether holographic complexity has features which mimic this parameter. In addition, it would be interesting to directly investigate how the circuit gets modified once the geodesics deviate from the τ = 0 plane. It is our hope that the present comprehensive analysis stimulates such further work. A

B Unitary decomposition of the TFD state
In this appendix we introduce a useful decomposition formula which allows us to express the TFD state as a unitary operation on the product of vacuum states. We begin with the following operators 58 which satisfy the algebra The decomposition formula taken from appendix 11.3.3 of ref. [121] reads The requirement that the operator in eq. (252) be unitary implies α + = −α * − and ω ∈ R. For the special case where ω = 0 and α + = −α − = α ∈ R we obtain Now acting on the vacuum state |0 L |0 R , we have Hence if we act with the expression in eq. (253), we obtain Hence, if we identify tanh α = exp(−βω/2) , which also implies cosh α = (1 − e −βω ) −1/2 , then eq. (255) reproduces the desired thermofield double state (17)- (18). This formula also extends to the complex case with where z = α e −iωt . As before, tanh α is given by eq. (256) so that eq. (257) reproduces the time-dependent thermofield double state (30)- (31). Note however, that this expression does not capture the overall phase factor exp(−i ωt/2) in eq. (30).

C Matrix generators for Sp(4, R)
To illustrate the ideas of subsection 3.5, let us consider Sp(4, R) as it will be useful for the discussion in section 4. We will be using the L(eft)-R(ight) basis associated with the canonical variables ξ = (q L , q R , p L , p R ) The matrix generators for the Sp(4, R) are given by eq. (95). We can split them to the Sp(2, R) subalgebra acting on the left oscillator only the Sp(2, R) subalgebra acting on the right oscillator only

D Comments on bases
In this appendix, we comment on a few details with respect to our description of complexity in the scalar field theory. First recall from eq. (175) that the Hamiltonian describing the lattice regulated field theory can be written as which takes the form of a chain of N coupled harmonic oscillators. Implicitly, we are choosing periodic boundary conditions here with Q N +1 = Q 1 and P N +1 = P 1 . We wish to note that because the initial theory involved a real scalar field, the P a and Q a describe real degrees of freedom, i.e., P † a = P a and Q † a = Q a . However, the next step in our analysis in section 5.1 was to pass from the position basis along the chain to the normal mode basis by Fourier transforming as in eq. (178). At this point, the variablesP k andQ k are no longer real. 59 Rather, we haveP † k =P −k andQ † k =Q −k . This constraint indicates that the positive and negative momentum modes are mixed, e.g., with the canonical commutation relations [Q k ,P −l ] = i δ k,l . In other words, the two complex degrees of freedom labeled by (Q k ,P k ) and (Q −k ,P −k ) only contain two real degrees of freedom. This mixing is also evident in the corresponding Hamiltonian (180), which takes the form Of course, we can also define Hermitian combinations of these operators. If we assume that N = 2n + 1, we would havẽ where 1 ≤ k ≤ n in all of these expressions. Now (Q R,k ,P R,k ) and (Q I,k ,P I,k ) obey canonical commutation relations, but otherwise these operators commute with each other, as they correspond to different real modes. 60 Further, with these new variables, the Hamiltonian (263) becomes which is a sum of N = 2n + 1 independent real harmonic oscillators. A similar line of arguments can be used in order to express the ultralocal Hamiltonian (182) in terms of real harmonic oscillators. It is really when the regulated scalar field theory is expressed in terms of these (Hermitian) operators that we can easily evaluate the complexity by drawing upon the results in section 4 for the TFD state comprised of two (real) harmonic oscillators. Recall that in general, the complexity of a single mode depends on the frequency of the oscillator. Now with the assumption (which we proved for the L 2 norm with λ R = 1) that the total complexity is given by the sum of the complexities of the individual modes, we have where in the second equality we have used the fact that ω k = ω N −k . Using these arguments, eq. (192) follows immediately by substituting the different frequencies in eqs. (150) and (152), and a similar expression is obtained for the C (LR) 1 complexity using eq. (151). The relation between the real and complex variables can also be used to express the covariance matrix directly in terms of the complex variablesG a,b k (t) = TFD(t)|(ξ a kξ †b k + ξ †b kξ a k )|TFD(t) , whereξ k = (q L k ,q R k ,p L k ,p R k ), see eq. (236). To obtain the explicit expression for this expectation value, we should consider the covariance matrix (104) for the real modes 59 With the exception of k = 0, and for even N , also k = N/2. 60 These modes could be constructed by replacing the complex exponentials in the Fourier transform (178) with cos(2πka/N ) and sin(2πka/N ). together in terms of an 8 × 8 matrix ordered according to V := (Q L R,k ,Q L I,k ,P L R,k ,P L I,k , . . . R ) T , where the ellipsis indicates the same ordering of the right (R) degrees of freedom. We then rotate this matrix to the complex coordinatesṼ := (Q L k ,Q L −k ,P L k ,P L −k , . . . R ) T using eq. (264), namelyṼ = AV where The new covariance matrix will be related to the real mode expression according toG = AGA T which results in a block diagonal form for the k and −k modes where for instance the 4 × 4 block for the k modes reads where as defined in eq. (187), we have substituted λ = ω k /µ g into eq. (104). This also agrees with eq. (236) where implicitly, we made the choice µ g = L, but this choice is immaterial to the evaluation of the entanglement entropy in section 6.
Next, we would like to consider expressing the generators of our quantum circuits in terms of annihilation and creation operators. Following [23], in eq. (91), we expressed these generators in terms of a basis of Hermitian operators (Q i , P i ) with the standard canonical commutation relations, Here, the indices i, j = 1, . . . , N may specify either the positions along the one-dimensional lattice, or the real and imaginary parts of the Fourier modes, described above. This choice will not be important for the following discussion. Now in working with QFTs, it is also natural to express the quantum circuits in terms of annihilation and creation operators, e.g., see ref. [51][52][53]. Therefore let us introduce the operators satisfying [a i , a † j ] = δ i,j as usual. Here, to produce dimensionless operators, we need to introduce scales µ i but we have left these scales unspecified for the moment. If operators were chosen in the momentum basis, 61 then it would be natural to choose µ i = ω i /δ, as was done in refs. [51][52][53], in which case the a i annihilate the usual vacuum state of the (regulated) scalar field theory. Alternatively, in the present context of our study of complexity, it would also be natural to fix µ i = µ/δ in which case the a i instead annihilate the unentangled reference state, i.e., the ground state of the ultralocal Hamiltonian introduced in eq. (181). Of course, changing these scales can be accomplished with a subgroup of Bogoliubov transformations, R N ∈ Sp(2N, R), e.g., see ref. [56].
Turning now to the generators of our unitary circuits, we can write all of the quadratic combinations of these annihilation and creation operators as follows While we have not constructed Hermitian generators above, this is easily done by considering the combinations T + T , i( T − T ), A + A † and i( A − A † ). We note that the generators A span the u(N ) subalgebra of the full sp(2N, R) algebra. These are the compact generators while the remaining generators encoded in T and T are all noncompact [77].
In comparing the expressions above with the sp(2N, R) generators in eq. (91), we conclude that it is rather unnatural to choose different scales µ i for each of the (a i , a † i ) pairs. Following the more logical approach of setting all of these scales to be the same, we denote this common scale as the gate scale, i.e., µ i = ω g , and then we find where W , V and Z are the generators in eq. (91). This gives us a new insight into the role of the gate scale which was simply introduced to ensure that the sp(2N, R) generators were dimensionless. Here, we see that these generators are then naturally written in terms of annihilation and creation operators which annihilate a state where the degrees of freedom are all unentangled and the variance of all these Gaussians is set by ω g . Choosing this state to be the reference state may again seem like the natural choice, i.e., it seems natural to choose µ i = ω g = µ/δ. These observations may give us some additional confidence in focusing on the case of λ R = 1 in our analysis of the complexity in the main text.

E Complexity of formation in the diagonal basis
We have described in the main text the general behaviour of the complexity of formation in the diagonal basis given by eq. (209) if the reference scale is in the UV, with µ > Λ. We will now analyze two more examples: when the reference scale is exactly equal to the cutoff µ = Λ and when the reference scale is in the IR, with µ < Λ. For small temperatures, the profile of the ratio of complexity of formation over temperature becomes small. For higher temperatures, it develops a dependence on the temperature and the cutoff scale Λ. The fact that the complexity of formation is either zero or not proportional to the entropy contrasts with the holographic results of ref. [27].
Let us start by the choice µ = Λ. For instance, in cMERA circuits this is suggested as the natural reference scale for an ultralocal scalar theory. We will find that there are two important regimes to consider in this case. We focus on the only argument of the absolute value in eq. (209) that may change sign which is 1 2 log k µ + α k . The sign flips when there is some k f in the range [0, Λ] that satisfies the transcendental equation There are then two interesting ranges of temperatures namely, 0 < T < Λ/4 and T > Λ/4. Focusing first on the case 0 < T < Λ/4, we find that there is a k f in the range of integrated momenta that satisfies eq. (273). In this case for k < k f , we have that 1 2 log k µ + α k is always negative, and the integrand in eq. (209) vanishes, while for k > k f , we have that 1 2 log k µ + α k is always positive. The complexity of formation is then given by For T > Λ/4, we have that 1 2 log k µ +α k is always positive and does not change sign for k in the range [0, Λ]. All the other arguments of the absolute values are negative, and the complexity of formation then reads We show in fig. 26 the integrated complexity of formation in the diagonal basis divided by the entropy for several dimensions and µ = Λ. For small temperatures with respect to the cutoff scale Λ, the ratio of the complexity of formation and entropy goes to zero, and for large temperatures it develops a nontrivial profile. The sign flips in the arguments of the absolute values are slightly more complicated if the reference scale is in the IR, or µ < Λ. There are many possible sign flips in the arguments of the absolute values in eq. (209). We denote by k sum values of the momentum related to the sign flip of the absolute value of 1 2 log k µ + α k and by k sub values of momentum related to the sign flip of absolute value of 1 2 log k µ − α k . There are sign flips if k sum and k sub satisfy the transcendental equations Therefore, there are two critical temperatures associated with sign flips given in eq. (276), For the term in eq. (209) with 1 2 log k µ + α k , if the temperature satisfies T < T IR c1 , then one has to solve the transcendental equation (276) for k sum to find the momentum at which there is a sign flip. In this regime, the sum in question is negative for k < k sum and positive for k > k sum . If T > T IR c1 , then the sum 1 2 log k µ + α k is always positive. For the argument of the absolute value in eq. (276) given by 1 2 log k µ − α k , there are the following cases to analyze. If T < T IR c2 , one has to solve eq. (276) for k sub , and the subtraction is negative for k < k sub , and positive for k > k sub . If T > T IR c2 , the subtraction is always negative in the range of momentum [0, Λ]. In order to study the overall behaviour of eq. (276) when µ < Λ, we should break down the different cases depending on whether T IR c1 is smaller or bigger than T IR c2 . Without loss of generality, we will consider T IR c1 < T IR c2 , which from eq. (277) is satisfied if µ < 0.833557Λ. There are three separate regimes to consider in this case. If T < T IR c1 , one has to solve for k sum and k sub in eq. (276) in order to find the momentum where the sign flips for all the individual terms. The complexity of formation in this regime reads Next, in the range T c1 < T < T c2 , the sum 1 2 log k µ + α k is always positive, and the subtraction 1 2 log k µ − α k changes sign at k sub given by eq. (276). The expression for the complexity of formation reads Finally, if T > T c2 , then the sum 1 2 log k µ +α k is always positive and the subtraction 1 2 log k µ −α k We show the profile of the complexity of formation for µ = Λ/2 in fig. 27. Analogously to the previous cases, the ratio of the complexity of formation and the entropy has a nontrivial profile that in principle depends on T and the cutoff scale Λ. In fig. 28, we compare the ratio of the complexity of formation over the entropy for different reference state scales for d = 4. We notice that for smaller reference state scale, the peak of the profile is larger and happens at smaller temperatures. F Minimal geodesics for N degrees of freedom with λ R = 1 In section 5, we have proposed that the shortest geodesic (with respect to the unpenalized metric) connecting a Gaussian reference state |G R with a Gaussian target state |G T is just given by a combination of squeezing operations on independent normal modes. In this appendix, we prove this result for the full group Sp(2N, R) but only for the F 2 or D κ=2 norms with the choice λ R = 1. While the latter choice does not modify the underlying geometry on the space of unitaries, it does modify the boundary conditions describing a particular target state. Recall that the latter is not achieved by a single unitary transformation U (acting on the reference state) but rather by a family of unitaries U U R , where U R belongs to the stabilizer group (83) of the reference state. As we will see below, with λ R = 1, the stabilizer group takes a particularly simple form and allows us to make use of a particular fiber bundle structure. Note that a very similar discussion of the fermionic case can be found in ref. [56]. Our proof requires some Lie group techniques together with some well-known decompositions of symplectic group elements. We also rely on the conventions established in section 3. We assemble the N degrees of freedom together as ξ := (q 1 , · · · , q N , p 1 , · · · , p N ). Of course, with this choice of the standard basis of conjugate positions and momenta, the canonical commutation relations are given by [ξ a , ξ b ] = iΩ a,b where Ω a,b is the canonical symplectic form in eq. (42). The Gaussian states |G are then completely characterized by the symmetric part of their covariance matrix G a,b in eq. (43). 62 Finally recall that as described in section 3.2, we are naturally led to a representation of the symplectic group acting on the covariance matrix where the generators K ∈ sp(2N, R) are given by where k a,b define the quadratic operators appearing in the gates, as in eq. (48).

Lie group geometry
In the Nielsen approach to complexity, we equip the Lie group Sp(2N, R) with a right invariant and positive metric. Such a metric is completely characterized by its value at the identity where we identify the tangent space T 1 Sp(2N, R) with its Lie algebra sp(2N, R). We represent generators A ∈ sp(2N, R) as matrices, namely linear maps A a b on phase space. In order to define our right invariant metric, we need to specify its value at the identity by giving a prescription how to compute the inner product A, B 1 between the two generators A, B ∈ sp (2N, R). For a given Gaussian reference state |G R , there is a natural metric that we call the unpenalized metric because it is the generalization of the unpenalized metric (109) on Sp(2, R) to systems with more degrees of freedom. This matrix is defined as Given two tangent vectors X, Y ∈ T U Sp(2N, R) represented as matrices at a point U ∈ Sp(2N, R), we can find their inner product from the right-invariance by multiplying with U −1 from the right, namely where we can apply eq. (282).

Fiber bundle structure
The choice of the reference state |G R equips the Lie group Sp(2N, R) with a fiber bundle structure. There exist symplectic group elements U that leave the reference state invariant, namely the stabilizer subgroup Sta G R ⊂ Sp(2N, R) in eq. (83). We explicitly introduce the choice λ R = 1 because this simplifies the covariance matrix for the reference state to being a 2N × 2N identity matrix, i.e., G R = 1 2N for λ R = 1 (compare to eq. (47) for N = 1). Then the elements of the stabilizer group are both symplectic (with respect to Ω) and orthogonal (i.e., the reference state is left invariant: G R = U G R U with G R = 1 2N ), and so they form the subgroup It is well known that U(N ) is the largest subgroup of Sp(2N, R) and that the different choices of how to embed U(N ) into Sp(2N, R) are in one-to-one correspondence to the different choices of metrics G R . This is not surprising because every choice of a metric G R chooses a different subset Sta G R that are all isomorphic to U(N ). We define the equivalence relation U ∼Ũ if and only if U G R U =Ũ G RŨ . This means acting with U andŨ on G R will give the same target state. In particular, the stabilizer subgroup, i.e., U(N ), is equal to the equivalence class [1] of the identity. Moreover, for every pair U ∼Ũ , there exists a u ∈ U(N ), such that U u =Ũ . 63 Therefore, Sp(2N, R) becomes a fiber bundle where the fibers correspond to the different equivalence classes diffeomorphic to U(N ) and the base manifold is given by the quotient This space is diffeomorphic to R N (N +1) as a manifold and generally referred to as a symmetric space of type CI. We will refer to U as the space of pure Gaussian states and identify a point [U ] ∈ U with the Gaussian state |U G R U up to an overall complex phase.

Polar decomposition
Identifying the Lie algebra sp(2N, R) with the tangent space at the identity, we have a natural "vertical" subalgebra u(N ) ⊂ sp(2N, R) that is tangential to the fiber [1] = U(N ). A priori, there is no natural "horizontal" complement to write the Lie algebra as a direct sum of a vertical and a horizontal part. However, by equipping the Lie algebra with the inner product ·, · 1 in eq. (282), we can choose the orthogonal complement sym(N ) := A ∈ sp(2N, R) A, B 1 = 0 ∀ B ∈ u(N ) .
In contrast to u(N ), sym(N ) is not a subalgebra. Its name stems from the fact that the decomposition sp(2N, R) = sym(N ) ⊕ u(N ) is equivalent to splitting the set of generators into symmetric and antisymmetric matrices with respect to the metric G R of the reference state: • Vertical subspace u(N ) A generator B in the subspace u(N ) must preserve the reference state G R . It satisfies which is equivalent to B being an antisymmetric matrix in a basis where G R is the identity.
• Horizontal subspace sym(N ) = u ⊥ (N ) A generator A that is orthogonal to all elements B ∈ u(N ) satisfies 0 = A, B 1 = Tr(AG R B g R ) .
In a basis where G R and g R are just the identity, we are searching for a matrix A that has zero trace when multiplied with any antisymmetric matrix B. This condition is equivalent to stating that A is a symmetric matrix, namely satisfying We can refer to sym(N ) as orthogonal complement of u(N ), i.e., u ⊥ (N ).
We can exponentiate the space sym(N ) to define the N (N + 1)-dimensional submanifold Sym + (N ) = exp (sym(N )) = e A A ∈ sym(N ) consisting of all symplectic group elements that are symmetric with respect to G R and have positive eigenvalues. The fact that symmetric generators A are diagonalizable with real eigenvalues implies that the exponential map provides a diffeomorphism and Sym + (N ) is thus diffeomorphic to R N (N +1) . Let us emphasize that Sym + (N ) is not flat with respect to our right-invariant metric, but rather some complicated embedded curved surface. The polar decomposition of a symplectic group element U is given by U = T u with T = U G R U g R ∈ Sym + (N ) and u = T −1 U ∈ U(N ) .
It is unique and provides a diffeomorphism between the symplectic group and the Cartesian product Sym + (N ) × U(N ). In particular, it provides a trivialization of the fiber bundle on Sp (2N, R) where the base manifold is identified with surface Sym + (N ), from which we can move up and down along the fiber by multiplying with group elements u ∈ U(N ). Due to the fact that Sym + (N ) is diffeomorphic to sym(N ), we can use the pair (A, u) as generalized polar coordinates for any group element U (A, u) U (A, u) = e A u with A ∈ sym(N ) and u ∈ U(N ) .

Cylindrical foliation
Using the polar coordinates (293), we can foliate the symplectic group by generalized cylinders defined as C s = e A u A ∈ sym(N ), A = s, u ∈ U(N ) The Lie group can be represented as fiber bundle over its quotient given by the symmetric space Sp(2N, R)/U(N ). This base manifold can be interpreted as the space of Gaussian quantum states. The fiber over the reference state |G R is given by the subgroup U(N ) ⊂ Sp(2N, R), while the fiber e A U(N ) over any target state |G T is not a subgroup. We consider a path γ in the group that connects 1 to some other group element U = e A u. Such a point lies on the cylinder C s with s = A . Every curve γ in the group can be projected down to the curve π • γ in the base manifold. The vertical submanifold Sym + (N ) = exp sym(N ) is generated by exponentiating sym(N ) and it plays an important role because it contains the minimal geodesics. In particular, the straight line e tA connecting 1 with e A will turn out to be the minimal geodesic between 1 and the fiber e A U(N ). We do not show the vector field R consisting of radially outwards pointing unit vectors on the cylindrical surfaces C s , such that the curves e tA u are its integral curves. Note that a similar sketch was used in ref. [56] to describe the setting of fermionic Gaussian states.
with the topology S N (N +1)−1 × U(N ). 64 Moreover, we will define the radial vector field R at point U (A, u) ∈ Sp(2N, R) given by We will prove that this vector fields points radially outwards and is everywhere orthogonal to the cylindrical surfaces C s . We will show the orthogonality by considering different directions individually. Note that the normalization 1/ A is irrelevant here.
• Orthogonality to the U(N ) fiber: We show that R is orthogonal to any vector pointing along the U(N ) fiber. Let X ∈ u(N ), so that e A uX points in the direction of the U(N ) fiber at point U (A, u). We can compute the inner product We define Y = uXu −1 which lies in u(N ) because u(N ) is a subgroup. This implies uX = Y u. We can therefore compute e A uX, e A Au e A u = e A Y u, e A Au e A u = e A Y, e A A e A = e A T e −A , A 1 .
At this point, we can use the explicit form of the metric at the identity given by where we have used G R A g R = A for A ∈ sym(N ).
• Orthogonality to a generator A ∈ sym(N ) preserving C s : This second computation is slightly more involved. Let us look at a point U = e A u and ask what are the directions B ∈ T U Sp(2N, R) that are tangential to the surface C s , but also to the surface exp sym(N ) u. We can describe such elements by choosing a second generator B ∈ sym(N ) that is orthogonal to A with A = B . The circle γ(t) = e (cos (t) A+sin(t) B) u (299) lies in Sym + (N ) and on C s with s = A = B . This gives rise to the tangent vectoṙ We can compute the inner product with R U (A,u) using A = G R A g R as At this point, we write out the full exponential as where we have used the fact that trace is cyclic and that B was chosen orthogonal to A. Note that the sum just gives the identity.
This proves that we have indeed a vector field R that is everywhere orthogonal to the cylindrical surfaces C s . Furthermore, we can quickly confirm that this vector field indeed has a constant length equal to 1, by computing Given a trajectory γ : [0, 1] → Sp(2, R) : t → γ(t), we can compute how the coordinate s(γ(t)) changes. Due to the fact that the vector field R is orthogonal to the surface C s of constant s and correctly normalized, we have ds = R γ(t) ,γ(t) γ(t) . (305)

Inequality for the geodesic length
We will now use the cylindrical structure to bound the geodesic length from below. Given an arbitrary point U (A, u) = e A u on the cylinder C s , let us assume that we have already found the shortest path connecting the identity 1 with U (A, u). This path may be given by γ(t) with γ(0) = 1 and γ(1) = U (A, u). We can compute the change ds as the inner product Clearly, if we integrate this inner product we find how far we move in the s-direction. This follows directly from the fact that moving in the direction of R increases s with a constant rate, while moving along any orthogonal direction does not change s. Therefore, we have s = We can compare this with the actual length of the geodesic given by At this point, we should note that γ(t), R γ(t) γ(t) ≤ γ(t) for all t. This follows from the fact that we are projecting onto the unit vector R, so this projection is at most the length oḟ γ(t). We can combine these two equation to find the important inequality stating compactly that any path connecting 1 with U ∈ C s must have a length of s or more.
Shortest path to a fiber e A U(N ) At this point, we have not proven that for every U ∈ Sp(2N, R) there exists a path with length s and there certainly are points U where we cannot find such a shortest path. However, we are interested in the minimal geodesic that connects the identity 1 with an arbitrary point in the fiber [U ]. This means that if we find a single path that does this with length s, we have proven that this is indeed the optimal path and there is no shorter one.
We can do this for the fiber of e A U(N ) by checking that the path γ(t) = e tA (310) satisfies exactly these conditions and reaches the representative e A at t = 1. This path has length γ = A = s. At this point, we have proven that for the "unpenalized" inner product discussed at the beginning, the shortest path is indeed always given by e tA with A ∈ sym(N ).
We can now ask how A is related to the target state G T . We must have Now requiring that A ∈ sym(N ) implies that U = e A is symmetric with respect to the basis where G R is the identity. In an invariant language, we have With this in hand, we can claim that the linear map U = √ G T g R will do the job. Importantly, U satisfies U G R = G R U . We can explicitly verify that The algebra element that generates U is given by A = log U = 1 2 log G T g R . We have s = A = 1 2 G T g R . Let us note at this point that all expressions, such as log G T g R and √ G T g R are well defined, because G T g R is a positive symmetric, symplectic matrix in a basis where g R is the identity. This fact implies that G T g R is (a) diagonalizable and (b) has positive non-zero eigenvalues.
Of course, the linear map G T g R is precisely the relative covariance matrix (86) between our target state and reference state, This matrix encodes the invariant information about the relation between the reference state |G R and the target state |G T . The eigenvalues of ∆ come in conjugate pairs (e i , 1/e i ). We can compute the geodesic distance, which is equal to the norm A , directly 65 from ∆:

Normal mode decomposition
The result for the minimal geodesic is equivalent to stating that for any two Gaussian states |G R and |G T , there exists a set of preferred normal modes, such that the optimal geodesic just corresponds to a linear combination of single-mode squeezing operations on these independent normal modes. 65 Another useful quantity that can be directly computed from ∆ is the inner product Let us consider a reference state |G R and a target state |G T . We can choose a symplectic basis, such that the covariance matrix G R is simply the 2N ×2N identity matrix: This is always possible because G R is a positive, symmetric bilinear form. The covariance matrix G T will be another general symmetric matrix. However, we can still change basis by acting with a matrix u in the stabilizer group of G R , which leaves eq. (316) invariant. As in eq. (284), u is then an orthogonal matrix which acts by a similarity transformation on G T and by choosing u appropriately, we can put the latter in a diagonal form. Due to the fact that the covariance matrix of a pure Gaussian state is itself a symplectic matrix, the diagonal form of G T will consist of conjugate pairs of eigenvalues, i.e., We can refer to our final basis as ξ := (q 1 , · · · , q N , p 1 , · · · , p N ). In this basis, the matrix representation of ∆ = G T g R will be same as the one of G T , because G R and g R are represented by the identity. However, the eigenvalues e i are only matrix invariants of ∆, but not of G R nor of G T . The basis chosen above provides a normal mode decomposition, where each conjugate pair (q i , p i ) corresponds to a normal mode with a single degree of freedom. Both |G R and |G T can be written as tensor products over normal mode states, |G R = |0 ⊗ · · · ⊗ |0 , |G T = |e 1 ⊗ · · · ⊗ |e N , where the state |0 is the ground state of H = 1 2 (p 2 i + q 2 i ), while the |e i 's are the ground states of H = 1 2 (p 2 i + e 2 i q 2 i ). The generator A is diagonal in the same basis and given by The squeezing operator producing |G T = e −iÂ |G R is given byÂ = 1 2 ξ a ω a,c A c d ξ d (where ω a,c Ω c,b = δ b a ). This can be explicitly written aŝ log e i 4 (q ipi +p iqi ) where the last expression uses the notation in eq. (91). Hence the generator producing the optimal circuit simply consists of N commuting single-mode squeezing operatorsÂ i which squeeze each of the normal modes independently. Furthermore, it is now straightforward to evaluate the complexity using the above results. However, let us first note that a similarity transformation was needed to put G T in the diagonal form given in eq. (317). Hence we should focus on the F 2 or κ = 2 cost functions (given in eqs. (8) and (9), respectively) since they are invariant under rotations of the basis. From eq. (320) or by comparing the matrix generator A in eq. (319) with the sp(2N, R) generators in eq. (95), we see that the tangent vector to this circuit is simply given by and hence the F 2 complexity is given by where the extra factor of 1/ √ 2 appears in the second expression because the eigenvalues of the relative covariance matrix appear in conjugate pairs, e.g., see eq. (317). This last expression gives a simple covariant formula for the F 2 complexity in terms of the relative covariance matrix, which can be applied for any reference and target states without any further calculation. In addition, since the tangent vector (321) is constant, it follows that the κ = 2 complexity is simply related to the above expression with Application to Sp(2, R) In section 4.2, we examined the special case of Sp(2, R) using the coordinates (ρ, θ, τ ) given in eq. (111). In eq. (119), it was found that the final state only depends on ρ and the combination χ = θ + τ . Further, it was found that the optimal circuit preparing a particular target state in the relevant equivalence class was given by the simple straight-line geodesic in eq. (121). That is, ρ(σ) = ρ 1 σ , θ(σ) = θ 1 , τ (σ) = 0 , where ρ 1 and θ 1 characterize the target state G T . It is interesting to understand this result from the perspective presented in this appendix and so we consider here applying the previous analysis to the special case of N = 1. First, we can consider Sp(2, R) as a U(1) fiber bundle over the plane parameterized by (ρ, θ) = (ρ, χ) with fixed τ . The subgroup U(1) that preserves the reference state G R = 1 is generated by which is an antisymmetric matrix in accord with eq. (288). The U (1) fiber above the identity is then given by simple rotation matrices e τ K 3 = cos τ − sin τ sin τ cos τ , (1) Sym (1) e sK 1 1 Figure 30: This figure illustrates the geometry of Sp(2, R) in the coordinates (ρ, θ, τ ) in the left picture and (ρ, χ, τ ) in the right picture with χ = θ + τ . The identity element 1 is in the center on the surface Sym(1) corresponding to τ = 0. The subgroup U (1) corresponds to the vertical line with ρ = 0. The equivalence class e K U (1) of group elements that prepare the same target state are given by spiral (left picture) and vertical line (right picture). For a reference state with |G(ρ, χ) , all group elements preparing this state necessarily lie on the cylindrical surface C ρ . A general geodesic γ winds around, but the minimal geodesic corresponds to a straight line of the form e sK .
which we recognize to agree with U (0, 0, τ ) in eq. (111). This is simply the stabilizer group of the reference state in the case λ R = 1, see eq. (85). This means that K 3 spans the subalgebra u(N ) for N = 1 whose orthogonal complement u ⊥ (N ) with respect to the right-invariant inner product (283) is spanned by This subspace consists of symmetric matrices with respect to G R , in accord with eq. (290). These two generators are referred to as sym(N ) = u ⊥ (N ) with N = 1. Applying the exponential map to an arbitrary generator ρ(cos θ K 1 + sin θ K 2 ) ∈ sym(1) gives e ρ(cos θ K 1 +sin θ K 2 ) = cosh ρ − sin θ sinh ρ cos θ sinh ρ cos θ sinh ρ cosh ρ + sin θ sinh ρ , = U (ρ, θ, 0) , which we recognize in the last equality to be given by elements in the plane with τ = 0, again symmetric matrices. We therefore refer to this subset as Sym(N ) = exp(sym(N )) with N = 1. Let us emphasize that this is not a subgroup, because the product of two symmetric matrices will in general not be symmetric. In this notation, the straight-line geodesic in eq. (324) preparing the state G T (ρ = ρ 1 , χ = θ 1 ) is described by γ(σ) = U (ρ = σρ 1 , θ = θ 1 , τ = 0) = e σρ 1 (cos θ 1 K 1 +sin θ 1 K 2 ) . where following the framework of ref. [95]. This is obtained by acknowledging that the TFD in eq. (18) has the form of a pure vacuum state, time evolved under a quadratic Hamiltonian with the real number α encoding the inverse temperature formally taking the role of a time. This is a convenient form of the covariance matrix of the TFD of a single decoupled mode. One can easily verify that W (α) Ω W (α) = Ω, where Ω is the symplectic form in this convention, so that indeed W (α) ∈ Sp(4, R) for all values of α. Again following ref. [95], the time evolved TFD state (31) can be expressed as a matrix exponential G(α, t) = e tK G(α)e tK , with K having components K a b = Ω a,c (h ⊕ h) c,b . This is a concise form of the covariance matrix of the time-dependent TFD state of a single mode. Having paved the ground in the case of a single mode and its double, we now turn to the TFD state of the full quantum field theory with the Hamiltonian (175). This Hamiltonian can be written with respect to the following choice of coordinates ξ = (q L,1 , q L,2 , . . . , q L,N , p L,1 , p L,2 , . . . , p L,N , q R,1 , q R,2 , . . . , q R,N , p R,1 , p R,2 , . . . , p R,N ), (338) and x = δ −1 m 2 1 N + δ −3 Toeplitz(2, −1, . . . , −1).
The latter is a banded matrix with the entry 2 on the main diagonal and −1 on the first off diagonal. This is a concise way of expressing the Hamiltonian of the UV regularized quantum field theory. The matrix x can be diagonalized as with a suitable O ∈ O(N ), so that k can be diagonalized according to (4N, R). This is easy to see: The coupling is in the positions only, so the momenta will be left unchanged and captured by the diagonal matrix 1 N δ by this real orthogonal transformation, while the positions are diagonalized to D. This way presents an alternative to the direct complex Fourier transform. We now turn to a closed form of the covariance matrix of the entire time dependent thermofield double state of the quantum field theory. We will keep the expressions entirely in the form of matrix exponentials. Using eqs. (333) and (343), the covariance matrix of the initial vacuum state is found to be G = diag δ 1/2 x −1/2 , δ −1/2 x 1/2 , δ 1/2 x −1/2 , δ −1/2 x 1/2 , in terms of the matrix square root of x and its inverse [87]. The covariance matrix of the full N -mode thermofield double is G(α) = W (α)GW (α) , where now This is again a convenient form that reflects the fact that the transformation that diagonalizes the position part of the Hamiltonian commutes with the transformation that maps the vacuum onto its TFD state for single modes. Following eq. (30) we obtain the final expression for the covariance matrix of the time-dependent TFD state which reads G(α, t) = exp(tk/2)G(α) exp(tk /2), where the factor of 1/2 is a result of the convention taken in eq. (30). Here, the expression is concise, as the Fourier transform does not have to be explicitly performed, the covariance matrix being expressed directly as a matrix function of the coupling matrix.