Complexity of mixed Gaussian states from Fisher information geometry

We study the circuit complexity for mixed bosonic Gaussian states in harmonic lattices in any number of dimensions. By employing the Fisher information geometry for the covariance matrices, we consider the optimal circuit connecting two states with vanishing first moments, whose length is identified with the complexity to create a target state from a reference state through the optimal circuit. Explicit proposals to quantify the spectrum complexity and the basis complexity are discussed. The purification of the mixed states is also analysed. In the special case of harmonic chains on the circle or on the infinite line, we report numerical results for thermal states and reduced density matrices.


Introduction
The complexity of a quantum circuit is an insightful notion of quantum information theory [1][2][3][4][5][6]. During the last few years it has attracted increasing attention also because it has been proposed as a new quantity to explore within the (holographic) gauge/gravity correspondence between quantum (gauge) field theories and quantum gravity models from string theory. In this context, different proposals have been made to evaluate the complexity of a quantum state by considering different geometric constructions in the gravitational dual [7][8][9][10][11][12][13][14][15][16].
A quantum circuit constructs a target state by applying a specific sequence of gates to a reference state. The circuit complexity is given by the minimum number of allowed gates that is needed to construct the target state starting from the assigned reference state. This quantity depends on the target state, on the reference state, on the set of allowed gates and, eventually, on the specified tolerance for the target state. Notice that this definition of complexity does not require the introduction of ancillary degrees of freedom.
Remarkable results have been obtained over the past few years in the attempt to evaluate complexity in quantum field theories . Despite these advances, it remains an interesting open problem that deserves further investigations.
In order to understand the circuit complexity in continuum theories, it is worth exploring the complexity of a process that constructs a quantum state in lattice models whose continuum limit is well understood. The free scalar and the free fermion are the simplest models to consider. In particular, it is worth focussing on the Gaussian states because they provide an interesting arena that includes important states (e.g. the ground state and the thermal states) and that has been largely explored in the literature of quantum information [41][42][43][44][45]. The bosonic Gaussian states are particularly interesting because, despite the fact that the underlying Hilbert space is infinite dimensional, they can be studied through techniques of finite dimensional linear algebra.
Various studies have explored the complexity of quantum circuits made by pure Gaussian states in lattice models [17][18][19][20][21][22][23][24][25][26]. In these cases the gates implement only unitary transformations of the state. It is important to extend these analyses by considering quantum circuits that involve also mixed states; hence it is impossible to construct them by employing only unitary gates [6]. A natural way to construct mixed states consists in considering the system in a pure state and tracing out some degrees of freedom. This immediately leads to consider the entanglement entropy and other entanglement quantifiers (see [46][47][48][49] for reviews). The same consideration holds within the context of the holographic correspondence, where the gravitational dual of the entanglement entropy has been found in [50][51][52] (see [53][54][55] for recent reviews). Quantum circuits involving mixed states naturally lead to introduce the notions of spectrum complexity and basis complexity [28].
A method to quantify the complexity of circuits involving mixed states has been recently investigated in [23]. In this approach, the initial mixed state is purified by adding ancillary degrees of freedom and the resulting pure state is obtained by minimising the circuit complexity within the set of pure states. This procedure requires the choice of a fixed pure state to evaluated this circuit complexity for pure states.
In this manuscript we explore a way to evaluate the complexity of quantum circuits made by mixed states within the framework of the Information Geometry [56][57][58]. The method holds for bosonic Gaussian states and it does not require the introduction of ancillary degrees of freedom. It relies on the fact that, whenever the states provide a Riemannian manifold and the available gates allow to reach every point of the manifold, the standard tools of differential geometry can be employed to find the optimal circuit connecting two states. Since the pure states provide a submanifold of this manifold, this analysis also suggests natural quantum circuits to purify a given mixed state.
We focus only on the bosonic Gaussian states occurring in the Hilbert space of harmonic lattices in any number of dimensions. These are prototypical examples of continuous variable quantum systems; indeed, they can be described by the positions and the momenta, which are continuous variables. The bosonic Gaussian states are completely characterised by their covariance matrix, whose elements can be written in terms of the two point correlators, and by their first moments. The covariance matrices associated to these quantum states are real symmetric and positive definite matrices constrained by the validity of the uncertainty principle [41][42][43][44][45]. We mainly explore the bosonic Gaussian mixed states with vanishing first moments. This set can be described by a proper subset of the Riemann manifold defined by the symmetric and positive definite matrices [59][60][61][62][63] equipped with the metric provided by the Fisher information matrix [56,57,[64][65][66]. We remark that our analysis considers quantum circuits that are made by Gaussian states only. Despite this important simplifying assumption, the resulting quantum circuits are highly non trivial because non unitary gates are involved.
The manuscript is organised as follows. In Sec. 2 we introduce the quantities and the main results employed throughout the manuscript: the covariance matrix through the Gaussian Wigner function, the Fisher-Rao distance between covariance matrices and the corresponding geodesics, that provide the optimal circuits. The particular cases given by pure states, thermal states and coherent states (the latter ones need further results discussed in Appendix B) are explicitly considered. In Sec. 3 we provide explicit expressions to evaluate the spectrum complexity and the basis complexity, by employing also the first law of complexity [67,68]. The purification of a mixed state is explored in Sec. 4, where particular optimal circuits are mainly considered. In Sec. 5 we discuss some lower and upper bounds on the complexity. In Sec. 6 we focus on the circuits that do not contain pure states because they can be also parameterised through the entanglement hamiltonian matrices. The Gaussian channels underlying the optimal circuits are briefly discussed in Sec. 7. In Sec. 8 we describe the approach to the complexity of mixed states based on the purification of a mixed state through ancillary degrees of freedom. The last analysis reported in Sec. 9 focuses on the periodic harmonic chain in one spatial dimension and on its limiting regime given by the harmonic chain on the infinite line. Numerical results are reported both for some quantities introduced in the other sections and for other quantities like the mutual complexity for the thermofield double states and for the reduced density matrices. Finally, in Sec. 10 we summarise our results and discuss future directions.
Some appendices (A, E and D) contain the derivation of selected results reported in the main text and related technical details. Other appendices, instead, provide complementary analyses that expand the discussion of the main text, adding further results. In particular, in Appendix B we explore Gaussian states with non vanishing first moments, in Appendix C the Bures and the Hilbert-Schmidt distances are discussed, in Appendix F the complexity of the thermofield double states is explored and in Appendix G we describe the two particular bases employed in [23] to study the complexity of mixed states through the F 1 cost function.

Complexity as Fisher-Rao distance and the optimal path
In Sec. 2.1 we introduce Gaussian Wigner functions (defined in terms of the covariance matrix and of the first moment) to characterise a generic Gaussian state. The Fisher-Rao distance and other distances are defined in Sec. 2.2. In Sec. 2.3 we discuss the Williamson's decomposition of the covariance matrix, a crucial tool largely employed throughout the manuscript. The optimal circuit in the Fisher information geometry is analysed in Sec. 2.4. The special cases given by pure states and thermal states are explored in Sec. 2.5 and Sec. 2.6 respectively. Finally, in Sec. 2.7 some results about the complexity for the coherent states are discussed.

Gaussian states in harmonic lattices
The hamiltonian of a spatially homogeneous harmonic lattice made by N sites with nearest neighbour spring-like interaction with spring constant κ reads where the second sum is performed over the nearest neighbour sites. The position and momentum operatorsq i andp i are hermitian and satisfy the canonical commutation relations [q i ,q j ] = [p i ,p j ] = 0 and [q i ,p j ] = iδ ij (we set = 1 throughout this manuscript). The boundary conditions do not change the following discussion, although they are crucial to determine the explicit expressions of the correlators. Collecting the position and momentum operators into the vectorr ≡ (q 1 , . . . ,q N ,p 1 , . . . ,p N ) t , the canonical commutation relations can be written in the form [r i ,r j ] = iJ ij , where J is the standard symplectic matrix being 1 the N ×N identity matrix and 0 the matrix with the proper size having all its elements equal to zero. Notice that J 2 = −1 and J t = J −1 = −J.
The density matrixρ, that characterises a state of the quantum system described by the hamiltonian (2.1), is a positive definite, hermitean operator whose trace is normalised to one. When the state is pure, the operatorρ is a projector.
A useful way to characterise a density matrix is based on the Wigner function w(r), that depends on the vector r made by 2N real components. The Wigner function is defined through the Wigner characteristic function associated toρ, that is [43][44][45] χ(ξ) ≡ Tr ρ e ir t J ξ = Tr ρ D ξ ξ ∈ R 2N (2.3) where in the last step we have introduced the displacement operator as The Fourier transform of the Wigner characteristic function provides the Wigner function w(r) ≡ 1 (2π) 2N χ(ξ) e −i r t J ξ dξ (2.5) where dξ = 2N i=1 dξ i denotes the integration over the 2N real components of ξ.
In this manuscript we focus on the Gaussian states of the harmonic lattices, which are the states whose Wigner function is Gaussian [42,43,46,[74][75][76][77] w G (r; γ, r ) ≡ e − 1 2 (r− r ) t γ −1 (r− r ) (2π) N det(γ) . (2.6) The 2N × 2N real, symmetric and positive definite matrix γ is the covariance matrix of the Gaussian state, whose elements can be defined in terms of the anticommutator of the operatorŝ r i as follows The covariance matrix γ is determined by N (2N + 1) real parameters. The expressions (2.6) and (2.7) tell us that the Gaussian states are completely characterised by the onepoint correlators (first moments) and by the two-points correlators (second moments) of the position and momentum operators collected into the vectorr. It is important to remark that the validity of the uncertainty principle imposes the following condition on the covariance matrix [43,70] γ + i 2 J 0 . (2.8) In [63] a real, positive matrix with an even size and satisfying (2.8) is called Gaussian matrix. Thus, every symmetric Gaussian matrix provides the covariance matrix of a Gaussian state.
A change of baser →r = Sr characterised by S ∈ Sp(2N, R) induces the transformation γ → γ = S γ S t on the covariance matrix.
In this manuscript we mainly consider Gaussian states with vanishing first moments, i.e. having r i = 0 (pure states that do not fulfil this condition are discussed in Sec. 2.7). In this case the generic element of covariance matrix (2.7) becomes γ i,j = 1 2 {r i ,r j } = Re r irj (2.9) and the Wigner function (2.6) slightly simplifies to w G (r; γ) = e − 1 2 r t γ −1 r (2π) N det(γ) . (2.10) where we have lightened the notation with respect to (2.6) by setting w G (r; γ) ≡ w G (r; γ, 0). The quantities introduced above characterise generic mixed Gaussian states. The subclass made by the pure states is discussed in Sec. 2.3.1.
The most familiar way to describe the Hilbert space is the Schrödinger representation, which employs the wave functions ψ(q) = q|ψ on R N (elements of L 2 (R N ) depending on q ≡ (q 1 , . . . , q N ) t ) for the vectors of the Hilbert space and the kernels O(q,q) = q| O |q for the linear operators O acting on the Hilbert space [69]. In the Appendix A.1 we relate the kernel ρ(q,q) = q|ρ |q of the density matrix to the corresponding Gaussian Wigner function (2.10). In the Appendix A.2 we express the kernel ρ A (q A ,q A ) for the reduced density matrix of a spatial subsystem A in terms of the parameters defining the wave function of the pure state describing the entire bipartite system.

Fisher-Rao distance
The set made by the probability density functions (PDF's) parameterised by the quantities γ is a manifold. In information geometry, the distinguishability between PDF's characterised by two different sets of parameters γ 1 and γ 2 is described through a scalar quantity D(γ 1 , γ 2 ) called divergence [57,58], a function such that D(γ 1 , γ 2 ) 0 and D(γ 1 , γ 2 ) = 0 if and only if γ 1 = γ 2 and D(γ, γ + dγ) = 1 2 i,j g ij dy i dy j + O (dy) 3 (2.11) where g ij is symmetric and positive definite and y denotes the vector collecting the independent parameters that determine γ = γ(y). In general D(γ 1 , γ 2 ) = D(γ 2 , γ 1 ); nonetheless, notice that the terms that could lead to the loss of this symmetry are subleading in the expansion (2.11). Thus, every divergence D introduces a metric tensor g ij that makes M a Riemannian manifold.
A natural requirement for a measure of distinguishability between states is the information monotonicity [57,58]. Let us denote by s = s(r) a change of variables in the PDF's and byD(γ 1 , γ 2 ) the result obtained from D(γ 1 , γ 2 ) after this change of variables. If s(r) is not invertible, a loss of information occurs because we cannot reconstruct r from s. This information loss leads to a less distinguishability between PDF's, namelyD(γ 1 , γ 2 ) < D(γ 1 , γ 2 ). Instead, when s(r) is invertible, information is not lost and the distinguishability of the two functions is preserved, i.e.D(γ 1 , γ 2 ) = D(γ 1 , γ 2 ). Thus, it is naturally to require that any change of variables must lead to [57,58] D(γ 1 , γ 2 ) D(γ 1 , γ 2 ) . (2.12) This property is called information monotonicity for the divergence D.
Let us consider a geometric structure on M induced by a metric tensor g ij associated to a divergence satisfying (2.12). An important theorem in information geometry due to Chentsov claims that, considering any set of the PDF's, a unique metric satisfying (2.12) exists up to multiplicative constants [56,57].
The Wigner functions of the bosonic Gaussian states (2.6) with vanishing first moments are PDF's that provide a manifold M G parameterised by the covariance matrices γ. The Chentsov's theorem for these PDF's leads to introduce the Fisher information matrix [56,57,64,66,78] g ij = w G (r, γ) ∂ log[w G (r; γ)] ∂y i ∂ log[w G (r; γ)] ∂y j dr (2.13) which provides the Fisher-Rao distance between two bosonic Gaussian states with vanishing first moments. Denoting by γ 1 and γ 2 the covariance matrices of these states, their Fisher-Rao distance reads [60-63, 65, 79] d(γ 1 , γ 2 ) ≡ Tr (log ∆) 2 ≡ log γ −1/2 1 (2.14) This is the main formula employed throughout this manuscript to study the complexity of Gaussian mixed states.
In Appendix B we report known results about the Fisher-Rao distance between Gaussian PDF's with non vanishing first moments [64,65,78,[80][81][82]. We remark that (2.14) is the Fisher-Rao distance also when the reference state and the target state have the same first moments, that can be non vanishing [65,79,83]. Although an explicit expression for the Fisher-Rao distance in the most general case of different covariance matrices and different first moments is not available in the literature, interesting classes of Gaussian PDF's have been identified where explicit expressions for this distance have been found [79,[83][84][85].
The distance between two states can be evaluated also through the distance between the corresponding density matrices. Various expressions for distances have been constructed and it is natural to ask whether they satisfy a property equivalent to the information monotonicity (2.12), that is known as contractivity [86][87][88]. A quantum operation Θ is realised by a completely positive operator which acts on the density matrixρ, providing another quantum state Θ(ρ) [43,86,87] (see also Sec. 7). A distance d between two states characterised by their density matricesρ 1 andρ 2 is contractive when the action of a quantum operation Θ reduces the distance between any two given states [86,88], namely 1 This is a crucial property imposed to a distance in quantum information theory.
The main contractive distances are the Bures distance, defined in terms of the fidelity F as follows the Hellinger distance d 2 H (ρ 1 ,ρ 2 ) = 2 1 − Tr ρ 1 ρ 2 (2.17) and the trace distance d L 1 (ρ 1 ,ρ 2 ) ≡ Tr ρ 1 −ρ 2 . (2.18) The trace distance is the L p -distance with p = 1 and it is the only contractive distance among the L p -distances. For p = 2 we have the Hilbert-Schmidt distance [87] d HS (ρ 1 ,ρ 2 ) = Tr ρ 1 −ρ 2 2 (2. 19) which is non contractive. In Appendix C we further discuss the Bures distance and the Hilbert-Schmidt distance specialised to the bosonic Gaussian states.

Williamson's decomposition
The Williamson's theorem is a very important tool to study Gaussian states [89]: it provides a decomposition for the covariance matrix γ that is crucial throughout our analysis.
The Williamson's theorem holds for any real, symmetric and positive matrix with even size; hence also for the covariance matrices. Given a covariance matrix γ, the Williamson's theorem guarantees that a symplectic matrix W ∈ Sp(2N, R) can be constructed such that where D ≡ diag(σ 1 , . . . , σ N ) ⊕ diag(σ 1 , . . . , σ N ) and σ k > 0. The set {σ k } is the symplectic spectrum of γ and its elements are the symplectic eigenvalues (we often call D the symplectic spectrum throughout this manuscript, with a slight abuse of notation). The symplectic spectrum is uniquely determined up to permutations of the symplectic eigenvalues and it is invariant under symplectic transformations. Throughout this manuscript we refer to (2.20) as the Williamson's decomposition 3 of γ, choosing a decreasing ordering for the symplectic eigenvalues. The real dimension of the set made by the covariance matrices is N (2N + 1) [45].
Combining (2.8) and (2.20), it can be shown that σ k 1 2 [43]. A diagonal matrix is symplectic when it has the form Υ ⊕ Υ −1 . This implies that a generic covariance matrix is not symplectic because of the occurrence of the diagonal matrix D in the Williamson's decomposition (2.20).
Another important tool for our analysis is the Euler decomposition of a symplectic matrix S (also known as Bloch-Messiah decomposition) [71]. It reads where Λ = diag(Λ 1 , . . . , Λ N ) with Λ j 0. The non-uniqueness of the decomposition (2.21) is due only to the freedom to order the elements along the diagonal of Λ. By employing the Euler decomposition (2.21) and that the real dimension of K(N ) is N 2 , it is straightforward to realise that the real dimension of the symplectic group Sp(2N, R) is 2N 2 + N , as already mentioned in Sec. 2.1. The simplest case corresponds to the one-mode case, i.e. N = 1, where a 2 × 2 real symplectic matrix can be parameterised by two rotation angles and a squeezing parameter Λ 1 .
The quantities explored in this manuscript provide important tools to study the entanglement quantifiers in harmonic lattices. For instance, the symplectic spectrum in (2.20) for the reduced density matrix allows to evaluate the entanglement spectrum and therefore the entanglement entropies [49, 74-76, 90, 91] and the Euler decomposition (2.21) applied to the symplectic matrix occurring in the Williamson's decomposition of the covariance matrix of a subsystem has been employed in [92] to construct a contour function for the entanglement entropies [93,94]. The Williamson's decomposition is also crucial to study the entanglement negativity [74,[95][96][97][98] a measure of the bipartite entanglement for mixed states.

Covariance matrix of a pure state
A Gaussian state is pure if and only if all the symplectic eigenvalues equal to 1 2 , i.e. D = 1 2 1. Thus, the Williamson's decomposition of the covariance matrix characterising a pure state reads The last expression, which has been found by employing the Euler decomposition (2.21) for the symplectic matrix W , tells us that the covariance matrix of a pure state can be determined by fixing N 2 + N real parameters.
The covariance matrix of a pure state satisfies the following constraint [99] After a change of basis characterised by the symplectic matrix S, the covariance matrix where K ∈ K(N ), the covariance matrix drastically simplifies to γ = 1 2 1. In the Schrödinger representation, the wave function of a pure Gaussian state reads [71] where E and F are N × N real symmetric matrices and E is also positive definite; hence the pure state is parameterised by N (N + 1) real coefficients, in agreement with the counting of the real parameters discussed above. The L 2 norm of (2.24) is equal to one.
The covariance matrix corresponding to the pure state (2.24) can be written in terms of the matrices E and F introduced in the wave function (2.24) as follows [71] γ = 1 2 where the symplectic matrix W and its inverse are given respectively by The expression (2.24) is employed in the Appendix A.2 to provide the kernel of a reduced density matrix in the Schrödinger representation.

Mixed states
Considering the set P(N ) made by the 2N × 2N real and positive definite matrices, the covariance matrices provide the proper subset of P(N ) made by those matrices that also satisfy the inequality (2.8).
The set P(N ) equipped with the Fisher-Rao distance is a Riemannian manifold where the length of a generic path γ : [a, b] → P(N ) is given by 4 [60][61][62][63]66] The unique geodesic connecting two matrices in the manifold P(N ) has been constructed [62]. In our analysis we restrict to the subset made by the covariance matrices γ. Considering the covariance matrix γ R and the covariance matrix γ T , that correspond to the reference state and to the target state respectively, the unique geodesic that connects γ R to γ T is [62] where s parameterises the generic matrix along the geodesic (we always assume 0 s 1 throughout this manuscript) and it is straightforward to verify that The geodesic (2.28) provides the optimal circuit connecting γ R to γ T . In the mathematical literature, the matrix (2.28) is also known as the s-geometric mean of γ R and γ T . The matrix associated to s = 1/2 provides the geometric mean of γ R and γ T . We remark that being γ R and γ T symmetric Gaussian matrices, it can be shown that also the matrices belonging to the geodesic (2.28) are symmetric and Gaussian [63].
By employing (D.1), we find that the geodesic (2.28) can be written in the following form The Fisher-Rao distance between γ R and γ T is the length of the geodesic (2.28) evaluated through (2.27). It is given by This distance provides the following definition of complexity It is straightforward to realise that, in the special case where both γ R and γ T correspond to pure states, the complexity (2.33) becomes the result obtained in [22] for the F 2 complexity, based on the F 2 cost function; hence we refer to (2.33) also as F 2 complexity in the following. The matching with [22] justifies the introduction of the numerical factor 1 2 √ 2 in (2.33) with respect to the distance (2.31). Equivalently, also the κ = 2 complexity given by C κ=2 ≡ C 2 2 can be considered.
We remark that the complexity (2.33) and the optimal circuit (2.28) can be applied also for circuits where the reference state and the target state have the same first moments [65,79,83].
Performing a change of basis characterised by the symplectic matrix S, the matrix ∆ TR changes as follows From this expression it is straightforward to observe that the Fisher-Rao distance (2.31), and therefore the complexity (2.33) as well, is invariant under a change of basis. We remark that (2.33) is invariant under any transformation that induces on ∆ TR the transformation (2.35) for any matrix S (even complex and not necessarily symplectic).
From the expression (2.28) of the geodesic connecting γ R to γ T , one can show that the change s → 1 − s provides the geodesic connecting γ T to γ R ; indeed, we have that 6 Another interesting result is the Fisher-Rao distance between the initial matrix γ R and the generic symmetric Gaussian matrix along the geodesic (2.28) reads [62] d γ R , G s (γ R , γ T ) = s d(γ R , γ T ) . (2.37) The derivation of some results reported in the forthcoming sections are based on the geodesic (2.28) written in the following form 7 This expression is interesting because the generic matrix of the optimal circuit is written in a form that reminde a symplectic transformation of γ R through the U s . Nonetheless, we remark that in general U s is not symplectic because the covariance matrices are not symplectic matrices. The steps performed to obtain (2.39) lead to write (2.36) as follows and then employing (D.1) in both the expressions within the square brackets of (2.38) with f (x) = x s/2 . D R < l a t e x i t s h a 1 _ b a s e 6 4 = " / R 4 E N c y 5 b m w G 0 D T H e 8 1 0 k 6 2 H B 3 I = " > A A A C C H i c b V C 7 T g J B F J 3 F F + I L t D I U T i Q m V m Q X C y 0 s S L S w B C O P h C V k d h h g w u z s Z u a u k W y 2 t P F X b C w 0 x t Y P s L D T r 3 F 4 F A q e Z D I n 5 9 y b e + / x Q s E 1 2 P a X l V p a X l l d S 6 9 n N j a 3 t n e y u d 2 6 D i J F W Y 0 G I l B N j 2 g m u G Q 1 4 C B Y M 1 S M + J 5 g D W 9 4 M f Y b t 0 x p H s g b G I W s 7 Z O + 5 D 1 O C R i p k z 3 A r k 9 g Q I m I L 5 N O 7 A K 7 A + W b n 8 s R v k 6 S T r Z g F + 0 J 8 C J x Z q R Q z r n 5 7 + r + R 6 W T / X S 7 A Y 1 8 J o E K o n X L s U N o x 0 Q B p 4 I l G T f S L C R 0 S P q s Z a g k P t P t e H J I g o + M 0 s W 9 Q J k n A U / U 3 x 0 x 8 b U e + Z 6 p H G + t 5 7 2 x + J / X i q B 3 1 o 6 5 D C N g k k 4 H 9 S K B I c D j V H C X K 0 Z B j A w h V H G z K 6 Y D o g g F k 1 3 G h O D M n 7 x I 6 q W i c 1 I s V U 0 a 5 2 i K N M q j Q 3 S M H H S K y u g K V V A N U X S P H t E z e r E e r C f r 1 X q b l q a s W c 8 e + g P r / Q d o i J 0 o < / l a t e x i t > D T < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 4 u B H Q S p p o b q 2 9 Z 0 O y k 5 G D S w H 5 H z z c z n G t S T p Z g t 2 0 Z 4 C L x N n T g r l n J v / r h 5 + V L r Z T 7 c X 0 M h n E q g g W r c d O 4 R O T B R w K l i S c S P N Q k J H Z M D a h k r i M 9 2 J p 4 c k + M Q o P d w P l H k S 8 F T 9 3 R E T X + u x 7 5 n K y d Z 6 0 Z u I / 3 n t C P o X n Z j L M A I m 6 W x Q P x I Y A j x J B f e 4 Y h T E 2 B B C F T e 7 Y j o k i l A w 2 W V M C M 7 i y c u k U S o 6 Z 8 V S 1 a R x i W Z I o z w 6 R q f I Q e e o j G 5 Q B d U R R f f o E T 2 j F + v B e r J e r b d Z a c q a 9 x y g P 7 D e f w B r l J 0 q < / l a t e x i t > T < l a t e x i t s h a 1 _ b a s e 6 4 = " Y a v w c I I S O 2 r g n X t 1 t p g n b 2 H M v 2 g = " > A A A C A n i c b V C 7 S g N B F J 3 1 G e N r V R D E Z j A I V m E 3 F l p Y B G w s I + Q F S Q i z k 0 k y Z G Z 2 m b m r L k u w 8 T + s b C w U s f U r 7 P w A / 8 P J o 9 D E A 5 d 7 O O d e Z u 4 J I s E N e N 6 X s 7 C 4 t L y y m l n L r m 9 s b m 2 7 O 7 t V E 8 a a s g o N R a j r A T F M c M U q w E G w e q Q Z k Y F g t W B w O f J r N 0 w b H q o y J B F r S d J T v M s p A S u 1 3 Y N m j 0 h J 2 m k T 2 B 1 o a T t X C S 4 P h 2 0 3 5 + W 9 M f A 8 8 a c k V 9 y v 3 k b J 9 2 O p 7 X 4 2 O y G N J V N A B T G m 4 X s R t F K i g V P B h t l m b F h E 6 I D 0 W M N S R S Q z r X R 8 w h A f W 6 W D u 6 G 2 p Q C P 1 d 8 b K Z H G J D K w k 5 J A 3 8 x 6 I / E / r x F D 9 7 y V c h X F w B S d P N S N B Y Y Q j / L A H a 4 Z B Z F Y Q q j m 9 q + Y 9 o k m F G x q W R u C P 3 v y P K k W 8 v 5 p v n B t 0 7 h A E 2 T Q I T p C J 8 h H Z 6 i I r l A J V R B F 9 + g J v a B X 5 8 F 5 d t 6 c 9 8 n o g j P d 2 U N / 4 H z 8 A O 8 f m 3 E = < / l a t e x i t > R < l a t e x i t s h a 1 _ b a s e 6 4 = " c 1 0 P n Y + 5 / m + T 2 R H T r w w X A K 6 s B r 4 = " o Q s G m l r U h + N M n z 5 J K I e 8 f 5 w u X N o 0 z N E Y G 7 a M D d I R 8 d I K K 6 A K V U B l R d I + e 0 A t 6 d R 6 c Z + f N e R + P z j m T n R 3 0 B 8 7 H D + w T m 2 8 = < / l a t e x i t > Figure 1: Pictorial representation of the optimal circuit (2.28) connecting γ R to γ T (solid black curve). Coloured solid curves represent the sets made by symmetric Gaussian matrices having the same symplectic spectrum. The red curve corresponds to D R and the blue curve to D T .
It is enlightening to exploit the Williamson's decomposition of the covariance matrices discussed in Sec. 2.3 in the expressions for the complexity and for the optimal circuit. The Williamson's decomposition (2.20) allows to write γ R and γ T as follows where D R and D T contain the symplectic spectra of γ R and γ T respectively. Let us introduce also the Williamson's decomposition of the generic matrix along the geodesic (2.28), namely It would be insightful to find analytic expressions for W s and D s in terms of γ R and γ T . This has been done later in the manuscript for some particular optimal circuits.
In Fig. 1 we show a pictorial representation of the optimal circuit (2.28), which corresponds to the solid black curve. The figure displays that the symplectic spectrum changes along the geodesic because the black curve crosses solid curves having different colours, which correspond to the sets of matrices having the same symplectic spectrum.
In order to write the complexity (2.33) in a convenient form depending on the symplectic spectra and on the symplectic matrices W R and W T , let us employ that, after a canonical transformation characterised by the symplectic matrix S, the covariance matrices in (2.41) become where K D R is symplectic and such that K D R D R K t D R = D R (the set of matrices made by K D R is a subgroup of Sp(2N, R) called stabilizer [22]), we have that (2.43) become respectively where we have introduced the symplectic matrix W TR defined as follows For later convenience, let us consider the Euler decomposition (defined in Sec. 2.3) of the symplectic matrix W TR , namely being Λ TR a diagonal matrix with positive entries. By specifying (2.35) to (2.44), we find that which allows to write the F 2 complexity (2.33) as 8 This expression is independent of K D R and tells us that, in order to evaluate the F 2 complexity (2.33) we need the symplectic spectra D R and D T and the symplectic matrix (2.45).
By employing the Euler decomposition (2.46), the second covariance matrix in (2.44) can be decomposed as follows which is simpler to compute than (2.28) because γ R is diagonal. Let us remark that G s (γ R , γ T ) is different from G s (γ R , γ T ) but they have the same length given by (2.49). Furthermore, while the optimal circuit (2.51) depends on the matrix K D R , its length (2.49) does not.
For pure states, both (2.49) and (2.51) simplify in a significant way, as discussed in Sec. 2.5.

One-mode mixed states
For mixed states defined by a single mode (i.e. N = 1), the results discussed above significantly simplify because the diagonal matrices D R and D T are proportional to the 2×2 identity matrix; hence the covariance matrices of the reference state and of the target state become respectively where σ R 1/2 and σ T 1/2.
In this case the Williamson's decomposition for the optimal circuit (2.28) can be explicitly written. Indeed, from (2.52) one finds that and this leads to write the expression (2.39) for the optimal circuit as follows which provide the Williamson's decomposition of the generic matrix along the optimal circuit.
By specialising the complexity (2.49) to the one-mode mixed states in (2.52) we get Thus, the formal expression for the complexity does not simplify significantly for the one-mode mixed states with respect to the general case with N 1.

Pure states
It is very insightful to specialise the results presented in Sec. 2.4 to pure states.
When both the reference state |ψ R and the target state |ψ T are pure states, the corresponding density matrices are the projectorsρ R = |ψ R ψ R | andρ T = |ψ T ψ T | respectively. In this case the symplectic spectra drastically simplify to where 1 is the 2N × 2N identity matrix. This implies that the Williamson's decompositions in (2.41) become respectively The complexity of pure states can be easily found by specialising (2.49) to (2.56). The resulting expression can be further simplified by employing (2.46), (2.47) and the cyclic property of the trace. This gives the result obtained in [22] which can be also obtained through the proper choice of the base described below.
Since we are considering pure states, (2.26) can be employed to write W R and W T in terms of the pairs of symmetric matrices (E R , F R ) and (E T , F T ) occurring in the wave functions (2.24) of the reference state and of the target state respectively. The matrix W TR in (2.45), that provides the complexity (2.58), can be written as follows which becomes block diagonal for real wave functions (i.e. when F R = F T = 0).
As for the optimal circuit (2.28), by specialising the form (2.30) to the covariance matrices of pure states in (2.57), we obtain We find it instructive also to specialise the expression (2.39) for the optimal circuit to pure states. Indeed, in this case ∆ TR is symplectic and the result reads This expression provides the Williamson's decomposition of the optimal circuit made by pure states, being W s ∈ Sp(2N, R).
A proper choice of the basis leads to a simple expression for the optimal circuit. Since D R = 1 2 1, we have that K D R introduced in the text below (2.43) is an orthogonal matrix. For pure states the convenient choice is K D R = R TR . Indeed, by specifying (2.44) to this case we obtain that in this basis the covariance matrices γ R and γ T become the following diagonal matrices We remark that this result has been obtained by exploiting the peculiarity of the pure states mentioned in Sec. 2.3, namely that, after a change of basis that brings the covariance matrix into the diagonal form 1 2 1, another change of basis characterised by a symplectic matrix that is also orthogonal leaves the covariance matrix invariant. The occurrence of non trivial symplectic spectra considerably complicates this analysis (see (2.49) and (2.51)).
Specialising the form (2.30) of the optimal circuit to the covariance matrices in (2.62), the following simple expression is obtained [22] This expression tells us that, for pure states, this basis is very convenient because the optimal circuit is determined by the diagonal matrix X TR .

Thermal states
The thermal states provide an important class of Gaussian mixed states. The density matrix of a thermal state at temperature T ≡ 1/β isρ th = e −βĤ /Z, where H is the hamiltonian (2.1) for the harmonic lattices that we are considering and the constant Z = Tr e −βĤ guarantees the normalisation condition Trρ th = 1.
In order to study the Williamson's decomposition of the covariance matrix associated to a thermal state, let us observe the matrix H phys in (2.1) can be written as where P phys = 1 m 1 and Q phys is a N × N real, symmetric and positive definite matrix whose explicit expression is not important for the subsequent discussion.
Denoting by V the real orthogonal matrix that diagonalises Q phys (for the special case of the harmonic chain with periodic boundary conditions, V has been written in (9.6) e (9.7)), it is straightforward to notice that (2.64) can be diagonalised as follows where Ω 2 k are the real eigenvalues of Q phys /m. It is worth remarking that the 2N × 2N matrix V is symplectic and orthogonal, being V orthogonal. By employing the argument that leads to (D.10), the r.h.s. of (2.65) can be written as where we have introduced the following symplectic and diagonal matrix X phys = diag (mΩ 1 ) 1/2 , . . . , (mΩ N ) 1/2 , (mΩ 1 ) −1/2 , . . . , (mΩ N ) −1/2 .
Following the standard quantisation procedure, one introduces the annihilation operatorŝ b k and the creation operatorsb † k aŝ which satisfy the well known algebra given by [b i ,b j ] = J ij . In terms of these operators, the hamiltonian (2.70) assumes the standard form Thus, the symplectic spectrum in (2.69) provides the dispersion relation of the model.
The operator (2.72) leads us to introduce the eigenstates |n k of the occupation number operatorb † kb k , whose eigenvalues are given by non negative integers n k , and the states |n ≡ N k=1 |n k . The expectation value of an operator O on the thermal state reads Considering the covariance matrix of a Gaussian state defined in (2.9), by employing (2.70), where W phys is a real matrix, one finds that the covariance matrix of the thermal state can be written as in terms of the covariance matrix in the canonical variables collected intoŝ, whose elements are given by the correlators q kqk , p kpk and q kpk . These correlators can be evaluated by first using (2.71) to write Re ŝŝ t = Re Θ bb t Θ t , where we remark that Θ is not a symplectic matrix because it does not preserve the canonical commutation relations. Then, by exploiting (2.73) and the action ofb k andb † k onto the Fock states, one computes bb t . This leads to a diagonal matrix Re ŝŝ t whose non vanishing elements are given by [43,45] Re q kqk = Re p kpk = 1 2 coth(βΩ k /2) Re q kpk = 0 .
where the symplectic eigenvalues entering in the diagonal matrix D th and the symplectic matrix W th are given respectively by We remark that W th is independent of the temperature.
Taking the zero temperature limit β → +∞ of (2.76), one obtains the Williamson's decomposition of the covariance matrix of the ground state. This limit gives σ th,k → 1/2, as expected from the fact that the ground state is a pure state, while W th does not change, being independent of the temperature. Thus, the Williamson's decomposition of the covariance matrix of the ground state reads where W phys has been defined in (2.69).
It is worth considering the complexity when the reference state and the target state are thermal states having the same physical hamiltonian H but different temperatures (we denote respectively by β R and β T their inverse temperatures). From (2.76), we have that the Williamson's decomposition of the covariance matrices of the reference state and the target state read respectively where W th is independent of the temperature; hence W R = W T . This means that W TR = 1 in this case (see (2.45)); hence the expression (2.49) for the complexity significantly simplifies to (2.80) The optimal path connecting these particular thermal states is obtained by plugging (2.79) into (2.39). Furthermore, by exploiting (D.1) and some straightforward matrix manipulations, we find that the Williamson's decomposition of the generic covariance matrix belonging to this optimal path reads where the same symplectic matrix W th of the reference state and of the target state occurs and only the symplectic spectrum depends on the parameter s labelling the covariance matrices along the optimal path.
It is worth asking whether, for any given value of s, the covariance matrix G s (γ th,R , γ th,T ) in (2.81) can be associated to a thermal state of the system characterised by the same physical hamiltonian underlying the reference and the target states. Denoting by σ s,k the symplectic eigenvalues of (2.81) this means to find a temperature T s ≡ β −1 s such that σ s,k = 1 2 coth(β s Ω k /2). This equation can be written more explicitly as follows We checked numerically that a solution T s = T s (T R , T T ) for any 1 k N does not exist.
The quantities discussed above are further explored in Sec. 9.3, where the thermal states of the harmonic chain are considered.

Coherent states
The coherent states are pure states with non vanishing first moments [45]. They can be introduced through the displacement operator defined in (2.4), where the real vector a ∈ R 2N can be parameterised in terms of the complex vector α ∈ C N as a t = √ 2 (Re(α) t , Im(α) t ).
The displacement operator (2.4), which is unitary and satisfies D −1 a = D −a , shifts the position and the momentum operators as follows The distance (2.31), that is mainly used throughout this manuscript to study the circuit complexity of mixed states, is valid for states having the same first moments [62,65,79], as reported in the Appendix B.
In the Appendix B it is also mentioned that an explicit expression for the complexity between coherent states is available in the literature if we restrict to the coherent states having a diagonal covariance matrix and a = ( √ 2 α, 0, . . . , 0) [79,84,85]. These states provide a manifold parametrised by (α, v 1 , . . . , v 2N ), being v 2 k the k-th entry of the diagonal covariance matrix, and whose metric is given by (B.10) with n = 2N and µ = √ 2α, namely Let us remind that the covariance matrices that we are considering must satisfy the constraint (2.8), which is equivalent σ k 1/2 for the symplectic eigenvalues, being k = 1, . . . , N . By using (D.10), one finds that the symplectic eigenvalues of the diagonal covariance matrix diag(v 2 1 , . . . , v 2 2N ) are σ k = v k v k+N for k = 1, . . . , N . Thus, in our case the manifold equipped with the metric (2.88) must be constrained by the conditions v k v k+N 1/2 for k = 1, . . . , N .
The coherent states are pure states, hence their covariance matrices must satisfy the condition (2.23), which holds also when the first moments are non vanishing [45,99]. For the class of coherent states that we are considering, the constraint (2.23) leads to which saturate the constraints v k v k+N 1/2 introduced above. By imposing the conditions (2.89), the metric (2.88) becomes which is twice (B.10) with n = N and µ = α. The geometry given by (2.90) has been found also in [18]. The constraint (2.89) tells us that the metric (2.90) is defined on a set of pure states, but we are not guaranteed that the resulting manifold is totally geodesic. This is further discussed in the final part of this subsection.
Given a reference state and a target state parametrised 9 by φ R = (α R , v R,1 , . . . , v R,N ) and φ T = (α T , v T,1 , . . . , v T,N ) respectively, the square of the complexity of the circuit corresponding to the geodesic connecting these states in the manifold equipped with (2.90) is easily obtained from (B.11). The result reads (2.91) By adopting the normalisation in (2.33), which is consistent with [17,18], one can introduce the complexity between coherent states as follows Setting α R = 0 (or α T = 0, equivalently) in (2.92), one obtains the complexity between a coherent state in the particular set introduced above and the ground state. As consistency check, we observe that, by setting α R = α T in (2.92), the complexity (2.58) between pure states is recovered.
It is instructive to compare (2.92) with the results reported in [18], where the complexity C κ=2 = C 2 2 between the ground state and a bosonic coherent state has been studied through the Nielsen's approach [1][2][3]. The analytic expression for the complexity in [18] has been obtained for reference and target states with diagonal covariance matrices and first moments with at most one non vanishing component. Since these are the assumptions under which (2.92) has been obtained, we can compare the two final results for the complexity. The analysis of [18] allows to write the complexity C κ=2 in terms of a free parameter x 0 which does not occur neither in the reference state nor in the target state. We observe that the square of (2.92) with α R = 0 coincides with the result in [18] 10 with x 0 = 2v R,1 .
In the following we consider circuits in the space of the Gaussian states with non vanishing first moments such that the reference and the target states are given by two coherent states (2.84) originated from the same ground state, denoting their first moments by a R and a T respectively. These states have the same covariance matrix γ 0 (see (2.86)), whose symplectic 9 The vector φ corresponds to the vector θ used in Appendix B restricted by the condition (2.89). 10 In Eqs. (4.23) and (4.24) of [18], set i = 1, √ 2 ω R,there = 1/v R,k , √ 2 m there ω k,there = 1/v T,k and a i,there = a there /x0 = αT/( √ 2 vR,1).
eigenvalues are equal to 1/2, being the coherent states pure states. Parametrising the reference state by θ R = (a R , γ 0 ) and the target state by θ T = (a T , γ 0 ), a recent result obtained in [83] and discussed in Appendix B allows us to write the circuit complexity as follows where d FR has been defined in (B.13). We are not able to prove that σ k 1/2 for the symplectic eigenvalues of the symmetric and positive matrices making the geodesic whose length is (2.93).
It is worth remarking that the expressions (2.92) and (2.93) for the complexity are defined for different sets of Gaussian Wigner functions with a non vanishing intersection. Indeed, (2.93) holds between PDF's with the same covariance matrix (that can be also non diagonal), while (2.92) is valid for diagonal covariance matrices that can be different. Moreover in (2.92), a R and a T can have only one non vanishing components, while in (2.93) they are generic. Thus, in order to compare (2.92) and (2.93) we have to consider reference and target states which have the same diagonal covariance matrix and and whose first moments have only one non vanishing component. Setting v R,k = v T,k ≡ v k with k = 1, . . . , N in (2.92), we obtain (2.94) and a S = ( √ 2 α S , 0, . . . , 0) for S = R and S = T in (2.93), one finds (2.95) A simple numerical inspection shows that (2.95) is always smaller than (2.94). This example allows to conclude that the submanifold of pure states with diagonal covariance matrix equipped with the metric (2.90) is not totally geodesics.

Spectrum complexity and basis complexity
In this section we discuss the spectrum complexity and the basis complexity for mixed Gaussian states in harmonic lattices.
By exploiting the Williamson's decomposition we introduce the W path as the optimal circuit connecting two covariance matrices with W R = W T ≡ W and the D path as the optimal circuit connecting two covariance matrices having D R = D T ≡ D. In order to study these circuits, in Sec. 3.1 we discuss the first law of complexity for the Gaussian states that we are considering. The lengths of a W path and of a D path are employed to study the spectrum complexity (Sec. 3.3) and the basis complexity (Sec. 3.4) respectively. In Fig. 2 the dashed curves correspond to W paths (see (3.18)). H z z c z n G t S T p Z g t 2 0 Z 4 C L x N n T g r l n J v / r h 5 + V L r Z T 7 c X 0 M h n E q g g W r c d O 4 R O T B R w K l i S c S P N Q k J H Z M D a h k r i M 9 2 J p 4 c k + M Q o P d w P l H k S 8 F T 9 3 R E T X + u x 7 5 n K y d Z 6 0 Z u I / 3 n t C P o X n Z j L M A I m 6 W x Q P x I Y A j x J B f e 4 Y h T E 2 B B C F T e 7 Y j o k i l A w 2 W V M C M 7 i y c u k U S o 6 Z 8 V S 1 a R x i W Z I o z w 6 R q f I Q e e o j G 5 Q B d U R R f f o E T 2 j F + v B e r J e r b d Z a c q a 9 x y g P 7 D e f w B r l J 0 q < / l a t e x i t > R < l a t e x i t s h a 1 _ b a s e 6 4 = " Z D R d 7 I c w J o u 5 3 I h v 3 R H 6 8 9 l Q J a E = " > A A A C C n i c b V A 9 S w N B E N 2 L 3 / E r f n Q 2 q y J Y h b t Y K F a C j a W K U S E X w t 5 m E p f s 7 h 2 7 c 2 I 4 r r b x r 9 h Y K C J Y W d j Y 2 P l v 3 C Q W m v h g m M d 7 M + z O i x I p L P r + l 1 c Y G 5 + Y n J q e K c 7 O z S 8 s l p a W z 2 2 c G g 5 V H s v Y X E b M g h Q a q i h Q w m V i g K l I w k X U O e z 5 F 9 d g r I j 1 G X Y T q C v W 1 q I l O E M n N U r r I Q r Z h C x s M 6 V Y 3 s h C h B s 0 y n W h u / Q 0 z x u l T b / s 9 0 F H S f B D N g 9 W X 8 e S j / f 0 u F H 6 D J s x T x V o 5 J J Z W w v 8 B O s Z M y i 4 h L w Y p h Y S x j u s D T V H N V N g 6 1 n / l J x u O a V J W 7 F x p Z H 2 1 d 8 b G V P W d l X k J h X D K z v s 9 c T / v F q K r b 1 6 J n S S I m g + e K i V S o o x 7 e V C m 8 I A R 9 l 1 h H E j 3 F 8 p v 2 K G c X T p F V 0 I w f D J o + S 8 U g 5 2 y p U T l 8 Y + G W C a r J E N s k 0 C s k s O y B E 5 J l X C y S 2 5 J 4 / k y b v z H r x n 7 2 U w W v B + d l b I H 3 h v 3 5 t 9 n y I = < / l a t e x i t > T < l a t e x i t s h a 1 _ b a s e 6 4 = " 1 v I K C M k d W k S K X p y H c + p f I D R f W 2 g = " > A A A C C n i c b V A 9 S w N B E N 1 L / I j x K 3 5 0 N q s i W I U 7 L R Q r w c Z S w a i Q C 2 F v M 0 m W 7 O 4 d u 3 N i O K 6 2 8 a / Y W C g i W F n Y 2 N j 5 b 9 w k F m p 8 M M z j v R l 2 5 0 W J F B Z 9 / 9 M r F C c m p 6 Z L M + X Z u f m F x c r S 8 r m N U 8 O h x m M Z m 8 u I W Z B C Q w 0 F S r h M D D A V S b i I e k c D / + I K j B W x P s N + A g 3 F O l q 0 B W f o p G Z l P U Q h W 5 C F H a Y U y 5 t Z i H C N R r k u d J + e 5 X m z s u l X / S H o O A m + y e b h 6 k s x e X 9 L T 5 q V j 7 A V 8 1 S B R i 6 Z t f X A T 7 C R M Y O C S 8 j L Y W o h Y b z H O l B 3 V D M F t p E N T 8 n p l l N a t B 0 b V x r p U P 2 5 k T F l b V 9 F b l I x 7 N q / 3 k D 8 z 6 u n 2 N 5 v Z E I n K Y L m o 4 f a q a Q Y 0 0 E u t C U M c J R 9 R x g 3 w v 2 V 8 i 4 z j K N L r + x C C P 6 e P E 7 O d 6 r B b n X n 1 K V x Q E Y o k T W y Q b Z J Q P b I I T k m J 6 R G O L k h d + S B P H q 3 3 r 3 3 5 D 2 P R g v e 9 8 4 K + Q X v 9 Q u e i Z 8 k < / l a t e x i t > T < l a t e x i t s h a 1 _ b a s e 6 4 = " Y a v w c I I S O 2 r g n X t 1 t p g n b 2 H M v 2 g = " > A A A C A n i c b V C 7 S g N B F J 3 1 G e N r V R D E Z j A I V m E 3 F l p Y B G w s I + Q F S Q i z k 0 k y Z G Z 2 m b m r L k u w 8 T + s b C w U s f U r 7 P w A / 8 P J o 9 D E A 5 d 7 O O d e Z u 4 J I s E N e N 6 X s 7 C 4 t L y y m l n L r m 9 s b m 2 7 O 7 t V E 8 a a s g o N R a j r A T F M c M U q w E G w e q Q Z k Y F g t W B w O f J r N 0 w b H q o y J B F r S d J T v M s p A S u 1 3 Y N m j 0 h J 2 m k T 2 B 1 o a T t X C S 4 P h 2 0 3 5 + W 9 M f A 8 8 a c k V 9 y v 3 k b J 9 2 O p 7 X 4 2 O y G N J V N A B T G m 4 X s R t F K i g V P B h t l m b F h E 6 I D 0 W M N S R S Q z r X R 8 w h A f W 6 W D u 6 G 2 p Q C P 1 d 8 b K Z H G J D K w k 5 J A 3 8 x 6 I / E / r x F D 9 7 y V c h X F w B S d P N S N B Y Y Q j / L A H a 4 Z B Z F Y Q q j m 9 q + Y 9 o k m F G x q W R u C P 3 v y P K k W 8 v 5 p v n B t 0 7 h A E 2 T Q I T p C J 8 h H Z 6 i I r l A J V R B F 9 + g J v a B X  Eq. (3.18) tells us that the dashed black curves represent the W R path and the W T path that pass through γ R and γ T respectively (the auxiliary covariance matricesγ R andγ T have been defined in (5.1)). The arcs of the dashed curves that connect the blue curve to the red curve have the same length given by (3.19).

First law of complexity
It is worth investigating the first law of complexity [67,68] for the states described in Sec. 2. The derivations of the results reported below are given in the Appendix E.
Let us consider the following functional where q 0 = q(t 0 ) and q 1 = q(t 1 ) are the initial and final configurations respectively.
It is well known that the first variation of (3.1) under an infinitesimal change of the boundary conditions q i → q i + δq i for i = 0, 1 evaluated on a solution of the equations of motion is The functional we are interested in is the length functional (2.27) and the solution of its equations of motion is given by the optimal circuit (2.28), that satisfies the boundary conditions (2.29). In order to apply (3.2), one considers the infinitesimal variations γ T → γ T + δγ T and γ R → γ R + δγ R of the covariance matrices of the reference and of the target states that preserve the properties of these matrices. In other words, these variations are such that also the resulting matrices are covariance matrices.
The length functional (2.27) leads to introduce the following cost function By applying (3.2) to the length functional (2.27), one obtains the first law of complexity where the r.h.s. is evaluated on the geodesic (2.28).
Equivalent expressions for the variation (3.4) have been derived in Appendix E. For instance, it can be written as where G s is the geodesic (2.28). Another useful expression for (3.4), which is simpler to evaluate than (3.5), reads where ∆ TR has been defined in (2.32).
We find it worth providing also an expression for the variation (3.4) that is based on the Williamson's decompositions (2.41). It is given by where, by using (D.1), one can write the matrices within the square brackets in terms of the matrix W TR introduced in (2.45) as follows The form (3.7) for δd tells us that this variation can be written as the sum of four contributions: two terms from the variations δW R and δW T of the symplectic matrices in (2.41) and two terms from the diagonal and non negative variations δD R and δD T of the symplectic spectra.

Solving δd = 0
It is natural to look for relations between γ R and γ T that lead to δd = 0 and, in order to find them, let us consider the first law of complexity written in the form (3.7).
First we focus on the variations of W R and W T . When δD R = δD T = 0 in (3.7), the equation A trivial solution of this equation is given by Another solution of (3.10) is W R = M W T , where M is a constant symplectic matrix whose elements are just real numbers, i.e. it does not contain parameters to vary. Notice that these two simple solutions require that both W R and W T are allowed to vary.
In order to find solutions to (3.10) for any choice of the independent variations δW R and δW T , let us first observe that, by taking the first variation of the relation W t J W = J, that characterises a symplectic matrix W t , it is not difficult to realise that δX ≡ W t J δW must be a symmetric matrix. By using that W −1 = J −1 W t J, for the two terms in the r.h.s. of (3.10) one obtains Tr Since δX corresponds to a generic infinitesimal real and symmetric matrix, Tr(Y δX) = 0 is satisfied for every δX e.g. when Y is a real antisymmetric matrix. These observations lead us to write (3.10) as where δX R ≡ W t R J δW R and δX T ≡ W t T J δW T are real and symmetric matrices. Thus, from (3.12), we have that (3.10) is satisfied for where Y is a real antisymmetric matrix that can depend on γ R or γ T . We remark that (3.13) solves (3.10) for any choice of the independent variations δW R and δW T .
It is worth asking when the case W R = M W T mentioned above, where M is a symplectic matrix made of real numbers, becomes a special case of (3.13). The Williamson's decompositions (2.41) and . By inserting the identity 1 = J −1 J between all the factors within the argument of the logarithm occurring in this matrix, using (D.1) again and exploiting the properties of the symplectic matrices, one arrives to log(γ −1 Comparing this expression with the one reported above, we conclude that the matrix log(γ −1 . This is the case for a symplectic matrix that is also orthogonal, i.e. M ∈ K(N ). In particular, the special case given by M = 1 ∈ K(N ) corresponds to (3.11). Summarising, in the special case given by As for the terms corresponding to the variations of DR and DT in (3.7), let us observe that, for a diagonal matrix Λ > 0, we have that Λ −1 δΛ = δ log Λ and that Tr(H δ log Λ) = 0 holds for a generic δΛ when all the elements along the diagonal of H vanish. The matrices having vanishing elements on their main diagonal are called hollow matrices. By employing these observations in the equation δd = 0 with δd given by (3.7), where the variations of DR and DT are independent, we conclude that the main diagonals of the matrices within the square brackets in (3.7) must vanish. By introducing two non vanishing hollow matrices Z and Z, this gives log γ −1 which correspond to the terms containing δDT and δDR respectively in (3.7). By employing (3.8) and (3.9), one finds that the relations in (3.14) can be written respectively as These observations tell us that δd = 0 for generic variations of γ R and γ T when these covariance matrices are related by (3.13) with Y constrained by the conditions that the elements on the diagonals of A rough analysis shows that this problem is too constrained for N = 1 and N = 2. Indeed, the antisymmetric matrix Y depends on N (2N − 1) parameters and imposing that the diagonals of Z and Z vanish provides 4N constraints. In particular, when N = 1 the 2 × 2 antisymmetric matrix Y has only one non vanishing off diagonal element a in the top right position and it is straightforward to check that Y J = −a 1. By using also (3.13) and (3.14) specialised to this case, we obtain that the above procedure leads to impose that W S Y J W −1 S , with S = {R, T}, must have vanishing elements along the diagonal. This is possible only for a = 0, i.e. Y = 0. Thus, when N = 1 we cannot find a solution of the form (3.13) for the equation δd = 0 with δd given by (3.7).

Spectrum complexity
It is worth exploring the possibility to define the circuit complexity associated to the change of the symplectic spectrum.
Let us consider a reference state and a target state such that in the Williamson's decompositions of their covariance matrices γ R and γ T (see (2.41)) the same symplectic matrix occurs, namely We call W path the optimal circuit (2.28) connecting these two covariance matrices.
In order to study the Williamson's decomposition of a matrix belonging to a W path, we consider the expression (2.39) for the optimal circuit. When (3.16) holds, from (2.32) it is straightforward to find that ∆ TR = W t D T D −1 R W −t . Then, by employing (D.1) both in U s and in U t s occurring in (2.39), we obtain which tells us that the Williamson's decomposition of the matrix along the W path is (2.42) with It is remarkable that the symplectic matrix W s is independent of s. This means that in the Williamson's decomposition of a matrix belonging to a W path the same symplectic matrix W occurs. In Fig. 2 the dashed curves correspond to the W R path and to the W T path. Considering e.g. the W R path in Fig. 2, from (3.18) we have that the Williamson's decomposition of a generic matrix γ belonging to this W R path is given by the symplectic matrix W R and by the symplectic spectrum corresponding to the coloured line intersecting the dashed line of the W R path at γ.
An interesting example of W path is given by the thermal states of a given model at different temperatures (see Sec. 2.6). Indeed, in the Williamson's decomposition (2.76), the symplectic matrix W th is independent of the temperature.
For a W path we have δd = 0 (see (3.11)); hence the W paths provide a preferred way to connect the set of covariance matrices with symplectic spectrum D R to the set of covariance matrices with symplectic spectrum D T .
We find it natural to define the spectrum complexity as the length of a W path because this quantity is independent of the choice of W . In particular, from (3.16), we have that W TR = 1, hence (2.49) simplifies to Fig. 2 the arcs of the dashed curves that connect the blue curve to the red curve have the same length given by (3.19).
Another natural definition for the spectrum complexity is the distance between the set of covariance matrices whose symplectic spectrum is D R (red curve in Fig. 2) and the set of covariance matrices whose symplectic spectrum is D T (blue curve in Fig. 2). It reads where the minimisation over the symplectic matrices W and W is difficult to perform. It is straightforward to realise thatd spectrum (D R , D T ) d spectrum (D R , D T ).
In the simplest case of one-mode mixed states (i.e. when N = 1), the optimal circuit (3.17) simplifies to which tells us that the 2 × 2 matrix belonging to the W path is a proper rescaling of the covariance matrix of the reference state.

Basis complexity
In order to study the circuit complexity associated to a change of basis, let us consider the Williamson's decompositions of two covariance matrices γ R and γ T having the same symplectic spectrum, i.e.
that have been obtained by setting D R = D T = D in (2.41). An important example is given by states whose density matricesρ T andρ R are related through a unitary transformation U , namelyρ T = Uρ R U † . Indeed, this means that the corresponding covariance matrices are related through a symplectic matrix (that does not change the symplectic spectrum) [43,45].
We denote as D path the optimal circuit connecting the covariance matrices having the same symplectic spectrum, identifying its length as a basis complexity. This basis complexity can be found by specifying (2.49) to (3.22) and the result is 11 where W TR has been defined in (2.45). Notice that we have not required that all the matrices along a D path have the same symplectic spectrum.
We find it reasonable to introduce also another definition of basis complexity as the minimal length of an optimal circuit that connects a covariance matrix whose Williamson's decomposition contains the symplectic matrix W R (i.e. that lies on the dashed curve on the left in Fig. 2) to a covariance matrix having the symplectic matrix W T in its Williamson's decomposition (i.e. that belongs to the dashed curve on the right in Fig. 2). This basis complexity is defined as follows where the minimisation is performed over the set Diag(2N, R) made by the diagonal matrices of the form diag(σ) ⊕ diag(σ), being σ a vector of N real numbers Specifying the form (2.39) for the optimal circuit to (3.22), it is straightforward to find that the D path is given by where we remark that W s is not symplectic in general.
It is worth asking when W s is symplectic because in these cases (3.25) provides the Williamson's decomposition of the D path. The requirement U s ∈ Sp(2N, R) leads to When this condition holds, (3.23) simplifies to the following expression which is independent of D.
For pure states, which have D = 1 2 1, the condition (3.26) is trivially verified. Another interesting example where (3.26) holds is given by the one-mode states, where D is proportional to the 2 × 2 identity matrix. In this case we can always connect two covariance matrices with the same symplectic spectrum through the optimal circuit (2.53), that can be written as Writing W TR as a block matrix made by four N × N matrices, it is straightforward to find that the condition (3.26) holds whenever every block of W TR commutes with diag(σ 1 , . . . , σ N ). Then, we can exploit the fact that a diagonal matrix with distinct elements commutes with another matrix only when the latter one is diagonal 12 . Thus, if the symplectic spectrum is non degenerate, all the blocks of W TR must be diagonal to fulfil the condition (3.26). We remark that the non-degeneracy condition for the symplectic spectrum is not guaranteed; indeed, the symplectic spectrum has some degeneracy in several interesting cases. For instance, for pure states all the symplectic eigenvalues are equal to 1 2 . Another important example is the reduced covariance matrix of an interval in an infinite harmonic chain with non vanishing mass [90].
We find it worth discussing the relation between the optimal circuits considered above to study the basis complexity and the solutions of the equation δd = 0 described in Sec. 3.2. For the set of paths occurring in (3.24), which includes the D paths, we have δW R = δW T = 0 in (3.7). In this case, in Sec. 3.2 we found that a solution of δd = 0 is given by (3.15), where Z and Z are non vanishing hollow matrices. Restricting to the cases of D paths that satisfy also (3.26), these relations simplify respectively to W TR W t TR = e Z and W t TR W TR = e Z , whose solution is non trivial because a matrix does not commute with its transpose in general (the matrices satisfying this property are called normal matrices). Notice that W TR ∈ K(N ) are not admissible solutions because Z and Z are non vanishing.
A slightly more general solution can be obtained by restricting to the D paths (see (3.22)). In this case, from (3.24) with δW R = δW T = 0 and δD R = δD T ≡ δD, we have that δd = 0 becomes By using (3.8), (3.9) and the discussion made in Sec. 3.2, one finds that (3.29) is solved when is a non vanishing hollow matrix. When (3.26) also holds, this condition simplifies to the requirement that log(W TR W t TR ) − log(W t TR W TR ) is a non vanishing hollow matrix, which is independent of D.

Purification through the W path
The purification of a mixed state is a process that provides a pure state starting from a mixed state. This procedure is not unique. Considering the context of the bosonic Gaussian states that we are exploring, in this section we discuss the purification of a mixed state by employing the results reported in Sec. 3.
Given a mixed state that is not pure and that is characterised by the covariance matrix γ R , any circuit connecting γ R to a pure state provides a purification path. A purification path connects the covariance matrices γ R to γ T whose Williamson's decompositions are given respectively by A q e Z D I n 5 9 y b e + / x Q s E 1 2 P a X l V p a X l l d S 6 9 n N j a 3 t n e y u d 2

o H z z c z n G t S T p Z g t 2 0 Z 4 C L x N n T g r l n J v / r h 5 + V L r Z T 7 c X 0 M h n E q g g W r c d O 4 R O T B R w K l i S c S P N Q k J H Z M D a h k r i M 9 2 J p 4 c k + M Q o P d w P l H k S 8 F T 9 3 R E T X + u x 7 5 n K y d Z 6 0 Z u I / 3 n t C P o X n Z j L M
h v 3 5 t 9 n y I = < / l a t e x i t > T < l a t e x i t s h a 1 _ b a s e 6 4 = " x n e m I w = < / l a t e x i t > W T < l a t e x i t s h a 1 _ b a s e 6 4 = " T 8 D V a J 2 Z C p 9 6 G 8 2 n A D 7 z u H t e Z T k = " k C e y L N 3 7 z 1 6 L 9 7 r c H T E + 9 p Z I j / g v X 0 C 1 n + e q g = = < / l a t e x i t >˜ Figure 3: The optimal purification paths for γ R and γ T correspond respectively to the W R path and to the W T path, that are represented through dashed lines. The straight black solid line represent the set of the pure states, whose symplectic spectrum is given by D = 1 2 1.
where W R ∈ Sp(2N, R) and D = 1 2 1 are assigned, while W T ∈ Sp(2N, R) is not. Among all the possible paths, the optimal circuit is obtained by specifying (2.28) to (4.1). The result is which depends on the symplectic matrix W T that determines the final pure state. The length of the purification path (4.2) can be found by evaluating (2.31) for the special case described by (4.1). It reads The optimal purification path is the purification path with minimal length, which can be found by minimising (4.3) as W T ∈ Sp(2N, R) varies within the symplectic group. This extremization procedure selects a symplectic matrix W 0 that determines the pure state through its covariance matrix 1 . This is a special case of the analysis performed in Sec. 3.2 corresponding to δD R = δD T = δW R = 0.
In Sec. 3.2 we have shown that a W path provides a solution to this equation, namely which is the trivial solution corresponding to δW R = 0. In the following we focus on the purification process based on the W paths. We cannot prove that, among all the solution of δd 0 = 0 (see Sec. 3.2), the W path corresponds to the one having minimal length.
The W R path connects the mixed state γ R = W t R D W R to the pure state γ 0 = 1 2 W t R W R . By specialising (3.17) to D T = 1 2 1, we find that this W R path is given by and its length can be easily obtained by setting D T = 1 2 1 in (3.19), finding an expression that depends only on D It is instructive to focus on the one-mode mixed states, when the covariance matrices (4.1) Any purification path corresponding to a geodesic can be written as in (2.53) with σ s = 1 2 (2σ) 1−s and W s defined in (2.54). In particular, for the W R path we have W s = W R and its length is given by The thermal states are interesting examples of mixed states to explore. The Williamson's decomposition of the covariance matrix of a thermal state is given by (2.76). By specialising (4.5) to this case, we obtain the W path that purifies a thermal state. It reads where it is worth reminding that the symplectic matrix W th , given in (2.77), does not depend on the temperature of the thermal state, but only on the parameters occurring in the hamiltonian.
It is natural to ask whether the W th path (4.7) is made by thermal states. This is the case if, for any given s ∈ [0, 1], the symplectic spectrum of (4.7) is thermal at some inverse temperature β s determined by the inverse temperature β of the thermal state that plays the role of the reference state in this purification path. Using (2.77), this requirement leads to the following condition coth(βΩ k /2) 1−s = coth(β s Ω k /2) (4.8) which corresponds to (2.82) when β T → ∞, as expected. This condition depends on the dispersion relation of the model. A straightforward numerical inspection for the periodic chains (see Sec. 9.1) shows that (4.8) cannot be solved by β s = β s (β) for any 1 k N ; hence we conclude that the purification path (4.7) is not entirely made by thermal states.
The W paths provide a natural alternative way to connect two generic mixed states γ R and γ T by using a path that passes through the set of pure states. In particular, by exploiting the Williamson's decompositions given in (2.41), one first considers the W R path that connects γ R to the pure stateγ R,0 and the W T path that connects γ T to the pure stateγ T,0 . From (4.5), these W paths are given respectively by Then, within the set of the pure states, we consider the geodesic connectingγ R,0 toγ T,0 . Our preferred path to connect γ R and γ T passing through the set of pure states is obtained by combining these three paths as follows The length d pur (γ R , γ T ) of this path can be found by summing the lengths of its three components. From (3.23) and (4.6), we get which can be written more explicitly as follows (4.13) This expression provides an upper bound d(γ R , γ T ) d pur (γ R , γ T ) on (2.31).

Bounding complexity
The complexity (2.33) is proportional to the length of the optimal circuit (2.28) connecting γ R to γ T ; hence it is straightforward to observe that the length of any other path connecting these two covariance matrices provides an upper bound on the complexity. The analysis reported in Sec. 3 and in Sec. 4 naturally leads to consider some particular paths.
The simplest choice is a path made by two geodesics that connect γ R and γ T to an auxiliary covariance matrix γ aux that does not belong to the optimal circuit (2.28) (i.e. that does not lie on the black solid curve in Fig. 2). Natural candidates for γ aux are the covariance matrices whose Williamson's decompositions contain either D R or D T or W R or W T . For instance, we can choose for γ aux a covariance matrix whose symplectic spectrum is D R or a covariance matrix whose symplectic spectrum is D T (that lie respectively on the red solid curve and on the blue solid curve in Fig. 2). Different choices for γ aux lead to different bounds; hence it is worth asking whether a particular choice provides the best bound. The answer depends on the set where γ aux is allowed to vary.
Let us consider some natural paths where only a single auxiliary covariance matrix γ aux is involved. In Sec. 3.2 we have shown that the W paths satisfy the first law of complexity with δd = 0. Thus, natural candidates to consider for γ aux arẽ that have been represented by black squares in Fig. 2 and Fig. 3.
By first applying the triangle inequality for the paths γ R →γ R → γ T and γ R →γ T → γ T , and then picking the one that provides the best constraint between the D paths, we obtain Denoting byd(γ R , γ T ) the r.h.s. of this inequality, by using (3.19) and (3.23) we find that The path γ R →γ T → γ T corresponds to an explicit realisation of the proposal made in Fig. 6 of [28] within the approach that we are considering, that does not require the addition of ancillary degrees of freedom.
Better bounds could be obtained by constructing paths that involve more auxiliary covariance matrices γ aux . For instance, one can consider paths γ R → γ aux,1 → γ aux,2 → γ T that involve two auxiliary covariance matrices γ aux,1 and γ aux,2 . Referring to Fig. 2, natural paths to consider within this class are the ones where γ aux,1 belongs to the W R path and γ aux,2 to the W T path, or the ones where γ aux,1 belongs to the red curve (i.e. its symplectic spectrum is D R ) and γ aux,2 belongs to the blue curve (i.e. its symplectic spectrum is D T ).
Another interesting path to consider is the one constructed in (4.11): it involves the two auxiliary matrices γ aux,1 =γ R,0 and γ aux,2 =γ T,0 and its length is (4.12) (see Fig. 3). It is straightforward to observe that d(γ R , γ T ) d pur (γ R , γ T ), but it is non trivial to find the best bound betweend(γ R , γ T ) and d pur (γ R , γ T ). Since we cannot provide a general solution to this problem, in the following we focus on simple special cases where we can show that . Another class of special cases that we find interesting to consider is given by the pairs (γ R , γ T ) such that all the matrices along the D R path connecting γ R toγ R have the same symplectic spectrum D R and, similarly, all the matrices along the D T path connecting γ T toγ T have the same symplectic spectrum D T . This means that (3.26) holds for both D = D R and D = D T ; hence (5.3) simplifies tõ The first square root in the r.h.s. can be bounded as follows where we have employed first that all the elements of 2D are larger than or equal to 1 (in order to discard a positive term under the square root) and then the inequality √ a + b √ a + √ b, that holds for any a and b. By employing (5.5) in (5.4) and comparing the result against (4.13), we can conclude thatd(γ R , γ T ) d pur (γ R , γ T ).
In the most general case where the first moments are non vanishing r ≡ a = 0, a closed expression for the Fisher-Rao distance is not known, as also remarked in the Appendix B, where the notation µ = a and Σ = γ has been adopted. Nonetheless, lower and upper bounds on the complexity can be written by employing some known results about the Gaussian PDF's [79,83,85,101,102].
Given a reference state and a target state, that can be parameterised by θ R = (a R , γ R ) and θ T = (a T , γ T ) respectively, let us introduce the following (2N + 1) × (2N + 1) matrix A lower bound for the Fisher-Rao distance, first obtained in [101], is given by Upper bounds for the Fisher-Rao distance have been also found for non vanishing first moments [79,83,85,102]. An upper bound can be written through d , and O is the orthogonal matrix whose columns are the eigenvectors of γ Another upper bound has been found in [85]. It has been written by introducing the 2N × 2N orthogonal matrix O such that O a T,R = (|a T,R |, 0, . . . , 0) and the following 2N × 2N matrices These matrices are employed to identify the states corresponding to the following vectors where d a T,R is defined in (B.6) and d diag in (B.9). Since an inequality between the two upper bounds in (5.8) and (5.11) cannot be found for any value of θ R and θ T [83], we pick the minimum between them.
Combining the above results, one obtains In order to provide a consistency check for these bounds, let us consider the case a R = a T = a. From (5.6) we obtain By employing a formula for the determinant of a block matrix reported below (see (8.14)), one finds that the first 2N eigenvalues of (5.13) are the eigenvalues of γ T γ −1 R , while the last eigenvalue is equal to 1. Thus, d lower in (5.7) becomes (2.31) in this case, saturating the lower bound. As for the upper bound in (5.12), we have (5.10) simplify to θ 0 = θ O = θ γ = (0, 1) in this case. This implies that d upper,2 in (5.11) becomes (2.31); hence also the upper bound is saturated.

Optimal path for entanglement hamiltonians
The density matrix of a mixed state can be written as followŝ where the proportionality constant determines the normalisation ofρ. We denote the operator K as entanglement hamiltonian, with a slight abuse of notation. Indeed, the operator K is the entanglement hamiltonian whenρ =ρ A = Tr H Bρ 0 is the reduced density matrix obtained by tracing out the part H B of a bipartite Hilbert space H = H A ⊗ H B . Instead, for instance, the thermal states are mixed states that do not correspond to a bipartition of the Hilbert space. For these states K = β H, being H the hamiltonian of the system and β the inverse temperature.
The entanglement hamiltonians associated to some particular reduced density matrices have been largely studied for simple models, both in quantum field theories [47,[103][104][105][106][107][108][109][110] and on the lattice [49,91,[111][112][113][114][115][116][117]. The spectrum of the entanglement hamiltonian, that is usually called entanglement spectrum [118], is rich in information. For instance, in conformal field theories the entanglement spectrum provides both the central charge [119] and the conformal spectrum of the underlying model [106,109,110,116,117,[120][121][122] For the bosonic Gaussian states that we are considering, the entanglement hamiltonians are quadratic operators in terms of the position and momentum operators; hence they can be written as follows where H is a 2N ×2N symmetric and positive definite matrix that characterises the underlying mixed state. We denote H as the entanglement hamiltonian matrix. It can be written in terms of the corresponding covariance matrix γ as follows [47,91,[111][112][113]117] where J is the standard symplectic matrix (2.2). The expression (6.3) holds for covariance matrices that are not associated to pure states. Thus, in particular, the purification procedure reported in Sec. 4 cannot be described through the entanglement hamiltonian matrices H defined by (6.2).
Since the matrix H is symmetric and positive definite, we can adapt to the entanglement hamiltonian matrices many results reported in the previous discussions for the covariance matrices.
Given the matrices H R and H T corresponding to the reference state γ R and to the target state γ T respectively, we can consider the optimal path connecting H R to H T , namely The symplectic spectrum of H can be determined from the symplectic spectrum of γ as follows [91,111,112] This formula cannot be applied for pure states, which have D = 1 2 1. The symplectic matrices W and W , introduced in (2.20) and (6.6) respectively, are related as follows [113,117] W ≡ J t W J = W −t . (6.8) We find it worth expressing the distance (6.5) in terms of the matrices occurring in the Williamson's decompositions of H R and H T , as done in Sec. 2.4 for the covariance matrices. These decompositions read where W R and W T are related respectively to W R and W T through (6.8). By using (2.49) and the following relation we can write the distance (6.5) as The expression (6.3) (or equivalently (6.7) and (6.8)) provides a highly non trivial relation between the set made by the covariance matrices γ that are associated to the mixed states that are not pure states and the set of the entanglement hamiltonian matrices H. The map h in (6.3) is not an isometry, hence the distances are not preserved and geodesics are not sent into geodesics. Thus, we find it worth comparing the distance d(γ R , γ T ) = d(h −1 (H R ), h −1 (H T )) from (2.49) and the distance d(H R , H T ) in (6.11).
For the sake of simplicity, let us explore the case of one-mode mixed states, where D = σ 1 and E = ε 1 are proportional to the 2 × 2 identity matrix. In this simple case the expressions for d(γ R , γ T ) 2 from (2.49) and for d(H R , H T ) 2 from (6.11) take the form 13 with a = log(σ T /σ R ) ≡ a σ and a = log(ε T /ε R ) ≡ a ε respectively, which can take any real value.
Being d(γ R , γ T ) symmetric under the exchange γ R ↔ γ T , we can assume σ R σ T without loss of generality. Then, since the function 2arccoth(2x) x is a properly decreasing function when i.e. σ T /σ R ε T /ε R , once (6.7) has been used; hence a σ a ε . This does not provide a relation between d(γ R , γ T ) and d(H R , H T ) because the r.h.s. of (6.12) does not have a well defined monotonicity as function of a, being log(W TR W t TR ) non vanishing in general. Thus, the one-mode case teaches us that W TR plays a major role in finding a possible relation between d(H R , H T ) and d(γ R , γ T ). In order to find this relation in some simple cases, the expression (6.12) naturally leads us to consider the special cases of one-mode mixed states such that log(W TR W t TR ) = 0. In this cases (6.12) a 2 ε is equivalent to (a σ − a ε )(a σ + a ε ) 0, we observe that the latter inequality is satisfied because a σ a ε and 14 a σ + a ε = log σ T arccoth(2σ T ) When N 1 and W T = W R , i.e. W TR = 1 (this includes the thermal states originated from the same physical hamiltonian), the distance (6.11) simplifies to while d(γ R , γ T ) is given by (3.19). By applying the above analysis made for the one-mode case to the k-th mode, we can conclude that [log(σ T,k /σ R,k )] 2 [log(ε T,k /ε R,k )] 2 for any given k; hence d(γ R , γ T ) 2 d(H R , H T ) 2 is obtained after summing over the modes.
By using the decompositions (6.9), one can draw a pictorial representation similar to Fig. 1 and Fig. 2 also for the entanglement hamiltonian matrices H, just by replacing each γ with the corresponding H, each W with the corresponding W and where the solid coloured lines are labelled by the corresponding symplectic spectra E.
We find it worth discussing further the set of thermal states through the approach based on the entanglement hamiltonian matrices because the simplicity of these matrices in this case allows to write analytic results. For a thermal state H = βH phys , where H phys is the matrix characterising the physical hamiltonian (2.1) and β is the inverse temperature. This implies that the symplectic eigenvalues of H are ε th,k = β σ phys,k , being σ phys,k the symplectic eigenvalues of H phys .
We denote by β R and β T the inverse temperatures of the reference state and of the target state respectively. An interesting special case is given by thermal states of the same system, which have the same H phys . In this case H T H −1 where V is the number of sites in the harmonic lattice and the last expression has been obtained by specialising (6.11) to this case, where W TR = 1. Furthermore, from (6.4) it is straightforward to observe that in this case the entire optimal circuit is made by thermal states having the same H phys . The optimal circuit (6.4) significantly simplifies to By employing (2.68), one finds that the Williamson's decomposition of this optimal circuit reads where W phys is independent of s. Thus, (6.15) tells us that β s is the inverse temperature of the thermal state labelled by s along this optimal circuit.
In Sec. 9.3.3 the above results are applied to the thermal states of the harmonic chain with periodic boundary conditions. 7 Gaussian channels Quantum operations are described by completely positive operators acting on a quantum state, which can be either pure or mixed, and they are classified in quantum channels and quantum measurements [87,124]. The quantum channels are trace preserving quantum operations, while quantum measurements are not trace preserving [125].
The output Θ(ρ) of a quantum channel applied to the density matrixρ of a system is obtained by first extending the system through an ancillary system (the environment) in a pure state |Φ E , then by allowing an interaction characterised by a unitary transformation U and finally by tracing out the degrees of freedom of the environment [43,86,126], namely While within the set of the pure states the unitary transformations are the only operations that allow to pass from a state to another, within the general set of mixed states also non unitary operations must be taken into account.
In this manuscript we consider circuits in the space made by quantum Gaussian states; hence only quantum operations between Gaussian states (also called Gaussian operations) can be considered [125]. The quantum channels and the quantum measurements restricted to the set of the Gaussian states are often called Gaussian channels [43,99] and Gaussian measurements [44] respectively.
In the following we focus only on the Gaussian channels. A Gaussian state with vanishing first moments is completely described by its covariance matrix; hence the action of a Gaussian channel on a Gaussian state can be defined through its effect on the covariance matrix of the Gaussian state. This effect can be studied by introducing two real matrices T and N as [99] γ where T is unconstrained and the last inequality corresponds to the complete positivity condition. The Gaussian unitary transformations are the Gaussian channels with N = 0 and symplectic T . In this case the inequality in (7.2) is saturated. Further interesting results for Gaussian operations have been reported e.g. in [44,127].
We find it worth asking whether a matrix along the optimal circuit (2.28) can be obtained by acting with a Gaussian channel on the reference state. This means finding T s and N s that fulfil (7.2) for any 0 s 1 and such that where U s is defined in (2.39). Unfortunately, we are not able to determine T s and N s as functions of U s in full generality. In the following we provide some simple particular solutions.
A simple possibility reads which satisfies the inequality in (7.2), being G s a symmetric Gaussian matrix (see Sec. 2.4).
Another, less trivial, solution is given by where the complete positivity condition in (7.2) becomes i J 2 − i T s J 2 T t s 0. We have considered numerically some cases, finding that U s satisfies the complete positivity condition only when it is symplectic (in this case the complete positivity inequality is saturated).
An explicit example belonging to the class identified by (7.5) can be constructed by considering a particular D path where W T = X TR W R (see Sec. 3.4). In this case (3.25) holds, hence (7.5) is realised with 15 T A more general solution where both T s and N s can be non vanishing is obtained by imposing the following relation which solves (7.3) for any symmetric N s . The solution (7.4) is recovered from (7.7) with T s = 0. When N s = 0, the relation (7.7) gives T s = G s γ −1 R 1/2 = U s , where the last equality 15 The last expression in (7.6) is obtained by observing that WT = XTRWR and DT = DR into (2.39) give is obtained from (2.30) and (2.39). Plugging (7.7) into the complete positivity condition in (7.2), we obtain Thus, for any N s fulfilling this inequality, by using (7.7) we can implement our optimal circuit (2.28) through Gaussian channels.
An interesting class of N s that saturates (7.8) has been constructed in [128]. It is given by By plugging (7.9) into (7.7) first and then employing (2.30), we find the following equation whose solutions provide realisations of the optimal circuit (2.28) through Gaussian channels.
Plugging the definition of K s given in (7.9) into (7.10) we find The real part of this relation tells us that T s = ∆ s/2 TR = U s , while from the imaginary part we find that T t s J T s = J, i.e. that T s is symplectic. The latter result and (7.9) lead to K s = N s = 0.
Let us conclude by emphasising that all the explicit expressions for the Gaussian channels given above saturate the complete positivity condition in (7.2) (more details can be found in [128]). It would be interesting to explore also Gaussian channels where this inequality is not saturated, as done e.g. in (7.7) and (7.8).

Complexity of mixed states through ancillae
In this section we discuss the approach to the complexity of mixed states explored in [23], which is based on the introduction of ancillary degrees of freedom. Consider a quantum system in a mixed state characterised by the density matrixρ. A pure state can be always constructed fromρ by adding ancillary degrees of freedom. This purification procedure consists in first extending the Hilbert space of the system to a larger Hilbert space H extended ≡ H ⊗ H anc through an auxiliary Hilbert space H anc , and then finding a pure state |Ω ∈ H extended such that the original mixed state is obtained as the reduced density matrix given byρ = Tr Hanc |Ω Ω| (8.1) where the ancillary degrees of freedom have been traced out. We remark that the purifications discussed in Sec. 4 do not involve ancillary degrees of freedom.
There are infinitely many ways to construct H extended and |Ω such that (8.1) is satisfied; hence a purification criterion must be introduced. Different purtification criteria have been considered in the literature to study different quantities. An important example is the entanglement of purification [129][130][131][132]. In this manuscript we are interested in the purification complexity [28], that has been employed in [23] to study the complexity of mixed states.

Covariance matrix of the extended system
We are interested in a generic harmonic lattice made by N sites in the Gaussian mixed state characterised by the covariance matrix γ and by vanishing first moments. The covariance matrix γ can be decomposed in blocks as follows where Q and P are N × N symmetric matrices, while M is a generic N × N real matrix; hence N (2N + 1) real parameters must be fixed to determine γ.
We consider a simplification of the purification process by focussing only on Gaussian purifications. This means that a mixed state characterised by the covariance matrix (8.2) is purified by introducing ancillary degrees of freedom and construcing a 2N ext ×2N ext covariance matrix γ ext that corresponds to a Gaussian pure state |Ω for the extended lattice having N ext ≡ N + N anc sites. For the sake of simplicity, we assume also that |Ω has vanishing first moments, i.e. Ω|r ext |Ω = 0, beingr t ext ≡ (q t ,q t anc ,p t ,p t anc ), where we have separated the ancillary degrees of freedom from the ones associated to the physical system.
By writing also γ ext through the block decomposition (8.2) we have where Q ext and P ext are N ext × N ext symmetric matrices. Since the covariance matrix (8.3) corresponds to a pure state, the condition (2.23) must hold. This tells us that the blocks occurring in (8.3) are related by the following constraints The first relation tells us that M ext is determined by the product Q ext P ext , while the remaining two relations mean that M ext Q ext and P ext M ext are symmetric. Thus, (8.3) is determined by the symmetric matrices Q ext and P ext , that depend on 2 Next(Next+1) 2 real parameters, as expected also from Sec. 2.3.1 (see (2.24)).
We can impose that γ ext is the covariance matrix of a pure state also by using (2.22), i.e. by requiring that the Williamson's decomposition of (8.3) reads where W ext ∈ Sp(2N ext , R) and the last expression has been obtained from the Euler decomposition of W ext , that includes R ext ∈ K(N ext ). The symplectic matrix W ext can be partitioned through N ext × N ext matrices as follows The relation (8.5) provides the blocks in (8.3) through to the ones in (8.6). The result reads Another useful way to impose the purity condition on the final state of this purification process exploits the general form (2.24) for the wave function of a pure state and the corresponding covariance matrix (2.25). This allows us to write the covariance matrix of the extended system as 16 where the N ext × N ext symmetric matrices E ext and F ext are related to the blocks occurring in (8.3) as follows The second relation in (8.4) ensures that F ext is symmetric. We remark that (8.9) also tell us that the relation P ext = 1 2 (E ext + F ext E −1 ext F ext ) coming from the second block on the diagonal in (8.3) becomes the first relation in (8.4).
In order to relate γ in (8.2) to γ ext in (8.3), one observes that, beingr t ext ≡ (q t ,q t anc ,p t ,p t anc ), we have that the N ext × N ext blocks occurring in (8.3) can be partitioned in blocks as follows where Q, P and M are the N × N blocks of γ in (8.2), while Q anc and P anc are N anc × N anc symmetric matrices. Instead, M anc is a generic N anc × N anc real matrix. Indeed, by plugging (8.10) into (8.3), it is straightforward to observe that the covariance matrix (8.2) is obtained by restricting γ ext to the sites corresponding to the original degrees of freedom. Instead, by restricting γ ext to the ancillary sites, one gets the following 2N anc × 2N anc symmetric matrix By changing the order of the rows and the columns, the matrix in (8.3) becomes γ Γ Γ t γ anc (8.12) where γ is the covariance matrix (8.2), γ anc is the symmetric matrix defined in (8.11) and By using that (8.12) is positive definite, it can be shown that also γ anc is positive definite 17 ; hence γ anc can be interpreted as the covariance matrix of the ancillary system made by N anc sites.
An alternative approach exploits the expressions in the Schrödinger representation discussed in Appendix A. In particular, given the covariance matrix γ in the block matrix form (8.2), we can construct the N × N complex matrices Θ and Φ by using (A.11) and (A.12). Then, (A.19) provide the constraints for the blocks of E ext and F ext in terms of the complex matrices Θ and Φ.
There are many ways to construct the pure state |Ω . They correspond to the freedom to fix N anc first and then to choose e.g. the blocks in (8.10) that are different from Q, P and M , provided that the constraints (8.4) are satisfied.

One-mode mixed states
We find it instructive to consider explicitly the simplest case of a one-mode mixed state, i.e. N = 1. The minimal choice for the number of ancillae is N anc = 1.
When N = 1, only a non trivial symplectic eigenvalue σ occurs; hence the Williamson's decomposition (2.20) and the Euler decomposition (2.21) of a symplectic matrix provide the 2 × 2 covariance matrix given by where λ is the squeezing parameter and R is a 2 × 2 rotation matrix, which is completely fixed by the rotation angle θ. Notice that Q, P and M are real parameters in (8.15). Le us remark that the pure state condition (2.23) for (8.15) gives 1 − 4d γ = 0, where we have introduced d γ ≡ det(γ) = QP − M 2 . This implies that 1 − 4d γ = 0 for the covariance matrices (8.15) that correspond to the mixed states that are not pure.
When N anc = 1, the covariance matrix (8.5) of the pure state for the extended system reads This 4×4 covariance matrix corresponds to a pure state, hence it depends on N ext (N ext +1) = 6 real parameters (2 2 from R ext and two squeezing parameters λ i ), being N ext = 2. Writing the 4 × 4 covariance matrix γ ext in the form (8.3), it is straightforward to realise that 3 elements are given by the real parameters Q, P and M . Thus, we are left with three real parameters to construct the pure state for the extended system. 17 By employing the following formula for the determinant of a block matrix where it is assumed that D is invertible, one finds that the eigenvalues of γanc are also eigenvalues of γext. If A is invertible, a formula similar to (8.14) can be written where det(A) is factorised and this result can be used to show that the eigenvalues of γ are eigenvalues of γext as well.
We find it instructive to write explicit expressions for the elements of the covariance matrix γ ext . The constraints (8.4) for the 2 × 2 matrices Q ext , P ext and M ext provide six equations: four from the first relation and one from each one of the other two relations (that can be written in the form X = 0, where X is a 2 × 2 antisymmetric matrix).
When Γ M = 0, the solution of this system can be written in terms of Γ Q , Γ P and Γ M as When Γ M = 0 and Γ P = 0, we find while a solution does not exist for Γ M = Γ P = 0. Notice that 1 − 4d γ = 0 in these expressions because γ does not corresponds to a pure state.
We remark that also the analysis based on the Schrödinger representation reported in the Appendix A.2 allows to conclude that the purification of a one-mode mixed state can be realised through a pure state in an extended lattice with N ext = 2 that depends on three real parameters.

Block diagonal covariance matrices
Many interesting mixed states are described by a block diagonal covariance matrix γ = Q⊕P . In this cases M = 0 in (8.2).
It is worth considering a pure state for the extended system such that M ext = 0 in the corresponding covariance matrix (8.3). In this case (8.4) reduce to where Γ Q Γ t P = 0, being γ = Q ⊕ P the covariance matrix of a state that is not pure.
A common choice consists in considering purifications where the extended system has twice the degrees of freedom occurring in the original one, namely N anc = N . In these cases Γ Q , Γ P and Γ M are N × N matrices.
Considering the purifications with N anc = N , a drastic simplification corresponds to require that γ = γ anc , which is equivalent to impose that Q = Q anc and P = P anc . In this case, a solution is given by symmetric and commuting matrices Γ Q and Γ P that can be related through the last two equations in (8.19), which give Setting Γ P = α Γ −1 Q with α ∈ R, the last equality is solved while the remaining relation Q = − 1 α Γ Q P Γ Q , whose validity is not guaranteed, provides Γ Q . A different solution for the matrix equations in (8.20) can be written when Q and P can be decomposed through three real matrices A, B and Λ as follows In this case, we can construct Γ Q and Γ P as where a new matrixΛ that commutes with Λ has been introduced.
It is straightforward to check that (8.21) and (8.22) satisfy the matrix equations in (8.20). Notice that Γ P is not proportional to Γ −1 Q in (8.22). An important example where N anc = N and γ = γ anc is the thermofield double state (TFD). In Appendix F a detailed analysis for this pure state for harmonic lattices is reported. The relations (F.25) and (F.26) tell us that the TFD corresponds to a special case 18 of (8.21) and (8.22).
The simplest case corresponds to N = N anc = 1, which has been discussed in Sec. 8.1.1 in the most general setting. Solving the system (8.19) for this case, one finds When M = 0, the observation in the text below (8.15) tells us that 4QP − 1 = 0 in order to have a mixed state that is not pure to purify. Notice also that, by setting M = Γ M = 0 in (8.18), that holds for Γ M = 0, one finds (8.23) and M anc = 0. Thus, when M ext = 0 and N = N anc = 1 we can parameterise the pure state of the extended system through a single parameter. This is consistent with the analysis reported in [23]. As final remark about the purifications having N = N anc = 1, let us observe that the second equation in (8.20) is trivially satisfied, while the first one is obtained by setting Q = Q anc and P = P anc in (8.23).

Selection criterion for the pure state
In the previous discussion we have explored the constraints guaranteeing that the covariance matrix γ ext corresponds to a pure state under the condition that γ ext provides the covariance matrix γ of the given mixed state once the ancillary degrees of freedom have been traced out.
These constraints identify the parameter space of the pure states allowed by γ for a given value of N anc . Within this space of parameters, it is natural to introduce a quantity F whose minimisation provides a particular pure state with certain properties. Thus, F characterises the criterion to select the pure state provided by the purification procedure as follows where γ ext is the covariance matrix for the extended system, that is constrained as described in Sec. 8.1, and F denotes the minimal value of F as γ ext spans all the pure states allowed by γ. For the bosonic Gaussian states that we are considering, the calculations can be performed by employing either the wave functions or the covariance matrices.
For instance, the entanglement of purification for a bipartite mixed state [129][130][131][132] is (8.24), with F given by the entanglement entropy of a particular bipartition of H extended .
In [23,28], the purification complexity has been introduced to quantify the complexity of a mixed state. The definition of purification complexity is given by (8.24) in the special case where F is the complexity of the pure state corresponding to γ ext with respect to a given fixed pure state in H extended , whose covariance matrix is denoted by γ ext,0 . This definition of purification complexity requires the choice of a cost function. The purification complexity explored in [23] reads C r (γ) ≡ min γext C r (γ ext , γ ext,0 ) (8.25) where either r = 1 or r = 2, depending on whether the F 1 cost function or the F 2 cost function is adopted. In [23] the purification complexity based on the F 1 cost function has been mainly studied because, for the pure states, the divergence structure of the complexity evaluated through the F 1 cost function is closer to the one obtained from holographic calculations [17,32]. The complexity defined through the F 1 cost function depends on the choice of the underlying basis, while the F 2 cost function leads to a complexity that is independent of this choice.
This approach to the complexity of mixed states is different from the one considered in this manuscript. The main difference is due to the fact that in the purification procedure described in Sec. 4 ancillary degrees of freedom have not been introduced. Moreover, the purification complexity defined in (8.25) depends on the choice of the pure state corresponding to γ ext,0 (in [23] this pure state has been fixed to the one whose wave function (2.24) has E ext ∝ 1 and F ext = 0). Furthermore, in the evaluation of the complexity of a mixed state through (8.25), no cost is assigned to the purification process of extending the system through ancillary degrees of freedom, being the circuit considered in (8.25) entirely made by pure states in H extended .
Explicit computations through (8.25) are technically involved and discussing them is beyond the scope of this manuscript. We refer the interested reader to the detailed analysis performed in [23]. Focussing on the simple case of one-mode thermal states, in Sec. 9.7 we compare the complexity evaluated through the Fisher-Rao distance with the results found in [23] for the C 1 complexity of mixed states based on the purification complexity. The latter quantity depends on the basis: in Appendix G we discuss the diagonal basis and the physical basis, that have been introduced in [23] to evaluate this C 1 complexity.

Harmonic chains
In this section we further study some of the quantities discussed in the previous sections by focussing on the one-dimensional case of the harmonic chain, either on the circle (i.e. with periodic boundary conditions) or on the infinite line. In this case we obtain analytic expressions in terms of the parameters of the circuit for some quantities and provide numerical results for the quantities that are more difficult to address analytically. After a brief discussion of the model in Sec. 9.1, circuits whose reference and target states are either pure or thermal are considered in Sec. 9.2 and Sec. 9.3 respectively. In Sec. 9.4 we study the mutual complexity for the thermofield double states (TFD's). Numerical results for the complexity and the mutual complexity associated to subregions are presented in Sec. 9.5 and Sec. 9.6. Finally, in Sec. 9.7 we consider a simple comparison between the complexity for mixed states discussed in this manuscript and the one based on the purification complexity recently proposed in [23].
For the sake of simplicity, in this section we consider only examples that involve states whose covariance matrices are block diagonal. We remark that the results discussed in the previous sections hold also for states characterised by covariance matrices that are not block diagonal. For instance, these states typically occur in the out-of-equilibrium dynamics of the harmonic lattices [117,[133][134][135].

Hamiltonian
The hamiltonian of the periodic harmonic chain made by L sites, with frequency ω, mass m and elastic constant κ reads wherer ≡ (q 1 , . . . ,q L ,p 1 , . . . ,p L ) t collects the position and momentum operators and the periodic boundary conditionq L+1 =q 1 is imposed. The canonical transformation given bŷ q i →q i / 4 √ mκ andp i → 4 √ mκp i allows to write (9.1) as and we are naturally led to introduceω The non vanishing elements of the symmetric matrix T in (9.3) are T i,i+1 = T i+1,i = 1 with 1 i L − 1 and T 1,L = T L,1 = 1.
In order to find the Williamson's decomposition (2.68) for (9.3), first one observes that the matrix T in (9.3) is diagonalised by the following unitary matrix U r,s ≡ e 2πi r s/L √ L (9.5) that implements the discrete Fourier transform and it is independent of the parameters ω, m and κ. This implies that H phys in (9.3) is diagonalised by U ≡ U ⊕ U .
Since the symplectic matrix entering in the Williamson's decomposition (2.68) is real, let us consider the proper combinations of the eigenvectors entering in (9.5) that correspond to the same eigenvalue. This leads to introduce the L × L real and orthogonal matrix V , whose generic element for even L is given by  The matrix V diagonalises both T and 1 in (9.3); hence, by introducing the orthogonal matrix V ≡ V ⊕ V , that is also symplectic, we have [45] H phys = V κ/m diag Ω 2 1 , . . . , Ω 2 L , 1, . . . , 1 V −1 (9.8) where Ω k provides the dispersion relation, which depends on the parameterω defined in (9.4) as follows Ω k ≡ ω 2 + 4 (sin[πk/L]) 2 k = 1, . . . , L . In these expressions the zero mode corresponds to k = L and its occurrence is due to the invariance of the system under translations. The comparison between the expressions reported throughout this section and the corresponding ones in Sec. 2.6 can be done once the canonical transformation above (9.2) has been taken into account 19 .

Pure states
In this subsection we study the circuit complexity for pure states that are the ground states of periodic harmonic chains having different frequencies [17].

Covariance matrix
The two-point correlators in the ground state of the hamiltonian (9.2), where periodic boundary conditions are imposed, read being Ω k the dispersion relation (9.9). The periodic boundary conditions make this system invariant under translations. The expressions in (9.12) define the elements of the correlation matrices Q gs and P gs respectively. These matrices provide the block diagonal covariance matrix γ gs = Q gs ⊕ P gs .
By introducing the discrete Fourier transform of the operatorsq j andp j in the standard way [136], the matrices Q gs and P gs can be written as follows Q gs = U Q gs U −1 P gs = U P gs U −1 (9.13) where the matrix U have been defined in (9.3), while Q gs and P gs are diagonal matrices whose elements read (Q gs ) k,k = 1 2 Ω k (P gs ) k,k = 1 2 Ω k k = 1, . . . , L . (9.14) In order to find the Williamson's decomposition of γ gs , we have to consider the symplectic matrix V introduced in Sec. 9.1. Then, the observation made in the final part of the Appendix D specified to γ gs leads us to introduce the following symplectic and diagonal matrix where X phys has been introduced in (9.11). This matrix is related to (9.14) as follows By introducing the symplectic matrix we have that the Williamson's decomposition of γ gs reads Notice that the symplectic matrices W phys and W C , defined in (9.11) and (9.17) respectively, are related as follows where we have also used (9.15), (9.17), the property S −t = J t SJ of the symplectic matrices S and the fact that symplectic V is also orthogonal (see also (2.77)).
As consistency check, we can plug (9.17) into (9.18) first and then use (9.16), finding that This tells us that V is also the orthogonal matrix that diagonalises the symmetric matrix γ gs . Let us remark that V depends only on the number of sites L of the harmonic chain.

Complexity
We consider the circuit complexity where the reference state and the target state are the ground states of periodic harmonic chains whose hamiltonians are characterised by the parameters (ω R , κ R , m R ) and (ω T , κ T , m T ) respectively.
For the sake of simplicity, in our analysis we set κ R = κ T ≡ κ and m R = m T ≡ m; hence only the parameterω distinguishes the reference and the target states. In this case, from (9.20) and the fact that V is independent of the parameters ω, κ and m, it is straightforward to find that (2.32) becomes ∆ TR = V Q gs,T Q −1 gs,R ⊕ P gs,T P −1 gs,R V −1 (9.21) where the diagonal matrices Q gs,R , Q gs,T , P gs,R and P gs,T can be easily obtained by writing (9.14) for the reference and for the target state. By employing (9.21) and (9.14), it is straightforward to find that, in this case, the complexity given by (2.31) and (2.33) simplifies to [17] where Ω S,k ≡ ω 2 S + 4 (sin[πk/L]) 2 k = 1, . . . , L S ∈ R, T . It is instructive to obtain (9.22) also as a special case of (2.58), that is written in terms of the matrix W TR introduced in (2.45). For the pure states that have been chosen, whose covariance matrices have the form (9.18), W R and W T can be obtained by specialising (9.17) to the reference and the target states considered. Since the matrix V is the same for both of them, (2.45) simplifies to , . . . , Ω T,L Ω R,L 1/2 (9.24) where the last expression has been found by using (9.15). Thus, in this case W TR = X TR = X C,T X −1 C,R and this leads to obtain (9.22) from (2.58). In the thermodynamic limit L → ∞, the expression (9.22) for the complexity becomes [17] where the subleading terms have been neglected and the coefficient of the leading term reads being Ω S,θ ≡ ω 2 S + 4 (sin θ) 2 θ ∈ (0, π) (9.27) which can be easily obtained from (9.23). As consistency check, notice that a(ω R ,ω R ) = 0, as expected. For large ω T , the leading term of (9.26) is We find it interesting to observe that, once the limit L → ∞ has been taken, eitherω R or ω T can be set to zero. For instance, settingω T = 0 in (9.26) gives the following finite result a(ω R ,ω T = 0) 2 = 1 16π On the other hand, it is well known that the correlators q iqj in (9.12) diverge when the frequency of the chain vanishes because of the occurrence of the zero mode; hence we cannot evaluate C 2 for a finite chain when eitherω T = 0 orω R = 0. This tells us that the limits L → ∞ andω T → 0 do not commute.
In Fig. 4 we show the complexity C 2 as function of the size L of the periodic chain. The numerical results discussed in this manuscript have been obtained for κ = 1 and m = 1, unless otherwise specified; henceω R = ω R andω T = ω T . In the left and right panels of Fig. 4 we have ω T > ω R and ω T < ω R respectively (notice that ω R = 1 for all the panels). In the top panels the numerical data are compared against the expression (9.25) (solid lines) obtained in the thermodynamic limit: while in the top left panel the agreement is very good at large L, from the top right panel we conclude that larger values of L are needed to observe a reasonable agreement asω T → 0. In the bottom panels of Fig. 4 we consider the subleading term in (9.25): while in the bottom left panel the data agree with the horizontal lines corresponding to a(ω T ,ω R ) given by (9.26), in the bottom right panel the agreement gets worse asω T → 0. Notice that the solid lines in the right panels of Fig. 4 accumulate on a limiting line as ω T → 0. This line can be found by plugging (9.29) with ω R = 1 into (9.25).
We find it worth considering a perturbative expansion of the complexity of these pure states when the target state is infinitesimally close to the reference state. This means that ω T =ω R + δω with |δω| ω R . Assumingω R = 0, we expand (9.26) as δω/ω R → 0. The first order of this expansion gives Including also the O((δω) 2 ) term in the expansion of (9.26), we find In Fig. 5 we compare the exact formula (9.22) for finite L against the first order result (9.30) (dashed lines) and against (9.31), that includes also the second order correction (solid lines).  Figure 5: The complexity C 2 between the ground states of harmonic chains with κ = m = 1 and different frequencies, for a given ω R and as function of ω T . The data points come from (9.22). The dashed lines correspond to the first order approximation (9.30) and the solid lines to the second order approximation (9.31), in the thermodynamic limit and in the expansion whereω T =ω R + δω.

Thermal states
The thermal states are the most natural mixed states to consider. In the following we evaluate the complexity (2.33) when both the reference and the target states are thermal states of the harmonic chain.

Covariance matrix
The two-point correlators of a periodic chain in a thermal state at temperature T read where Ω k is given by (9.9) and we have introduced The correlators (9.32) and (9.33) provide the generic elements of the correlation matrices Q th and P th respectively, which are the non vanishing blocks of the covariance matrix γ th = Q th ⊕P th of the thermal state.
Following the standard procedure, also for the thermal state one first performs the discrete Fourier transform through the matrix U in (9.5), finding that (9.13) can be written also for the correlation matrices Q th and P th , i.e.
By employing the results obtained in Sec. 9.2 for the covariance matrix of the ground state, it is not difficult to find that the Williamson's decomposition of the covariance matrix γ th of the thermal state reads The matrix W C is the symplectic matrix (9.17) occurring in the Williamson's decomposition (9.18) of the ground state and D th ≡ diag(σ th,1 , . . . , σ th,L ) ⊕ diag(σ th,1 , . . . , σ th,L ), with the symplectic eigenvalues given by [44] σ th,k = 1 2 coth Ω k /(2 T ) = 1 2 coth σ phys,k /(2T ) k = 1, . . . , L . (9.38) From these observations, we find that where V is the same symplectic and orthogonal matrix introduced through (9.6) and (9.7) for the ground state. It is straightforward to check that (9.39) becomes (9.20) as T → 0, as expected.

Complexity
In order to explore the complexity of two thermal states of a periodic chain, let us consider a reference state characterised by frequency ω R and temperature T R and a target state characterised by frequency ω T and temperature T T , assuming again that κ R = κ T = κ and m R = m T = m. Like (9.20), we have that also (9.21) can be adapted to this case, simply by replacing Q gs,M with Q th,M and P gs,M with P th,M taken from (9.36), with M ∈ {R, T}. Thus, the complexity given by (2.33) and (2.31) for these thermals states becomes where we have introduced being Ω S,k the dispersion relation in (9.23). Notice that Ω M,N,k → Ω M,k as T N → 0; hence in the limit given by T R → 0 and T T → 0, the expected expression (9.22) for pure states is recovered.
In the special case whereω R =ω T ≡ω, the expression (9.41) simplifies to This result is consistent with the general expression (2.49) for the complexity in the special case where W TR = 1. This relation can be verified by settingω R =ω T in (9.24).
Another interesting special case to explore is given by a pure reference state, i.e. T R = 0. In this limit (9.41) becomes (9.44) In the low temperature limit (i.e. when T T ω R ,ω T ), by using that coth x 1 + 2e −2x + . . . as x → ∞ and that log(1 + x) x + . . . as x → 0, we find that (9.44) simplifies to where the dots denote subleading terms. This expansion tells us that the first correction to the pure state result (9.22) as T T → 0 is exponentially small. The high temperature regime corresponds to T R ω R and T T ω T . In this limit, by using that coth x 1/x when x → 0, we find that (9.41) becomes (9.46) In the thermodynamic limit L → ∞, for the complexity (9.41) we find being Ω M,N,θ ≡ Ω M,θ coth Ω N,θ /(2 T N ) , with M, N ∈ {R , T} (see (9.42)), written in terms of the dispersion relation Ω S,θ given by (9.27). Notice that a(ω R ,ω T , T R , T T ) → a(ω R ,ω T ) when T R → 0 and T T → 0, where a(ω R ,ω T ) has been defined in (9.26).
In Sec. 9.2 we have observed that the massless limit of the coefficient of the leading term of the complexity of pure states is finite (see (9.29)). This happens also for thermal states. Indeed, by settingω R =ω T = 0 in (9.48), we find Figure 6: The complexity C 2 for thermal states as function of the size L of the periodic harmonic chain. Here κ = 1 and m = 1; henceω R = ω R ,ω T = ω T , T R = T R and T T = T T . In all the panels T R = 0 and various values of T T are considered. Top panels: ω R = ω T = 1 (left) and ω R = ω T = 10 −2 (right). The solid lines correspond to (9.47). Bottom panels: subleading term log C 2 − 1 2 log L as function of log L, for T R = 0 and various values of T T . In the left panel ω R = ω T = 1 and ω R = ω T = 10 −2 in the right panel. The horizontal solid lines correspond to the constant values obtained from (9.48).
In Fig. 6 and Fig. 7 we report some numerical results for the complexity (2.33) between thermal states with different temperatures and ω R = ω T ≡ ω. The data have been taken for κ = m = 1, hence T R = T R and T T = T T .
In Fig. 6 we consider the complexity as function of the length L of the periodic harmonic chain. For the sake of simplicity, the reference state is the ground state (i.e. T R = 0) and the target state is a thermal state with temperature T T . In the left panels ω = 1, while in the right panels ω = 10 −2 . In the top panels the data are compared against the expression (9.47) (solid lines), while in the bottom panels the subleading term of the same expression is considered (horizontal solid lines). The data having ω = 1 agree very well with the predictions, while for the ones with ω 1 the agreement is worse because in these cases the values of L considered are not large enough.
In Fig. 7 the same quantity considered in Fig. 6 is shown as function of T T . The increasing behaviour of the curves tells us that the distance between the states increases with T T , as  expected. In the left panel we test numerically the analytic expression (9.41) for L = 50 and different values of ω. Instead, in the right panel we test numerically the formula (9.48), obtained in the thermodynamic limit L → ∞: the agreement is very good when ω 1, while it gets worse when ω 1. Thus, when ω is very small, larger values of L should be explored to observe the expected agreement between the numerical data and the curve (9.48). When T R = 0, in the latter case the curves for (9.48) collapse onto the limiting curve (9.49), obtained by setting ω R = ω T = 0.

Optimal path for entanglement hamiltonians and its complexity
In Sec. 6 we have discussed the map that provides the entanglement hamiltonian in terms of the covariance matrix of a mixed state. In the following we explore further the optimal path of entanglement hamiltonian matrices for the periodic harmonic chain in the special case where both the reference and the target states are thermal states.
The entanglement hamiltonian matrices H R and H T of a reference state and of a target state that are both thermal can be obtained by applying the map (6.3) to the covariance matrix γ th = Q th ⊕ P th introduced in Sec. 9.3.1, whose Williamson's decomposition is (9.37). The symplectic spectrum of the entanglement hamiltonian matrix of a thermal state can be easily obtained by plugging (9.38) into (6.7), finding This provides the elements of the diagonal matrix E th entering in the Williamson's decomposition (6.6) for the thermal state. Comparing (9.50) with (9.34) and (9.38), we get ε th,k = βσ phys,k , as discussed in Sec. 6 for the thermal states in any number of dimensions.
The distance (6.11) between H R and H T can be evaluated by employing (9.24) and (9.50). The result reads which can be obtained also by replacing D th with E th in (9.43).
In the special case of ω R = ω T , the summation in (9.51) can be easily performed, finding which corresponds to (6.14) specified to the one-dimensional harmonic chain. In this case we can employ the discussion made in Sec. 6 for the cases where W TR = 1 (see (9.24)) to conclude that d(γ R , γ T ) d(H R , H T ).
In the periodic harmonic chain the Williamson's decomposition of the optimal circuit connecting H R and H T is given by (6.16), with the symplectic eigenvalues (9.10) and the symplectic matrix (9.11). Thus, the symplectic eigenvalues for the matrix labeled by s ∈ [0, 1] along this optimal circuit are where Ω k is the dispersion relation (9.9). This means that the optimal circuit is made by the entanglement hamiltonian matrices of thermal states, as also discussed in Sec. 6. This is not a feature of the optimal circuit connecting the covariance matrices of two thermal states, as discussed in Sec. 2.6. This discrepancy is consistent with the fact that the map (6.3) does not send geodesics into geodesics.

Mutual complexity of TFD's
The thermofield double state (TFD) is a pure state obtained by entangling two equal copies of the harmonic lattice and such that a thermal state of the original system is obtained after the partial trace over one of the two copies. A detailed analysis of the TFD and of the circuit complexity between two TFD's is reported in the Appendix F.
It is worth comparing the circuit complexity of two thermal states with the one obtained from the corresponding TFD's. Following [23], we introduce the mutual complexity for the TFD's as where C th and C TFD are given by (9.41) and (F.38) respectively. More explicitly, (9.54) reads log Ω R,k coth β T Ω T,k /4 Ω T,k coth β R Ω R,k /4 2 which can be written also in terms of Ω M,N,k defined in (9.42) as follows After expanding the squares and a bit of manipulation, one obtains For fixed k, the argument of the sum in (9.57) only depends onβ T Ω T,k andβ R Ω R,k and it is symmetric under the exchange of T and R; hence we can fixβ T Ω T,k >β R Ω R,k for every k without loss of generality. This allows to show that every term of the sum (9.57) is negative 20 and therefore M TFD ω R ,ω T ,β R ,β T is always negative.
The thermodynamic limit L → ∞ of (9.55) gives M TFD ω R ,ω T ,β R ,β T = a TFD ω R ,ω T ,β R ,β T L + . . . (9.59) where the coefficient of the linear divergence can be written in terms Ω S,θ in (9.27) as a TFD ω R ,ω T ,β R ,β T = (9.60) We remark that the massless limit of M TFD /L diverges when L < ∞, while it is finite once L → ∞ is considered. Indeed, by settingω R =ω T = 0 in (9.60) we find This feature has been observed also for the complexity of pure states (Sec. 9.2) and for the complexity of thermal states (Sec. 9.3).
In the limitβ R → ∞ both the TFD in (F.8) and the thermal reference state become the product of two ground states. In this regime, (9.61) slightly simplifies to which depends only onβ T and can be easily studied. This function is negative for every value ofβ T and it vanishes whenβ T → ∞, as expected. Whenβ T → 0 in (9.62), we find the following logarithmic divergence In Fig. 8 we compare the mutual complexity for the TFD in (9.55) with its thermodynamic limit in (9.59) for various values of the parameters. In the left panel we show M TFD /L (dashed lines) as function of ω T for fixed ω R = 1 and for two values of β R = β T ≡ β. As L increases, the dashed curves approach the solid curves representing a TFD given in (9.60).
When ω T → 0, M TFD /L evalueted for finite L diverges, while its thermodynamic limit is finite, as observed above. In the right panel we show M TFD /L as function of β T when β R → ∞ and ω R = ω T ≡ ω = 1. Remarkably, the curves obtained for finite number of sites coincide with their thermodynamic limit already for L = 5. In the same panel we plot a TFD ω R = 0, ω T = 0, β R → ∞, β T in (9.62) (red solid curve), checking also that its behaviour for β T 1 is well reproduced by (9.63) (green dot-dashed curve).

Reduced density matrices
Important mixed states to explore are the reduced density matrices of a subsystem A.
Consider the density matrixρ R andρ T of the reference and of the target states respectively and introduce a spatial bipartition A ∪ B of the system that induces a factorisation of the Hilbert space, as already discussed in Sec. 6. For the Gaussian states that we are interested in, let us denote by γ R,A and γ T,A the reduced covariance matrices corresponding to the subsystem A, that characterise the reduced density matricesρ R,A ≡ Tr BρR andρ T,A ≡ Tr BρT respectively. We remark that, whenever B = ∅, the reduced density matricesρ R,A andρ T,A are mixed states, even whenρ R andρ T are pure states. The reduced covariance matrix γ A is obtained by just restricting the indices of the covariance matrix of the entire system to the ranges identifying the subsystem A.
By applying (2.33) to these mixed states, one obtains the subregion complexity In the context of the gauge/gravity correspondence, the subregion complexity has been studied e.g. in [10,28,137,138]. . We also report a TFD in (9.60) for the same values of the parameters (solid curves). In the right panel the dependence on β T is investigated by plotting (9.55) for three different values of L (dashed curves) and (9.60) (black solid curve) both for β R → ∞ and ω R = ω T ≡ ω = 1. The massless limit in (9.62) is also reported (red solid curve) and its small β T behaviour in (9.63) is checked (green dot-dashed curve).
where the additive constant depends on ω T . Comparing the left panels and the right panels, we observe that larger values for are needed to reach the behaviour (9.72) for these choices of ω R > ω T . We checked numerically that, when ω R = 1, (9.72) holds with a subleading term that depends also on ω R .
In Fig. 10 we have reported the subregion complexity for a block made by consecutive sites in an infinite harmonic chain when the entire system is in a thermal state and ω R = ω T . In particular, we have chosen T R = 0 and various values of T T > 0. In the left panels we have considered ω T = ω R = 1, finding a reasonable agreement with (9.72), where the subleading constant term depends on T T . In the right panels we have fixed ω T = ω R = 10 −6 , finding that the behaviour (9.72) is more difficult to observe as T T → 0 because larger values for are needed.

Mutual complexity of reduced density matrices
The complexity of the ground states and of the thermal states, considered in Sec. 9.2 and Sec. 9.3 respectively, grow like √ L as L → ∞, being L the number of sites composing the entire periodic chain. Furthermore, considering an interval made by sites in an infinite harmonic chain, the numerical results reported in Sec. 9.5 tell us that the subregion complexity for this interval grows like √ as → ∞.
Given a spatial subregion A and the density matricesρ R andρ T , which can correspond to pure or mixed states, let us denote by C R,T (A) the subregion complexity between the reduced density matricesρ R,A andρ T,A introduced in Sec. 9.5.
The subregions A 1 and A 2 can be either disjoint or have a non vanishing intersection. Since C R,T (A) 2 grows with the volume of A as the number of sites in A diverges, we are naturally led to introduce the mutual complexity for subregions as follows [23,137] which is finite as the number of sites in A 1 and A 2 diverges.
In an infinite chain, let us consider the mutual complexity when A 1 and A 2 are two equal and disjoint intervals made by sites and separated by d sites. In Fig. 11 we report the numerical results of the mutual complexity for this configuration as function of , while d/ is kept fixed (d/ = 1/2 for the data in the figure).
In the top and middle panels of Fig. 11, the reference state and the target state are the ground states of the chains characterised by ω R and ω T respectively. The mutual complexity is shown as function of : for the data shown in each panel ω R is fixed and the different curves are associated to different values of ω T . When ω T < ω R (top panels) the numerical curves for small values of are increasing until they reach a maximum at a value of that depends on ω T (top left panel). After the maximum, the mutual complexity decreases with , but for many values of ω T we cannot appreciate the finite asymptotic limit as → ∞ because larger values of are needed. In the top right panel, for small enough values of ω T , the values of that we consider are too small to appreciate the occurrence of a maximum.
When ω T > ω R (middle panels) a similar behaviour is observed: also in these cases the position of the maximum of the curve depends on ω T . In these cases we observe that, as → ∞, the mutual complexity decreases until the zero value is reached. By comparing the two panels in the middle, one observes that the value of where the data vanish increases when ω R decreases. Furthermore, from the middle panels we can appreciate also the fact that the sign of M R,T is not definite: it is mainly positive, but for some values of the parameters (ω T close to ω R and sufficiently small) the curve is negative.
In the bottom panel of Fig. 11 the reference state is a ground state again, whileρ T is a + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +  Figure 11: The mutual complexity for two disjoint equal intervals made by sites and separated by d sites in an infinite harmonic chain, as function of . Since we set κ = m = 1, we haveω R = ω R ,ω T = ω T , T R = T R and T T = T T . In all the panels d/ = 0.5. In the top and middle panels the chain is in its ground state and ω R = ω T : both ω T < ω R (top panels) and ω T > ω R (middle panels) are considered. In the bottom panel a case involving thermal states with T R = 0 and various T T > 0 is considered.
thermal state at temperature T T > 0. The curves corresponding to different values of T T > 0 decrease with and the asymptotic value depends on T T .

A comparison with the approach based on the purification complexity
We find it worth comparing our results in Sec. 9.3 for the complexity of thermal states with the corresponding ones obtained in [23] through the approach based on the purification complexity, that has been discussed in Sec. 8.2.
The results from [23] that we consider have been obtained using the F 1 cost function (see Appendix G for details), while the complexity (2.33) is based on the F 2 cost function. These different cost functions lead to a different scaling with the total size L of the chain. In particular, the F 1 cost function provides a complexity that diverges with L, while the F 2 cost function leads to the milder divergence given by √ L. This feature, which has been observed already in [17] for pure states, holds also for thermal states, as remarked in [23] for the F 1 cost function and in (9.47) for our approach, that is based on the F 2 cost function. Because of this different scaling, a meaningful comparison between these two approaches can be done only for one-mode mixed states, where L = 1. When both the reference and the target states are pure and L = 1, both the F 1 cost function and the F 2 cost function provide the same complexity [17].
Let us consider the circuit made by one-mode mixed states where the reference state is the ground state with frequency ω R and the target state is a thermal state at inverse temperature β with frequency ω T .
Specialising the complexity in (9.44) to this case and for m = κ = 1 we obtain In this simple case, analytical results have been found in [23] also through the approach based on the purification complexity. The results for the C 1 complexity, which is defined through the F 1 cost function, are basis dependent. In [23] two particular bases have been considered, which have been called physical basis and diagonal basis (see Appendix G for their definitions). In the diagonal basis, the analytic result found in [23] for this case reads (9.75) For the physical basis analytic results are not available. In the regime βω T 1, the following perturbative expansion has been found [23] C phys The expressions for the complexity in (9.74), (9.75) and (9.76) depend only on βω T and on the ratio ω R /ω T . As consistency check, we notice that in the limit βω T → ∞, where the circuit is made by pure states, all the expressions in (9.74), (9.75) and (9.76) become 1 2 | log(ω R /ω T )|, as expected from [17].
In Fig. 12 we show the expressions for the complexity in (9.74), (9.75) and (9.76) in terms of βω T for a fixed value of βω R (we choose βω R = 1 in the left panel and βω R = 10 in the right panel). The curve for C phys 1 occurs only in the right panel because it exists only in the regime of βω T 1. We find it worth remarking that curves for C 2 always lie below the curves corresponding to the complexity evaluated through the purification complexity. Furthermore, as βω T growths, all the curves collapse on the same curve, as expected from the above observation, since βω T → ∞ corresponds to the limit where the circuit is made by pure states. In the right panel one also notices that the complexity in the diagonal basis is smaller than the complexity in the physical basis, as already remarked in [23].

Conclusions
In this manuscript we have studied the circuit complexity of the mixed bosonic Gaussian states occurring in the Hilbert space of harmonic lattices in any number of dimensions by employing the Fisher-Rao distance between Gaussian Wigner functions.
Considering mixed states with vanishing first moments, by applying a known result for the symmetric and positive definite matrices [59][60][61][62][63] to the covariance matrices of the model we provide the optimal circuit (2.28) when the set of allowed gates includes all the Gaussian states that can be constructed through a covariance matrix. The length (2.31) of this optimal circuit in the geometry determined by the Fisher information matrix is identified with the circuit complexity (2.33) to obtain a target state from a given reference state (the tolerance is zero for these circuits). In the special case of pure states, the known results of [17,22] for the C 2 complexity have been recovered. For thermal states originated from the same hamiltonian, the expression (2.80) has been obtained. The Williamson's decomposition of the covariance matrix (see Sec. 2.3) is the main tool employed throughout our analysis. This decomposition identifies the symplectic spectrum, that is invariant under changes of basis that preserve the canonical commutation relations. The role of the symplectic spectra and of the basis in the computation of the C 2 complexity is made manifest in the expression (2.49). Furthermore, the Williamson's decomposition leads to natural ways to introduce the spectrum complexity and the basis complexity for mixed bosonic Gaussian states (see Sec. 3.3 and Sec. 3.4). This provides an explicit realisation of the proposal made in [28].
The optimal circuits described in this manuscript allow us to study the purification of mixed states without ancillae. Motivated by the first law of complexity, in Sec. 4 we have mainly considered the purification of a given mixed state through the W path (4.5). Further future analyses could lead to find the optimal purification path.
The Gaussian mixed states that are not pure in harmonic lattices can be characterised also through their entanglement hamiltonian matrices. The optimal circuit and the corresponding complexity in terms of the entanglement hamiltonian matrices have been investigated in Sec. 6.
It is important to understand how to construct the optimal circuits. A preliminary analysis has been carried out in Sec. 7, where the possibility to express the optimal circuit in terms of Gaussian channels has been explored. We have not been able to find a general solution to this interesting problem, hence further future investigations are needed.
It is instructive to compare alternative quantitative approaches to the complexity of mixed states. The approach described in this manuscript holds only for the bosonic Gaussian states occurring in harmonic lattices and it provides computable expressions for a generic number of degrees of freedom. The method discussed in [23] (see Sec. 8), that is based on the introduction of ancillary degrees of freedom, can be formulated for every model but it leads to expressions that are more difficult to evaluate.
A detailed analysis has been carried out in the simplest case of the harmonic chain either on the circle or on the infinite line (see Sec. 9). Analytic or numerical results have been reported. For the thermal states we have explored the optimal path, the corresponding circuit complexity (see e.g. (9.41)) and the purification. Analytic and numerical results have been found for the mutual complexity of thermofield double states (see Sec. 9.4). Finally, for the mixed states given by reduced density matrices, we have studied the circuit complexity for an interval in the infinite line and the mutual complexity of two disjoint intervals (see Sec. 9.6). Interestingly, in Fig. 11 we observe that the mutual complexity of two disjoint intervals vanishes as the size of the two equal intervals increases, while the ratio d/ is kept fixed, being d the number of sites separating them.
Our analysis mainly focusses on bosonic Gaussian states with vanishing first moments. It is very interesting to explore also the complexity of mixed Gaussian states whose first moments are non vanishing. The expression (2.33) for the circuit complexity holds also when the reference state and the target state have the same first moments, that can be non vanishing [65,79,83]. In Sec. 2.7 we have provided results for the coherent states (pure states with non vanishing first moments) and the complexity (2.92) has been discussed [79,85], finding agreement with [18]. We emphasise that an explicit expression of the Fisher-Rao distance in the most general case of mixed states with non vanishing first moments is not available in the literature [79,83]. Upper and lower bounds for the complexity have been discussed in Sec. 5.
The circuit complexity of mixed states is a challenging task deserving many future studies.
The analysis reported in this manuscript in the simple setup of bosonic Gaussian states can be extended in various directions. For instance, it is a straightforward application to study the C 2 circuit complexity in harmonic lattices in the presence of boundary and of defects [24,36,39] or in time dependent scenarios [35,[154][155][156].
One of the main motivations of our work is to provide some tools to study complexity in quantum field theories. Evaluating complexity of mixed states in quantum field theories remains an important challenge. The complexity of pure states in quantum field theories has been explored in various studies [21,[29][30][31][32][33][34][35][36][37][38][39][40] and it would be instructive to extend these analyses to mixed states. The tools of Information Geometry, that we have largely employed in our analysis, could provide further tools to handle this interesting problem [157].
Let us remark that our analysis has been performed by assuming that all the states of the quantum circuits are Gaussian. It is important to go beyond this limitation by exploring the complexity of circuits involving mixed states that are not Gaussian.
Finally, we remind that the results reported in this manuscript have been obtained in the ideal situation where the maximal freedom is allowed in the choice of the gates. Typically, only a limited number of gates can be employed in the construction of quantum circuits. It is worth trying to adapt our analysis to more realistic cases by introducing a tolerance and various kinds of restrictions in the set of the allowed gates [4,5].
In this Appendix we briefly discuss two aspects of the Gaussian mixed states described in Sec. 2.1 in the Schrödinger representation. In Sec. A.1 we report the kernel ρ(q,q) = q|ρ |q of the density matrix corresponding to the Gaussian Wigner function (2.10). In Sec. A.2 we consider the spatial bipartition A∪B of a system in a pure state, writing the kernel ρ A (q A ,q A ) for the reduced density matrix of the spatial subsystem A in terms of the parameters occurring in the wave function of the pure state of the entire system.

A.1 Wigner-Weyl transform
The density matrixρ, that fully characterises a mixed state, is hermitean, positive definite and its trace is equal to one. Being a linear operator on the Hilbert space,ρ is completely determined by its kernel ρ(q,q) = q|ρ |q in the Schrödinger representation.
In this manuscript we are interested in the states whose kernel ρ(q,q) is Gaussian [69]. This means that where Θ are Φ are N × N complex matrices. Since the argument of the exponential in (A.1) must be invariant under transposition, we have Θ = Θ t and Φ = Φ † . This implies that (A.1) is fixed by choosing N (2N + 1) real parameters: N (N + 1) real parameters in Θ and N 2 real parameters in Φ. The normalisation condition R N ρ(q, q) dq = 1 for (A.1) gives which is well defined when Re(Θ) − Re(Φ) is strictly positive. We remark thatρ =ρ † is equivalent to ρ(q,q) * = ρ(q, q) [69]. This condition is satisfied by (A.1).
The Wigner-Weyl transform (also called Moyal transform) of ρ(q,q) is defined as By using (A.1) and (A.2) into (A.3) and performing a Gaussian integral 21 , one finds that the Wigner-Weyl transform of (A.1) is Gaussian as well. In particular, it reads where The following Gaussian integral has been employed, where dx ≡ n i=1 dxi Since Re(Θ) and Re(Φ) are symmetric, T and C are N × N real and symmetric matrices, while I is a generic N × N real matrix, being Im(Θ) symmetric and Im(Φ) antisymmetric; hence (A.5) is determined by N (2N + 1) real parameters, as expected. The complex matrices Θ and Φ can be written in terms of the real matrices T , C and I by inverting (A.6). The result is The expression (A.5) can be written in the form (2.10) with γ given by (8.2) with This is obtained by noticing that The matrices T , C and I can be expressed in terms of the blocks of γ in (8.2) by inverting (A.8). The result is Thus, in terms of the blocks of γ in (8.2), the complex matrices in (A.1) and (A.7) read These matrices are real when γ = Q ⊕ P in (8.2) is block diagonal.
We remark that the Wigner characteristic function in (2.3) is related to the kernel ρ(q,q) through the following relation Indeed, the Wigner function (2.5) is recovered by plugging (A.13) into (A.3).

A.2 Reduced density matrix
In the Schrödinger representation, the kernel ρ A (q A ,q A ) corresponding to the reduced density matrixρ A of the subsystem A of a bipartite harmonic lattice in a pure state can be computed as follows.
Considering the wavefunction (2.24) for the pure state of the entire system, the spatial bipartition A ∪ B of the harmonic lattice naturally leads to write the real and symmetric matrices E and F in (2.24) as the following block matrices (A.14) where .16) and The kernel ρ A (q A ,q A ) corresponding to the reduced density matrixρ A is obtained by tracing out the degrees of freedom corresponding to the part B of the harmonic lattice. By employing (A.15) and the Gaussian integral (A.4), one obtains Notice that Θ A is symmetric and Φ A is hermitean, as expected from the general expression in (A.1). We find it worth remarking that F B does not occur in (A. 19). The real and the imaginary parts of Θ A read respectively and for Φ A we have respectively Imposing that the trace of (A.19) is one leads to If Ω AB is left invertible i.e. the N B × N A matrix Ω −1 AB exists such that Ω −1 AB Ω AB = 1, we have that Ω † AB is right invertible with (Ω † AB ) −1 = (Ω −1 AB ) † . Assuming also that Φ A is invertible, we can isolate Ω A and E B in (A.19), finding In the special case F = 0, that has been considered e.g. in [23,132], the matrices in (A.16) are real. Furthermore, from (A.21) and (A.23), we have that also Θ A and Φ A are real.
In Sec. 8 we have discussed the purification of a mixed state with N sites through the introduction of an auxiliary lattice with N anc sites. The results reported in this appendix can be employed in Sec. 8 by setting N = N A and N anc = N B . In particular, in the simplest case, which is given by N A = N B = 1, the above counting tells us that we have three parameters to choose. This result has been found also in Sec. 8.1.1 by using (8.17).
The above results provide a lower bound for the number N B of ancillary degrees of freedom that must be introduced to purify a mixed state. A theorem of matrix algebra [158] guarantees that, given two matrices M and N , the rank of their product is bounded by rank ( N B , we can conclude that N B rank(Φ A ). In [23] this argument has been applied for real matrices.

B On the Fisher-Rao distance between Gaussian PDF's
In this appendix we report some known results about the Fisher-Rao distance between Gaussian probability distribution functions [57,64,65,[78][79][80][81][82][83][84] in order to apply them to the analysis of the complexity of mixed bosonic Gaussian states.
A Gaussian probability distribution function (PDF) in one variable (also called univariate PDF) reads where µ ∈ R is the mean and v > 0 is the standard deviation.
These Gaussian PDF's provide a manifold M 1 once the metric is introduced through the Fisher information matrix [57,64,78] where i, j ∈ {1, 2}. Plugging (B.1) into (B.2), one obtains the diagonal matrix diag(1/v 2 , 2/v 2 ), that provides the following infinitesimal distance [64,65] which characterises the hyperbolic upper half plane H 2 after the rescaling µ → √ 2 µ. Thus, by equipping the space of the univariate Gaussian PDF's parameterised by the pair (µ, v) with the metric characterised by the Fisher information matrix (B.2), the geodesics are either the lines with constant µ or the half-ellipses with eccentricity 1/ √ 2 ending on the axis v = 0.
We are interested in the manifold M n made by the Gaussian PDF's in n variables x t ∈ R n , which are (see (2.6) with n = 2N ) where µ t ∈ R n is the mean vector and Σ is a n × n positive definite symmetric matrix called covariance matrix. The parameter space for θ has n+n(n+1)/2 real dimensions: n parameters for µ and n(n + 1)/2 for Σ. In this space, it would be interesting to have a closed form for the Fisher-Rao distance that generalises (B.4) to n 1. Nonetheless, important explicit results have been obtained for some interesting submanifolds of M n .
In 1976, S. T. Jensen [65] found that the n(n + 1)/2 dimensional submanifold M µ 0 defined by the Gaussian PDF's with the same µ = µ 0 is totally geodesic 23 and that the Fisher-Rao distance in this case becomes where λ i are the eigenvalues of Σ −1/2 1 The distance (B.6) is employed throughout this manuscript to evaluate the complexity of mixed bosonic Gaussian states (see (2.31)).
Another interesting submanifold M Σ 0 to consider is given by the Gaussian PDF's with the same covariance matrix Σ = Σ 0 . The Fisher-Rao distance on this submanifold becomes the Mahalanobis distance [79,81] We remark that M Σ 0 is not a totally geodesic submanifold of M n [79,84].
It is worth considering also the submanifold M diag made by the Gaussian PDF's whose covariance matrix is diagonal, namely Σ = diag(v 2 1 , . . . , v 2 n ), with v i > 0. In this submanifold the infinitesimal distance becomes [84,160] The normalisation of (B.4) is different from the one used in [79] 23 A submanifold M ⊂ M is totally geodesics if any geodesics in M is also geodesics in M [159].
which suggests that it is convenient to arrange the parameters as θ = (θ n ), with θ (1) i ≡ (µ i , v i ) in this case. The infinitesimal distance (B.8) leads to write the distance between two PDF's in M diag in terms of (B.4) as follows From (B.8) one concludes that the geodesics in M diag are the curves θ(s) such that θ (1) i (s) is a geodesic in hyperbolic upper half plane equipped with the metric (B.3), for all 1 i n. Notice that we are not guaranteed that a geodesic in M diag is also a geodesic in M n because M diag is not a totally geodesic submanifold of M n . Instead, a totally geodesics submanifold of M n is M diag ⊂ M diag , which is made by the elements of M diag such that µ is a given eigenvector of Σ (see e.g. Proposition II.1 in [85] 24 ). For instance, the Gaussian PDF's whose covariance matrices are proportional to the identity are contained in M diag and in this case µ is the generic element of R n .
Consider a diagonal Σ and the eigenvector µ t = (µ, 0, . . . , 0) [85]. In this case the metric (B.8) becomes and the corresponding geodesics can be found as discussed above [84,85]. By specialising (B.9) to this case and employing (B.4), one obtains (B.11) Notice that, when Σ has a degenerate spectrum, its eigenvectors can have more than one non vanishing components.
The Mahalanobis distance (B.7) can be applied on the submanifold M Σ 0 , which is not totally geodesic. Very recently, a closed form for the distance d FR (θ 1 , θ 2 ) between PDF's in M n having the same covariance matrix Σ 0 has been found [83]. Since M Σ 0 is not a totally geodesic submanifold of M n , the Mahalanobis distance (B.7) does not necessarily correspond to the length of a geodesic connecting two PDF's with the same covariance matrix in M n . Instead, the distance d FR (θ 1 , θ 2 ) provides the length of the shortest path in M n between two PDF's with the same Σ 0 . Since we are not restricting to a submanifold of M n , this is the proper Fisher-Rao distance in M n between two PDF's with the same covariance matrix. Thus, given two Gaussians PDF's with the same covariance matrix Σ 0 , we have that Given two Gaussian PDF's in M n parametrised by θ 1 ≡ (µ 1 , Σ 0 ) and θ 2 ≡ (µ 2 , Σ 0 ), let us consider the orthogonal matrix Π such that Π(µ 2 − µ 1 ) = (|µ 2 − µ 1 |, 0, . . . , 0) ≡ |µ 2 − µ 1 |e 1 . Being Σ 0 symmetric and positive definite and Π orthogonal, also the matrix Π Σ 0 Π t is symmetric and positive definite, hence it can be decomposed as [83] where U is an upper triangular matrix with all the diagonal entries equal to one and S Σ 0 is a diagonal matrix with positive entries. The Fisher-Rao distance between θ 1 = (µ 1 , Σ 0 ) and θ 2 = (µ 2 , Σ 0 ) in M n found in [83] reads where d diag is defined in (B.11).
In order to construct the matrices Π and S Σ 0 , let us introduce the unit vector m ≡ (µ 2 − µ 1 )/|µ 2 − µ 1 |, observing that the orthogonal matrix Π satsfies Π m = e 1 . This matrix can be constructed by considering the basis of R n given by B = m, e 1 , . . . , e k−1 , e k+1 , . . . , e n , where m k = 0 is a non vanishing component of m and e i is the unit vector having only the i-th component equal to one. The standard Gram-Schmidt procedure [158] allows to construct an orthonormal basisB = m, u 1 , . . . , u n−1 from B. Then, the orthogonal matrix Π in (B.12) is the matrix whose columns are the vectors ofB.
The Cholesky decomposition [161] allows to write a symmetric and positive definite matrix M in a unique way as M = L c L t c , where L c is a lower triangular matrix. This result can be related to (B.12) by considering the matrix I having 1 on the antidiagonal and 0 elsewhere, which satisfies I = I t = I −1 . The matrix I Π Σ 0 Π t I is symmetric and positive definite, hence its Cholesky decomposition tells us that it can be written as I Π Σ 0 Π t I = L c L t c in term of a lower triangular matrix L c . This gives Π Σ 0 Π t = I L c L t c I = I L c I (I L c I) t . Being L c a lower triangular matrix, we have that U c ≡ I L c I is an upper triangular matrix and it satisfies Π Σ 0 Π t = U c U t c . (B.14) This decomposition agrees with (B.12), provided that U c = U S 1/2 Σ 0 . For any upper triangle matrix U , we have that 25 where U has 1 along the diagonal. Applying this to U c gives S 1/2 Σ 0 = diag(U c ) and U = U c . When n = 1, the distance d FR (θ 1 , θ 2 ) in (B.13) becomes d The above discussion can be employed to define the complexity for coherent states, which are pure states described by Gaussian Wigner functions with non vanishing first moments (see Sec. 2.1 and Sec. 2.7) [45]. Let us restrict to the coherent states with diagonal covariance matrices and first moments with a single non vanishing component. Since the coherent states are pure states, their covariance matrices are constrained by (2.23) [45,99]. Applying the constraints to (B.10), one obtains the metric (2.90). This metric and the distance (B.11) lead 25 Writing (B.15) in components we have U j,k = l U j,l δ l,k U k,k = U j,k U k,k . When j = k, the identity is verified because Uj,j = 1. When j > k, we have that U j,k = 0 implies U j,k = 0, being U k,k > 0 (which comes from the Cholesky decomposition).
to the expression (2.92) for the complexity for coherent states. This is consistent with the results found in [18], as discussed in Sec. 2.7 in a more detailed way. In Sec. 2.7 we have also exploited the distance (B.13) to compute the complexity (2.93) between two coherent states defined by (2.84) from the same ground state. These states have the same covariance matrix, but different first moments.

C Bures distance and Hilbert-Schmidt distance
In the literature of quantum information, different distances have been constructed for mixed states, even in the simple case of the bosonic Gaussian states. In this appendix we discuss the Bures distance and the Hilbert-Schmidt distance [87], that have been introduced in Sec. 2.2. In particular, we report their expressions in terms of the covariance matrices and then consider the special case of thermal states. An application of the Bures metric in the context of the complexity is discussed in [162].
The Bures distance between quantum states (defined in (2.16) from the fidelity) is Riemannian and contractive [88] (see Sec. 2.2). An explicit expression for the fidelity between two bosonic Gaussian states in terms of the corresponding covariance matrices γ 1 and γ 2 has been found in [113]. For vanishing first moments, it reads 26 where F tot is defined as (C. 2) The Bures distance in terms of the covariance matrices can be easily obtained by plugging (C.1) into (2.16). A canonical transformation characterised by the symplectic matrix S induces the change γ i → Sγ i S t for the covariance matrices γ i , with i = 1, 2. Simple matrix algebra based on the property of the symplectic matrices leads to conclude that the auxiliary covariance matrix γ aux in (C.2) changes as γ aux → Sγ aux S t and also that γ aux J → S(γ aux J)S −1 . Thus, both F 4 tot in (C.1) and F(γ 1 , γ 2 ) in (C.2) are left invariant by a canonical transformation. We refer to [113,163,164] for the Bures distance between two density matrices that are infinitesimally close.
Let us focus on γ 1 and γ 2 corresponding to thermal states having temperatures T i in harmonic chains with frequencies ω i , elastic constants κ i and masses m i , being i = 1, 2.
By using (9.39) and exploiting the fact that V depends only on the size of the chain, we can easily diagonalise γ 1 + γ 2 as follows In [113] the fidelity (C.1) between two Gaussian states with non vanishing first moments is also provided.
where the elements of Q i and P i with i = 1, 2 are defined in (9.36) and the matrix V has been introduced in Sec. 9.1. By employing (9.20) and the fact that V is orthogonal and symplectic, we observe that for γ aux in (C.2) we have Notice that M 1,2 is not diagonal, while M 1,2 J is diagonal. From (C.3) and (C.4) one realises that V cancels in (C.1) and (C.2), leaving the diagonal matrices Q i and P i . After some algebra, we find that the fidelity (C.1) for the thermal states γ i with i = 1, 2 becomes The Bures distance is easily obtained by substituting (C.6) into (2.16). The result reads As consistency check of this expression, we can consider the limit T i → 0, which provides the Bures distance between pure states. In this limit all the symplectic eigenvalues are 1 2 ; hence, from (C.7) we get B k = 1 4 . Then, the fidelity (C.6) simplifies to (Ω 1,k + Ω 2,k ) 2 4 Ω 1,k Ω 2,k −1/4 = L k=1 cosh 1 2 log Ω 2,k Ω 1,k −1/2 (C.9) and the Bures distance (C.8) becomes which is equal to the Fubini-Study distance between the two states, as expected.
The other distance that we consider is the Hilbert-Schmidt distance, which has been defined in (2.19) for two generic density operators. When the two density matrices are infinitesimally close to each other (i.e.ρ =ρ + dρ), this definition gives Focussing on Gaussian states, the Hilbert-Schmidt distance (2.19) between two mixed states can be written in terms of their covariance matrices as follows [165] d HS (γ 1 , γ 2 ) = 1 Since a canonical transformation characterised by the symplectic matrix S induces the transformation γ i → γ i = S γ i S t on the covariance matrices and det(S) = 1, it is straightforward to check that d HS is invariant under canonical transformations. The infinitesimal distance for (C.12) reads [165] The Hilbert-Schmidt distance (C.12) between the thermal states introduced in the text above (C.3) can be evaluated by employing (9.39) and (C.3), where the diagonal matrices are given by (9.36). Thus, for the determinants involved in (C.12), one finds where i = 1, 2 and Plugging (C.14) and (C. 15) into (C.12), in terms of the notation in (C.7) we get 27 d HS (γ 1 , γ 2 ) = (C.16) In the special case of pure states, all the symplectic eigenvalues are equal to 1 2 ; hence (C.16) simplifies to (Ω 1,k + Ω 2,k ) 2 4 Ω 1,k Ω 2,k It is worth comparing the Bures and the Hilbert-Schmidt distances in the case of pure states. From (C.10) and (C.17), one obtains The occurrence of this relation should be related to the fact that the Fubini-Study distance is the natural distance between pure states [87,88].

D Comments on some matrix identities
In this appendix we discuss some matrix identities employed throughout this manuscript.
In many matrix computations we have used the following property It is straightforward to prove these matrix identities when f (x) = x n and n is an integer number. Nonetheless, (D.1) has been often employed for f (x) = log x or for f (x) = x s with 0 s 1; hence in the following we show that (D.1) holds also for these functions.
The logarithm of a matrix M is defined through the series expansion [123] log This definition gives Since for the k-th term of this series we have Another remark deserving more detailed comments concerns (2.31), where we introduced the Fisher-Rao distance d(γ R , γ T ) as || log(γ where in the last equality the cyclic property of the trace has been used. On the other hand, which tells us that d(γ R , γ T ) cannot be written as || log(γ T γ −1 R )|| 2 . We find it worth providing further details about the construction of the symplectic matrices occurring in the Williamson's decompositions of H phys and of some covariance matrices in Sec. 9.

E Details on the first law of complexity
In this appendix we report some technical details concerning the first law of complexity and the derivations of the results reported in Sec. 3.
The variation of the square of the distance (2.31) under the independent variations γ T → γ T + δγ T and γ R → γ R + δγ R of the covariance matrices of the reference and of the target state reads δd 2 = 2 Tr log ∆ TR δ log ∆ TR = 2 Tr log ∆ TR ∆ −1 TR δ∆ TR . (E.1) The last expression can be found by first observing that, since M and δM do not commute in general, we can exploit the following formula [166] δ log M = Finally, δd 2 = 2d δd and (D.1) with f (x) = log x provide the expression (3.6).
In the following we compute separately the two sides of (3.4). Considering the r.h.s. of (3.4) first, from (3.3) one obtains ij ∂F ∂γ ij δγ ij = Tr γ −1γ γ −1 δγ From the expression (2.28) for the geodesic, it is not difficult to find that (E.10) For the subsequent discussion, let us remark that, by specifying (E.10) to s = 0 and s = 1, and using (D.1), we find The denominator in the r.h.s of (E.7) along the geodesic (2.28) reads Tr G −1 s ∂ s G s 2 (E.13) Combining (E.8) and (E.9), we observe that, for any 0 s 1, this expression is equal to d defined in (2.31). Furthermore, the numerator in the r.h.s. of (E.7) at the endpoints of the geodesics can be expressed by using (E.8), (E.9) and (E.10). Thus, from (E.7) we get ij ∂F ∂γ ij δγ ij (E.14) As for the l.h.s. of (3.4), let us consider d 2 from (2.31). First, one notices that (D.1) gives The expressions obtained by specialising this result to f (x) = log x and (log x) x −1 allow to write (E.6) as follows Then, by using (E.10) in (E. 16), we obtain that δd = (2d) −1 δd 2 can be written as while ∂F ∂γ can be read from (E.7). The expression (E.19) and ∂F ∂γ along the geodesics (2.28) can be found by using (E.8), (E.9) and (E.10). Then, some algebra leads to (E.18).
We find it instructive specialising the above results to pure states. Considering the geodesic given by (2.63), whose initial and final covariance matrices are given in (2.62), we have that both ∆ TR and ∆ TR + δ∆ TR are diagonal, hence they commute. This implies that (E.6) can be obtained directly from (E.1). Indeed, since in this basis ∆ TR = X 2 TR , we have that (E.1) becomes δd 2 = 2 Tr log X 2 TR X −2 TR This leads to write (E.6) for pure states as follows where we have also used that δγ T = 1 2 δX 2 TR in the basis that we are considering. Since δγ R = 0, one immediately realises that (E.23) corresponds to (E.17) for pure states.
Another natural value for s to choose in (E.21) is s = 1/2, that corresponds to the middle point of the geodesic. Comparing (E.21) with (E.20), it is natural to consider s a = −1, i.e. a = −2, finding that (E.20) can be written as We remark that X TR is the diagonal matrix providing the complexity in the case of pure states.
Another useful expression for δd comes from the Williamson's decompositions (2.41). Considering variations δγ such that γ + δγ is also a covariance matrix (in particular, δγ is symmetric). Given the Williamson's decomposition (2.20) for γ, let us express δγ in terms of the variations δD and δW of the symplectic spectrum and of the symplectic matrix W respectively. By using δW t = (δW ) t , Tr(M N t ) = Tr(M t N ) and the fact that ∂ s G −1 s is symmetric 28 , we can write (E.17) as in terms of the four contributions coming from the basis variations δW R and δW T and from the spectra variations δD R and δD T . The expression (3.7) for δd can be easily obtained by plugging (E.11) and (E.12) into (E.25). element ofb l andb r respectively. In terms of the components ofb d , the hamiltonian (F.1) becomes Ω k b † l,kb l,k +b † r,kb r,k + 1 . (F.7) The standard quantisation procedure leads to introduce the eigenstates |n r , n l ≡ |n l l |n r r of the number operator, that can be factorised through the eigenstates |n l l and |n r r of the number operators corresponding to the two parts. The eigenstates with n r = n l ≡ n allow to define the thermofield double state (TFD) as follows [167] |TFD = When β → ∞, the TFD becomes the product state of the two ground states |0 l |0 r .
Tracing out the degrees of freedom corresponding to one of the two parts, e.g. the right part, in (F.8) one obtains Tr Hr |TFD TFD| = N k=1 1 − e −βΩ k n e −β N k=1 Ω k n k |n l l n| (F.9) which is the thermal density matrix for the left half system at temperature 1/β.

F.1 Covariance matrix
The covariance matrix of the TFD can be found through a slight generalisation of the procedure described in Sec. 2.6 for the thermal states. From (2.9) and (F.6), the covariance matrix of this pure state can be written as Re q l,kql,k = Re q r,kqr,k = Re p l,kpl,k = Re p r,kpr,k = 1 2 coth(βΩ k /2) Re q r,kql,k = Re q l,kqr,k = −Re p r,kpl,k = −Re p l,kpr,k = 1 2 1 sinh(βΩ k /2) (F.11) where the notation O ≡ TFD| O |TFD has been adopted. By using (F.11), the covariance matrix Re TFD|ŝ dŝ Plugging (F.12) into (F.10) and employing the definition of W d in (F.5), for the covariance matrix of the TFD we find 16) where (F.3) and V d ≡ V ⊕ V have been employed to write the last expression, which is given in terms of the following 2N × 2N symmetric matrices as expected, being the TFD a pure state whose covariance matrix have non vanishing blocks only along the diagonal.
In order to write the Williamson's decompositions of γ TFD , let us observe that the matrices Υ where the 4N × 4N symplectic matrix W TFD is It is instructive to express the fact that the TFD is a particular purification of a thermal state (see (F.9)) by identifying it within the analysis reported in Sec. 8. This can be done by setting N anc = N and by rewriting the covariance matrix of the TFD in terms of the matrices occurring in (8.10).
Comparing (F.16) with (8.3), we easily conclude that in this case Q TFD and P TFD correspond to Q ext and P ext respectively, while M ext = 0. Then, by employing the block diagonal matrices (F.3), (F.12) and V = V ⊕ V , where V is the N × N orthogonal matrix (see (2.65) and (9.6)-(9.7) for the periodic harmonic chain), we can write Q TFD and P TFD as the partitioned matrices in (8.10) with and We remark that (F.25) and (F.26) satisfy the conditions in (8.20). Furthermore, Q ⊕ P = Q anc ⊕ P anc constructed from (F.25) provides the covariance matrix of a thermal state given in (2.76), as expected. Thus, the TFD is a purification of the thermal state and its covariance matrix satisfies (8.19).

F.2 Complexity
The TFD are pure states; hence the complexity of a target TFD with respect to a reference TFD can be computed by employing (2.58). In the most general case where the target TFD and the reference TFD originate from different hamiltonians, complicated expressions occur because W TFD depends on the physical hamiltonian through S and V in a non trivial way.
For the sake of simplicity, let us focus on the special case where the same hamiltonian underlies both the target TFD and the reference TFD, which are only distinguished by their inverse temperatures β R and β T . In this case both the reference state and the target state have the same S and V . Moreover, since (F.20) tells us that O does not contain parameters, the reference and target states that we are considering can be distinguished only through their matrices X TFD,R and X TFD,T . In this case, by employing (F.24) we find that the matrix defined in (2.45) crucially simplifies to the following diagonal matrix W TFD,TR = X TFD,T X −1 TFD,R ⊕ X −1 TFD,T X TFD,R . (F.27) The circuit complexity corresponding to this choice of TFD's can be obtained by plugging (F.27) into (2.58). The result reads Tr log X 2 TFD,T X −2 TFD,R 2 ⊕ log X −2 TFD,T X 2 TFD,R 2 (F. 28) which can be written more explicitly by employing (F.22) for these TFD's, finding log coth β T Ω k /4 coth β R Ω k /4 2 . (F.29) 1. In this limit the reference state is the product of the ground states of the two parts because only n = 0 contributes in (F.8). In this regime the complexity (F.29) simplifies to log coth β T Ω k /4 2 (F. 30) which is consistent with the results reported in [22].
We find it worth generalising (F.29) by considering a circuit where the reference state and the target state correspond to different hamiltonians that have the same matrix V d in their decompositions (F.2). This is the case e.g. for the periodic harmonic chain explored in Sec. 9.1, where V defined in (9.6) and (9.7) depends only on the number of sites of the chain, hence it is independent of the parameters occurring in the hamiltonian of the chain. From (F. 16 where Λ TR andΛ TR are N × N diagonal matrices whose entries read Λ TR k,k = Ω R,k Ω T,k coth β R Ω R,k /2 coth β T Ω T,k /2 − 1 sinh β R Ω R,k /2 sinh β T Ω T,k /2 Λ TR k,k = Ω R,k Ω T,k coth β R Ω R,k /2 sinh β T Ω T,k /2 − coth β T Ω T,k /2 sinh β R Ω R,k /2 1 k N . where the entries of the diagonal matrices within the square brackets are given by Λ TR +Λ TR k,k = Ω R,k coth β T Ω T,k /4 Ω T,k coth β R Ω R,k /4 Λ TR −Λ TR k,k = Ω R,k coth β R Ω R,k /4 Ω T,k coth β T Ω T,k /4 . (F.37) Plugging (F.36) into (F.32), we find that the orthogonal matrices V and O do not contribute to the complexity. By employing also (F.37), one obtains log Ω R,k coth β R Ω R,k /4 Ω T,k coth β T Ω T,k /4 2 + log Ω R,k coth β T Ω T,k /4 Ω T,k coth β R Ω R,k /4 2 . (F. 38) In the regime β R Ω R,k 1 this expression simplifies to log Ω T,k coth β T Ω T,k /4 Ω R,k 2 + log Ω R,k coth β T Ω T,k /4 Ω T,k 2 (F. 39) which is consistent with the results reported in [22] 29 .
G Diagonal and physical bases for the C 1 complexity In this appendix we briefly discuss the definition of the C 1 complexity, which is based on the F 1 cost function, hence it is a base dependent quantity. We also introduce the diagonal basis and the physical basis, slightly extending the definition given in [23]. Some results reported in Sec. 9.7 have been obtained by employing these bases.
In the Nielsen's geometric approach to complexity between pure states [1][2][3], the circuit connecting the reference and the target states is made by the unitary matrices U N (s), with s ∈ [0, 1], which are written as follows where ← − P is the path-ordered exponential indicating that the circuit is constructed from right to left as s increases, K I are the hermitian generators of the unitary transformation and the functions Y I (σ), that are called control functions, characterise the gates at a given point of the circuit. The circuit depth is defined through cost function F as follows The complexity corresponds to the minimal circuit depth, obtained by comparing all the possible unitary circuits connecting the reference state to the target state. The allowed cost functions must satisfy some properties that have been discussed e.g. in [17].
In this manuscript we consider only the F 1 cost function and the F 2 cost function. These cost functions are defined respectively as and, through (G.2), they provide the C 1 complexity and the C 2 complexity respectively.
Consider the harmonic lattice and the corresponding covariance matrix introduced in Sec. 2.
In [22], the complexity of pure states has been studied by employing the fact that, since a unitary circuit can be represented as a circuit in Sp (2N, R), instead of (G.1), for this model we can equivalently consider U (s) = P e where K I are the generators of Sp (2N, R), hence the index 1 I N (2N + 1).
The symplectic matrix U (s) in (G.4) has been discussed also in Sec. 2.5. In particular, from (2.61) and Eq. (57) in [22], we have being γ R and γ T the covariance matrices of the reference pure state and the target pure state respectively. In [22] the symplectic matrix U (s) has been written in terms of the matrix ∆ TR After a change of basis, the generators K I and the control functions in (G.4) change respectively as follows [17,22] where O IJ are the entries of an orthogonal N (2N + 1) × N (2N + 1) real matrix O. Being O orthogonal, in (G.7) the C 2 complexity is invariant while the C 1 complexity is not.
In [23] the complexity of mixed states based on the purification complexity (see Sec. 8) is mainly studied by employing the F 1 cost function. In particular, this C 1 complexity is investigated in two different bases: the diagonal basis and the physical basis.
The diagonal basis in the extended system, which is in a pure state, is defined by the change of basis corresponding to the symplectic and orthogonal matrix R introduced in (2.22).
In order to introduce the physical basis, let us consider the wave function (2.24) of the pure state characterising the extended system. This wave function is completely described by N ext × N ext complex symmetric matrix E ext + iF ext , that can be written as follows where E and F are N ×N real symmetric matrices, E anc and F anc are N anc ×N anc real symmetric matrices, while Γ E and Γ F are N × N anc real matrices.
In the physical basis both E + iF and E anc + iF anc are diagonal matrices. By employing a result of matrix algebra (see Corollary 4.4.4 of [168]), the complex and symmetric matrices E + iF and E anc + iF anc can be diagonalised as follows D = X(E + iF )X t D anc = X anc (E anc + iF anc )X t anc (G. 10) where D and D anc are real diagonal matrices with non negative entries and the matrices X and X anc are unitary. The physical basis is defined through the change of basis characterised by the matrix X phys ≡ X ⊕ X anc , that brings the blocks on the diagonal of E ext + iF ext in (G.9) in their diagonal forms. In the special case of F ext = 0, the definition of physical basis given in [23] is recovered.