Entanglement of Purification in Many Body Systems and Symmetry Breaking

We study the entanglement of purification (EoP), which measures total correlation between two subsystems $A$ and $B,$ for the free scalar field theory on a lattice and the transverse-field Ising model. We numerically compute the EoP when the subsystems $A$ and $B$ are of the same size and find interesting properties which are common to both of these models. First, we find that the EoP becomes a non-monotonic function of the distance $d$ between $A$ and $B$ when the total number of lattice sites is small. When it is large, the EoP becomes monotonic and shows a plateau-like behavior at small $d$. Next, we find that the $Z_2$ reflection symmetry, which exchanges $A$ and $B$, can get broken in a purified system in optimal purifications. Moreover, in the ferromagnetic phase of the Ising model with transverse magnetic field, the above-mentioned $Z_2$ reflection symmetry is also spontaneously broken. Finally, we interpret our results in terms of the interplay between classical and quantum correlations.


Introduction
For pure states, entanglement entropy (EE) is essentially a unique measure of quantum entanglement [1]. When we decompose the total quantum system into two subsystems A and B, the EE is defined as S A = −Tr[ρ A log ρ A ]. The reduced density matrix is ρ A ≡ Tr B |Ψ AB Ψ| AB , where |Ψ AB describes a pure state in the total system. The EE helps us to extract essential properties of quantum field theories [2,3], especially conformal field theories (CFTs) [4]. Since it also possesses a simple geometrical interpretation in gravity [5,6], it has recently played an important role in the context of the holographic Anti de-Sitter space/conformal field theory (AdS/CFT) correspondence [7].
If the total state is not pure, EE is neither a good measure of quantum entanglement nor classical correlations. There are several excellent measures of quantum entanglement for mixed states, such as entanglement of formation and squashed entanglement, which coincide with the EE when the total system is pure. Common to the definitions of these quantities is a required minimization over infinitely many quantum states. Such minimizations are very challenging in quantum field theories and almost no results have been known.
The purpose of this letter is to analyze a quantity called the entanglement of purification [8,9], which is a close cousin of such entanglement measures and also involves a minimization procedure. The entanglement of purification (EoP), written as E P (ρ AB ), is a measure of total correlation between two subsystems A and B for a mixed state ρ AB and is defined as follows. Consider a purification of a given ρ AB , denoted as |Ψ AÃBB , which lives in an enlarged Hilbert space Among infinitely many different choices of the purifications |Ψ AÃBB , the entanglement of purification is defined by minimizing the EE S AÃ as The EoP is a measure of total correlation, not only quantum entanglement, in the sense that it vanishes only for product states and monotonically decreases under local operations. Remarkably, in spite of that, its regularization has an operational meaning as a number of EPR pairs needed to create a given state under LOq procedure [8]. Moreover, a geometric formula for EoP, based on AdS/CFT, was conjectured in [10,11] and has been actively studied . This formula also motivates us to study field theoretic computations of EoP. An earlier preliminary analysis of EoP in a free scalar field theory was done in [34] for very small subsystems. Refer also to [35] for a conformal field theoretic approach, which matches with the holographic formula for specific examples. In this letter, we will present numerical calculations of EoP both in a free scalar field theory for a wider range of subsystem sizes assuming a Gaussian wave functional ansatz and in the transverse-field Ising spin chain. In our Ising model analysis, when the sizes of the subsystems are minimal, our calculation of EoP gives an exact answer without assuming any ansatz. In both models, we will observe intriguing non-monotonic behavior of EoP as a function of the distance d between A and B. Moreover, we will find a spontaneous breaking of the Z 2 reflection symmetry which exchangesÃ andB (for a system symmetric under A, B reflection), which was not taken into account in the earlier paper [34].
2. EoP in free scalar field theory Consider a free massive scalar field theory in 1 + 1 dimensions, defined by the Hamiltonian As in [2,34,36], the ground state wave function Ψ 0 of this lattice scalar theory is given by the Gaussian function The matrix W is defined by where N is the number of total lattice sites. The parameter a is the lattice spacing and we set a = 1 below by rescaling m. Notice that W is a symmetric and realvalued matrix. We consider masses between m = 10 −1 and m = 10 −4 near the conformal (massless) limit (refer to Fig. 1 for a sketch of our lattice setup for N = 20). We divide the total Hilbert space into three parts We denote the number of lattice sites in A, B, · · · by |A|, |B|, · · · and the distance between A and B by d (see Fig. 1). We rewrite (4) as with the sub-matrices P, Q, R following from (5) after determining the (relative) position of the subsystems A and B.
From this wave functional, we can directly calculate the mutual information I(A : B) = S A + S B − S AB and the logarithmic negativity E N (ρ AB ), both of which are shown in Fig. 2. The latter is a useful probe of quantum entanglement between A and B, defined as E N (ρ AB ) = log Tr|ρ Γ B AB | [37], where ρ Γ B AB is the partial transposition with respect to B. Refer to appendix A for the details of computing E N (ρ AB ). We observe that E(ρ AB ) takes the largest value at d = 0 and for d ≥ 1 shows exponential decay. On the other hand, the mutual information slowly decreases as function of d, following a power-law. To calculate the EoP, we purify the system by adding auxiliary subsystemsÃ andB. Assuming the purified wave functional is Gaussian, we restrict it to the form where we have decomposed the Gaussian coupling matrix V into three sub-matrices J, K, L. The condition (1) requires J = P . Furthermore, assuming subsystems of equal width w = |A| = |B|, L becomes a w × w square matrix and is related to K by the equation Use of a symmetry transformation [34] allows the simplification of the K to the form: Thus, all free parameters are contained in the w × w matrices K A,B and K B,Ã . If one assumes a Z 2 symmetry which reflects AÃ and BB, we will have K A,B = K R B,Ã , where we define M R of a matrix M as the inverse ordering of all rows and columns, i.e.
To quantify the Z 2 symmetry breaking, we also define the Z 2 asymmetry where ||M || 2 is the 2-norm over all entries of M . In Appendix B we show that E P is Z 2 -invariant. The entanglement entropy S AÃ can be computed from the eigenvalue spectrum {λ k } of the matrix Λ ≡ −V −1 AÃ,BB · V BB,AÃ [2] as follows: (12) The EoP, E P (ρ AB ), is the minimum of S AÃ over all purifications Ψ AÃBB [φ], i.e. under varying the components of K A,B and K B,Ã . We calculated this minimum for subsystem sizes w = 1, 2, 3, 4 and studied its dependence on the distance d. The results were computed using a numerical L-BFGS optimization implemented with the C++ package dlib [38].
As the Z 2 reflection symmetry is a property of the original system (ρ AB ) and leaves the entanglement of purification invariant, it might be natural to assume, as in [34], that the optimal purification is Z 2 -symmetric. Surprisingly, we observe that finding the true minimum of S AÃ requires us to enlarge the parameter space by breaking the Z 2 exchange symmetry betweenÃ andB. The results for N = 60 are shown in Fig. 3 (for the result assuming Z 2 symmetry, refer to Fig. 13 in Appendix B). From d = 0 to d = 1, we observe a plateau-like behavior of E P at large w whose width appears indendent of w, suggesting a finite-size effect. This notion is supported by the form of K for minimal S AÃ , shown in Fig. 4 for m = 10 −4 and w = 4. At d = 0 and d = 1, the couplings between A 4 ↔B 1 andÃ 4 ↔ B 1 are enhanced, implying additional short-range entanglement. The Z 2 symmetry breaking, while hard to discern from Fig. 4, clearly appears at d = 1 when considering the Z 2 asymmetry parameter A we defined in (11), shown in Fig. 6. Within numerical accuracy, A = 0 for any d = 1. Note that the symmetry breaking becomes more pronounced with increasing w. At w = 1, it is not observable within numerical accuracy, while is clearly visible for the w = 4 data in Fig. 6.
Furthermore, we find that the EoP does not need to be a monotonically decreasing function of d at small N , as shown in Fig. 5 (left) for w = 2. As we increase N , this non-monotonicity gradually disappears and we find a plateau at large w (Fig. 5, right). While one may be inclined to expect the EoP, being a correlation measure, to decrease monotonically with the distance between A and B, in the next section we will conclusively show that it is not a monotonic function. The same behavior can also appear in spin chains.
For small masses m and block widths 1 w N , we observe the scale invariance E P (w, d, m) = E P (nw, nd, m) with n ∈ N. As evident from Fig. 7, we find E P (d, w, m 1) = a(m) (d/w) −p(m) in the continuous limit, with positive scaling coefficients a(m) and and m = 10 −4 (right), no Z2 symmetry being assumed. (9)) for minimal entanglement of purification between physical sites AB and auxiliary sitesÃB for mass m = 10 −4 , block width w = |A| = |B| = 4, N = 60 and various block distances d.
p(m). We find that a(m) diverges logarithmically as m → 0, while p(m) converges sublinearly to zero.

EoP in the transverse-field Ising model
We will now compute the EoP for spin systems. First we review a method to compute the EoP in a finite dimensional system. Let us denote Hilbert space dimension by D such that D A = dim H A etc. In general, the dimension of auxiliary Hilbert space DÃDB should be at least as large as rank ρ AB to purify a given mixed state ρ AB , and there is no upper bound for this in general. Fortunately, however, the true minimum of S AÃ can be found for DÃ, DB ≤ rank ρ AB in a system with finite-dimensional Hilbert space [39], which enables us to compute EoP in practice. For convenience, we call the  purification minimal when DÃDB = rank ρ AB , and maximal when DÃDB = (rank ρ AB ) 2 . A simple example of purification is the thermofield double state (TFD) where we diagonalized the density matrix such that ρ AB = i p i |i i| AB with i p i = 1, p i ≥ 0, and the ancilla state ρÃB is identical to the original one ρ AB . The terminology thermofield double arises from the fact that the modular Hamiltonian K AB = − log ρ AB can be identified with the thermal Hamiltonian with inverse temperature β = 1.
In general, all possible purifications of a fixed dimension can be obtained by acting with unitary operators on the auxiliary systems |ψ(U ) ABÃB = I AB ⊗ UÃB |ψ 0 ABÃB , where |ψ 0 ABÃB is an initial state. We also vary the dimensions of the auxiliary systems DÃ, DB to achieve both the maximal and minimal purification. In principle, the maximal purifications are needed to obtain the EoP. However, we will find that often the minimal purification is sufficient to find the true minimum of S AÃ .
We should remark that our numerical minimization relies on a method which is a variation of the steepest descent method. Results obtained from this are only guaranteed to converge to a local minimum of the objective function. To obtain the global minimum, we start from several (random) initial purifications and ensure that the same point of convergence is reached. Because of the nature of this technique, the existence of narrow local optimal valleys cannot be excluded, in which case the numerical results only provide an upper bound to the true EoP. A similar caveat exists for the scalar field theory results.
We deal with a 1D transverse-field Ising model where i, j denotes the summation over nearest neighbors with periodic boundary condition and N is the number of total spin sites. First, we focus on the ground state of the critical Ising model (h = 1) on N sites. We compute the EoP for the subsystems |A| = |B| = 1 as a function of the distance d (see Fig 8). While the optimization is performed for the maximal set of purifications, the optimal purification always corresponds to the minimal purification for this case.
Remarkably, when N = 4, one can rigorously show that the EoP is a non-monotonic function of the physical distance d between A and B. Namely, E P must coincide with S A at d = 1 because ρ AB has support only on a symmetric subspace [40], while E P < S A at d = 0 (and d = 2) clearly follows from the numerical computation (refer to the appendix C for details). This provides us with an exact example in which EoP behaves nonmonotonically (Fig. 8, right), which puts it in contrast to other correlation measures.
The Z 2 reflection symmetry, which exchanges AÃ and BB, is explicitly broken at d = 1. This is obvious as SÃ = SB as show in the Fig. 9. As in the scalar case, the Z 2 symmetry breaking leads to two degenerate configurations for E p , related by an AÃ ↔ BB flip. Moreover, S A = SÃ for any d indicates that the optimal purifications are different from TFD purifications. Note that, when the subsystem size is minimal, the Z 2 symmetry is gradually recovered as N gets larger.
We also compute the EoP with larger subsystem size |A| = |B| = 2. In this case, the optimization was done within minimal purifications to expedite the computation. Interestingly, we observe a non-monotonic behavior of EoP, which weakens as N increases (Fig. 10). Note that the Z 2 symmetry breaking is also found at d = 1, which remains even in the large N limit.
Both a plateau and a Z 2 symmetry breaking occur also for a class of two-qubit states called Werner states, which coincide with the ground state of the Heisenberg spin chain. For details, refer to the appendix C.

Phase transition and Z 2 symmetry breaking in the Ising model
Furthermore, we compute the EoP as a function of the magnetic field h for the nearest-neighbor minimum subsystems in the thermodynamic limit N → ∞. We consider the whole system to be in the thermal ground state for which the analytic form of the reduced density matrices is obtained [41,42]. The result in Fig. 11 shows that the EoP has an inflection point at h = 1. It indicates that the EoP correctly captures the phase transition of the original system. Remarkably, the Z 2 reflection symmetry of AÃ and BB gets broken only in the ferromagnetic phase h < 1. However, the thermal ground state maintains a Z 2 flip symmetry in any phase.

Conclusions and Discussion
Finally, we seek to provide an interpretation of our EoP results. For both free scalar theory and the critical Ising model we observed a plateau-like behavior of the FIG. 12. A simplified toy model for EoP for 1D many body systems, assuming only short-range quantum entanglement (zigzag lines). The original system |Ψ ABC is shown above, the purification |Ψ AÃBB below. We set w = |Ã| = |B| = 3.
EoP for small d. When the size of the total system N is small, the EoP even has a peak at d = 1. Though this non-monotonic behavior disappears when N is large, we still find its remnant as a plateau-like behavior. These behaviors are very special to the EoP and cannot be found in the mutual information (MI) I(A : B), which is another measure of total correlation. To interpret the difference between EoP and MI, we consider a split of total correlation into quantum entanglement and classical correlations. We assume (half of) MI evaluates them equivalently. In contrast, we assert that the EoP captures them differently: the EoP enhances classical correlations, while treating the quantum entanglement just as the MI does. Indeed, for purely quantum correlations (i.e. pure states), they coincide with each other: E P = I/2. On the other hand, for purely classical correlations (i.e. separable states), the EoP actually enhances this correlation at least twofold: E P ≥ 2(I(A : B)/2) [8]. Based on the above observation, we provide an explanation for the non-monotonic behavior of the EoP as follows: Since the quantum entanglement can be estimated by the log negativity which falls off quickly, we can say that the correlations at d = 0 mainly come from quantum entanglement, while those at d ≥ 1 arise from classical correlations and thus are enhanced in the EoP. Indeed, for d 1, the EoP is at least twice as large as half of the MI, and both MI and EoP monotonically decrease following a similar power law. At d = 0, however, the EoP may lose this enhancement of classical correlation, which results in a plateau or non-monotonicity of EoP. Possible connections to analogous quantities such as discord [43,44] may be an interesting future work.
These behaviors can be approximated by a toy model of a purified system, depicted in Fig. 12, which also illustrates the Z 2 symmetry breaking at d = 1. We take into account only nearest-neighbor quantum entanglement. At d = 0, two boundary sites of A and B are strongly entangled, and this entanglement remains after the purification, unlike classical correlations at d > 0. At d = 1, A and B are separated by one extra site C. As C is strongly entangled with both A and B, tracing it out will lead to a highly mixed state and the direct quantum entanglement between A and B is very small. Therefore, the optimal purification requires strong entanglement for A ↔Ã, B ↔B,Ã ↔B,Ã ↔ B and A ↔B to become pure, while the entanglement A ↔ B should be negligible. The result of this complicated competition is a Z 2 reflection symmetry breaking, where only eitherÃ ↔ B or A ↔B exhibit strong entanglement (Fig. 12, center), which we confirmed both for the free scalar field theory and the Ising model. This implies that the Z 2 symmetry breaking occurs when classical correlations of ρ AB are strong compared to quantum entanglement, which is what we observe at d = 1. For d ≥ 2, however, A and B are separated by multiple sites and quantum entanglement between A and B is negligible. The EoP decreases along with the small remaining classical correlations as d increases.
In our analysis of the EoP for the transverse-field Ising model, we found that the Z 2 -broken region coincides with the ferromagnetic phase. This suggests an interesting connection between a symmetry breaking in the optimal purification for the EoP and a quantum phase transition. This deserves future studies.
Acknowledgements We thank Pawel Caputa, Horacio Casini, Jens Eisert, Masamichi Miyaji, Masahiro Nozaki, Kazuma Shimizu, and Brian Swingle for useful conversations. We are very grateful to Yoshifumi Nakata for valuable comments on the draft. AB and KU are supported by JSPS fellowships. AJ is supported by a Studienstiftung fellowship. TT is supported by the Simons Foundation through the "It from Qubit" collaboration. One simple characterization of quantum entanglement between subsystems A and B for a mixed state ρ AB is the logarithmic negativity [37]. For this we introduce the socalled partial transposition Γ B , which is the transposition acting only for the subsystem B. It is well-known that for separable states, the partially transposed density matrix ρ Γ B AB is still positive, while for non-separable (=entangled) states , this positivity is not preserved in general. We would like to note that even if ρ Γ B AB is positive, we cannot say AB is separable, while the converse statement is true.
The logarithmic negativity is defined by where we introduced If we write the eigenvalues of ρ Γ B AB as λ i , then we can write Note here that since ρ AB and ρ Γ B AB are both normalized, we have i λ = 1 and thus i |λ| > 1. The logarithmic negativity is vanishing if and only if all the eigenvalues λ i are non-negative. This quantity is known to be monotonic under LOCC and satisfies at least the minimal property of an entanglement measure for mixed states. Also note that when the total state ρ AB is pure, E N (ρ AB ) is not equal to the (von Neumann) entanglement entropy but is equal to the n = 1/2 Rényi entropy, defined by S (1/2) A = 2 log Tr(ρ A ) 1/2 . Now we compute logarithmic negativity for the ground state Ψ 0 for our free scalar lattice model. We divide the total lattice system into subregions A, B and C such that H tot = H A ⊗ H B ⊗ H C . We define their lattice sizes to be |A|, |B| and |C|. In this setup, we wish to compute the logarithmic negativity which measures the quantum entanglement between A and B. First remember that the ground state wave functional is given by (6). Then the reduced density matrix ρ AB = Tr C [|Ψ ABC Ψ ABC |] is obtained by integrating out C: where M and N are symmetric real matrices. They are defined as For later purpose, it is useful to decompose M and N , which are (|A| + |B|) × (|A| + |B|) matrices, into |A| × |A|,|A| × |B| and |B| × |B| matrices as follows: where T is the standard transposition. Given this density matrix we now proceed to compute negativity. For that, we have to first perform the partial transpose Γ B , which is equivalent to interchanging φ B and φ B . After we rearrange this as a matrix whose arguments are of the form (φ A , φ B ), (φ A , φ B ), we obtain: Now we can perform a field redefinition: where V is a orthogonal matrix and M D is a diagonal matrix. We choose them such that we have We apply the same transformation on φ AB . Then we find To diagonalizeÑ we perform another transformation, where S is another orthogonal matrix. Finally, up to a normalization factor we have Here µ i are the eigenvalues of the matrixÑ , equivalently the eigenvalues of the matrix M −1Ñ . This is because we can writeÑ asÑ = ( Once we numerically obtain these eigenvalues µ i we can calculate the logarithmic negativity in a similar way to the entanglement entropy in [2,36]. As a toy model, consider a scalar φ in quantum mechanics with the density matrix We can diagonalize (30) and find the eigenvalues where λ is defined by Thus we obtain which is non-negative and is positive when µ is negative. Now notice that our ρ Γ B AB given by (28) can be regarded as |A| + |B| copies of this kind of quantum mechanics. Thus, finally, we can evaluate the logarithmic negativity as follows: where We first show that the entanglement of purification E P is invariant under a Z 2 symmetry transformation. In terms of the w × w matrix K determining the coupling between the AB andÃB systems, this symmetry is expressed as: For the w×w matrix L determining the coupling withinÃB, we use (8) to find where we used S 2w QR −1 Q T S 2w = QR −1 Q T , as the initial system A and B is symmetric under the Z 2 symmetry. For the same reason, J → S 2w JS 2w . As a result, the matrix Λ from whose eigenvalue spectrum we compute E P becomes This is merely a similarity transform (since S = S −1 ) and a transpose, which do not affect the eigenvalues. Hence, for any purification in ABÃB that is not itself Z 2invariant, there exists another purification with identical EoP that is produced by acting with the Z 2 symmetry. A study of Z 2 -invariant purifications for w = 1, 2 at m 1 was already considered in [34]. Extending to w = 3, 4 yields the data shown in Fig. 13. We observe a peculiar peak of E P at d = 1 appearing at larger w and a power-law decay at d > 1. As discussed in the main text, the peak disappears after relaxing the Z 2 constraint.

Details of Scaling Properties
In the conformal limit m 1 and at distances d 0 larger than the lattice spacing, we observe a power law scaling of the entanglement of purification E P :  shown along with their EoP counterparts in Fig. 14. We observe that EoP decays slower than the mutual information in this range.

Details of Numerical Calculation of EoP
Let us briefly review the EoP in a finite-dimensional system. Given a bipartite state ρ AB , any purification of ρ AB can be created from an arbitrary initial purification |ψ 0 by acting with local unitary operations on the auxiliary system: As an initial purification, one may use the standard purification [8] or the TFD purification (13). Note that we can regard |ψ 0 as a vector in higher-dimensional purification space, especially for the maximal purification DÃB = (rankρ AB ) 2 , without loss of generality. Therefore, the minimization of EoP over all possible purifications can be equivalently expressed in terms of unitary operators on auxiliary systems: (42) where the minimization is also taken over all possible divisions of the ancilla Hilbert space into HÃ and HB, imposing DÃDB ≤ (rankρ AB ) 2 . Note that the optimal purification has a trivial redundancy since S AÃ is invariant under any local unitary UÃB = UÃ ⊗ UB. In other words, SÃ = SB indicates a non-trivial degeneracy of the optimal purification.

Werner state
An interesting type of quantum state is the Werner state on 2 qubits system where p ∈ [0, 1] is a parameter of states, P sym and P asym are the projections onto the (anti-)symmetric subspace in H AB (45) in {|00 , |01 , |10 , |11 } basis, I is the 4 × 4 identity matrix, and |Bell := 1 √ 2 (|01 − |10 ). The Werner state is also related to an isotropic state, which can appear as the ground state of the antiferromagnetic Heisenberg model. The EoP of the Werner state has already been calculated numerically in [8,45] (using a slightly different definition |Bell := 1 √ 2 (|00 + |11 ), which does not change the correlation). Here, we computed the EoP for Werner state again and found more fine-grained phase structures related to the Z 2 symmetry breaking.
The result is shown in Fig. 15. There are 4 different regimes classified by a configuration of optimal purifications: In the phase (a), we found that there are two choices (non-equivalent) for the optimal purifications: DÃ = 2, DB = 3 and DÃ = DB = 3 which produces the same results for EoP up to certain the numerical accuracy. Indeed, they give the same value for S AÃ (after minimization) up to 10 digits around the transition point. For the results shown in the figure, we have used DÃ = 2, DB = 3. In the phase (d), it was observed that the EoP is strictly smaller than S A (except p = 1), which was missed in the previous works. Now let us focus on the optimal purifications around the transition points p 1 0.319 and p 2 0.401. SÃ, SB and I(Ã :B) around this region are shown in Fig. 17. It is clear that the optimal purification breaks the Z 2 reflection symmetry which exchangesÃ andB in the phase (a). It is analogous to what we found in the free scalar field theory and in the critical Ising spin chain. Note that SÃ = SB holds for phase (b) and (c), contrary to the Ising model. We also find that Z 2 reflection symmetry is similarly broken in the phase (d). The entanglement of purification for the antiferromagnetic isotropic Heisenberg model (left) and for a chaotic spin chain (right).

Heisenberg model
Let us consider the anti-ferromagnetic isotropic Heisenberg model H Heisenberg = i,j (σ x i ⊗ σ x j + σ y i ⊗ σ y j + σ z i ⊗ σ z j ). (47) For even N , the reduced density matrix ρ AB of size |A| = |B| = 1 constructed from the ground state is equivalent to the Werner state (46) [47]. We set the total number of spins to N = 12. Interestingly, the resulting EoP exhibits a slight peak at the farthest distance d = 5 while MI monotonically decreases (Fig. 18). This peak also shows that the EoP does not necessarily monotonically decrease along with other correlation measures.

A Chaotic Spin Chain
Finally we consider a non-integrable model by adding a parallel magnetic field to the Ising model, We set the parameters h = 1.05 and g = −0.5 following [48]. We use the same setup as the Heisenberg model and find that the long range correlations are almost vanishing in the vacuum (Fig. 18).
Non-monotonicity of EoP for spin chain with N = 4 The non-monotonicity of EoP for N = 4 ( Fig. 8) is common in any homogeneous spin chain. The key observation is that, when N = 4, the symmetric and antisymmetric projectors P sym , P asym acting on A and B located on the diagonal position (d = 1) are symmetries of the system, i.e. they commute with the Hamiltonian [H, P sym ] = [H, P asym ] = 0. Indeed, each term <ij> σ l i ⊗ σ l j , i σ l i (l = x, y, z) commutes with these projectors, and thence they are the symmetries of the system regardless of the values of coupling parameters. Since P sym and P asym are orthogonal, its unique vacuum (and any other non-degenerate excited state) always belongs to either symmetric or anti-symmetric subspace of H AB (for example, in the vacuum of the antiferromagnetic isotropic Heisenberg model, ρ AB is P sym /3 itself).
On the other hand, the EoP coincides with S A when ρ AB has support either on the symmetric or on the antisymmetric subspace of H AB [40]. Thus we have E P = S A at d = 1, while E P < S A at d = 0, 2 in general, leading a non-monotonic behavior of the EoP.