Binding Complexity and Multiparty Entanglement

We introduce"binding complexity", a new notion of circuit complexity which quantifies the difficulty of distributing entanglement among multiple parties, each consisting of many local degrees of freedom. We define binding complexity of a given state as the minimal number of quantum gates that must act between parties to prepare it. To illustrate the new notion we compute it in a toy model for a scalar field theory, using certain multiparty entangled states which are analogous to configurations that are known in AdS/CFT to correspond to multiboundary wormholes. Pursuing this analogy, we show that our states can be prepared by the Euclidean path integral in $(0+1)$-dimensional quantum mechanics on graphs with wormhole-like structure. We compute the binding complexity of our states by adapting the Euler-Arnold approach to Nielsen's geometrization of gate counting, and find a scaling with entropy that resembles a result for the interior volume of holographic multiboundary wormholes. We also compute the binding complexity of general coherent states in perturbation theory, and show that for"double-trace deformations"of the Hamiltonian the effects resemble expansion of a wormhole interior in holographic theories.


Introduction
The importance of quantum computational complexity in computer science became apparent after Shor [1] proved that the quantum circuit model could solve integer factorization in polynomial time. The typical notion of quantum computational complexity counts the minimal number of simple unitary operations needed to reach some target state from a specific initial state. For instance, one may be interested in how hard it is to prepare the (generically entangled) ground state of a given Hamiltonian starting from an initial state which is factorized across all degrees of freedom. This prompts the related question of whether there is a relation between the strength and structure of the entanglement between degrees of freedom in a quantum state and the complexity of preparing that state. In this work, we answer this question in the affirmative for a type of complexity we call binding complexity that counts the number of quantum gates acting on multiple parties simultaneously.
A motivating example that the binding complexity might be connected to the strength of entanglement comes from examination of the two inequivalent classes of multiparty entanglement between three qubits [2], the GHZ and W states, and their n-party generalizations: √ n (|00 . . . 01 + |00 . . . 10 + . . . + |10 . . . 00 ) . (1. 2) The GHZ states are separable upon tracing out any subset of the parties, whereas the W states are not. In this sense, the W states can be thought of as possessing more robust entanglement. We can understand the structure of these states better by computing the entanglement entropy of one qubit with the rest, as the number of qubits n grows large. We would normally understand this quantity as a diagnostic of the strength of entanglement between parties.
In more detail, the entanglement entropy corresponding to a partition (A,Ā) of degrees of freedom in a quantum state is defined as the von Neumann entropy of the reduced density matrix on A: S A = −Tr(ρ A ln ρ A ). For the GHZ states, we find that entanglement entropy of a single qubit with the rest of the system is S 1,GHZ = ln 2, (1.3) which is constant, nonzero, and independent of n. By contrast, for the W states the single qubit entropy is S 1,W = − n − 1 n ln n − 1 n + 1 n ln 1 n → 0 as n → ∞. (1.4) It is tempting to conclude from (1.3) and (1.4) that the GHZ states possess "stronger" entanglement, at least for large n, since there is always maximal entanglement between even a single party and the rest. This seems to be in qualitative tension with our conclusion above that the W states have a more robust pattern of entanglement.
However, one should reinterpret these equations using the principle of monogamy of entanglement [3]: although S 1,GHZ is constant, tracing out one party removes all of the entanglement as the remaining state is separable. Conversely, S 1,W is small because as the number of parties grows large, tracing out one party only removes a very small amount of entanglement: nearly all of the entanglement remains tied up between the remaining n − 1 qubits which are still approximately in a W state. We would like to define a quantity that captures this sort of "robustness" of entanglement: it is distributed between several parties and is difficult to destroy.
Correspondingly, let us consider quantum circuits preparing the GHZ and W states, and the binding complexities associated to them. To compute complexity we must fix a set of allowed gates that we may use to prepare states. Here, we take the gate set to be the set of all one-qubit or two-qubit unitary operators, although typically we will want further restrictions on which unitaries are allowed.
In the case of the GHZ states, it is very easy to explicitly write down a circuit that prepares a generalized GHZ state from the factorized state |0 ⊗n (Fig. 1). This circuit uses n gates, n−1 of which act on multiple parties. The binding complexity, i.e. the number of gates acting on multiple parties at once, is simply n − 1. It is easy to understand that one cannot write a more efficient circuit to construct a GHZ state because a minimum of n − 1 two-party gates are required simply to couple all of the qubits; otherwise, the state will factorize across some partition. In the case of the W states, it is not simple to write down a circuit, and there is no proof of minimality. However, [4] gives a deterministic construction of arbitrary W states that requires 1 2 n(n + 1) − 2 ∼ O(n 2 ) two-qubit gates. To our knowledge no asymptotically more efficient construction has been found. In fact, we should expect that none exists -intuitively, since the W state is not separable upon tracing out any number of parties, it is as if n 2 ∼ O(n 2 ) gates have been used to entangle all pairs of qubits. Consequently, at least in the qubit context, we see that the binding complexity is a natural diagnostic of the robustness of entanglement -the minimal number of gates required to entangle different parties naturally controls how entangled the parties become in the final state. Indeed, we will demonstrate bounds relating binding complexity to other measures of robustness such as entanglement negativity, which quantifies non-separability of quantum states.
We will study binding complexity in a toy model of a free scalar field [5], which reduces to a system of harmonic oscillators. Binding complexity is defined as the minimum number of gates acting on multiple parties that is needed to prepare the state starting from a specified reference. In Nielsen's geometric approach [6][7][8] to complexity, one places a Riemannian metric on the space of unitaries, so that complexity is measured by the geodesic distance between the identity and the unitary operator that makes the state of interest. We choose a metric that is infinitesimal in directions that act only on a single party, so that the geodesic length measures the binding complexity that we want to study. 1 A key step in the computation of circuit complexity is the choice of the gate set. The vacuum wavefunction for a coupled oscillator system is a Gaussian of the schematic form e − x T Ω x . In [5], the gate set acting on such states was chosen to change the components of Ω. We will divide the oscillators into "parties" defined by block structure in Ω. We want to compute the binding complexity of states that are entangled between these parties. Essentially, this involves only counting the gates from [5] that act across parties -we will call these the relevant gates. To calculate binding complexity we employ the Euler-Arnold approach to simplify the geodesic equation using the Lie algebra of the gate set. 2 It has been suggested that entanglement in quantum field theory can be holographically realized by wormholes between otherwise disconnected regions of spacetime [32,33]. In these contexts complexity in field theory has been conjectured to be dual to the volume or action of an interior region of the wormhole [34][35][36][37]. It is also possible in 2+1 dimensions to construct wormholes that connect multiple asymptotic regions [38][39][40][41][42][43]. Recently these geometries were 1 The Nielsen approach was previously extended to free fermion fields in [9,10], coherent states of free scalar fields in [11], states in φ 4 theory in [12], applied to the study of complexity growth following a quench in [13,14], and used to study the complexity of Hamiltonians and quantum phase transitions in [15]. An axiomatic study of the Finsler geometry in the Nielsen approach and comparisons to the holographic expectation in thermofield double states and their time evolutions was undertaken in [16][17][18][19]. Other approaches to field theory definitions of complexity include complexity from the Fubini-Study metric using momentum space [20] and complexity from optimization of the Euclidean path integral [21][22][23][24][25][26]. Current work has taken first steps towards understanding the Nielsen complexity in CFT and connecting it to the path-integral complexity [27,28]. Most recently, it was argued that the Nielsen complexity is superior to several of the other methods as only the Nielsen complexity displays the correct behavior under certain forward-and backward-time evolutions [29]. 2 Nielsen suggested applying the Euler-Arnold equation in [8], which was originally explained in [30]. A nice review can be found in [31].
used to study multipartite holographic entanglement [33,[44][45][46]. Since binding complexity measures the difficulty of entangling the wavefunctions of multiple otherwise disconnected parties, we conjecture that it is related to the interior volume of multiboundary wormholes, i.e., Binding Complexity = Volume of Wormhole Interior.
We address this conjecture by computing the binding complexity for a natural class of multiparty entangled states in our toy model, and showing that it has a linear dependence on entanglement entropy like the interior volume of the multiboundary wormholes of [33,41,42]. The CFT states dual to these wormholes were prepared by the Euclidean path integral on a branched bulk topology [33]. Consequently, we consider states in our toy model which are prepared by the Euclidean path integral on certain branched graphs with wormhole-like structure. 3 We find that such states have binding complexity and entanglement structure that are (a) similar to properties of the wormhole interior volume, and (b) reminiscent of the bit thread perspective on holographic states. 4 As a further check on the conjecture we test that adding small double-trace interactions between distinct parties causes binding complexity to increase linearly in the expansion parameter, as expected from the volume increase of holographic wormholes in the Gao-Jafferis-Wall approach [47].

Lower bounds
To begin, we will demonstrate some elementary lower bounds on binding complexity in terms of other measures of the entanglement structure of a state, such as the entanglement entropy, separability, etc. For simplicity, we will focus on a gate set G consisting only of one and two-party gates, although our arguments can be generalized to k-local gates.
We begin with the simplest example. Imagine that our Hilbert space can be decomposed into two tensor factors: where A consists of N A parties and B consists of the remaining N B parties. Let ψ be a state in this Hilbert space, and consider a unitary quantum circuit which builds ψ from the reference state |00 . . . 0 where the U i are one and two party gates which are allowed within our gate set. Of these gates U i , those which act within A or B do not contribute to the entanglement between A and B; only the two-party gates which act across this partition will contribute to the entanglement. Let n AB be the number of such gates which act across the partition. As discussed in the introduction, the binding complexity of the state ψ with respect to the partition H A ⊗ H B is equal to the minimum value n AB in the set M ψ,G of all the quantum circuits which construct ψ using the gate set G: Figure 2: We can "cut" a two-party gate (denoted by the red dashed line) by using its operator Schmidt decomposition into a sum of products of one-party operators.
In order to study the entanglement structure of ψ given such a quantum circuit in M ψ,G , we introduce the concept of "cutting a gate" (see Fig. 2). Any two-party gate can always be written in the form where s a are positive real numbers, the {O (a) 1/2 } are a basis of (not necessarily unitary) operators on the first/second party, and J is called the operator Schmidt rank of U . This is referred to as operator Schmidt decomposition [48]. Some examples of the operator Schmidt decomposition of two-qubit gates are: Turning to our original state ψ in (2.2), we repeatedly employ operator Schmidt decomposition to cut all the two-party gates which act across the partition H A ⊗ H B , while leaving all other gates untouched. This allows us to rewrite the state in the form (see Fig. 3 where a = (a 1 , · · · a n AB ). If we denote by J G the maximum operator Schmidt rank of any gate in the gate set G, then the above formula shows that the rank of the reduced density matrix a 1 a 2 a 3 Figure 3: A sample piece of a unitary quantum circuit. The red lines denote the subsystem A and the blue lines denote the subsystem B. We have "cut" all the two-party gates acting across the bipartition by using their operator Schmidt decomposition into sums of products of one-party operators.
on A (or B) will be upper bounded by J n AB G . Therefore, the entanglement entropy between A and B satisfies the upper bound 5 (2.8) While this upper bound is satisfied by every quantum circuit which constructs ψ from the given gate set G, the bound will be the tightest for the circuit which minimizes n AB . Therefore, we conclude that 9) or equivalently This bound shows that the binding complexity of the state with respect to a bipartition is lower bounded by the entanglement entropy. Intuitively this is clear, because if we are to build a state with a certain amount of entanglement, then we will need sufficiently many gates to achieve this.
We can easily generalize this bound to multipartite systems. Consider for example a tripartite system H A ⊗ H B ⊗ H C consisting of N A , N B , and N C qubits respectively. Then by cutting arguments similar to those used above, we obtain For the tripartite system, the binding complexity is defined as the minimum value of (n AB + n BC + n CA ) across all circuits in M ψ,G , and therefore we obtain Similarly, the n-partite generalization of this result is (2.14) So far, we have focused on bounds involving the entanglement entropy. However, as we discussed in the introduction, the entanglement entropy is not always sufficient to probe the fine-grained multiparty entanglement structure of the state. For this purpose, it is useful to consider other information theoretic concepts such as separability. Let us consider a tripartite quantum system H = H A ⊗ H B ⊗ H C . If we trace out A, then the reduced density matrix on BC is called separable if and only if it can be written in the form where ρ i B/C are density matrices on B/C, and p i are positive real numbers which sum up to 1. In this case, we interpret ρ BC as having no quantum entanglement, i.e., tracing out the subsystem A has destroyed the quantum entanglement between B and C. On the other hand, if ρ BC is not separable, the state retains quantum entanglement despite tracing out A. A necessary (but not sufficient) criterion for separability is the Peres-Horodecki positivity of partial transpose [49][50][51]. Here, we are instructed to construct the partial transpose ρ Γ BC of the density matrix, which is defined as: where |j B , j C and |j B ,j C denote basis vectors for H B ⊗ H C . If ρ Γ BC has any negative eigenvalues, then this necessarily implies that the density matrix ρ BC is not separable. We can therefore quantify the amount of non-separability of ρ BC by the number of negative eigenvalues of the partial transpose ρ Γ BC . We will denote as E A|BC the logarithm of one plus the number of negative eigenvalues of ρ Γ BC . Another measure of the non-separability is the entanglement negativity N A|BC , which is defined as where ||A|| = Tr √ A † A is the trace norm.
Going back to our state ψ ∈ H A ⊗ H B ⊗ H C , consider once again some unitary quantum circuit in M ψ,G which builds the state from the chosen gate set. By cutting all the two-party gates which act across the tripartition, we can now express the state in the form where as before a = (a 1 , · · · a n AB ), b = (b 1 , · · · b n BC ) and c = (c 1 , · · · c n CA ). It is clear from this expression that if we trace out A, then the number of negative eigenvalues of ρ Γ BC will be upper bounded by the maximum allowed rank of ρ BC minus one (there needs to be at least one positive eigenvalue so the trace can be one), i.e., (J n AB +2n BC +n CA G − 1). Therefore, Using the same argument by in turn tracing out B and C, we obtain Once again, the tightest bound is obtained for the circuit which minimizes the right hand side, from which we conclude At least in the case of qubit systems, the same bound is also true for the (logarithmic) entanglement negativity, i.e., if we replace E → ln(1+2N ) in all the terms above. This follows from the fact that the magnitude of all negative eigenvalues is always upper bounded by 1/2 [52]. However, the bound is tighter when stated in terms of E. The bound in equation (2.21) shows that the binding complexity is a much more fine grained probe of the entanglement structure than the entanglement entropy, and in particular is sensitive to multiparty entanglement measures such as separability.

Computation of the binding complexity
To begin, we review the definition of complexity as a geodesic length in the space of unitary operators and explain how the group structure of this space simplifies the geodesic equation.
We will keep all sums explicit below as some repeated indices will not be summed. Let us consider a general quantum system with Hilbert space H. We start by fixing some base state ψ 0 , such as a completely factorized state. Now consider some other pure state ψ of the entire system, which we wish to study -for instance, ψ could be the ground state of some interesting Hamiltonian. Let U be the space of all unitary maps on H, and let {O I } be a basis for its Lie algebra u: We may think of O I as generators of the elementary unitary gates at our disposal (thus e iO I are the elementary gates). Let U ∈ U be an operator such that In order to define the complexity of U , we need a notion of distance on the group manifold U. One possibility is the standard bi-invariant metric which is obtained from the Cartan-Killing form K IJ on the Lie algebra u, defined in terms of the structure constants as 6 7 If we allow gates that can act on any number of qubits at the same time, arbitrarily complex operations could be done in a single step. Thus, it is necessary to restrict the gate set to be "local" in some sense. We will require gates to be "bilocal", acting on no more than two qubits at the same time.
So far our discussion has been general, but now we wish to specialize the notion of complexity to study multiparty entanglement. To this end, let us consider a system which has a natural tensor factorization of the form In order to study the multiparty entanglement structure (with respect to the above partition) of a state in this Hilbert space, we define the binding complexity as the minimal number of gates, required in a quantum circuit construction of U , which act on more than one factor at a time, i.e., they act across the chosen partition. Gates which act within a tensor factor do not add to entanglement, and as such are treated as irrelevant. However, gates which act on two or more factors do contribute to the entanglement between various parties, and as such will be regarded as relevant. We wish to optimize over the number of relevant gates in building the unitary U . 6 The Cartan-Killing form satisfies N fIJ N KNM = − N fIM N KNJ , which is simply the statement that it is invariant under adjoint action of the group, i.e., Tr ad ([Z, X]Y ) + Tr ad (X[Z, Y ]) = 0 for any three elements X, Y, Z of the Lie-algebra. 7 We have defined the Cartan-Killing form up to overall sign and normalization here, since our main goal is using it to construct a right-invariant Riemannian metric whose normalization is fixed by "cost factors" (see (3.5)).
To accomplish this we can define a different inner product G IJ on the Lie algebra, which assigns a different "cost" for gates acting on one vs. multiple parties. We define the inner product by the metric where the c I are the cost factors for the operators O I . We can then construct a right-invariant metric g on U as follows: if X = dU dt is a tangent vector to U at some point U , then we can define a corresponding Lie algebra element XU −1 . Then the metric is defined by To define the cost factors, let us split our generators O I into O α ∈ R and Oᾱ ∈R, where O α are the relevant generators which simultaneously act on multiple factors, while Oᾱ are irrelevant and act within individual factors. Then we can simply take the metric G IJ to be of the form in (3.5) with the cost factors given by where is a small parameter that will be taken to zero at the end. We now define the Nielsen binding complexity (or simply binding complexity, for brevity) C b of a unitary U as the minimal distance between U and the identity with respect to the above metric, in the limit → 0. Many different unitary operators may prepare the same state -e.g., we can always multiply one such unitary by others that rotate the part of the Hilbert space that is orthogonal to the reference state. Consequently, the complexity (binding or otherwise) of a state as opposed to an operator is defined as the complexity of the simplest unitary operator preparing that state.
In the examples we study it will turn out that there is a unique operator preparing each state, so we can avoid this subtlety.
From this perspective the binding complexity C b of a unitary operator is its minimal geodesic distance from the identity in the metric (3.6) [6][7][8]. 8 For group manifolds with rightinvariant metrics of the form discussed above, the geodesic equation takes a simple form, often referred to as the Euler-Arnold equation (perhaps familiar from rigid-body dynamics). Let U (s) be a geodesic on U, and let v(s) = dU ds U −1 ∈ u be the velocity vector pulled back to the identity. Then the Euler-Arnold equation is where we have used the invariance property of the Cartan-Killing form, explained in footnote 6. Alternatively, if we define L = I, Note that it is crucial that the structure constants mix generators with different cost factors for the term on the right to survive. In order to obtain the geodesics, we must solve equation (3.9) for the velocity v I . We then use this to obtain the geodesic, which satisfies and implement the boundary conditions U (0) = 1 and U (1) = U , where U is the unitary whose complexity we wish to study.

Complexity of Gaussian states
Our starting point is the toy model of [5], which takes a system of harmonic oscillators as an approximation to a free scalar field theory on an n-point lattice. Since we are interested in using this as a setting for the study of multiparty entanglement, we partition the oscillators into m groups of N oscillators each, so that n = N m. We will refer to each group of oscillators as a "party". The operator content of the theory are the "position" and "momentum" operatorsφ i andπ i at each site, with i = 1, 2, . . . , n and canonical commutation relation [φ i ,π j ] = iδ ij . We consider Gaussian states of the form: where d ϕ = dϕ 1 . . . dϕ n , | ϕ = |ϕ 1 ⊗ . . . ⊗ |ϕ n with |ϕ i an eigenstate ofφ i , Ω a symmetric matrix with positive-definite eigenvalues, and the coefficient out front is for normalization. The vacuum state has a specific Ω.
We will determine the binding complexity of such states with respect to a reference in which Ω is diagonal. The gate set for measuring complexity will consist of the Hermitian operators: 9Ô , so theÔ(A) operators form a representation of gl(n, R). Choosing as generators of gl(n, R) the elementary matrices (M ij ) k = δ ik δ j , we correspondingly define the generators of the gate set:Ô (That is, e i i,j ijÔ ij are the gates we use.) A short computation gives the structure constants To be clear, theÔ ij are the O I in the discussion above (3.1), where now I = ij is a doubleindex since we are working with a matrix Lie group. A unitary operator that prepares the general Gaussian state (3.12) from the reference state can then be reached from the identity by a continuous sequence of unitary operators, described by the path-ordered exponential where s parameterizes a trajectory in the space of unitary operators and the V ij (s) describe the instantaneous direction in the tangent space gl(n, R), i.e., "velocity" in the space of unitary operators. We pick the boundary condition so thatÛ (1) is the unitary operator that prepares the desired state.
To define binding complexity we follow the geodesic formalism described above. In terms of the non-degenerate metric on the space of generators (3.5), G IJ ≡ G ij,k , operator complexity is defined as the length of the geodesic trajectory connectingÛ (s = 1) to the identity, If there are multiple such geodesics, complexity is defined as the minimum of their lengths. The relevant and irrelevant operator directions are defined by the "costs" in the metric (3.7), so that G ij,k = (c ij + c k )K ij,k /2. Here c ij = 1 ifÔ ij ∈ R and c ij = 2 ifÔ ij ∈R where R is the set of operatorsÔ ij such that oscillators i and j are located in different parties. We take K ij,k to be the Cartan-Killing form for gl(n, R), where we have included an additional normalization factor of 1 2n for convenience as compared to (3.3). In the end, will be taken to zero and is included to make sure that G is nondegenerate.
A subtlety here is that the Cartan-Killing form for gl(n, R) has a degenerate direction, which in our notation reads 10 iÔ ii . This leads to a degeneracy in the metric, which we had wanted to avoid. Fortunately, the direction with a vanishing line element is irrelevant (i.e. it represents a gate acting within parties, as opposed to between them). So the degeneracy does not affect the binding complexity. However it can potentially lead to an ambiguity in the geodesic equation (3.8), because in the degenerate directions the equation becomes 0 = 0. Fortunately, in the rigid-body form (3.9), degeneracies arising from the Cartan-Killing form drop out allowing us to avoid this subtlety. Other than that, our metric is block diagonal (i.e., does not mix relevant and irrelevant directions), permutation-symmetric between parties, and only the relevant operators creating entanglement between parties contribute to binding complexity.
In the → 0 limit, the binding complexity is then To compute the velocities V ij (s) on geodesics we use the Euler-Arnold equation (3.9). In the present case this equation takes the form k, where the structure constants are given in (3.15), and the matrix I ij k = c ij δ i k δ j . 11 To solve (3.20), we must consider two cases: either i, j are in the same party, or they are in different parties. The resulting equations are These are in general solved by: where the v ij are integration constants. We are going to choose final states that are symmetric between the parties just like the initial states. Thus we expect to find a geodesic that is permutation-symmetric between the parties, and also between the oscillators within each party. Enforcing this permutation symmetry, we take all v ii = a to be identical, as a consequence of which , and all v ij = c when i = j andÔ ij ∈ R (relevant operators). Therefore, by requiring total permutation symmetry, we have restricted the matrix of velocities to three independent parameters that determine the final unitary operatorÛ (s). Of course, permutation symmetry between parties as opposed to oscillators is not essential; for example, we could consider final states that are not symmetric in this way. In Appendix A we demonstrate how to compute binding complexity for a less symmetric case and conjecture a solution for the completely general case.
Since the parameters a, b, and c determine the operatorÛ (s) that evolves from initial state to final state, we fix them by specifying these boundary conditions on the wavefunction. Namely, we take the initial wavefunction to be determined by the matrix Ω (i) = diag(ω 0 , ω 0 , . . . , ω 0 ) and the final wavefunction to be determined by (3.25) Thus, the initial wavefunction is the product of Gaussians in every oscillator; it contains no entanglement. The final state wavefunction contains "couplings" ω of each oscillator to itself, couplings λ 1 between different oscillators in the same party, and couplings λ 2 between oscillators in different parties. The structure of the final wavefunction above is meant to be a permutation-symmetric toy model to mimic the structure of entanglement in a generic quantum field theory state, where if we partition our system into m parties (i.e., either subregions or boundaries in the multiboundary case), then the state will have some internal entanglement within each party, in addition to entanglement between different parties.
In Ω (f ) , the couplings λ 1 create the internal entanglement between the oscillators inside each block/party, while the couplings λ 2 create entanglement between different blocks/parties. Although the wavefunction does not have the expected "spatial locality" of a quantum field theory state within each party, this locality can be added to the wavefunction by further acting on it with local unitary transformations which act only within each block; since such unitaries do not change the binding complexity, they will not affect our result below. For illustration, in the N = 3 case, the matrix Ω (f ) takes the form Importantly, there are three independent couplings, matching the number of independent parameters of V mn : a, b, and c. We determine the velocities a, b, and c in terms of these couplings by examining how the matrix Ω flows under the infinitesimal action of the unitarŷ U (s). SinceÛ (s) does not take the wavefunction out of the set of Gaussian wavefunctions, we can label the state at an arbitrary time s as Over an infinitesimal parameter length ds, the state changes according to This follows because (3.28) is a Schrödinger equation, the solution of which for the operator U (s) is well-known to be the path-ordered exponential (3.16).
Using the expression (3.14) for theÔ ij operators andπ i = −i ∂ ∂ϕ i in the |ϕ i basis, the right-hand side becomes in this basis where V is the matrix of velocities V ij . The symmetry of both Ω and V has been used in deriving (3.29). The trace term can be absorbed into the wavefunction normalization, so the action of theÔ ij operators induces the following flow of the matrix Ω: Since Ω(s) has only three independent components ω(s), λ 1 (s), and λ 2 (s) by the ansatz (3.25), the matrix equation (3.30) reduces to the three independent equations The coefficients above have been derived by expanding (3.30) and counting the number of terms of each type. Solving with the boundary conditions Ω (i) and Ω (f ) specified earlier by taking ω(1) = ω, λ 1 (1) = λ 1 , and λ 2 (1) = λ 2 , we determine the constants a, b, and c. In terms of the three independent eigenvalues (λ + , λ 0 , λ − ) of Ω, Plugging into (3.19) and counting the number of relevant operatorsÔ ij ∈ R, the binding complexity of the general Gaussian wavefunction is therefore Unlike conventional circuit complexity [5], the binding complexity as computed here is finite in the N → ∞ continuum limit of a large number of oscillators.
We can also write the binding complexity in terms of the dimensionless, UV-finite parameter µ = N λ 2 ω+(N −1)λ 1 as This parameterization is convenient because the entanglement entropy of a single party of oscillators relative to the rest is also controlled by µ. Using the method of Srednicki [53], it is straightforward to compute that this entanglement entropy is For S to be finite, we must have 1 1−m < µ < 1; if we require S to remain finite in the large N limit, this similarly constrains λ 2 λ 1 . At the points µ = 1 and µ = 1 1−m , the entanglement entropy associated with a single party as well as the binding complexity both diverge. Expanding about either point, where the argument of the logarithm in (3.41) becomes large, as does the macroscopic entanglement entropy (i.e. we are at high temperature), we find that the binding complexity and entanglement entropy are related as (also see Fig. 4): where S i refers to the entanglement entropy associated to the ith party. (For our symmetric wavefunctions all S i = S are equal). That is, the binding complexity scales linearly with the entanglement entropy, up to a constant term and corrections exponentially small in the entropy. As we will discuss below, this scaling of binding complexity with entropy resembles expectations from holographic duality.
That the binding complexity scales linearly with the entanglement entropy with both controlled by the same parameter µ is remarkable. Nevertheless, as discussed in the introduction, the single-party entanglement entropy may yield a misleading characterization of the robustness of entanglement in quantum states. To diagnose this robustness in the Gaussian states (3.12), we use the Peres-Horodecki separability criterion as written by Simon [54]. This criterion is a necessary and sufficient condition for separability of a two-oscillator Gaussian state. Therefore, for the remainder of this section we work in the special case N = 1, so that there are m total oscillators with a single oscillator in each of the m parties. We will check the separability of the reduced density matrix upon tracing out m − 2 parties.
We will call N g the Gaussian negativity. Taking N = 1 and tracing out m − 2 parties yields a state ρ on two oscillators for which N g evaluates to where µ = λ 2 ω is the N → 1 limit of the same parameter µ previously defined above (3.41). Recalling that 1 1−m < µ < 1 for the entropy and binding complexity to be finite, we see that this same condition leads to N g > 0. We conclude that the Gaussian states (3.12) for N = 1 are robustly entangled like the W states. When N > 1, the condition N g ≤ 0 is no longer equivalent to separability [55]. 12 However, the similar structure of the wavefunction leads us to expect that the states will remain robustly entangled when N > 1.
Since N g is also controlled by the parameter µ, we may again expand about the point where the binding complexity becomes large to find that the binding complexity scales linearly with the logarithm of the Gaussian negativity up to exponential corrections (see Fig. 5),   Since the single-party entanglement entropy is controlled by the same parameter µ as the binding complexity C b and the Gaussian negativity N g , it is not obvious if a large binding complexity ultimately stems from a robust entanglement structure rather than merely a large entanglement entropy. To address this question, Fig. 6 shows that even at fixed entropy S, the binding complexity per party C b /m increases with ln N g . Consequently, binding complexity does diagnose robustness of entanglement.

The interior volume of multiboundary wormholes
Multiboundary wormholes [38][39][40][41][42][43] are vacuum solutions of Einstein's equations in 2+1 dimensions that have multiple asymptotic regions (Fig. 7). Recently, properties of these geometries were used in [33,44,45] to investigate the entanglement structure and complexity of the boundary CFT state. Tensor network models for multiboundary wormholes were presented in [46]. Like the two-sided BTZ black hole [57], the multiboundary wormholes are constructed as quotients of AdS 3 space. On the t = 0 slice, AdS 3 is just hyperbolic space H 2 , which has an isometry group PSL(2, R). The t = 0 slice of the wormhole is obtained by quotienting this H 2 by a discrete diagonal subgroup Γ ⊂ PSL(2, R) with hyperbolic generators. The action of Γ will identify pairs of boundary-anchored geodesics in H 2 , so M = H 2 /Γ will be a Riemann surface with m boundaries (each topologically S 1 ), where m − 1 is the number of generators of Γ. Since any two disjoint boundary-anchored geodesics in H 2 are joined by a unique minimal length geodesic, the endpoints of the latter join to form causal horizons for the newly disjoint conformal boundary. The set of causal horizons bounds the interior of a wormhole that connects all the asymptotic regions together. A holographic observer with access to observables on just a single boundary cannot access physics in the wormhole interior. It was shown in [33] that the CFT state dual to these wormholes contains multipartite entanglement between degrees of freedom localized on the different boundaries.
Following [46], we can think of the complexity of the quantum state dual to a wormhole in holographic terms by imagining a tensor network that tiles the bulk Cauchy slice. Such a tensor network will prepare a state with the necessary pattern of entanglement (Fig. 8a). The complexity of the state is then proposed to be related to the number of gates in the tensor network [21,44,45,58,59], an idea which correlates nicely with the proposal that complexity is holographically dual to the volume of spatial slices [58]. In this tensor network construction of the boundary state dual to the wormhole, the tensors outside the horizons correspond to unitary operations acting within each boundary (Fig. 8b). On the other hand, tensors enclosed within the wormhole interior may be thought of as corresponding to unitary quantum gates acting simultaneously on multiple boundaries. Thus, we might expect that binding complexity corresponds to the interior volume of the wormhole.
To make this comparison, we compute the interior volume of the wormhole. Since all of our calculations pertain to an equal-time slice of the 2 + 1-dimensional spacetime, the volume of the interior is really an area. It is easy to compute this area using the Gauss-Bonnet theorem in terms of the number of asymptotic boundary regions m, the genus of the interior g, and the geodesic curvature of each causal horizon. The interior W g,m is topologically a Riemann surface of genus g with m punctures, and the area of this surface, with the constant curvature metric inherited from H 2 , is given by where the second term on the right hand side is the integral of the geodesic curvature on the boundary of the interior.
We will set g = 0 for simplicity, i.e., the wormhole has a spherical internal topology. 13 If we choose the interior region to end strictly at the causal horizons (which are geodesic), then the geodesic curvature term vanishes. In this case the area (4.1) vanishes for the BTZ black hole (which has m = 2). 14 Nevertheless, we know that there is bipartite entanglement between the two boundaries of BTZ, and there will be an associated binding complexity. Thus the interior volume on the t = 0 slice cannot be literally equal to complexity.
In view of this, we are led to consider "stretched horizons", which are non-geodesic curves pushed slightly away from the true horizons in the full wormhole geometry toward the asymptotic boundaries (see [60] and references therein). In the tensor network picture of complexity, we interpret this procedure as including tensors just outside the horizons which still contribute to the entanglement between multiple boundary CFTs, c.f. Fig. 8b. This interpretation is substantiated by [46], which showed that for tensor network models built by quotienting the networks preparing vacuum AdS 3 states, it is possible for an "bipartite residual region" to remain after entanglement distillation in the m = 2 case. We are thinking 13 In principle, we could extend our toy model construction to higher genus by considering states of more complicated entanglement structure (see Sec. 5). For example, the (m = 4, g = 1) case might correspond to removing the y2-y4 and y1-y3 edges in Fig. 9c. The Euler-Arnold equation in the general form of this case becomes very difficult to solve, but such a calculation could serve as another check of our proposal that the (stretched) interior volume equals binding complexity.
14 This is because the causal horizons of the two asymptotic regions of the eternal BTZ black hole coincide on the t = 0 surface at the bifurcation point of the horizon, so that, unlike the multiboundary case, the internal volume vanishes.
of these residual tensors as living inside the stretched horizon. We will take the stretched horizon to be a surface of constant geodesic curvature k g .
In sum, we obtain a contribution to the area that is proportional to the length of each stretched causal horizon Here we used the fact that the horizon lengths are equal to 4G times the entropies of entanglement of the CFT on the i th boundary with all the other boundaries. 15 The O( −1 AdS ) constants a i are given in terms of the horizon lengths by where L i is the horizon length and ∂ i W 0,m is the i th boundary of the interior.
The formula (4.2) for the volume of the (stretched) wormhole is structurally similar to the formula (3.43) for binding complexity. Both expressions have a constant piece, and a part that is linear in the entanglement entropies of each disconnected party. Thus it is tempting to propose the correspondence 16 Binding Complexity = Volume of Stretched Wormhole Interior and vanishes for m = 2. The origin of this discrepancy may lie in the simplicity of the toy model of Sec. 3 and might be resolved with an appropriate generalization of the framework for computing binding complexity of states in a nontrivial conformal field theory with a semiclassical bulk dual. However, it also might simply be that the toy model states whose complexity we considered were not structured in the same way as in holographic theories. In Sec. 5 we will provide evidence that the latter is indeed the case by using the Euclidean path integral to construct a natural class of states in our toy model whose binding complexity reproduces the form of the stretched wormhole volume. Indeed in AdS 3 /CFT 2 [33] precisely such a Euclidean procedure constructs the CFT states dual to the multiboundary wormhole. 15 Note that assumes that we are in a region of the moduli space for the interior geometry of the wormhole where the entropy of each boundary is holographically given by the causal horizon separating it from the other asymptotic regions. Remarkably, there are other regions of the moduli space where the entropy of boundary i is actually give by the sum of areas of the causal horizons of all the other boundaries. This surprising fact is explained in [33]. 16 The volume here is being expressed in units of 1 AdS G N , as is usual in discussions of complexity.

Euclidean path integrals
In the previous section we argued that the binding complexity of Gaussian states that we calculated in Sec. 3 resembles the volume of the interior of multiboundary wormholes in AdS 3 /CFT 2 . However, there was a discrepancy in the scaling of the complexity with the number of entangled parties m which could arise if the permutation-symmetric states of Sec. 3 do not have the same entanglement structure as the states in AdS/CFT. In the AdS setting, the states dual to multiboundary wormholes can be constructed within the CFT by performing the Euclidean path integral on 2-manifolds with the topology of the bulk wormhole (i.e., the time-reflection symmetric Cauchy surface in the bulk) [33]. 17 To compare with the wormhole it would therefore be natural to compute the complexity of states in our toy model constructed in terms of similar Euclidean path integrals. In our case we have a collection of n harmonic oscillators. So, we should perform a path integral on a (0 + 1)-dimensional graph with n external legs. As will see, the binding complexity depends on the topology of the Euclidean graph.
A general 1D Euclidean path integral for a system of n = N m harmonic oscillators is computed on a graph G consisting of a set of vertices V G , n of which are external, and a set of edges E G each of different lengths. Such a graph may contain internal vertices. The value of the oscillator field at these vertices is a boundary condition which must be matched in the propagators at all incoming edges and integrated over. Each edge (v 1 , v 2 , β) of length β between vertices v 1 , v 2 at positions x 1 , x 2 respectively in the graph corresponds to a factor of the propagator K(x 1 , x 2 , β) in the integrand: where β is the length of the edge in the graph and M is oscillator mass. The Euclidean propagator for the harmonic oscillator can be computed exactly; it is a Gaussian function known as the Mehler kernel: Let us label all external vertices by the vector x and internal vertices by the vector y. The wavefunction of a state prepared by the Euclidean path integral on the graph G is therefore Since the propagator is Gaussian, the end result of the integrals over the internal vertices is also a Gaussian wavefunction, which can always be written where Ω is a real symmetric matrix and N is a normalization constant. Consequently, we may bring to bear the technology of Sec. 3 in computing the binding complexity.

Permutation-symmetric graphs
We are interested in the complexity of states in which the different parties are multiparty entangled. It is natural to imagine that such entanglement is produced in the Euclidean path integral if the graph is branched so as to connect between the parties. In Sec. 3 we considered states (3.25) in which the oscillators within parties were entangled with one strength, while the parties as a whole were entangled block-wise with other parties and with a different strength. We will first see how to construct such permutation-symmetric states through a Euclidean path integral.
In the Euclidean path integral, oscillators become entangled if their propagators meet at a vertex where a shared boundary condition is integrated over. This suggests that to construct the states in the previous section we need a graph with m groups of N external lines that each meet at a vertex to create the internal entanglement within parties. These vertices can then be connected by further propagators to create entanglement between the parties. Three such graphs are shown in Fig. 9. We label the vertices at the end of the external lines as x i j for the ith oscillator in the jth party. The internal vertices can have any number of lines ending on them -the internal structure of the graph can be completely arbitrary up to the permutation symmetry of the state we are trying to construct. In analogy with the holographic setting, we might refer to the internal part of the graphs in Fig. 9 as a "wormhole" connecting the exterior legs.
First consider the simplest graph Fig. 9a. The internal vertices on the ith branch are labeled y i , and the central vertex is labeled y c . We integrate over the boundary condition of the field at each vertex to perform the path integral. The lengths of the edges are moduli of the graph, and the wavefunction generated by the path integral is a function of these parameters. Permutation symmetry of the states (3.25) dictates that the external lines have the same length (β 1 ) and the internal lines have the same length (β 2 ). Similarly, Figs. 9b, 9c have three moduli. Performing the path integral on the family of graphs of Fig. 9a according to the procedure of (5.3), one obtains a Gaussian state (5.4) in the permutation-symmetric form (3.25) with parameters ω, λ 1 , λ 2 where ω and λ 1 quantify entanglement within a party and λ 2 quantifies entanglement between parties. We find that (see Appendix B for details) Since this is a permutation-symmetric Gaussian state, the binding complexity is given by (3.43) and the entanglement entropy of a single party is given by (3.42). Both quantities vanish in the limit N → ∞ with β 1 , β 2 , M fixed since λ 2 (which quantifies entanglement between parties) scales as 1/N 2 at large N . This disentangling at large N can be understood as a manifestation of the principle of entanglement monogamy: when the number of oscillators within a party grows large, most oscillators are entangled within their party rather than with other parties. We can compensate by taking a kind of 't Hooft limit in which β 1 N is held fixed as N → ∞, in which case both the binding complexity and entanglement entropy will be finite and nonzero since λ 2 approaches a finite value in the large N limit. The latter scaling limit can also be thought of as a rescaling of the couplings with the lattice scale so that the couplings remain finite in the continuum for a lattice quantization of scalar field theory. Fig. 10 shows the moduli dependence of the binding complexity for the graph Fig. 9a as computed in (5.5)-(5.7). The complexity increases as β 1 , β 2 become small. This is because as β → 0 the propagator in (5.2) becomes the identity, thus more closely coupling the values of the fields at either end of a line in the graph. In the other limit, as β → ∞, the propagator projects onto the ground state, essentially decoupling the external oscillators from the internal structure of the graph. Finally, consider wavefunctions associated with the graphs Fig. 9b and Fig. 9c. Because all the integrals are Gaussian, we will again get Gaussian wavefunctions and because the graphs are permutation-symmetric, the wavefunctions will be as well. Of course, the coefficients in the wavefunctions will contain different functions of the moduli in each case because the detailed integrals are different. However, all of these wavefunctions are necessarily of the form (3.25), and therefore the constant term in the binding complexity will scale as the logarithm of the number of entangled parties, unlike the linear scaling with parties of the interior volume of multiboundary wormholes.

Bipartite entanglement graphs
We would like to find graphs that generate states with complexity-entropy scaling relations that match the holographic form. First note that the scaling relation (3.43) between complexity and entropy holds in the large β limit in which the entropy associated with any single party is large. It was shown in [44] in the holographic setting that in this regime, the entanglement structure of the multiboundary wormhole is dominated by bipartite entanglement between boundaries. 18 Consequently, to better match the holographic expectations, we seek graphs on which the path integral will produce strongly bipartite entanglement. There is an independent reason to be interested in such graphs: in the "bit thread" interpretation of holographic entanglement entropy [61,62] one expects the correlations between independent tensor factors of a CFT to be dominated by bipartite entanglement (i.e., between the two qubits connected by a bit thread). In our setup, the mixing of terms in the wavefunction is dictated by topological connectedness in the graph on which we perform the path integral. Therefore, we engineer multipartite entangled states with locally bipartite entanglement structure by using graphs which factorize so that a given connected component of the graph connects only two parties.
In Fig. 11a, we display such a "bipartite entanglement graph", in which the oscillators in each party have been partitioned into groups that are only entangled with oscillators in one other party. The overall graph factorizes into a collection of the two-party permutationsymmetric graphs of Sec. 5.1. Let N = (m − 1)k be the number of oscillators in each party, where m is the total number of parties and k is the number of oscillators per grouping, so that each of the k groups connects to a different one of the other m − 1 parties (see Fig. 11 for details). We again choose β 2 to be the length of internal lines and β 1 to be the length of external lines. As drawn in Fig. 11a it appears that only part of each party is connected to part of another party. However, as before, one may always mix the oscillators in a single party The graph corresponding to ψ branch,k , with k oscillators in each of the two parties. via local unitaries which will not affect the binding complexity or the entanglement entropy associated with that party. Therefore, we may think of Fig. 11a as encoding locally bipartite entanglement between parties without restricting the entanglement to reside in some subsystem of each party. In other words, in Fig. 11a we have used local unitary transformations to "diagonalize" the entanglement structure in each party.
Since the manifold on which we are performing the Euclidean path integral is topologically disconnected, the path integral factorizes over the connected components, as does the resulting wavefunction. Consequently, ψ( x) is the product of m 2 = m(m−1) 2 permutation-symmetric wavefunctions. Let ψ branch,k be the wavefunction of the graph in Fig. 11b, which has k oscillators in each party. This is a permutation-symmetric graph as described in Sec. 5.1, so ψ branch,k is a two-party permutation-symmetric wavefunction. Then the wavefunction of the full bipartite entanglement graph can be explicitly written as where the product includes m(m − 1)/2 such terms corresponding to the bipartite connection between each pair of parties. The total wavefunction is still Gaussian and takes the form of (5.4), but Ω is no longer permutation symmetric within each party. Since we have the freedom to relabel oscillators so that topologically connected vertices are ordered adjacently in the matrix, Ω takes a block-diagonal form, consisting of m(m − 1)/2 identical permutationsymmetric subblocks each of size 2k × 2k. Each subblock is of the form (3.25) with the couplings ω, λ 1 , and λ 2 given by (5.5) -(5.7) with the replacement m → 2 and N → k. The structure of the matrix Ω in the special case m = 3, N = 4, k = 2 is shown below, where the solid lines demarcate parameters corresponding to the same party and dashed lines demarcate parameters corresponding to topologically connected oscillators: We have not yet relabeled oscillators above to bring Ω into block diagonal form, so that the grouping of oscillators in each party is clearer.
Although Ω is not permutation-symmetric, each of its subblocks is permutation-symmetric. Consequently, the entanglement entropy associated with a single party is where S branch,k refers to the entanglement entropy associated with a single party of the wavefunction ψ branch,k . Equation (5.12) follows automatically from the factorized form of the graph as shown in Fig. 11a: the wavefunction splits over each component in the graph, so the total entanglement entropy is the sum of the entropies of each component. 19 In other words, the entanglement entropy associated to a single party essentially counts the minimal number (m − 1) of edges which are "cut" in separating the oscillators in that party from the rest. This continues to hold for the entropy associated with other partitions: the prefactor m − 1 in (5.12) changes to the minimal number of edges cut in separating those parties from the rest. It is tempting to compare this result to bit threads and to the tensor network picture of holographic entanglement entropy, in that the entropy associated to a given party is directly proportional to the number of "threads" leaving that party. This counting property of entropy is thought to underlie the Ryu-Takayanagi formula for holographic entanglement entropy.
Upon tracing out m − 2 parties, the reduced density matrix ρ associated with two parties of a bipartite entanglement graph has a robust W-like entanglement structure. From the product structure of the wavefunction (5.10), it follows that ρ takes the schematic form ρ = ρ mixed 1 ⊗ ρ 12 ⊗ ρ mixed 2 , where the subscripts refer to the first and second party. Here ρ 12 is a pure state corresponding to the two-party permutation-symmetric graph that connects a single group of oscillators in each of the two parties, while ρ mixed 1 refers to the complicated mixed state of the remaining oscillators in the first party and similarly for 1 ↔ 2. In Sec. 3 we argued that a permutation-symmetric state like ρ 12 has a robust entanglement structure, so ρ will demonstrate this structure as well. In Fig. 11a and in (5.10), we have picked an adapted basis that has separated the oscillators in such a way that upon doing partial traces, the degrees of freedom that remain entangled are distinct from the degrees of freedom that are in a mixed state. In general, we can act with local unitary transformations so that all the degrees of freedom retain both entanglement and mixedness.
The binding complexity of these graphs is The complexity-entropy scaling relation that follows from (5.12) and (5.13) is Comparing to (4.2), one sees that the constant term now scales with the number of boundaries in the same O(m) fashion as the holographic expectation, at least in the large m limit, adding support for the idea that Binding Complexity = Wormhole Volume.

Complexity for coherent states in perturbation theory
The Nielsen formalism also allows us to compute how much the complexity of a state changes when it is perturbed. For example, suppose we want to compute the Nielsen complexity of a state of the form |ψ = e it I h I O I |ψ 0 , (6.1) relative to the base state ψ 0 , where t will be treated as a small parameter in which we do perturbation theory. In other words, we are interested in studying the complexity of the unitary operator U = e it I h I O I perturbatively in t. This situation can arise in several contexts; for example if we treat t as time and H = I h I O I as a Hamiltonian, then we obtain the small-time behavior of the complexity of time evolution. Alternatively, we may treat U as creating a coherent state on top of some base state ψ 0 , and t may be a small parameter which controls the size of the background deformation, as in the next subsection.
As before, we will take the operators O I to form the Lie algebra As discussed previously, we need to define a positive-definite, bilinear form G IJ on this algebra, which fixes the complexity of individual gates. We will be interested in where v is the local velocity and P is path-ordering. The geodesic equation in terms of the velocity is given by the Euler-Arnold equation The boundary conditions are where λ is a small parameter. For states of this form, we can solve the equations in perturbation theory with respect to λ. So let us take v I (s) = tv I (1) (s) + t 2 v I (2) (s) + t 3 v I (3) (s) + · · · . (6.6) First order in t: At leading order, the equation can be solved trivially: Therefore, the unitary U (1) is given by Comparing this with (6.5) at first order, we deduce that v I (1) = h I . (6.9) We can now compute the binding complexity of this state as the geodesic distance:  So, now the unitary becomes (6.14) Once again, comparing with equation (6.5), we find 15) and therefore to this order the velocity is then given by We can now use this result to compute the O(t 2 ) correction to the complexity, and we find that the O(t 2 ) contribution vanishes after performing the s-integral. Therefore, we obtain We can proceed in a similar fashion to obtain higher order corrections, for instance, the O(t 3 ) correction is shown in Appendix C. We see that for small t the binding complexity of the unitary U = e it I h I O I increases linearly in t, with the proportionality constant being the norm of the Hamiltonian H = I h I O I in the multipartite sector, that is, only the relevant operators which act simultaneously on multiple factors are included in the norm.

Double-trace deformations: towards creating wormholes
We can now how ask how the binding complexity changes if we perturb a state by acting with an operator that locally couples degrees of freedom in two distinct parties. In the holographic context this sort of "double-trace deformation" was shown in [47] to create or expand a wormhole in the geometric description of disconnected but entangled CFTs. Once again we consider the toy model with n free, decoupled harmonic oscillators φ i , with the Hamiltonian We can straightforwardly diagonalize H 0 by introducing the creation and annihilation op- The vacuum state for this Hamiltonian, which satisfies is a completely decoupled product state, and as such it will have no binding complexity. We now deform the Hamiltonian by a small bilinear coupling with The coupling clearly introduces some entanglement and binding complexity in the new vacuum; our aim here is to compute this binding complexity perturbatively in g. In order to diagonalize the new Hamiltonian H, let us introduce the orthogonal matrix V ij which diagonalizes M ij = δ ij + gC ij : Here ω 2 i are the eigenvalues of M . Then, we define the new operators which also satisfy the appropriate bosonic commutation relations. In terms of these new variables the full Hamiltonian becomes

Now diagonalize this Hamiltonian by introducing the new creation and annihilation operators
We can express these new creation and annihilation operators in terms of the old creation and annihilation operators as We can represent this Bogoliubov transformation in terms of conjugation by a unitary operator: where the real, anti-symmetric matrix v ij is defined as V = e v . Therefore, the new vacuum ψ in presence of the bilinear interaction can be related to the old vacuum ψ 0 as where We can also re-express this state in terms of the gl(n, R) generatorsÔ ij = 1 2 φ iπj +π jφi , which were discussed in Sec. 3.1: The binding complexity of the state can now be computed perturbatively in g, following our discussion in the previous section. The leading order contribution is where B = i,j B ijÔij and the · · · indicate higher order corrections which enter at O(g 3 ) (as discussed in the previous section).
This result shows that adding "double-trace deformations" to the Hamiltonian creates binding complexity in the vacuum. If binding complexity measures the interior volume of wormholes, our result implies that the deformation has created a wormhole where none previously existed. This is in analogy with the holographic results of [47] where double-trace deformations of a product of CFTs enlarged a wormhole between the corresponding geometric asymptotic regions. We computed our results above in a toy model of oscillators, but we expect that a similar calculation will go through in the case of generalized free fields describing the large N limit of CFTs, which is the limit in which field theories are holographically described by classical geometry.

Discussion
We have suggested an interpretation for the volume of multiboundary wormhole interiors in AdS/CFT in terms of the binding complexity of the dual state. However, our discussion was limited to the interior volume of the time-reflection symmetric Cauchy slice in the bulk. If we consider a generic Cauchy surface ending at the times (t 1 , · · · t n ) on the boundaries, then the volume of the wormhole interior will, in general, be larger. However, the binding complexity should be independent of the times t i because changing these times simply corresponds to local Hamiltonian evolution in the different CFTs, and does not add any entanglement. This observation suggests that the covariant version of the bulk dual to binding complexity should be given by minimizing the interior volume over all the bulk Cauchy surfaces and over the different boundary times {t i }. Note that if we consider the maximum volume slice in the bulk ending at the times t i , then its volume is expected to be dual to the total complexity of the boundary state, which indeed depends on the t i because local Hamiltonian evolution adds to the total complexity. However, the corresponding circuit is not the minimal one from the point of view of binding complexity. Fig. 12 illustrates that the maximal volume Cauchy slice in the two-sided wormhole corresponding to the BTZ black hole can have a large interior volume, but it is always possible to find a different Cauchy slice that passes through the bifurcation surface. Figure 12: Cauchy slices of maximal volume (blue) and of minimal interior volume (red) in the BTZ geometry, both anchored at boundary times t 1 in the left CFT and t 2 in the right CFT. The volume of the interior of the maximal slice (dark blue) increases over time, but the corresponding circuit does not minimize binding complexity.
The relation between binding complexity and wormhole interiors was most concrete for certain states created by performing the Euclidean path integral on a graph with locally bipartite connections between parties, but which can nevertheless have multipartite entanglement. This occurs if some local degrees of freedom in each party have bipartite entanglement with local degrees of freedom in different parties. This is a structure resembling the W-state on qubits (1.2). However, we know that states with holographic duals satisfy the additional condition that mutual information is monogamous [63], implying that it is of the perfect tensor type [64]. In the bit-thread picture of entanglement, it seems necessary to sum over different bit-thread configurations to achieve this constraint [61,65]. In our picture this would mean summing over multiple (perhaps all) Euclidean graphs that produce states on a given partition of external variables. It would be interesting to consider the binding complexity for these kinds of states -it is not obvious that the complexity will simply be a weighted sum of the complexities of the individual graph states.
Our notion of binding complexity has similarities to the idea of quantum communication complexity, where several independent parties attempt to collaborate on some particular computation. 20 We can define the quantum communication complexity of a task to be the minimum number of qubits that must be exchanged between all the parties in order to complete the computation. Binding complexity measures a similar quantity, namely the number of gates that affect more than one party's qubits. In this way, both binding and quantum communication complexity increase as the computation requires more cooperation or interaction between the parties. In fact, we can obtain a strict relationship between the two quantities. Suppose all the gates in an n-qubit quantum circuit U are k-qubit gates. Then the quantum communication complexity of applying U to some distributed set of qubits is bounded by the binding complexity, since we may always transmit qubits across party lines in order to apply one of our gates. If the distributed parties run into a gate that contributes to the binding complexity during the application of U , they may simply communicate all the qubits to one of the involved parties, apply the unitary locally, and then send the qubits back to their proper owners (the bound is improved by a factor of 2 if we drop this last requirement). Each cross-boundary gate therefore contributes a maximum of 2k to the quantum communication complexity, and we obtain the upper bound Note that if we described this in a holographically dual geometry, the required multiboundary wormhole need not be traversable -there is no wormhole-based "quantum FedEx" that would allow qubit transfers between the different boundaries, which we are treating as the distributed parties attempting to build the unitary U . However, we can obtain a bound on the communication complexity of the problem by studying this geometry, assuming our conjecture holds. It would be interesting to make this analogy between binding and quantum communication complexity more precise in holography.

A Binding complexity for more general states
In this appendix, we compute the binding complexity for a state with less symmetry than that of (3.25). This will be an educational exercise that suggests a solution procedure for a totally arbitrary state. Consider a wavefunction for a four-party state taking the general Gaussian form (5.4) with the matrix Ω taking a block structure like As in (3.26), each entry above is an N ×N matrix, where N is the number of oscillators on each boundary. The elements λ (i) 2 are the matrices all of whose elements are couplings similarly labeled λ 2 refers to the coupling, not the full matrix). The elements ωλ 1 are matrices that are ω on the diagonal and λ 1 on all off-diagonals. This Ω is not completely general: in the language of Sec. 5, it corresponds to the path integral on a graph with Z 2 × Z 2 symmetry between the four parties.
The solution of the Euler-Arnold equation (3.20) is independent of the structure of the wavefunction, so the velocity matrix V again is constant. In general, one can show that choosing the structure of the velocity matrix V to have the same form as Ω will allow for solution of the flow equation (3.30). Consequently, we choose V to take the same form as (A.1) with a replacing ω, b replacing λ 1 , and three cross-party velocities c 1 , c 2 , c 3 replacing λ where Ω = ω λ 1 λ T arranges the s-dependent couplings of the matrix Ω into a vector and For comparison, note that the equations (3.31) -(3.33) can be written as a similar 3 × 3 matrix equation. The matrix M has five distinct eigenvalues: Solving (A.2) with the usual boundary conditions of Ω (i) = diag(ω 0 , ω 0 , . . . , ω 0 ) at s = 0 and Ω (f ) as given by (A.1) at s = 1, we find ω = ω 0 4N (4(N − 1)e κ 0 + e κ 1 + e κ 2 + e κ 3 + e κ + ) (A.9) We remark that the five distinct eigenvalues of Ω are given by: 2 ) (A.15) ρ 2 = ω + (N − 1)λ 1 + N (−λ 2 ), (A. 18) closely related to the eigenvalues of M . In the permutation-symmetric limit, ρ + corresponds to λ + , ρ 0 corresponds to λ 0 , and ρ 1 , ρ 2 , ρ 3 all approach λ − . Solving the system (A.9) -(A.13) for the eigenvalues κ we obtain Rewriting the c k that determine the binding complexity in terms of the eigenvalues ρ i , Lastly, a short combinatorial computation determines the binding complexity where |c| = (c 1 ) 2 + (c 2 ) 2 + (c 3 ) 2 . Notice that again the prefactor of N above cancels the N dependence of the c k so that the binding complexity is finite in the large N limit.
Unfortunately, the binding complexity does not arrange nicely in terms of a parameter µ as in Sec. 3, and it is prohibitively diffcult to evaluate the entanglement entropy associated with a single party of the state specified by (A.1) to obtain a complexity-entropy scaling. Nevertheless, this computation is instructive to understand how to compute the binding complexity for a (more) general Gaussian state. In general, we expect that if we arrange the couplings in Ω into a vector Ω, the eigenvalues ρ of the matrix Ω will be some linear combination of the couplings: ρ = A Ω. In this case, choosing V to have the same matrix structure as Ω gives rise to a lower-dimensional matrix equation for the couplings in terms of a matrix M = 2A. The eigenvalues κ of M will be κ = M V , and the solution for the velocities will looks like V i = 1 2 j (A −1 ) ij ln ρ j ω 0 . Note that V is the vector of velocities analogous to Ω. Once the velocities are obtained, it is straightforward to compute the binding complexity based on the particular combinatorics of a given setup.

B Wavefunctions of permutation-symmetric graphs
In this appendix, we compute the wavefunctions of the branched graphs presented in Fig. 9a, for an arbitrary number of parties m and number of oscillators per party N , working in the M → 0 limit for simplicity. In this limit, the propagator (5.2) remains Gaussian and takes the simple form K(x 1 , x 2 , β) ∝ e whereÑ is a normalization constant. Performing the Gaussian integral over y c , we obtain whereÑ is a new normalization constant. The remaining integral (B.5) is also Gaussian over the internal vertices y i , although it has a linear term. That is, it takes the form where the matrix A and vector B T are given by: That is, A takes value α on the diagonal and β on all off-diagonals. The constants α and β are given in terms of T 1 , T 2 , m, and N by The exact solution of the general Gaussian matrix integral of the form of (B.6) is well-known. Evaluating it gives whereÑ is another new normalization constant. The inverse of A has the same symmetry as A, with A −1 ij = P δ ij + Q(1 − δ ij ) and P = α + (m − 2)β α 2 + (m − 2)αβ − (m − 1)β 2 , Q = −β α 2 + (m − 2)αβ − (m − 1)β 2 .
(B.10) Therefore, in terms of the oscillator variables x i j , the wavefunction is Despite the cumbersome sum notation for the general case, one can check that this is indeed Gaussian and can be written in the standard Gaussian form ψ( x) =Ñ exp(− 1 2 x T Ω x) with Ω in the form of (3.25). To be completely explicit, Ω has the general permutation-symmetric form (3.25) with . (B.14) in agreement with the M → 0 limit of (5.5) -(5.7). This computation was entirely in the M → 0 limit, but the trick employed herein of rewriting the product over propagators as matrix Gaussian integrals works very generally. For any permutation-symmetric graph the computation goes through identically with possibly different values of α and β, even when M = 0.

C Perturbation theory to O(t 3 )
For completeness, we will show how to proceed at O(t 3 ) in the small-time perturbation theory in this appendix. At third order, we find Comparing with equation (6.5), this implies We can now compute the complexity at O(t 3 ), if we so desire.