Entanglement and Complexity of Purification in (1+1)-dimensional free Conformal Field Theories

Finding pure states in an enlarged Hilbert space that encode the mixed state of a quantum field theory as a partial trace is necessarily a challenging task. Nevertheless, such purifications play the key role in characterizing quantum information-theoretic properties of mixed states via entanglement and complexity of purifications. In this article, we analyze these quantities for two intervals in the vacuum of free bosonic and Ising conformal field theories using, for the first time, the~most general Gaussian purifications. We provide a comprehensive comparison with existing results and identify universal properties. We further discuss important subtleties in our setup: the massless limit of the free bosonic theory and the corresponding behaviour of the mutual information, as well as the Hilbert space structure under the Jordan-Wigner mapping in the spin chain model of the Ising conformal field theory.

Understanding quantum information theoretic properties of quantum field theories (QFTs) and, via holography, also of quantum gravity has been an enormously fruitful research front of the past two decades (as seen, for example, in [1][2][3][4][5]).
The main player in this endeavour has been the notion of entanglement and its entropy S. Starting with a pure state |Ψ and a subsystem A (its complement denoted by Ā), the entanglement entropy is defined as the von Neumann entropy of the reduced density matrix 1 ρ A = Tr Ā |Ψ Ψ| associated with A, specifically

S(A) ≡ − Tr
While entanglement entropy is very hard to calculate in a generic QFT, by now many results exist for free quantum fields, conformal field theories (primarily in two spatial dimensions) and strongly coupled QFTs with a holographic description.In the latter case, the entanglement entropy acquires a natural geometric description in terms of the Bekenstein-Hawking entropy of certain codimension-2 surfaces penetrating anti-de Sitter (AdS) geometries [11][12][13][14] and led to a wealth of results on quantum gravity in this setting.
Complexity is another quantum information-theoretic notion that made its appearance in the context of QFTs only recently and is directly motivated by holography.To this end, it was observed in [15][16][17][18][19][20] that codimension-one boundary-anchored maximal volumes and codimensionzero boundary-anchored causal diamonds have properties expected from the hardness of preparing states using tensor networks [21] in chaotic quantum many-body systems.
Subsequent articles starting with [22,23] saw in this conjecture a strong motivation to define the notion of complexity in the realm of QFTs in a similar spirit in which pioneering works [24,25] introduced the notion of entanglement entropy in the same context.The articles [22,23] were largely inspired by the continuous tensor network of cMERA [26] and viewed preparation of a pure target state |ψ T in QFT as a unitary transformation from some pure reference state |ψ R where the unitary U is obtained as a sequence of layers constructed by exponentiation of more elementary Hermitian operators O I U = Pe −i 1 0 dτ Following the approach of [27], which was originally devised to bound complexity of discrete quantum circuits, one can associate the cost of invocations of different gates generated by O I as related to the infinitesimal parameter Y I (τ ) dτ in the exponent.Translating this literally into a mathematical formula would lead to which is an integral over the circuit of a L 1 norm of a formal vector formed from the parameters Y I .Complexity C arises then as the minimum of (4) subject to the condition (2) As anticipated already in [23], cost functions based on L 1 norms such as (5) lead to challenging minimization problems.In the present work our rigorous results on complexity will be based on a particular choice of a L2 norm where, following [28,29], η IJ is going to be a particular non-negative definite constant matrix and O I are going to be normalized accordingly.This preferred by us form of η IJ is naturally induced from the reference state |ψ R (bosonic systems) or from Lie algebra (fermionic systems).
The essence of recent progress on defining complexity in QFT using broadly defined approach of [27] lies in making educated choices for O I and |ψ R , which allow one to perform minimization encapsulated by Eq. (5).In the vast majority of cases, it was achieved by focusing on free QFTs and utilizing powerful toolkit of Gaussian states and transformations [30].
The discussion so far concerned pure states, i.e., von Neumann entropy as an entanglement measure between a subregion and a complement in pure states and complexity as a way of quantifying hardness of preparation of pure states.Much less understood in the QFT context are quantum information properties of mixed states and the present paper concerns precisely this important subject.Of interest to us will be entanglement of purification (EoP) [31] and complexity of purification (CoP) [32].We will introduce these quantities in more detail in, respectively, sections V and VI.Here want to stress instead that the key motivating feature behind our work stems from both of these quantities involving in their definition scanning over purifications of mixed many-body states 2 .
Such purifications, i.e., embedding a mixed state in an enlarged Hilbert space in which it arises as a reduced density matrix, in the context we are interested in, i.e., QFT physics, are clearly challenging to operate with.Earlier works on EoP and CoP in high-energy physics include respectively [33,34] and [35] 3 and focus on free QFTs in which mixed states of interest, such as vacuum reduced density matrices or thermal states, are Gaussian.Gaussian mixed states allow for purification to pure Gaussian states, which underlied strategies employed in the aforementioned references.However, even purifications within the Gaussian manifold of states for large subsystems can be challenging to operate with and the above works made additional choices in this respect.This is where the key novelty of our present work appears, which is to consider the most general Gaussian purifications.To this end, we will consider free QFTs on a lattice and, whenever possible, encode reduced density matrices in terms of corresponding quadratic correlations represented by covariance matrices.Considering the most general Gaussian purifications amounts then to embedding mixed state covariance matrices as parts of larger covariance matrices corresponding to pure states.Utilizing efficient Gaussian techniques allows us to minimize the two quantities of interest, EoP and CoP for a judicious choice of a definition of pure state complexity [29], over general purifications to a given number of bosonic or fermionic modes.
Our primary focus is on a particularly simple yet telling setup of two-intervals vacuum reduced density matrices in free QFTs with a vanishing or very small mass.In quantum information context, such setup arose first in studies of mutual information (MI) defined using two subsystems A and B as in (1+1)-dimensional conformal field theories (CFTs), where A and B are two disjoint or adjacent intervals on a flat spatial slice as depicted by figure 1. MI will play an important role in our studies providing us with a guidance regarding both the behaviour of EoP, as in [33], as well as will help us to understand subtleties underlying our models.Our studies will mostly concern scaling of MI, EoP and MI with control parameters such as interval size, separation and, when present, system size and the mass.While EoP turns out be such a ultraviolet finite quantity by itself, for CoP we will consider a combination of single and two interval CoP results akin to (7) for which the leading ultraviolet divergences cancel and only milder divergences remain.Our paper is structured as follows.In section II, we review the two models we consider, the Klein-Gordon field in the massless limit and the critical transverse field Ising model, on a lattice paying a particular attention to description of their ground states in terms of covariance 3 One should also mention in this context [36], which, motivated by holographic complexity proposals, explored properties of CoP in the setting of a single harmonic oscillator.The subsystem that defines reduced density matrices for our discretized bosonic and fermionic models in their vacuum state consists of two intervals of a width of wA/δ and wB/δ sites and separated by a distance of d/δ sites, where δ is the lattice spacing.When d = 0, we will keep wA and wB generic.When d > 0, we will set for simplicity wA = wB ≡ w and the natural continuum combination is w/d.We will see that numerically determined MI and EoP approach in the continuum limit functions of w/d.With CoP the situation is more complicated, as it turns out to be ultraviolet divergent and brings in an additional dimensionful scale through the class of reference states of interest |ψR .
matrices.In section III, we benchmark our abilities to reach continuum limit in lattice calculations by comparing the results of our numerics with existing analytic formulas for MI in the aforementioned two interval case.In section IV, we discuss briefly the mathematics of purifications of Gaussian states as seen by covariance matrices, which is the working horse behind most of the results reported in the present article.Subsequently, we use this machinery to study EoP and CoP in the two-interval case of figure 1, respectively, in sections V and VI.In section VII, we comment on two subtleties relevant for our model, namely the zero mode when taking the massless limit for a bosonic theory and the different notions of locality in the spin vs. fermion picture of the Ising model.We summarize our results and present an outlook in section VIII.We also provide an extensive appendix that provides further details regarding our methods.

II. SETUP
In the present work we focus on two paradigmatic models: the Klein-Gordon field in the massless limit and the critical transverse field Ising model in 1+1 dimensions.For our numerical calculations, we discretize both theories either on a lattice with N sites and periodic boundary conditions(i.e., we identify the sites N + j ≡ j) or on an infinite lattice.Both theories describe CFTs in the respective limits with central charge c = 1 (Klein-Gordon) and c = 1 2 (Ising model).We will review the Hamiltonians of both models and their ground states with a particular focus on the covariance matrix formulation.The latter for free bosons will allow for an efficient calculation of EoP and CoP using Gaussian techniques.For the Ising model, we will discuss in detail how there are two distinct notions of locality associated to the spin and fermion formulation, respectively.

Analytical predictions
Gaussian numerics (fermions) Hologr.CFTs: subregion-CA Ising CFT (c = 1 2 ) 0.103 w δ + 0.0544 log w δ + 0.0894  [35,[42][43][44] and our numerical CoP results in free CFTs.The latter are obtained using Gaussian bosonic and fermionic states, the L 2 norm circuit complexity encapsulated by (59) and spatially disentangled reference states.The mutual complexity is defined differently for the holographic complexity proposals, see (49), and for our implementation of CoP, see (64).In the case of the CA proposal, L is the AdS curvature radius and CT is an arbitrary length scale arising from counter-terms [45,46].In the case of the bosonic calculation, µ is the reference state scale and functions f0, f1 and f2 are defined in (63).

A. Klein-Gordon field
We consider the well-known Klein-Gordon scalar field with a mass m that we will later take to zero.Its discretized Hamiltonian on a lattice with N sites is where δ represents the lattice spacing.We thus have a circumference We define canonical variables where a = 1, 2. It is well-known that the Hamiltonian can be diagonalized via Fourier transformations leading to N decoupled harmonic oscillators with frequencies The ground state |0 is Gaussian and fully characterized by its covariance matrix where a and b label the entries of the matrix.Continuum limit on a circle requires taking N → ∞ keeping product of meaningful continuum quantities m L = m δ N fixed.Each value of this product corresponds to a different QFT as a continuum limit within the class of free Klein-Gordon theories.Furthermore, when considering subsystems, as depicted in figure 1, continuum limit requires keeping ratios of w δ and d δ to L fixed as N → ∞.In practice, one takes N to be large but finite and requires that as N is increased well-defined quantities, for example the MI (7), stop changing significantly with N and stabilize in the vicinity of their QFT values.When w δ L 1, d δ L 1, then the results of the calculations should be effectively indistinguishable from the situation in which the spatial direction is a line.The mass m 1 δ becomes then the only dimensionful parameter in the continuum theory.Also, in this case the number k associated with discrete momenta in (11) gets incorporated into a continuum variable and a sum in (12) needs to be replaced by an appropriate integral, see for example [29].
We are particularly interested in the massless limit m → 0, where the Klein-Gordon field describes the CFT with central charge c = 1.More precisely, the c = 1 CFT with the periodic boundary conditions we imposed should be regarded as a 1-parameter family of theories arising in the path integral language the compactification of the bosonic field ϕ (i.e., periodically identified): The dimensionless parameter R is the radius of compactification in the field space and plays the role of a moduli specifying a particular c = 1 CFT.The scaling dimension of the lightest operator is then given by The above formula is a hint of an underlying duality between theories with field complactification radia of R and 2 R [47].The massless limit of (8) corresponds to the decompactification limit of compact free boson CFTs (R → ∞), which is a subtle limit since in light of (14) the gap in the operator spectrum approaches 0. While this limit leads to correct correlation functions of vertex operators or a single interval entanglement entropy, for other quantities the situation is more complicated.In particular, the modular invariant thermal partition function for the free boson reads [47] whereas the free massive boson calculation for m L 1 and upon keeping the regularized zero point energy is In both expressions η is the Dedekind eta function defined as The mismatch between the two calculations can be understood using the representation of the partition function on a circle as an Euclidean path integral on a torus.
In the case of ( 15), the zero mode contribution is neglected, as its inclusion would lead to an infinite volume term coming from the integration over the field space.
In the case of ( 15), the zero mode φ contribution to the path integral is included and is finite, as it originates in the path integral language from where the product β L is the torus spatial volume.Multiplying the partition function (15) by the factor (18) leads to (16), which explains the relation between the two partition functions.We will come back to these calculations in section VII A, where we discuss the influence of the zero mode on MIdecay with separation between two intervals.
In our studies, we will be using the free massive boson setup to extract the properties of the modular invariant c = 1 free boson CFT in the decompactification limit R → ∞.From this perspective, the partition function of interest, i.e., (15), can be indeed recovered from the massive boson Gaussian calculation ( 16) by dividing it by the known zero mode contribution (18).However, in the case of other quantities calculated using Gaussian techniques at non-vanishing mass the effect of the zero mode is not straightforward to isolate.As we already mentioned, numerical studies showed that Gaussian calculations with a small mass reproduce the universal entanglement entropy result for a single interval [1].Furthermore, one may expect the two interval case at small separations to be reliably described by the massive free boson calculation, as the zero mode affects primarily the long distance physics.As a result, these will be the mixed state setups that we will consider in our EoP and CoP explorations.On the other hand, the two interval case at large separations is tricky and we will return to it in the case of MI in section VII A.
Another subtlety that originates in the massless limit is that the ground state is only defined distributionally.The issue is best understood by diagonalizing the Hamiltonian (8) by transforming to momentum modes πk where we find N decoupled harmonic oscillators.For the oscillator with k = 0 (zero momentum mode), we have ω 0 = m, which vanishes in the massless limit.Consequently, the ground state of this mode approaches a delta distribution, which does not lie in Hilbert space.This leads to the divergence of certain terms in the covariance matrix (12).However, we are still able to define expectation values of observables and entanglement measures, such as the entanglement entropy, by computing those quantities for finite m and generating numerically results for values of m gradually approaching 0. In section VII A, we will discuss the role of the zero mode for such calculations in more detail using MI as an example.

B. Critical transverse field Ising model
We consider the transverse field Ising model in the critical limit J = h, where Ŝx,z i are spin-1 2 operators in the direction x or z at position i in the chain, i.e., related Ŝx,z i = 1 2 σ x,z i to the well-known Pauli matrices.The system consists of N spin- 1  2 degrees of freedom arranged in a circle, i.e., we assume periodic boundary conditions with N + i ≡ i.
The transverse field Ising model can be solved analytically by employing the Jordan-Wigner transform [48], i.e., eigenvalues and eigenvectors of the Hamiltonian Ĥ can be constructed in closed form.The transformation is based on introducing fermionic creation and annihilation operators f † i and fi .For the transformation, we write S ± i = S x i ± S y i as which leads to the almost quadratic Hamiltonian where h.c. stands form Hermitian conjugation and P = exp (iπ f † j fj ) is the parity operator.In this picture, the operators f † i and fi are fermionic creation and annihilation operators, but with a different notion of locality than the spin operators appearing in (20).From the Jordan-Wigner transformation (21), it is clear that the fermionic operator on site i are local to the whole region from site 1 to i in the fermionic picture, and vice versa.This ensures that bipartite entanglement of a connected region of sites is equivalent in the spin and the fermionic picture, because we can use translational invariance to identify this region with the sites {1, . . ., w}, for which the spin and the fermionic pictures are isomorphic, i.e., the density operators are unitarily equivalent leading to the same entanglement entropy.
It is well-known [49] that the fermionic Hamiltonian (22) can be written as a sum of two quadratic Hamiltonians Ĥ± of the form where P ± represent orthogonal projectors onto the Hilbert subspaces H ± of even and odd number of excitations, respectively, i.e., with f † i fi |n 1 , . . ., n n = n i |n 1 , . . ., n n .We can equivalently describe the Hamiltonian in terms of Majorona modes ξa i ≡ (q i , pi ) with which leads to the Hamiltonian given by If we define Majorana operators we can write the critical (h = J) model as This demonstrates explicitly that the two quadratic Hamiltonians Ĥ± in the fermionic picture only differ by the type of boundary conditions in the two sectors, namely anti-periodic boundary conditions for Ĥ+ and periodic boundary conditions for Ĥ− .We are particularly interested in the ground state |0 of the critical model, which is completely characterized by its covariance matrix where the canonical variables ξa i are now defined using (27) as κ = π N (2k + 1) with k = 0, . . ., N − 1 and the entries c κ (j) are given by for N being even 4 .An important subtlety arises if we compute bipartite entanglement of disconnected regions, because in this case the entanglement entropies are different for the spin vs. fermion picture.This subtle fact has been recognized numerous times in the literature [50][51][52] and plays an important role when relating the lattice model with the continuum CFT.The key observation is that the canonical anticommutation relations induce a different notion of tensor product and partial trace for fermions [53].Interestingly, this different notion only affects the bipartite entanglement entropy of disjoint regions, i.e., the reduced state in a subsystem consisting of two non-adjacent intervals on the circle will be different if we compute it using the spin vs. fermion picture.We comment on this in section VII B and review the respective literature in appendix B. In practical terms, this fact will lead us to apply our Gaussian numerics based on purifications only to the case when the two intervals are adjacent, i.e., d = 0 in figure 1.
Finally, note that the c = 1 free Dirac fermion CFT can be obtained from two copies of Ising model (i.e., Majorana fermion CFT) by imposing a different GSO projection [47].As a result, the discussion in the present section about spatial locality and Gaussianity applies also to the free Dirac fermion CFT.In effect, our fermionic Gaussian methods reproduce the properties of the free Dirac fermion CFT only for a single subregion or adjacent subregions and the answers in these cases are simply given by twice the answers for the corresponding Ising calculation.It is well-known that the Dirac fermion CFT is equivalent to the free compactified boson CFT at the compactification radius R = 1 (or equally R = 2) via the bosonization procedure [54,55].Let us also emphasize on this occasion that modular invariance is a property that is not always imposed in free fermion calculation available in the literature (see [56] for a discussion of the modular invariance in the context of entanglement entropy in CFTs).This sometimes leads to apparent tensions between CFT expectations and free fermion results.We will come back to this point in the next section.

III. MUTUAL INFORMATION
MI defined in (7) provides an important correlation measure between two subsystems A and B and below we summarize some of its properties.One reason to do this is to test our ability to reproduce them using our numerics before we apply it to a much less understood case of EoP and CoP.Another is to explore what kind of behaviour to expect from EoP and CoP.
MI is generically a non-universal quantity in CFTs, as it is related to a four-point function of twist operators and the latter is spectrum-dependent [39].
At large distances between the intervals, i.e. for d w in the notation of figure 1, the operator product expansion analysis predicts the following behaviour of MI where O min is the operator with lowest (but non-zero) conformal dimension and ∆ min = h min + hmin [40,41].At short separations, d w, one expects the following universal result [39,57] For d = 0, one can use a universal, i.e., only c-dependent, formula for a single interval entanglement entropy in the vacuum to arrive at a variant of (33) with d = δ.In table I we provided the form of (33) when the two intervals have arbitrary lengths.
Moving on to the two models we consider, for the free massless scalar QFT one has continuous and gapless spectrum of primary operators.As a result, the formula (32) does not apply as such and in section VII A we comment on a possible generalization.However, for a compactified free boson CFT, see (13), the spectrum of operators develops a gap (14) and calculations of entanglement entropy in [58] do reproduce this behaviour.
For c = 1/2 Ising CFT, in the limit d w, we expect I(A : B) ∼ (w/d) 1/2 , because O min is the spin operator σ ∼ Ŝx i , see (20), which has h = h = 1 16 , i.e., ∆ min = 1 8 .Furthermore, as the free Dirac fermion CFT is dual to the free compact boson theory at R = 1 (or, equivalently, R = 2), according to (14) one expect MI to decay as 1/d governed by the h = h = 1  8 operator.Such an operator would be natural to interpret as a product of two σ operators from each underlying Majorana fermion model.
Let us also note that existence of the following formula for MI for Dirac fermions [59] While this formula agrees at short distances with (33), at large distances it falls off as 1/d 2 rather than the aforementioned 1/d predicted by bosonization.This is related to the fact that the calculation in [59] utilizes torus partition function with the anti-periodic boundary condition for fermions.However, the modular invariant partition function leading to (32) includes also contributions from sectors in which fermions satisfy period boundary conditions.This is directly related to the discussion about reduced density matrices in the fermionic formulation of the Ising model mentioned in section II B and expanded later in section VII B and appendix B.
Having discussed analytic expectations, let us show how our bosonic and fermionic Gaussian states numerics reproduces it.This should be regarded as a cross-check of both our numerical lattice setup and its ability to reproduce features of the continuum limit.Also, it will illustrate to what extent considering a decompactified free boson with the zero mode regulated via non-vanishing mass capture long (d w) and short distance (d w) CFT expectations.
First, we consider the behavior of MI for a free bosonic field with central charge c = 1, shown in the first row of figure 2. As anticipated in section II A, employing Gaussian methods restricts us to the decompactified free boson with non-vanishing mass.
In the limit of small d/w, we find a logarithmic dependence similar to (33), specifically Note the logarithmic Lm dependence already observed in [34].The coefficient a 1 , expected to be c/3 ≡ 1/3 in the continuum limit, converges very slowly as the block widths w and the size L of the periodic system are increased.However, we can bound it by considering the behavior of I(A : B) and S(A ∪ B) separately.Estimating a 1 by a discrete derivative with respect to log(w/d), we find that this estimate approaches 1/3 from above for the I(A : B) data and from below for S(A ∪ B), as shown in the top-left corner of figure 2, suggesting an asymptotic ∝ 1/3 log w behavior identical to that of S(A).Our data, extending until N = 32000 and w = 2000 δ, yields a bound 0.27 a 1 0.40, consistent with our expectations.All shown data uses a distance d δ = 1 of one lattice site, as lattice effects on the value of the a 1 estimate were found to be negligible.
The behaviour at large d/w can be anticipated to be subtle in light of (32) as in the present case ∆ min → 0. In particular, assuming a power law behavior one would expect b 2 to vanish.The power b 2 can be estimated by discretizing the derivatives in the expression Indeed, we find its estimated value to gradually decrease at large d/w as N is increased, though MI can be well approximated by a power law with b 2 ≈ 0.2 in the range 1 < d/w < 100, consistent with earlier numerical studies in the d/w < 50 range [34].In the d/w → ∞ limit, we can bound b 2 0.15.The coefficient of a potential logarithmic growth I(A : B) ∼ b 1 log(w/d) in this limit can be bounded as b 1 0.06.While these results are obtained on a circle and extrapolated to a line, in section VII A we will discuss the large-distance behavior of free boson directly on a line, where we will also consider possible sub-logarithmic decay functions at large d/w.Such functional dependencies, resulting from subtle large distance behavior for free, nearly massless bosons, will reappear in the context of EoP studies in section V.
As we can only study the Ising CFT via a Gaussian fermionic model under the Jordan-Wigner transformation for adjacent intervals, MI computed in this approach is only relevant for the d = 0 case, which which follows directly from the entanglement entropy formula for a single interval.These formulas are included in table I for completeness.

IV. GAUSSIAN PURIFICATIONS
As discussed in section II, we focus on free theories as their ground states are Gaussian states, so we have a powerful machinery at our disposal to analytically compute the entanglement entropy and other quantities analytically from the covariance matrix of pure Gaussian states (see appendix C).Similar analytical formulas exist also for the L 2 circuit complexity of interest.The primary goal of this paper is to use this machinery to define and compute similar quantities, such as entanglement entropy and complexity, for mixed states leading to the notions of EoP [60] and CoP [32].They are defined in the following way: 1. We start with a function f (|ψ ) that is defined for arbitrary pure states |ψ .
2. For a mixed state ρ A ∈ H A , we construct the purification |ψ on a larger Hilbert space Of course, there is quite some freedom of how large the purifying Hilbert space H A can be.
3. The purification |ψ is not unique, but if we have found one purification |ψ , we can construct any other purification by acting with a unitary U = 1 A ⊗ U A , where U A is an arbitrary unitary on the purifying Hilbert space H A .
4. We then define a new function F (ρ A ) for the mixed state to be given by i.e., we minimize the original quantity f (|ψ ) defined for pure states over all purifications of the mixed state ρ A .
Note that there are some subtleties related to the fact that the purifying Hilbert space H A may not have a direct physical interpretation, e.g., if H A represents a local subsystem (region in real space) of a QFT, it is a priori not clear what the physical meaning of H A is. Consequently, the function f needs to be defined in an appropriate way that it can be meaningfully applied to arbitrary extended Hilbert spaces H = H A ⊗H A .While this is relatively straightforward for entanglement entropy, one needs to be a bit careful about circuit complexity that is usually defined with respect to a reference state that is chosen as spatially unenentagled with respect to a physical notion of locality.As explained in (51), one can show that this can also be done for the purifying Hilbert space H A in such a way that the resulting CoP is actually independent of the notion of locality or, put differently, the outcome of the minimization procedure can even be understood as equipping H A with a notion of locality.
problem for practical applications.The reason is that the required minimization procedure must be generally performed numerically, while the dimension of the respective manifold over which one needs to optimize grows quickly with the number of degrees of freedom.Therefore, EoP has been only studied for small systems and often only with respect to certain subfamilies of states, while CoP was exclusively studied via purifying individual degrees of freedom [35] rather than directly larger subsystems.
The key ingredient that enables the progress of the present paper is that we can efficiently compute EoP and CoP for the family of Gaussian states.For this, we start with a Gaussian mixed state ρ G A and compute a Gaussian purification |ψ G .When performing our minimization algorithm, we only sample over Gaussian states, i.e., we define the new function Clearly, we must have ) is an upper bound for the true minimum.Moreover, it is reasonable to assume that for many quantities, such as EoP and CoP, we actually have the equality . This was already conjectured in [33] and is further supported by [61].In the case of CoP, there is still limited progress in even defining circuit complexity for non-Gaussian states, which means that it is natural to only consider F G (ρ G A ) to start with.In both cases, it is therefore a meaningful restriction to only consider Gaussian purifications of Gaussian states.
For Gaussian states, we can use the covariance matrix and linear complex structure formalism as explained in appendix C. Rather than working with Hilbert space vectors, which would live in an infinite dimensional Hilbert space for bosons and a 2 N A -dimensional Hilbert space for fermions, we can fully characterize the Gaussian state by a 2N A -by-2N A dimensional matrix, where N A represents the number of bosonic or fermionic degrees of freedom.We restrict to Gaussian states with z a = tr( ξa ρ G A ), for which all relevant information is encoded in the so called restricted complex structure J A defined in (C7).For a mixed Gaussian state J A has purely imaginary eigenvalues ±ic i , where c i ∈ [1, ∞) for bosons and c i ∈ [0, 1] for fermions.Only if all c i = 1, we have a pure state.For every mixed state ρ G We can always purify such a state using a Hilbert space H A with the same number of degrees of freedom as H A , i.e., N A = N A .Then, there always exists a basis in the system A , such that the complex structure J of the purified state |ψ G takes the form [62] where (+) applies to bosons and (−) to fermions with and s i = c 2 i − 1 for bosons and s i = 1 − c 2 i for fermions 5 .From the perspective of Gaussian states, different purifications of ρ G A only differ in the choice of basis of the purifying system B, for which J takes the above standard form.Consequently, we can use the action of the respective Lie group G B (symplectic group Sp(2N A , R) for bosons, orthogonal group O(2N A , R) for fermions) to transform As reviewed in appendix D, we optimize over all Gaussian purification by taking the natural geometry (Fubini-Study metric) of the state manifold into account.Using the fact that this geometry is compatible with the group action, i.e., the Fubini-Study metric on the manifold of purifications is left-invariant under the group action of G A , we do not need to recompute the metric at every step, but can fix an orthonormal basis of Lie algebra generators equal to the dimension of the manifold.This enables us to efficiently perform a gradient descent search attuned the geometry of states, which scales polynomially in the number of degrees of freedom and enables us probe the field theory regime of our discretized models, which has not been possible previously in this setting.In particular, previous studies [33,34] of EoP restricted to special classes of Gaussian states (namely, real Gaussian wave functions generated by the subgroup GL(N A , R) ⊂ Sp(2N A , R)) for a small number of degrees of freedom.Similarly, CoP has been almost exclusively studied by purifying individual degrees of freedom (mode-by-mode purifications [35]) rather than whole subsystems for larger N A .
For purifications of small subsystems, e.g. of 1 + 1 or 2 + 2 sites, this optimization only takes a few seconds on a desktop computer, and is still feasible within a few hours for 10 + 10 sites, with efficiency depending on the optimization function, the accuracy threshold, and the hardware on which the computation is performed.For the particular case of CoP the optimization procedure for bosons was found to be an order of magnitude faster than the fermionic case for the same accuracy threshold, even for small subsystems.This implied that for larger subsystems, e.g., of order of 10+10 sites, the optimization parameters such as the gradient and function tolerance were lowered without compromising the results.For instance, lowering the gradient and function tolerance by a couple of orders of magnitude resulted in changes in the final value of the optimization in the third or fourth decimal.

V. ENTANGLEMENT OF PURIFICATION
We discuss our results for the EoP in bosonic and fermionic field theories using the purifications discussed in the previous section the algorithm described in appendix D.

A. Definition and existing results
EoP is a measure of correlations, which include both classical and quantum ones, and can be regarded as a mixed state generalization of entanglement entropy [60].When a mixed state ρ AB : H AB → H AB is given, we first purify it into a pure state |ψ ∈ H by extending the Hilbert space H AB according to (43) such that ρ AB = Tr A B (|ψ ψ|).The EoP E P (ρ AB ) is defined as the minimum of the entanglement entropy S(A ∪ A ) = − Tr(ρ AA log ρ AA ) for the reduced density matrix ρ AA = Tr BB (|ψ ψ|) over all possible purifications When ρ AB is pure, it simply reduces to the entanglement entropy as EoP is relatively new to the QFT setting and its understanding in this context is in development, which adds a strong motivation for our paper.Our knowledge about this subject is based on a conjecture in holography, results governed by local conformal transformations in CFTs and ab initio studies in free QFTs, which is the research direction the present work subscribes to.Below we briefly summarize the state of the art that sets the stage for the results of our research.
In strongly-coupled CFTs, a holographic formula which computes EoP was proposed in [31,63].Analytical calculations of EoP, based on the idea of path-integral optimization for CFTs [64], were given in [37].In particular, when the subsystem A and B are adjacent in a CFT, both holographic and path-integral result predict the universal formula where the widths of A and B are w A and w B , respectively and δ is the lattice spacing.Exploratory numerical calculations of EoP in a lattice regularization of (1+1)dimensional free scalar field theory have been performed in [33,34].Below we would like to extend such computations so that we can compare the result (45) with our discretized numerical calculations, as well as understand better the long distance physics (d w) in the QFT limit.The key technical difference on this front with respect to [33,34] is using bigger total system sizes, significantly bigger subsystems -both of which are desired to be closer to the QFT limit -and the most general Gaussian purifications discussed in section IV.

B. Numerical studies using the most general Gaussian purifications
Using the approach outlined in section IV and numerical techniques explained in appendix D, we can now compute both bosonic and fermionic EoP for purifications on the whole Gaussian manifold.In light of the discussion of our models in section II, for bosons we expect the Gaussian ansatz to describe well the CFT properties when the two intervals are adjacent (d = 0) or at small separation d w, as in these cases we do not expect the zero mode to be a significant contribution to the calculations we perform.For fermions, we expect the Gaussian ansatz to be appropriate for CFT calculations only when the two intervals are adjacent.Otherwise, the desired starting point of our calculations, spatially reduced density matrices for the Ising and Dirac fermion CFTs are non-Gaussian and our method is not applicable.In order to complete the picture, we will nevertheless provide results of our methods for bosons and fermions in the aforementioned regimes, however, they are not supposed to be seen as CFT predictions based on lattice calculations.
Starting with the adjacent intervals, our Gaussian lattice calculations perfectly reproduce the behavior (45) as shown in figure 3, in both our bosonic and fermionic (or equivalently, Ising spin) computations up to slight lattice effects at small w A or w B .
We move on to a more general case where the subsystem A and B are disjoint intervals in a free CFT.We again take the lengths of both intervals to be w and the distance between them to be d.When d w, both holographic [31,63] and path-integral approaches [37] predict the behavior which agrees with (45) under the replacement δ = d, w A = w B = w.On the other hand, no universal results have been known for d w and one possibility is a behavior similar to MI described in section III.
The numerical results for nearly a massless free scalar QFT are plotted in the second row of figure 2. As expected, we find a logarithmic dependence on w/d when it becomes large, given by with a convergence to c 1 ≈ 1 6 much faster than seen in MI.This result is consistent with the ∝ c 6 log w d behavior of (46).
In the case of small w/d, we observe that the bosonic EoP behaves extremely similar to bosonic MI.Such an observation for smaller subsystems and separations was already made in [34] and our results should be seen as a corroboration of this earlier finding.Given this similarity and our discussion of MI in section III, it should not come as a surprise that a power-law fit to the bosonic EoP in the regime of small w/d, is unstable as w/d → 0. The best we could do is to provide the upper bound on the power, d 2 0.15, which is consistent with the absence of a long-distance power behaviour.One should note that the power of such a quasipower-law for EoP agrees well with the one extracted for MI, as can be seen by comparing the two rows of figure 2 (right).

VI. COMPLEXITY OF PURIFICATION
In the present section we provide a comprehensive discussion of CoP in the single and two adjacent interval case (d = 0 in figure 1).We start by briefly reviewing the relevant results of the holographic complexity proposals, as well as the studies of pure state complexity in free QFTs.These results will guide us in the choice of a reference state and, also, in choosing the way to combine two and single interval CoPs to get a complexity analogue of MI.Subsequently, we discuss the CoP results obtained via optimization over the whole Gaussian manifold.We focus on the single and adjacent intervals to provide a clean message and we hope to report the d-dependence of CoP, which at least superficially seems involved, in further work.We also compare some of our results with a simplified version of a single mode purifications adopted in an earlier study of single interval CoP by [35] to avoid the technical problem our work addresses, i.e., optimizing over the full manifold of Gaussian purifications.Finally, we compare the properties of CoP with the notion of mixed state complexity discussed in [65].

A. Holographic predictions
Holographic complexity proposals relate novel gravitational observables associated with picking a time slice on the asymptotic boundary of solutions of AdS gravity with measures of hardness of preparing corresponding pure states in dual QFTs using limited resources.The first covariant notion is the spatial volume of the boundary-anchored extremal (codimensionone) bulk time slice (C V ) [17].The second covariant notion is the spacetime volume (i.e., a codimension-zero quantity) of the bulk causal development of such a time slice (C V 2.0 ) [20].The third covariant notion is also of a codimension-zero type and is the bulk action evaluated in the causally defined region (C A ) [18,19].The first two quantities are unique up to an overall normalization, whereas C A has an additional ambiguity related to the presence of null boundaries [45,46].
While there is also another evidence in support of the association of C V , C V 2.0 and C A with complexity, an important clue about the correctness of these conjectures comes from free CFT calculations of complexity of pure states along the lines of [22,23].In particular, such free CFT calculations are able to match the structure of leading divergences of holographic complexity [22,23,28,66] provided the reference state is taken to be a spatially disentangled state.Interestingly, in the case of bosonic calculations of [22,23] the scale entering the definition of a spatially disentangled reference state can be linked, via the leading divergence, both to the overall normalization freedom in the case of all three proposals, as well as to an additional ambiguity appearing in the C A case.Furthermore, the free boson CFT calculation in [29] explained qualitative features of the holographic complexity excess in thermofield double states as compared to the vacuum complexity reported in [67].
All three holographic complexity proposals acquire natural generalizations for mixed states represented as spatial subregions of globally pure states [68][69][70].Instead of considering extremal volumes or causal developments of a full Cauchy slice in the bulk, the mixed state version of holographic complexity proposals uses the corresponding notions applied to the relevant entanglement wedge [71][72][73].While there are certainly other possibilities regarding the kind of complexity the proposals [68][69][70] represent, see [32] for a discussion of some of the available options, we will treat their properties as a guiding principle to study CoP in free CFTs.
The results of these proposals applied to a single and two interval cases of interest can be found in table II.One can clearly see that the leading divergence of holographic complexity is in the volume of the combined subregions and there can be also subleading logarithmic divergences 6 .An earlier study of divergences encountered in the case of a single interval CoP in the vacuum of a free boson theory using restricted purifications is [35].In the present work we lift the restriction on purifications within the Gaussian manifold of states, include also the corresponding results for fermions and carefully resolve finite contributions to CoP including their dependence on the reference state scale and residual mass for bosons.The latter we achieve by considering the two adjacent intervals case.
For two intervals it is interesting to define a better behaved (less divergent) quantity in a manner similar to the definition of MI (7).Led by the form of leading divergences, as well as simplicity, the mutual complexity ∆C was defined in [74] as the sum of contributions for each individual intervals and subtract from it the holographic complexity of the union The results in the two intervals setup at a vanishing separations are included in table II and motivated us to seek for a logarithmic behaviour as a function of w A w B (w A +w B ) δ also in the analogous setting in free CFTs.This is also reminiscent of the behaviour of MI and EoP at d = 0, see table I.

B. Definition and implementation
CoP is defined in analogy with EoP as a measure of complexity for mixed states with the use of a definition of complexity for pure states minimized with respect to all purifications [32,75].This includes, in principle, purifications which contain an arbitrary number of ancilla greater or equal to number of the degrees of freedom in the subsystem.
Given a mixed state in a Hilbert Space H A characterized by a density matrix ρ A , we define a new Hilbert space with ancillary system A .There exist many purifications |ψ T ∈ H , such that ρ A = Tr H A (|ψ T ψ T |).In analogy to the EoP, see (44), we define the CoP C P as the minimum of the complexity function C with respect to a reference state |ψ R and to all purifications |ψ T : CoP inherits the richness of building blocks of complexity for pure states, such as dependence on the choice of a reference state |ψ R as well as on the cost function which evaluates the circuits built from the unitaries generated gence in the C A case by a logarithm of the cut-off and change subleading divergence.
by the Lie algebra of admissible gates.This is an additional complication with respect to EoP, which requires only minimization over purifications.For a generic definition of complexity underlying CoP one would not only need to optimize over purifications, but also for each purification one would need to solve an intricate optimization problem finding the optimal circuit.The idea behind the present work is to make use of a particularly natural definitions of cost function for Gaussian states defined in [28] (fermions) and [29] (bosons), which provide closed form and efficient to evaluate expressions for complexity.In this way we make the problem of calculating CoP in free CFTs much more manageable, as similarly to EoP, it requires now only one layer of numerical optimization.Of course, it would be very interesting to explore other cost functions in the CoP context and we leave this rather difficult problem for future studies.
To introduce the relevant cost function, let us recall from section IV (see also appendix C) that bosonic and fermionic Gaussian states |J can be efficiently characterized by their linear complex structure J a b .The latter can be constructed from their two-point function C ab 2 = ξa ξb , where we only consider Gaussian states with ξa = 0.As shown in [28,29] the geodesic distance between |J R and |J T within the Gaussian state manifold gives rise to a version of complexity based on a L 2 cost function To relate with the discussion in the introduction, the above definition of complexity corresponds to optimization with respect to the L 2 cost function (6) with where G ab = J R | ξa ξb + ξb ξa |J R for the reference state |ψ R .Due to the canonical anti-commutation relations, this normalization of the Lie algebra elements K I is independent of the reference state for fermions (see [28]), but for bosons (53) implies that K I is normalized based on the specific reference state, which was also referred to as equating reference and gate scale (as discussed in [22,23]).In the above expression K I are the respective Lie algebra elements (symplectic for bosons, orthogonal for fermions) associated to the quadratic operators O I in their fundamental representation acting on the classical phase space.
In the following, we focus on minimal purifications, i.e., , purifications whose ancilla have the same number of degrees of freedom as the reduced density matrix of the subsystem.Our focus on minimality comes as a result of a number of numerical computations for the cost function (52) which indicate that purifying the reduced density matrix with a larger number of ancilla does not lead to a lower CoP.It would be very interesting to explore if this feature is special to the cost function and the resulting complexity (52) we considered, however, we leave it for future investigations.
When applying the closed form complexity formula (52) As there is no a priori physical notion of locality in the ancillary system, we only require that |J R A is pure and Gaussian.We choose over spatially local sites i, where we only introduced a reference scale µ for bosons 7 .
The optimization over all purification in (51) could therefore be equivalently performed over reference or target state or even both.The minimum would always be the same, which can be seen as follows.By construction, the complexity function is invariant under the action of a single Gaussian unitary U acting on both states, i.e., we have where U is related to a group transformation M a b via U † ξa U = M a b ξb .In the case of Gaussian purifications, we optimize over all Gaussian purifications for the target state, i.e., if we have found such a purification |J T , any other purification is given by 1 A ⊗U A |J T .We thus find where the equalities follow from (55) and where both U A and V A are Gaussian unitarites on the system A .In practice, we can therefore start with a basis ξA , such that [J T ] A takes the mixed standard form (40).It can then be purified so that the purification takes the standard form with respect to the extended basis ξ = ( ξA , ξA ), as defined in (41).The reference state has the block diagonal form as it is a product state.We have (1 A ⊗ U A ) |J = |M JM −1 with M = 1 A ⊕M A , so the optimization gives which shows explicitly that we can think of the optimization as either applied to target or reference state.When we perform the optimization, it is actually advantageous to optimize over the reference state, as its stabilizer group is larger (i.e., there are more group elements that preserve J R than J T ) and so we can identify a fewer number of directions/parameters to optimize over.Our algorithm is described in more detail in appendix D and is one of subjects of the companion paper [61].

C. Single interval in the vacuum
The discussion in the previous section was very general and here and in the following we want to apply them to free CFTs on the lattice, as we did for EoP in section V.The first case to consider is a single interval in the vacuum.This case appeared earlier in the context of the aforementioned mode-by-mode purifications in [35].In the present section we readdress the same problem using the most general purifications, whereas in later section VI E we reconsider the same problem using our simplified take on the mode-by-mode purifications to make further contact with [35].In light of a general physics picture where circuits acting on a spatially disentangled state need to build entanglement at all scales to match features of CFT vacua, as well as explicit results in [35], one expects CoP to diverge in the continuum limit.
For a single interval on a line, fermionic CoP will be a function only of w δ as the system size N becomes large.Bosonic CoP, however, also contains two additional parameters, the reference state scale µ, see (54a), and the effective mass m δ.As changing µ → aµ is equivalent to rescaling the mass and lattice spacing according to m → m/a, δ → aδ, we can set µ = 1 in numerical calculations and restore it in analytical formulas containing the (now unitless and independent) m and δ.
We begin with the simpler case of fermionic CoP (at central charge c = 1  2 ).Here, we find a relationship of the form lim Note that we consider the squared CoP.We test this functional form by computing the discrete derivative with respect to w/δ, expected to be described by the expression e 2 + e 1 δ/w.As figure 4 (top) shows, it is indeed perfectly linear, allowing us to determine e 0 = 0.0894 , e 1 = 0.0544 , e 2 = 0.103 (61) with the given three significant digits corresponding to the numerical accuracy of the optimization algorithm.
The bosonic case (with c = 1) is somewhat more complicated, as we must subtract terms in m µ and δµ that diverge in the continuum limit m µ , δµ → 0. However, we still see in figure 4 (bottom) that the functional form of (60) still holds, but with dependencies lim This form accurately describes the w dependence over a large range of m/µ and µδ.The functions f 0 to f 2 are estimated as in the region m/µ, µ δ 1.The leading divergences in µ δ and m µ are visible in figure 5, where we plot C P (non-squared) and find linear divergences in log(µ δ) and log m µ , respectively.The dependence of the slope of the divergence on w, given by ≈ √ w 2 , is clearly visible in the    In both cases, the number of total sites is given by N = 100 w/δ.µ δ case, while the m µ term diverges with a constant slope ≈ 1 2 at m µ 1, consistent with the appearance of these terms in f 2 and f 0 , respectively.Note that f 0 and f 1 are estimated from the setup of two adjacent intervals, analyzed in the next section, where the linear term f 2 cancels.In particular, the square root term in f 0 can be seen in figure 6, where the leading divergence in m µ is subtracted.
To corroborate this discussion, let us also note that the structure of the leading divergences in the two cases matches the result of the vacuum complexity in free CFTs, see [22,23] for bosons and [28,66] for fermions.In this pure state case analogy, the role of w is played by the total system size measured in lattice units.The presence of the log 2 µ δ contribution in the vacuum case comes from the ratio of the highest momentum frequency of the order of the inverse lattice spacing to the reference state scale and the overall coefficient in front of the whole divergence is µ-independent.The logarithmic divergence is present because the symplectic group is non-compact.For fermions, the group of transformations is compact and there is no logarithmic enhancement of the leading divergence.

D. Two adjacent intervals in the vacuum
The next case to consider are two adjacent intervals in the vacuum.This is basically the application of the formulas from the previous section with the twist that it allows us to gain a better control over the finite term f 0 ( m µ , µ δ) in the bosonic single interval CoP (62).
To this end, we are interested in a better behaved combination of complexities akin to (49).We take it to be ∆C (2) where we put two on the LHS in the brackets to emphasize that it does not denote taking a square.The rational behind this expression is that when one keeps µ δ fixed, the whole power-law divergent part cancels between the three terms.Similar combinations to (64) in the aforementioned context of pure state complexity in thermofield double states appear in [29].Also, note the difference with respect to holographic mutual complexity (49).Simple manipulations lead to the following result for the Ising CFT (c = 1  2 ): lim For the free decompactified boson CFT (c = 1), we can again use the single-interval result (62) to arrive at lim N →∞

∆C
(2)  but where the logarithmic coefficient and constant term now depend on the divergences in m µ and µδ in precisely the same way as the single-interval expression.The form of ∆C (2) P for both bosons and fermions in shown in figure 7 in terms of the ratio w A w A +w B .The qualitative behavior is the same as for the EoP result shown in figure 3 and also the holographic complexity proposal results encapsulated in table II.However, one should bear in mind a different subtraction of complexity with the latter case related to the use of a L 2 norm in our Gaussian studies.Note also that ∆C (2) P is logarithmically ultraviolet divergent.Finally, it is interesting to observe that the mutual complexity depicted in figure 7 is positive, which indicates subadditivity of our CoP definition.This is in line with an earlier observation in the case of two coupled harmonic oscillations in [36].

E. Single mode optimization for bosons
A natural hope in the context of Gaussian states is to split the problem of finding the CoP for a system with N A modes into N A problems for a single mode.In [35], the authors use a formula for the L 1 norm complexity of a single bosonic degree of freedom for certain Gaussian states derived in [23], where they introduce two types of L 1 norm bases (called the physical and the diagonal one).Note that the authors use the geodesic with respect to the L 2 norm, but compute its length with respect to the L 1 , as it is difficult to prove that a path is minimal with respect to the L 1 norm (in particular, for several modes).When considering several modes of the free Klein-Gordon field, we need to distinguish two settings reviewed previously: • Thermal states.Here, we can choose a basis ξA , such that both [J T ] A of a mixed thermal Klein-Gordon state and [J R ] A of a spatially unentangled vacuum decompose into 2-by-2 blocks, so it is easy to argue that the Gaussian CoP results from optimizing over individual modes.
• Subsystems.If we consider a mixed Gaussian state ρ A resulting from restricting the Klein-Gordon vacuum to a region, it is typically not possible to bring both [J T ] A and [J R ] A of a spatially unentangled vacuum into block diagonal form, so it will not suffice to optimize over individual modes.
Let us emphasize that in both cases, we have a modeby-mode purification with respect to the standard form (41), but it is the reference state that will only be a tensor product over these modes in the case of thermal states, but not for local subsystems.However, we may encounter situations where the standard decomposition of the mixed target state also approximately decomposes the reference state into individual modes.
We consider a single bosonic mode with a pure Gaussian reference state and a mixed Gaussian target state that both do not have ϕπ-correlations, i.e., φi πj = 0. We extend this system to a system of two bosonic modes H = H A ⊗H A to have the extended reference state |J R and the purified target state |J T , such that the respective complex structures are given by where λ ∈ [1, ∞) is the same as c i for several degrees of freedom in (40), µ is the reference state frequency for the original single mode, and ν is a parameter in the reference state, for which we will minimize the complexity functional (52).The latter is given for us by where we defined the variables In order to find the minimum of C P = min ν C(λ, µ, ν), we need to solve a transcendental equation for ν.This can be done numerically for any given value of c and µ in a very efficient manner.Unfortunately there is no closed analytic expression for C P of a single mode, but we reduced the problem of a single mode (with vanishing ϕπcorrelations in reference and target states) to a problem that is much simpler than the optimization over a large manifold.Can we extend this to larger systems?As derived in [22,28] for the perspective of the L 2 norm, the optimal geodesic from reference to target state gives (52), so we only need to worry about what goes into this formula.We consider a mixed state of a region A with N A sites.With respect to the local basis ξA = ( φA 1 , . . ., φA the covariance matrix of reference and target state will have the following forms If the Gaussian target state were pure, it is well-known that we can find symplectic transformation where As soon as G T represents a mixed Gaussian state, there still exists a symplectic transformation M , such that GT = M G T M is diagonal, but now M will not be of the form (73) anymore and will thus not preserve G R , i.e., M G R M = G R .However, we could pretend to approximate the true M by only diagonalizing G ϕϕ with a respective O, such that Gϕϕ = OG ϕϕ O is diagonal.We then apply the respective M of the form (73) and pretend that also Gππ = OG ππ O is diagonal, i.e., we drop the off-diagonal terms which we hope to be sufficiently small.With this assumption, we can apply a mode-bymode optimization based on (69), such that where λ i is extracted from the diagonal entries of Gϕϕ and Gππ .For a pure target state, (74) becomes an equality, where both sides match the regular complexity (52).
We consider the restriction of the Klein-Gordon vacuum to the subsystem of a single interval as explored in section VI C. In this setup, we can compare the approximate single mode optimization with the full optimization.While the full optimization takes several hours on a regular desktop computer, our approximate scheme of optimizing (69) over individual modes only takes a few seconds.Figure 8 shows how the single mode optimization is almost indistinguishable from the full optimization for small µL, but the approximation becomes increasingly worse for larger µL.In [35], the authors perform a similar calculation for the L 1 norm, which they refer to as mode-by-mode purifications 8 .The difference between 8 As pointed out at the beginning of this section, any Gaussian purification is a mode-by-mode purifications, but what the authors of [35] mean is that they only optimize over the purifications in a specific way, as if reference and target state would decompose into the same individual modes as an approximation.our single mode optimization and what the authors in [35] do is two-fold: First, we optimize over the L 2 norm, while they consider the L 1 norm.Second, we change the target state by hand to decompose into a product over modes in the same basis as the reference state and then perform the optimization semi-analytically for individual modes, i.e., we optimize for each mode independently.In contrast, the authors [35] do not change reference or target state, but only consider a subset of possible purifications, i.e., they evaluate the full complexity function and optimize over a restricted subset of parameters (one parameter per mode).As they do not change the target state, they cannot evaluate the complexity for individual modes, so would need in principle to optimize over all parameters simultaneously, but find good convergence when optimizing over O(1) parameters at once.Clearly, the approximation in [35] and our single mode optimization work because the for small µ L reference and target states are close to decomposable over individual modes.This will not be the case for generic subsystems (such as two intervals), general states (such as those with ϕπcorrelations) and fermionic systems (which cannot be decomposed into single mode squeezings), in which case our full optimization algorithm is required.

F. Comparison with the Fisher-Rao distance proposal
In [65], the authors propose a measure of bosonic mixed state complexity based on the Fisher-Rao distance, which can be defined on the manifold P(N ) of 2N × 2N real and positive definite matrices, of which bosonic covariance matrices are a subset of.Without the need for any purifications, the proposal for the complexity of mixed states is formally equivalent to (52), where here G and G 0 are taken to be the covariance matrices of the mixed target and reference state, respectively.The motivation for this definition is that the Fisher-Rao distance function measures the geodesic distance in the manifold of covariance matrices.It is important to highlight that the authors in [65] focus on bosonic Gaussian states occurring in the Hilbert space of harmonic lattices, and hence the proposal of the Fisher-Rao distance function should be thought of as applicable in principle only to bosonic states, although one could conjecture that a simliar formula should be applicable to fermionic Gaussian states.
Nonetheless, it is interesting to compare the properties of such distance function with the bosonic CoP measure arising from the Gaussian optimization procedure developed in this paper.
For the single interval case, there is in fact a noteworthy qualitative and quantitative agreement between the two, as shown in figure 9.The Fisher-Rao distance function and the CoP measure offer a comparable measure of the complexity of the mixed state associated to the single interval, which is remarkable given the fact that the Fisher-Rao distance function is a geodesic distance being evaluated on the manifold of mixed states on the Hilbert space H A , whereas CoP is a geodesic distance on the manifold of pure states on the larger Hilbert space H = H A ⊗ H A and these two need not be the same or even comparable to one another, as explained on figure 12.This comparison seems to work better, the smaller the parameter µL is taken to be.As our studies indicate, for two adjacent intervals the distinction between the Fisher-Rao distance and CoP deviate significantly from each other, even though the qualitative behaviour remains comparable.

VII. COMMENTS
In our considerations, two important subtleties of free QFTs, known to the literature, played a key role and influenced the vacuum subregions we could consider to make genuine QFT predictions.They are the zero mode in the case of free boson QFTs and spatial locality of disjoint intervals under the Jordan-Wigner mapping between the Ising model and the free Majorana fermion theory.Below we provide an additional discussion of these FIG.10.MI on an infinite line in the small m δ limit.Shown are the general mass dependence (left) as well as the decay power in the d/w → ∞ limit compared to the periodic system in the large N limit (right).two important points.

A. Zero mode for free bosons
The presence of the zero mode is a known subtlety of the free boson theory in the massless limit, as we discussed in section II A. The simplest way to deal with it is to keep the mass term in the Hamiltonian (8) nonzero and try to numerically approach the limit m → 0, which is precisely the strategy we adopted in the calculations of the bosonic EoP and CoP.Given the scarceness of other methods to shed light on EoP and CoP in QFTs, it is important to learn about the role of zero mode in betterunderstood problems.
To some degree we already explored this issue in section II A when understanding what needs to be done in order to reproduce the modular invariant thermal partition function (15) from the Gaussian calculation in the massive theory (16).Here we will address another quantity, which is the two interval vacuum MI reviewed in section III.Fitting a power law for the data on a periodic chain of bosons does not lead to convergence of the estimated power to a nonzero value in the d/w → ∞ limit (see table I), implying a decay with slower asymptotic functional dependency.Indeed, the massless limit of our study can be thought of as a decompactification limit of the free boson theory on a circle in the field space (which can be seen as another way of dealing with the zero mode), as reviewed in section II A, and the predicted large distance behaviour in this case is not of a power-law type.In our periodic setup we considered the limit of a large number of sites N with the dimensionless scale mL = mN δ kept constant and small.In this limit, the mass dependence of both MI and EoP is accurately described by an additive − 1 2 log(Lm) term (as already noted in [34]), so that the dependence on w/d can be studied independently of the mass.
However, one may alternatively consider free bosons on an infinite line, i.e., taking the limit m δ → 0 only after the limit N → ∞.As the scale m δ N formally di-  77) and (78).
verges in the first limit, the resulting mass dependence is qualitatively different from the periodic setup.To investigate the behavior in this case, we adapted our numerical method to free bosons on an infinite line and computed MI at extended precision at small m δ.Comparable studies were performed earlier in [76] and the authors reported a power-law fall-off of the form w d 0.05 , see also (36).The value of I(A : B) on the line changes only very slowly with mδ, as shown in figure 10 (left).For d/w mδ, this mass dependence can be expressed by a constant offset where the factor of 1/2 is reproduced with a numerical accuracy of four significant digits.This double-logarithmic infrared divergence matches previous results for singleinterval entanglement entropies in a similar setup [77].
Note that the positivity of I(A : B) prevents this dependence persisting at finite m δ, and thus I(A : B) begins to decay exponentially as d/w m δ.Apart from the different mass dependence, the periodic and line setup in their respective limits yield equivalent results.As we show in figure 10 (right), the estimated decay power from both limits exactly matches, vanishing as d/w → ∞.
As the line setup is more efficient at probing large values of d/w due to an absence of finite-size effects, we use this setup to test for functional dependencies slower than the previously considered power law and logarithmic functions.As the functional dependence of MI and EoP match in this limit, we expect the results to extend to both measures.In particular, we consider the powerlogarithmic asymptotics as well as a power-double-logarithmic one, In figure 11 coefficients from both fits are shown.While the power e 2 of the power-logarithmic fit converges to a value e 2 0.1 that may be zero in the d/w → ∞, mδ → 0 limit, the double-logarithmic power f 2 clearly converges to a value f 2 1.3 that is visibly larger than zero.Surprisingly, this clear observation of sub-polynomial functional dependence both in the infinite line and periodic setups contradicts the aforementioned earlier numerical observations of a power law in [76].These earlier results may be the result of numerical estimates performed only at moderately large d/w; while the functional dependence in this range can be well approximated by a power law, as also shown in [34], the apparent power law is not stable as d/w → ∞.Curiously, the authors of [76] analyze in the same setup another measure of entanglementthe logarithmic negativity [78] -and our extended precision calculations in this case reproduce the exponential fall-off as reported in [76] for the same range of masses considered in figure 11.This shows that not all non-local quantities are affected by the zero mode problem.
An additional insight about the expected behaviour of MI in the case of the free boson CFT in the decompactification limit R → ∞ can be obtained from having another look at the modular invariant partition function on the circle (15).The partition function provides the information about the density of states, which, via the state-operator correspondence, describes also the density of operators in the spectrum.The latter quantity, in conjunction with (32), will provide an indication on what to expect from the two interval case in the large separation limit for the decompactified free boson CFT.To this end, the partition function of a CFT with a continuous spectrum of operators (for a decompactified free boson theory vertex operators are labelled with a continuous index) can be written as where the second exponent comes from the Casimir energy and ∆ is the scaling dimension of operators in the theory.Note that descendent operators appear in the sum only for ∆ > 1.The power law multiplying the Casimir contribution in the low temperature limit of (15) points to the density of operators behaving in a power law fashion in the vicinity of ∆ = 0.In particular, the behavior can be explained by Note that for a free decompactified boson In this case, the density of states can be easily understood to be given by ∆ −1/2 upon noting that the operators of interest are vertex operators : e i ν φ : specified by a real number ν.The scaling dimension of vertex operators is given by ∆ ∼ ν 2 .The density of operators is uniform when parametrized by ν and viewing it as a function of ∆ brings in the Jacobian ∼ ∆ −1/2 , which gives (82).Now we can come back to MI at large separation.Since the formula (32) incorporates the exchange of a single operator, in the absence of a gap in the spectrum one needs to sum over the continuum of light operators with their density given by (81).Following [41] we can schematically write where c T T O∆ is the three-point function coefficient between two twist fields and a primary with dimension ∆ and the ellipsis denote contributions with higher powers of w d 4 ∆ .At the present moment we do not have control over the two kinds of contributions.However, neglecting the additional contributions and assuming that c T T O∆ has a power-law dependence on ∆ for small scaling dimensions the long distance behaviour of MI becomes Let us re-stress that the above equation is based on unverified assumptions and the correct answer is likely to be more involved, yet in principle calculable.However, what (85) indicates is that MI may decay much more slowly with distance than a simple power law, as is also shown by our numerical results.In particular, the fitting ansatz corresponding to ( 85) is (77) consistent with −α − 2 κ 0.1.Finally, note also that while in the massive theory, the mass scale eventually triggers an exponential decay of MI to zero, the ansatzes ( 77) and ( 78) predict divergence when naively extrapolated to d w → ∞.It is unclear at the moment if this is a feature of regularizing the zero mode via introducing a non-vanishing mass, or if the behaviour (77) or ( 78) that we see using our Gaussian numerics persists in the free boson CFT in the decompactification limit.

B. Subsystems in Ising CFT vs. free fermions
There is a subtle difference for computations of MI and EoP between the XY spin model and the Majorana fermion model.The spin model in the continuum limit becomes the genuine c = 1/2 Ising CFT, which has a modular invariant partition function.However, the fermion model can only be identified with this CFT in the continuum limit if we impose the correct GSO projections, summing over the sectors of periodic and anti-periodic boundary conditions (also known as the R and NS sector, respectively).This projection, originating from anti-commutative nature of fermions as opposed to spins, maps entanglement entropy and related quantities like MI and EoP for standard choices of subsystems in the spin model to unconventional ones in the fermion theory.This difference due to the projections does not occur in the calculations of MI and EoP when the subsystem A and B are adjacent.However, when A and B are separated, the non-trivial topology in the replica computation of entanglement entropy leads to a difference between the spin and fermion calculation as the simple choice of subsystems does not respect the projections.This is consistent with our numerical results for Majorana fermions and adjacent intervals shown in figure 3. We also note that MI and EoP are smaller in the fermion model calculations without projects compared with those in spin model calculations.Refer to appendix B for more details.

VIII. DISCUSSION
In this manuscript, we have a presented a systematic and comprehensive analysis of EoP and CoP that characterize mixed states.We computed the EoP between two blocks of widths w A and w B at distance d in onedimensional periodic systems at large size for both critical bosons and fermions, the latter of which are equivalent to a discretized Ising CFT while d = 0. Furthermore, we compared these results with the well-studied MI.At d = 0, our data shows confirming previous expectations through analytical methods [37,38].For d > 0, we considered the symmetric setup w A = w B = w in the two limits d/w 1 and d/w 1.In the former limit, our data is consistent with for our bosonic model.The latter limit shows a subpolynomial (logarithmic or double-logarithmic) decay of both I(A : B) and E P at large d/w.In summary, MI and EoP show the same scaling at large distance d.This result is consistent with the observation in [34] that EoP appears to weight quantum and classical correlations differently from MI, leading to different qualitative behavior only when both become relevant, i.e., at small distances.Indeed, this is the regime in which our numerical results show new model-dependent features that distinguish both measures.For both the periodic and infinite line setup, two-interval MI and EoP are divergent in the zero-mass limit; this divergence can be regulated by a 1 2 log(mδN ) and 1 2 log log(mδ) term, respectively.Let us also emphasize that the large distance behaviour of MI in the free boson case at small masses is very subtle and is described by the fall-off slower than any power-law, contrary to earlier studies in the literature.We discussed this at length in section VII A.
In the studies of CoP, our only guidance were the predictions of holographic complexity proposals for subregions, as summarized in table II.It is interesting to note, that our studies reproduce qualitatively terms present in the holographic results.In particular, for appropriately defined mutual complexity, all complexity notions we considered lead to an analogous dependence on the sizes of two adjacent intervals to the one seen in for entanglement measures in (86).What is also worth a separate remark is an intricate dependence of the subleading divergent and the finite term in the bosonic CoP on the reference state scale µ and the mass m acting as the zero mode regulator, see (63).
The other important set of CoP results has to do with a comparison with an earlier study of CoP using a restricted set of Gaussian purifications in [35], as well as with the Fisher-Rao distance introduced in a related context in [65].As we saw respectively in sections VI E and VI F, both of these notions can reproduce qualitative behaviour of the CoP, but they can also exhibit significant deviations from the numerical answer predicted by the full Gaussian optimization.This indicates that in general one indeed needs to optimize over as large sets of states in the Hilbert state as possible to reproduce the true CoP.Of course, even if the reduced density matrix is Gaussian, this does not imply that the optimal circuit is necessarily such.Generalizing our study to non-Gaussian states is an outstanding challenge and in the last paragraph of this discussion, we sketch what we believe might be an interesting and workable example.
In the context of CoP, we also want to offer an intuitive, yet rigorous interpretation of the setup that we were using.Formula (52) was derived in [28,29] based on a certain metric on the group manifold, which coincided with the geodesic distance on the Gaussian state manifold (Fubini-Study metric), as studied in [22].Note that there is an ongoing debate what the most appropriate notion of distance is, both on the manifold of states and the Lie group.While the L 1 norm appears to have similar properties as the different holographic complexity proposals, only the L 2 norm induced by a certain Riemannian metric can be analytically minimized.This is the key reason why we focused on the L 2 norm in the present manuscript, as most other notions of distances are intractable in the setup we are considering (many degrees of freedom, most general Gaussian gates, analytical optimization over all trajectories from reference to target state, numerical optimization of all possible pu-pure 12. We sketch how the manifold of mixed states on Hilbert space HA is related to the manifold of pure states on the larger Hilbert space H = HA ⊗ H A .We indicate the manifold (red line) of all possible purifications |ψT AA related by Gaussian unitaries U A .The CoP CP is given by the geodesic distance (blue) between the purified reference state |ψR AA .rifications).Figure 12 provides a visualization of this interpretation: CoP becomes in essence another type of minimization, namely finding the minimal distance between the unique reference state |ψ R to the the set of all possible purifications (or all Gaussian purifications) that is fully determined by the single mixed state, whose CoP we are computing.
Let us emphasize that there have been several approaches to define and compute complexity for mixed states, which predominantly focus again on Gaussian states.On the one hand, there are approaches [32,35,36] based on purifications, to which a notion of pure state complexity is applied.This provides an elegant way to carry any definition of pure state complexity over to arbitrary mixed states.Our present draft uses the same philosophy, but performs the required optimization over all possible (Gaussian) purifications numerically.On the other hand, there are approaches [65] to define mixed state complexity directly on the set of mixed states.Here, one introduces additional non-unitary gates that allow to change the spectrum of the density operator, so that unitarily inequivalent mixed states can be reached.The resulting geodesic distance agrees with the L 2 norm when restricted to pure states and the procedure can be understood as measuring the geodesic distance on the manifold of (Gaussian) mixed states equipped with the Fisher information geometry.Remarkably, the resulting analytical formula in terms of covariance matrices for bosons agrees with the one (51) derived for pure states.One can expect the same result to also hold for fermions.In particular, we can use these gates to transform a pure reference state into a mixed target state.Finally, a complementary approach is based on path integral optimization [64,79], which uses gates being exponentials of the energy-momentum tensor operators with complex coeffi-cients and so both unitary, as well as Hermitian operators [80].In the case of free CFTs, these are Gaussian gates, however, in the interacting cases they are not.
An important feature of our works stems from the fact that the Ising model in the spin picture (rather than the fermionic picture) leads to genuinely non-Gaussian mixed states if we restrict to disconnected regions, i.e., the spectrum of the reduced density operator to such subsystem cannot be reproduced by bosonic or fermionic mixed Gaussian states.This may open the window to study non-Gaussian circuit complexity (of purification) in a genuine QFT limit.For this, we propose to consider two individual separated sites in the spin picture of the critical Ising model.This leads to a mixed state in a system of two qubits, i.e., in a four-dimensional complex Hilbert space.We can now study CoP in this setup as a function of the separation between the two sites 9 .While these are ideas for future work, the implemented algorithm to optimize over all possible purifications (in this cases non-Gaussian ones) can be used, once an appropriate notion of non-Gaussian circuit complexity (for generic gates) is defined.In this sense, we believe that our considerations of CoP in the context of the critical Ising model CFT presents a stepping stone to explore genuine non-Gaussian circuit complexity in the context of QFTs.k = (k + 1) mod 2N , so we find where |ψ + is a Gaussian state vector with Ptot |ψ + = |ψ + .In matrix form, this can be written as where Γ + [i,j],[k,l] is the sub matrix consisting of the rows from i, i + 1 up to row j and columns from k, k + 1 up to l.For a parity-odd Gaussian state vector |ψ − with P tot |ψ − = − |ψ − , however, the rows and columns corresponding to the first two Majorana modes are signflipped: .(B15) In consequence, only parity-even fermionic Gaussian states are equivalent under different Jordan-Wigner transformations related by cyclic permutations, while parity-odd ones acquire a sign flip in the two-point correlations related to modes that are moved from the end to the beginning of the fermionic ordering.Fortunately, this sign flip does not affect fermionic entanglement entropies S(A), which are computed from the spectrum of the submatrix of Γ corresponding to sites in A, denoted Γ A .
As the eigenvalue spectrum is not affected by change of sign across entire rows and columns of a matrix, Γ| Ã and Γ |A have the same eigenvalues, and hence, S(A) = S( Ã) independent of parity.Note that once we consider index transpositions more general than cyclic permutations, the relative fermionic ordering is broken and even Gaussian states are no longer equivalent under different Jordan-Wigner transformations.

Ising CFT representations
Here, we would like to discuss subtle differences between the spin model and Majorana fermion model when we calculate MI and EoP.Refer also to [83] for an earlier detailed analysis on this problem, which is essentially equivalent to our argument below.
First let us remember that the Ising spin CFT is defined from Majorana fermion CFT by the GSO projection [55].Explicitly the modular invariant torus partition function of the Ising CFT is schematically written as where the first trace is taken for the NS sector (H + sector), i.e., the anti-periodic boundary condition is imposed for the free fermion on a circle.Also 1+(−1) F 2 is the restriction to the even fermion number state.The second one is for the R sector (H − sector), i.e., the periodic boundary condition is imposed for the free fermion on a circle.Also 1−(−1) F 2 is the restriction to the odd fermion number state.Note that the spin operator σ is included in the R sector of the Majorana fermion model.
In the calculation of entanglement entropy S(A ∪ B) or its Renyi entropy in Ising CFT, we need to perform this GSO projection of the Majorana fermion along the interval C between A and B. This is illustrated in the upper pictures in figure 13 by having in mind the calculation of Tr[ρ 2 AB ] as an example, which is essentially the 2nd Renyi entropy and which is equivalent to a torus partition function (B16).If we write the fermion number on this interval C as F C (i.e., F C = i∈C f † i fi ), then the wave function for the calculation of Renyi entropy should be which is depicted in the lower pictures in figure 13 This procedure leads to the spin operator in the spectrum and leads to the expected behavior of MI I(A : B) ∼ (w/d) 1/2 when d w.However, this form (B17) is not included in our Gaussian ansatz of the numerical calculation.The Gaussian ansatz only takes into account the first term which ignores the twisted sector of Majorana fermion namely the spin operator sector (i.e., R sector).This causes the absence of the spin operator excitation and explains our numerical result I(A, B) ∼ (w/d) 2 for the Majorana fermion model when d w.Moreover, this difference between the XY spin model and the Majorana fermion model leads to a significantly different behavior of MI when d w.To see this, as a toy model which mimics the calculation of MI in XY spin model at d = 1, we would like to analyze a three qubit system, whose spins are called A, C and B, from the left to right.We consider the following spin state for this ABC system:  where G ab is a symmetric positive definite metric and Ω ab is non-degenerate antisymmetric symplectic form.Here, Ω ab is fixed by (C1) for bosons, while G ab is similarly fixed by (C2) for fermions.Consequently, only the respective other piece, i.e., G ab for bosons and Ω ab for fermions, which are often called the bosonic and fermionic covariance matrix, respectively, will depend on the state |ψ .We can define the linear map We refer to the state |ψ as Gaussian if and only if J squares to minus identity, i.e., |ψ is pure Gaussian state ⇐⇒ J 2 = −1 , (C6) in which case J is called linear complex structure.
The entanglement entropy S(A) of a Gaussian state |ψ with complex structure J can be efficiently computed from the restriction J A of J to the respective subsystem A. More precisely, if we have H = H A ⊗ H B with phase space decomposition V = A ⊕ B. The restriction represents then the 2N A -by-2N A sub matrix of J associated to the subspace A ⊂ V .While J associated to a pure Gaussian state has eigenvalues ±i, its restriction J A will have eigenvalues ±λ i with λ i ∈ [0, ∞) for bosons and λ i ∈ [0, 1] for fermions.The entanglement entropy can be computed from these eigenvalues using the formulas where we wrote covariant matrix equations, which are equivalent to replacing J A by its eigenvalues and the trace by a sum over them.Such analytical formulas for Gaussian entanglement entropy were first derived in [84,85] based on the covariance matrices and later also phrased in terms of linear complex structures [86][87][88].
Note that the definition of entanglement is a subtle issue.Usually, one decomposes the Hilbert space into a regular tensor product H = H A ⊗ H B , such that operators O A and O B that are only supported in one of the subsystems are given as O A = Â ⊗ 1 and O B = 1 ⊗ B, leading to [O A , O B ] = 0.For fermions, we would like to define a similar structure, but independent operators should anticommute, i.e., {O A , O B } = 0.This requires a slightly different definition of "fermionic tensor product", as recognized and discussed in [50][51][52]89].We discuss some of the resulting subtleties in the next section.

Appendix D: Algorithm's implementation
Our minimization algorithm has the goal to minimize a function f (J), where J is the linear complex structure of a pure Gaussian state.Intuitively, we would like to apply gradient descent, i.e., solve the differential equation where G µν is the matrix representation of the Riemannian metric on the manifold of states.Using coordinate has two disadvantages: First, in general it will be difficult or even impossible to find a global coordinate system depending on the topology of manifold (in particular, for fermions the topology is non-trivial).Second, the matrix G µν (x) will depend on the position and would need to be evaluated at every step, which will slow down the computation.Our approach based on the natural Lie group parametrization avoids both of these disadvantages.
Our parametrization is defined relative to an initial state |J 0 .We then choose a basis of Lie algebra gener-ators (Ξ µ ) a b satisfying the conditions {Ξ µ , J 0 } = 0 and One can show [90] that the dimension of this space is N (N + 1)-dimensional for bosons and N (N − 1)dimensional for fermions.Rather than using coordinates x i , we use the matrix M a b to parametrize all Gaussian states J = M J 0 M −1 connected to J 0 .The gradient descent equation for M reduces then to where the gradient vector is computed as [90] F µ (M ) = − d ds s=0 f M e sΞµ J 0 e −sΞµ M −1 .(D4) In our algorithm, we discretize (D3).Starting from the identity M 0 = 1, we perform individual steps as where Here, we approximate the exponential for small in such a way that M n+1 is ensured to lie in the symplectic or orthogonal group for bosons or fermions, respectively, which e Kn ≈ 1 + K n would not achieve.At each step, we choose 0 < < 1 in such a way that f (M n+1 ) < f (M n ), which is always possible to achieve for sufficiently small .
Our algorithm circumvents the disadvantages of a coordinate parametrization by distinguishing between state and tangent space.While we use M to parametrize our state |J with J = M J 0 M −1 , we construct an orthonormal basis of Lie algebra generators K i which can be identified with the respective tangent vector of the curve γ(s) = M e sΞµ at the point M .This allows us in particular that our Riemannian metric on the manifold of Gaussian states is left-invariant, such that the so constructed tangent vectors are orthonormal at every point.We therefore do not need to compute its matrix representation G µν (x), but can instead work with orthonormal frames at every point M .
FIG.1.The subsystem that defines reduced density matrices for our discretized bosonic and fermionic models in their vacuum state consists of two intervals of a width of wA/δ and wB/δ sites and separated by a distance of d/δ sites, where δ is the lattice spacing.When d = 0, we will keep wA and wB generic.When d > 0, we will set for simplicity wA = wB ≡ w and the natural continuum combination is w/d.We will see that numerically determined MI and EoP approach in the continuum limit functions of w/d.With CoP the situation is more complicated, as it turns out to be ultraviolet divergent and brings in an additional dimensionful scale through the class of reference states of interest |ψR .

1 6 1 12
log w A w B (w A +w B )δ log 2w A w B (w A +w B )δ

Lattice spacing δ = 10 - 3 FIG. 4 .
FIG.4.Discrete derivative of fermionic (top) and bosonic squared CoP (bottom) of a single interval of length /δ.m and δ are given in units where µ = 1.In both cases, the number of total sites is given by N = 100 w/δ.

FIG. 5 .
FIG.5.Divergences of bosonic CoP of a single interval of size w at fixed µ δ = 10 −3 and m µ → 0 (left) and at fixed m µ = 10 −3 and µ δ → 0 (right), in units of µ = 1.In both cases, the leading divergence is linear.While the µ δ divergence is wdependent, the m µ one is not.

10 - 4 , δ = 10 FIG. 7 .
FIG.7.Fermionic/Ising (left) and bosonic (right) CoP for two adjacent (d = 0) subsystems A and B, in units of µ = 1.The expected analytical forms,(65) for fermions and (66) for bosons, are plotted as dashed curves.We consider w A +w B δ = 14 total sites for fermions and 20 sites for bosons.The total system size is set to N = 100 w δ for each block width w.

μL = 1 FIG. 9 .
FIG. 9. Comparison of CoP for a single interval obtained using the Gaussian optimization algorithm (solid) and the Fisher-Rao distance function (dashed).The data are generated for mL = 10 −2 and N = 100.
FIG. 14.The first graph plots the value of I s (A : B) (blue) and I f (A : B) (orange) as a function of θ at x = y = 0.The second graph is a 3D plot of the ratio I f (A : B)/I s (A : B) as a function x (horizontal) and y (depth) at θ = π/8.The third graph is a 3D plot of the ratio I f (A : B)/I s (A : B) as a function x (horizontal) and y (depth) at θ = π/4.

TABLE I .
Overview over known analytical results and numerical fits with approximate coefficients for entanglement entropy S(A), MI I(A : B), and EoP EP , all for an infinite-size system.Analytical entries marked with a star (*) are guesses based on analogous behavior between MI and EoP.Twice the free fermion result in the d = 0 Ising is the c = 1 free Dirac fermion CFT quantity.

TABLE II
. Summary of mixed states complexity results consisting of predictions of holographic complexity proposals collected from to Gaussian purifications |J T of some mixed state ρ A , we need to think about what an appropriate reference state |J R can be.The two most immediate applications are thermal states and mixed states resulting from the reduction to spatial subsystem: • Thermal states.Our mixed state could be the thermal state ρ in a system H =: H A , which we purify to H = H A ⊗ H A .Here, we can always choose a spatially unentangled and pure reference state, which we can extend to the purifying system as |J R = |J R A ⊗ |J R A ∈ H . • Subsystems.We consider a pure Gaussian state |ψ ∈ H = H A ⊗H B , which we reduce to some local subsystem ρ A = Tr H B |ψ ψ|.In this subsystem, we have a pure and spatially unentangled Gaussian reference state |J R A , which we can extend to the purifying system as |J R= |J R A ⊗ |J R A ∈ H .The spatially unentangled character of |J R is a choice motivated by the fact that such a state is on one hand truly simple and, on the other, in the case of pure state complexity it reduces the kind of divergence encountered in the holographic complexity proposals.In both scenarios outlined above, only the target state |J T is entangled across H A ⊗ H A , while the reference state is a product state