Long-distance entanglement of purification and reflected entropy in conformal field theory

Quantifying entanglement properties of mixed states in quantum field theory via entanglement of purification and reflected entropy is a new and challenging subject. In this work, we study both quantities for two spherical subregions far away from each other in the vacuum of a conformal field theory in any number of dimensions. Using lattice techniques, we find an elementary proof that the decay of both, the entanglement of purification and reflected entropy, is enhanced with respect to the mutual information behaviour by a logarithm of the distance between the subregions. In the case of the Ising spin chain at criticality and the related free fermion conformal field theory, we compute also the overall coefficients numerically for the both quantities of interest.

In the present work, we will be concerned with CFTs in arbitrary number of dimensions emerging as a longdistance limit of lattice models regulated by a lattice spacing δ. We will be interested in entanglement for subsystems composed of two disjoint regions A and B, for which one often considers the mutual information (MI) defined as A quantity of significant recent interest in such a setup is also the entanglement of purification (EoP) [26], which can be regarded as a generalization of EE for bipartite mixed states. It requires purifying the reduced density FIG. 1. Illustration of our general setup for CFTs in two spacetime dimensions on a lattice: The mixed state ρAB on a subsystem of two disjoint regions AB separated by N d ≡ d δ sites is purified to a state with auxiliary factors A and B , taken to be of the same size NA ≡ wA/δ and NB ≡ wB/δ as A and B, respectively, where δ is the lattice spacing. Here we mostly consider NA = NB.
matrix ρ AB to a pure state |ψ in an enlarged Hilbert space on H AB → H AA BB such that ρ AB = tr A B |ψ ψ| (visualized in Fig. 1). The EoP is then defined as EoP is challenging to compute in QFT due to its inherent optimization procedure of finding a purification whose EE is minimal. Its current understanding in the intersection of quantum information and high-energy physics is based on Gaussian calculations [27][28][29], CFT techniques with a limited range of applicability [30][31][32], and on a conjectured realization in holography [33,34]. In the lat-ter case, EoP has been conjectured to be dual to the entanglement wedge cross section [35][36][37][38]. This led to many novel developments regarding the emergence of the gravitational hologram [39], see e.g., [40][41][42][43][44][45][46][47][48][49][50]. Another quantity closely related to EoP and also conjectured to be holographically dual to the entanglement wedge cross section is the reflected entropy RE [48,[51][52][53][54][55][56][57][58][59][60]. It is defined as EE  (3), which implies E P ≤ S R . RE is much easier to compute compared to EoP, as it does not require an optimization over all possible purifications. The aim of this letter is to to elucidate a particularly simple setting in which EoP and RE behave universally across CFTs, without relying on Gaussianity or Weyl rescalings. We achieve this by using spin chains and more general lattice models and focusing on universal inequalities satisfied by EE. We corroborate our studies using analytics and numerics in the Ising and free fermion CFTs [61], which allows us to extract prefactors in the asymptotic scaling of EoP and RE.
Setup. In our analysis, we will be concerned with CFTs on a lattice. Our general statements will be made in any number of dimensions, whereas our numerics will focus on CFTs in two spacetime dimensions.
The setting of interest will contain two spherical subregions of diameter w separated by a distance d. Fig. 1 illustrates it for CFTs in two spacetime dimensions in which case the subregions become intervals. At large distances d w 1, the decay of MI (2) in CFTs reads where and ∆ corresponds of the scaling dimension of the lowest non-trivial operator(s) in the theory, N denotes the possible degeneracy of such operators and the ellipsis denotes faster decaying terms [18,19,62,63]. The formula (5) assumes a gap in the spectrum of scaling dimensions and the lowest lying operator(s) being scalar(s). We will carry over this assumption in our studies of EoP and RE. Our aim is to find and prove an analogue of the scaling in (5) for EoP and RE. In the latter case, recent numerical studies in free CFTs in [56,57] led to the following fit where α is a positive model-dependent constant.
In this letter, we use quantum-many body techniques in conjunction with elementary EE inequalities to prove that the asymptotic form (7) holds both for EoP and RE in a general CFT amenable to a lattice realization.
Elementary proof of the large distance behaviour. To set up the general argument valid both for EoP and RE, we only need to assume that the density operator ρ AB of two subsystems A and B far away from each other takes the form (8) where the ellipsis denotes higher, not necessarily integer powers of ∆ and we do not make any assumptions about subsystem sizes. The ∆ term in (8) is needed to reproduce the power-law scaling of correlation functions involving insertions of the lowest lying scaling operator in both A and B. As we will show, the ρ (2) AB contribution turns out to not contribute to the leading order decay of EoP and RE.
In the following, we will regard (8) as originating from a perturbative purification where the product nature of the density matrix (8) for an infinite separation leads to Note that in our conventions, |ψ and, therefore, also |ψ (0) are normalized, which also leads to constraint for |ψ (j≥1) . The long distance behaviour of EoP and RE is determined by the small-∆ expansion of the eigenvalues µ j of ρ AA ≡ tr BB |ψ ψ| (11) via the definition of EE (1): The fact that ρ AB is a product state for ∆ = 0 implies that ρ AA ( ∆ = 0) is itself pure, see (10), and thus has eigenvalues µ 0 = 1 and µ j>0 = 0. This result gets modified at large but finite distances.
The linear correction to µ j vanishes, since we expect µ j to originate from a well-defined density matrix regardless of the sign of ∆ when viewed as a formal parameter. As a result, the possible leading behaviour of eigenvalues of ρ AA is given by (12) where Note that α j>0 ≥ 0 and if all of them vanished, the behaviour encapsulated by (12) would simply involve a higher-than-two power of ∆ . Let us consider now the asymptotics of EoP and RE resulting from (12). As we explained in the introduction, these quantities are given by S AA subject to additional conditions on purifications. For any purification leading to (12), S AA behaves as where and one sees as the leading order behaviour the structure (7) identified in fits to free CFTs RE numerics in [56,57] and the ellipsis denotes higher order terms in ∆ . Regardless of purification and long-distance limit, S AA is bounded from below [64] S AA (|ψ ) ≥ as was shown in equation (6) of [65]. Given (5), in order for (16) to be satisfied at large distances S AA cannot scale with a higher power than 2 ∆ . Since the eigenvalue analysis predicts this as the strongest possible power-law factor in the long-distance behaviour of S AA , α tot must be bigger than 0 and the behaviour predicted by (14) is necessarily the behaviour of both EoP and RE in any CFT with a gap in the operator spectrum and amenable to a lattice description.
As a corollary of this proof, from the definition of α tot in (13) we necessarily obtain that at least one of α j>0 > 0 and, as a result, the first subleading term encapsulated in (14) is also generically there. This is consistent with the findings of [56,57], which also identified such a contribution in RE for free CFTs on a lattice.
Finally, let us emphasize that our proof of the longdistance behaviour of EoP and RE did not rely on dimensionality of a CFT in question.
Properties of the overall coefficient. Our proof predicts only that the overall prefactor α tot is positive. It is possible, however, to extract more information about what ingredients affect the exact value of α tot using a rather general argument. To this end, notice that perhaps the easiest way to compute α tot is to extract it from where we suppressed higher order terms in ∆ . Starting with (9) and defining |ψ (i) we can write the reduction ρ AA as AA |), (19) which allows us to compute Tr(ρ 2 AA ) explicitly. Upon using the normalization condition ψ|ψ = 1 in the form of the following constraints where, in analogy with (18), |ψ (1) . Quite remarkably and as advertised below (8), the correction quadratic in ∆ to the density matrix ρ AB does not contribute to the leading order coefficient in the scaling of both EoP and RE. This looks like a potentially useful insight for any attempt to fix α tot for EoP and RE in terms of CFT data in an analogous manner to (5) for MI. Furthermore, obtaining the leading behaviour of the EoP amounts simply to minimizing a quadratic polynomial obtained from components of |ψ (1) subject to the constraint (20a) and the condition (with ρ (1) which generally leads to affine-linear constraints on |ψ (1) . The fact that the minimum exists follows from the argument presented in the previous section. While ρ AB does not affect the leading large distance behaviour encapsulated by α tot , individual α j>0 do depend on it and, via (14), so does the coefficient in front of the subleading term quadratic in ∆ .
Analysis in the critical Ising chain and free fermion CFT. So far, we have been completely general in our studies and in the following we will specialize to two closely related lattice models describing CFTs in two spacetime dimensions. This will allow us to obtain numerical values of the leading and first subleading coefficients in the behaviour of EoP and RE captured by (14) and, for RE, compare with earlier studies in [56].

Free Fermions (Gaussian)
Ising Spins (non-Gaussian) The Ising model realization of the c = 1 2 CFT on an infinite spatial line can be described by the critical lattice HamiltonianĤ more general forms of which we discuss in the Supplemental Material. TheŜ x,z i are spin operators defined by the Pauli matricesŜ x,z i = 1 2 σ x,z i . In the Ising CFT there is a non-degenerate (i.e., N = 1) lightest operator of scaling dimension ∆ = 1 /8, often denoted as the spin field σ and corresponding to aŜ x i lattice operator. The critical Ising model can be mapped to a free fermion theory, which is going to be another model in which we obtain numerical coefficients in EoP and RE. This formulation leads to two different notions of reduced density matrices for disjoint intervals, see, for example, [29,[66][67][68][69], and therefore provides in itself an independent example. For free fermion CFT there are two (N = 2) lowest lying operators with ∆ = 1 /2 and being simply the fermionic field operators.
Critical lattice models will describe CFT predictions for large enough sizes of subsystems at fixed w/d. Since we are dealing with purifications, which can lead to challenges when the relevant Hilbert space dimension becomes big, the key question is how big subsystems need to be to reproduce the continuum physics of interest. One hint comes from MI, for which we see that the continuum value of the prefactor in (5) is well attained at large distances already for w = 2 δ and 3 δ with the smallest subsystems of w = δ giving already reasonable predictions, see Fig. 2(a,b).
Given this encouraging result, we were in fact able to analytically compute the coefficients of the leading order MI, EoP and RE in the critical Ising model and for free fermions when w = δ. The results are summarized in Tab. I together with offsets and derivations can be found in Supplemental Material. Fig. 2(c-f) shows fits of our proven asymptotic formula (14) to fully numerical results, i.e., based on the full density matrix for disjoint intervals of both CFTs, we consider. We see strong indications of convergence to continuum values. In particular, for the critical Ising model EoP, the behaviour of α tot for w = δ corroborates our analytical prediction in Tab. I. Looking at the results for the largest atteinable values of w, we see that the leading fall-off coefficient changes from the analytic prediction at w = δ by only 2.6% and the subleading fall-off coefficient by only 3.9%. Generating data for w = 3 δ and above is numerically challenging as this requires computing very large matrices, slowing down calculations (see Supplemental Material).
For the EoP and RE of the Ising CFT, our results provide to the best of our knowledge new predictions, whereas for RE for the massless free fermions we find and display very good agreement with earlier studies in [56].
An interesting question regarding EoP concerns the dimensions of the enlarged Hilbert spaces. In our setup, when purifying the state of a system with N A + N B degrees of freedom by adding N A +N B additional ones (see also Fig. 1), there is a priori no constraint on N A , N B other than the basic requirement following from the definition of the Schmidt decomposition that N A + N B ≥ N A + N B . However, we show in the Supplemental Material that the choice of minimal purifications used so far yields the true minimum of EE as long as we choose N A = N A and N B = N B . This was already shown in [70] based on ideas of [71] for the EoP for Gaussian states (with Gaussian purifications).
Outlook. In this letter, we studied large distance behaviour of EoP and RE in a generic CFT with a gap in the operator spectrum for two spherical subsystems of diameter w in the large distance d limit. Using (8) in conjunction with elementary properties of EE we were able to prove that the large order behaviour of both EoP and RE is governed by (14). In comparison to the classic result (5) encapsulating large-distance behaviour of MI, EoP and RE get both enhanced by a logarithm of a separation. Subsequently, we explicitly calculated the large distance behaviour for both EoP and RE in one spatial dimension for the critical Ising model and massless fermions. This allowed us to establish the value of coefficients appearing in (14), see Tab. I and Fig. 2.
Our work opens a genuinely new avenue for studying EoP and RE in QFTs without restriction to free EoP results models. Perhaps the most interesting question concerns the dependence of the coefficients in the large order behaviour of EoP and RE on CFT data, akin to (5) for MI. An intermediate step could be to supplement our numerical code with large-distance reduced density matrices obtained with tensor networks for more complicated models, in particular determining model-dependent coefficients akin to (5). Optimizing over purifications outside the Gaussian realm inevitably leads to vast parameter spaces that quickly exhaust desktop-scale computational resources. However, the entanglement between reasonably-sized subsystems both mixed and purified is not large and it should be possible to represent purifications as manageable tensor networks, perhaps building on earlier works [34,72].
with spin operators represented by Pauli matrices σ α with α ∈ (x, y, z) bŷ We also use the identificationŜ α N +1 ≡Ŝ α 1 . This spin model can be converted to fermions by defining the 2N Majorana operators γ k via The Ising Hamiltonian then takes the form Here P is the total parity operator σ z ⊗N = k (−i γ 2k−1 γ 2k ). At the critical point J = h, the Hamiltonian thus simplifies tô Covariance matrix. For the critical ground state vector |0 which has a positive total parity, all correlations are encoded in the Majorana covariance matrix which in the infinite system size limit takes the form

Fermions Spins
FIG. S1. Subsystem setup of our analytical limits for fermions (top) with an inherent ordering and spins (bottom) without one. In both systems, we consider the subsystem AB consisting of two single sites A and B separated by d/δ sites.
The entropy of a Gaussian mixed state ρ with covariance matrix Ω is given by where ±iλ k are the eigenvalues of Ω. We will consider mixed states ρ AB or ρ AA , whose mixed state covariance matrices Ω AB or Ω AA result from restricting (S8) to the respective blocks.
Fermionic subsystem. We first compute reduced density matrices from the perspective of fermions, i.e., imposing an ordering between modes following from the anti-commuting variables γ i (see Fig. S1). A subsystem consisting of 1 + 1 sites (w = δ) separated by d/w = d/δ sites is then fully characterized by the restriction of the covariance matrix in (S8) and explicitly given by which corresponds to a lowest-dimension (Majorana) operator with scaling dimension ∆ = 1 /2. The associated fermionic density operator is then with respect to the basis (|↓↓ , |↑↓ , |↓↑ , |↑↑ ) and using D = 1 4 + 1 π + 1 π 2 , E = 1 4 − 1 π 2 , and F = 1 4 − 1 π + 1 π 2 . As in the main text, ∆ ≡ (w/d) 2∆ which here becomes 1 /2 = w/d. If we further restrict to a single site, we find the covariance matrix and density operator where ρ fer A is written with respect to the basis (|↓ , |↑ ).
Spin subsystem. We can perform a similar calculation in the original Ising spin system whose reduced density matrices can be constructed from the fermionic covariance matrix [68]. We need not repeat the single interval case, as entanglement entropies of connected regions are equivalent under a Jordan-Wigner transformation. However, we still need the reduced density matrix of a system of 1 + 1 sites in the large d limit, which we find to be with w/δ = 1 for the setup of 1+1 sites. As we are considering the spin Ising CFT, the lowest-dimension primary is the "order field" σ with scaling dimension ∆ = 1 /8. The constant C corresponds to the expectation value of an operator nonlocal in fermions, and can be computed from where M n is defined as the n × n matrix (S15) Using this construction, one finds [73] C = e 3ζ (−1) 2 23/12 ≈ 0.1612506 . (S16) MI for fermions. We now begin computing entanglement measures for the small subsystems whose reduced density matrices we just obtained explicit expressions for. The continuum limit corresponding to the Ising CFT is obtained by keeping d/w (or, equivalently, ∆ ) fixed and taking δ/w to 0. We will see that taking only a few lattice sites is sufficient to describe the qualitative and approximate quantitative behavior of the continuum limit. To demonstrate this, we now show that the large distance asymptotics of the MI (2) for the case of 1 + 1 sites yield results close to the continuum formula (5). In order to compute the MI, we need to determine the von Neumann entropy of a single site S A = S B and of both sites S AB . These entropies can be computed from the eigenvalues of the covariance matrix Ω associated to the respective Gaussian state ρ. As an antisymmetric matrix, Ω fer AB has pairs of purely imaginary eigenvalues ±i λ k , from which applying (S9) leads to where the eigenvalues of Ω fer AB are to leading order We can similarly expand S 1+1 ≡ S AB at large d, which results in a MI for w = δ of This reproduces the correct continuum power law of fermionic MI, but yields a coefficient lower than the continuum value (5) which also matches the large-distance expansion of earlier results for Dirac fermions [74] for two blocks of general width w.
MI for spins. We compute EE for the spin system directly from the eigenvalue spectrum of the reduced density matrix (S13). Its four eigenvalues µ j are where C is given by (S16) and from which we can directly compute the EE for AB via Note that in the following we will denote eigenvalues of any density matrix by µ j . This analysis leads to the Ising model prediction for spin MI at w = δ and large separations d of the form Comparing this formula with the CFT analytics (5) for ∆ = 1 /8, we see an exact match in the power-law behavior. Furthermore, the prefactor in (S24) is only 3.6% off from the continuum value ≈ 0.309 predicted by (5).
EoP for fermions. Analogous to the MI calculation for free fermions, we now calculate the EoP in the fermionic subsystem of two sites separated by d/δ sites, expressing all calculations in terms of covariance matrices. We purify Ω AB in the limit d/δ → ∞ as associated to systems (A, B, A , B ) with G = 2 π and L = √ 1 − G 2 , whose EE S AA is zero and we thus have lim d→∞ E P = 0, i.e., the EoP vanishes for large d/δ, as expected.
In order to find the asymptotic behavior of E P , we need to study the variation of the symplectic eigenvalues ±iλ i of Ω AA when perturbing Ω according to The requirement of Ω representing a purification implies Ω 2 = −1, which induces the constraints We further require that the restrictions Ω (1) AB and Ω (2) AB matches the ones of (S10) expanded in 1 /2 , i.e., The equations (S27) and (S28) can be solved iteratively up to some free variables. We first solve Ω (1) in terms of Ω (0) and then Ω (2) in terms of Ω (0) and Ω (1) .
In order to find asymptotics of the symplectic eigenvalues λ i , we can use the fact that Tr(Ω 2 AA ) = −2(λ 2 1 + λ 2 2 ) and Tr(Ω 4 AA ) = 2(λ 4 1 + λ 4 2 ) to solve for the asymptotics of λ i to be given by where α tot will depend on some of the free parameters contained in Ω (1) and Ω (2) . With this trick, one finds where the variables x ij represent unconstrained entries in the block Ω AB,A B . In order to find the asymptotics of EoP, we need to minimize α tot over these parameters to find the smallest possible EE S AA . Due to the fact that (S30) is quadratic in x ij , we can calculate this valua analytically as Expanding S AA ∼ i (log 2 − λi 2 ) through λ i up to second order in 1 /2 based on (S29) allows us to also find the offset analytically, namely we have Combining this with the result from (S31) gives which agrees with the form (14) in the main text. Note that the simplicity of Gaussian states allowed us to even find the analytical form of the constant offset. The accuracy of this analytical prediction was tested numerically, for which we presented the results in Fig. 2 in the main text.
EoP for spins. In the limit of an infinite distance between the two single site subsystems, we purify (S13) by the state |ψ (0) with Schmidt decomposition where the convention for factors ordering in the purification is ABA B . Note that in this analysis we assume that a minimal purification from two to four spin degrees of freedom suffices and we will subsequently provide supporting numerical evidence and an additional discussion. Moving on, we supplement this purification with finite distance corrections up to second order in 1 /8 as |ψ ∼ |ψ (0) + 1 /8 |ψ (1) + 1 2 2 1 /8 |ψ (2) . (S35) We will optimize over |ψ (1) and |ψ (2) subject to the normalization constraint ψ|ψ = 1 order by order in which follows from (S13). We expand the first order perturbation as where z i = x i + iy i and |φ i is the basis of H ABA B ordered as (|↓↓↓↓ , |↑↓↓↓ , |↓↑↓↓ , |↑↑↓↓ , . . . , |↑↑↑↑ ). We then need to implement the condition (22) in the main text based on (S36) together with the normalization constraint (20a) in the main text. We solve these affine linear constraints by the replacements We can then compute α tot according to (21) in the main text as quadratic polynomial in terms of the remaining free variables z i which leads to the rather involved expression In order to find the EoP, we need to minimize over the z i to find the smallest possible value α tot , which can be done analytically and leads to The non-vanishing α tot shows that the resulting EoP obtained from (14) again has the form which, as in the fermion case, exhibits a leading-order long-distance behavior enhanced with respect to that of MI (S24) by a logarithm of the distance.
When it comes to the subleading long-distance behavior encapsulated by j>0 α j (1 − log α j ) , we would need to extract the individual α j and optimize over the remaining parameters. While it is plausible this can be also done analytically, we determined the value quoted above numerically, as discussed in the main text. Note that for the free fermion case with w = δ considered above, we determined this term analytically in terms of α tot .

RE for fermions.
In the Gaussian case of free fermions, our starting point is the following perturbative expansion of the reduced density matrix ρ AB of a system of 1 + 1 fermions in the large d separation, AB given by (S11). We similarly construct the canonical purification of (S11) via Note that in contrast with fermionic MI and EoP, we do not need to phrase our computation of RE in terms of the covariance matrix formalism since we can construct the canonical purification | √ ρ AB exactly for the given form of the initial reduced density matrix ρ AB . In this case, the first-order perturbation |ψ (1) is simply given by with the same ordering of the basis |φ i as in the previous case. From the canonical purification's density matrix ρ := | √ ρ AB √ ρ AB | we consider a restriction to subsystems AA given by the reduced density matrix ρ AA = tr BB (ρ) which has the perturbative expansion AA explicitly given by 4π(π 2 −4) , andJ = 2 1/2 4π 2 . We once again compute the trace of the square of (S42) according to (21) in the main text from which we obtain which also shows that the reflected entropy S R (ρ AB ) = S AA (ρ) of the fermionic subsystem also exhibits a logarithmic enhancement of the power law decay for w = δ given by where we also computed the constant term in (S44) from the eigenvalues of (S42) according to (21) in the main text.
RE for spins. For the Ising spin case, we now describe the detailed computation of the reflected entropy RE for w = δ in the large d limit just as for fermions. The reduced density matrix for a spin system of 1 + 1 sites in the large d limit can again be computed according to (8), AB + . . ., yielding (S13). We now construct the canonical purification of (S13) via | √ ρ AB = i √ e i |e i ⊗ |e i = |ψ (0) + 1 /8 |ψ (1) for ρ AB |e i = e i |e i and where the eigenvalues e i are defined in (S22a). In this case, the first order perturbation |ψ (1) is given by where the states |φ i = |φ i ABA B form an orthonormal basis for the purified Hilbert space H ABA B ordered as (|↓↓↓↓ , |↑↓↓↓ , |↓↑↓↓ , |↑↑↓↓ , . . . , |↑↑↑↑ ). From the canonical purification's density matrix ρ := | √ ρ AB √ ρ AB | we consider a restriction to subsystems AA given by the reduced density matrix ρ AA = tr BB (ρ) which has the perturbative expansion ρ AA = tr BB (|ψ (0) ψ (0) |) + 2 1 /8 (2tr BB (|ψ (1) ψ (1) |))/2 explicitly given by , where the coefficient C is defined as in (S16). From here we follow the strategy of the main text and compute the trace of the the square of (S46) according to (21). In this case, we find a value of α tot computed via (21) to be As a consequence, the large d leading behaviour of the reflected entropy S R (ρ AB ) := S AA (ρ) exhibits a non trivial logarithmic enhancement of the power law decay according to (14) and where the constant contribution can be computed from the eigenvalues of (S46) leading to a reflected entropy S R of the Ising subsystem for w = δ of The constant term is again determined numerically in the main text. Numerical approach and asymmetric purifications. Our numerical methods are based on [29,70,71], which outline the construction of an efficient algorithm for local optimization over Gaussian states, based on a gradient descent approach exploiting the natural Lie group parametrization of the state manifolds. Our numerical results are obtained using an adaptation of this algorithm to the non-Gaussian case of interest.
To compute the EoP as given in (3), we minimise EE S over the manifold M of purified state density matrices. We first purify our initial mixed density matrix to a 2 N -dimensional pure ρ 1 via the Schmidt decomposition. Here, N = X N X with N X denoting the physical degrees of freedom in subsystem X. We parametrize elements ρ U ∈ M by transformations U = 1 ⊗ U with U ∈ U(2 N A +N B ), so that ρ U = U ρ 1 U −1 . The tensor product signifies that U only acts non-trivially on degrees of freedom in A and B . We then optimize by performing iterative steps along directions in M which locally minimize S AA [29,70], Here, K n = µ F µ (U n )Ξ µ /||F|| 2 and F µ : M → R is the gradient descent vector field with {Ξ µ } as basis of u(2 N A +N B ). We choose U 0 = 1 and we pick 0 < t < 1 in such a way that the value of S AA decreases with successive steps. The {Ξ µ } span the tangent space at U = 1 and, due to the left-invariance of the Riemannian metric on M, form orthonormal bases for the tangent spaces at all other points in M, too, where Ξ µ is identified with the tangent vector to the curve γ(s) = U e sΞµ at γ(0) [70].This saves us having to re-evaluate the matrix representation of the metric at each step, as we would have to if we had chosen a coordinate parametrisation of M. While this makes our algorithm more efficient than a naive gradient descent, the numerically accessible range is still highly limited: since N A +N B ≥ N A +N B , the dimension of M is at least dimu(2 N A +N B ) = 2 2N A +2N B − 1 and (S49) requires exponentiation of at least 2 (N A +N B ) × 2 (N A +N B ) matrices, with a typical step count of several hundred. This becomes extremely slow on a powerful desktop computer for N A +N B ≥ 5. For the symmetric purifications in the main text this corresponds with w > 2 δ, which explains the regime we were able to explore.
Given this limitation on our numerical capabilities, it is instructive to ask whether an optimization over minimal purifications corresponding with N A + N B = N A + N B yields the true minimum of EE -not least because for large systems this becomes the only numerically viable choice. A natural follow-up question is whether among the choices of minimal purifications, the intuitive choice of N A = N A and N B = N B suffices to reach the true minimum defined as EoP. More pertinently, we might ask whether it is even possible to reach the true minimum with a minimal purification for which N A = N A and N B = N B . In [70], a combination of numerical and analytical evidence was provided to show that the answer to this question is in affirmative for Gaussian states. While limited by the greater numerical challenge in the non-Gaussian case, we present similar numerical evidence in Table S1 to show that the same may be said for our model: the true minimum can only be reached if N A ≥ N A and N B ≥ N B , which indicates that the lowest-dimensional purification for which the EoP can be obtained is the minimal purification with N A = N A and N B = N B . * hugo.camargo@aei.mpg.de † lucas.hackl@unimelb.edu.au ‡ michal.p.heller@aei.mpg.de; On leave of absence from: National Centre for Nuclear Research, Pasteura 7, 02-093 Warsaw, Poland § a.jahn@fu-berlin.de ¶ bennet.windt17@imperial.ac.uk