Entanglement of puri(cid:12)cation in free scalar (cid:12)eld theories

: We compute the entanglement of puri(cid:12)cation (EoP) in a 2d free scalar (cid:12)eld theory with various masses. This quantity measures correlations between two subsystems and is reduced to the entanglement entropy when the total system is pure. We obtain explicit numerical values by assuming minimal gaussian wave functionals for the puri(cid:12)ed states. We (cid:12)nd that when the distance between the subsystems is large, the EoP behaves like the mutual information. However, when the distance is small, the EoP shows a charac-teristic behavior which qualitatively agrees with the conjectured holographic computation and which is di(cid:11)erent from that of the mutual information. We also study behaviors of mutual information in puri(cid:12)ed spaces and violations of monogamy/strong superadditivity


JHEP04(2018)132
sures for mixed states (see reviews [11,12]). However, they always involve minimization procedures, which are more complicated than the one for the EoP (refer to [30] for a Gaussian ansatz for entanglement of formation). For computational difficulity of entanglement measures refer to [31], where it has been established that they are NP-hard justifying that EoP might be a good starting point even if it is not an entanglement measure. In this sense, our analysis of EoP can be regarded as the first step toward computations of entanglement measures of mixed states in field theories.
This paper is organized as follows. In section two, we will briefly review the definition and properties of entanglement of purification (EoP) as well as its holographic counterpart. In section three, we present our general strategy to numerically calculate the EoP in a free scalar field theory. In section four, we will provide explicit numerical results of EoP. In section five, we will study the behaviors of mutual information between various subsystems. In section six, we will examine whether inequalities of monogamy and strong superadditivity are satisfied or not in our examples. In section seven we summarize our conclusions and discuss future problems.

Entanglement of purification and holographic dual
In this section, we briefly review the basics of the entanglement of purification, which include the definition and information-theoretic properties of EoP. We also give a summary of the recently conjectured holographic computation of EoP and its implications.

Definition of entanglement of purification and its properties
Let us consider a mixed state ρ AB in a bipartite system AB. We can always purify this mixed state by extending the Hilbert space from H A ⊗ H B to H A ⊗ H B ⊗ HÃ ⊗ HB such that the total state ρ AÃBB is pure and ρ AB is embedded in it: (2.1) Such a pure state |Ψ AÃBB is called a purification of ρ AB . Note that a purification of a given state ρ AB is not unique and in general there are infinitely many ways to purify it. The entanglement of purification (EoP) of ρ AB is defined by minimizing the entanglement entropy S AÃ (= S BB ) over all possible purifications of ρ AB [13]: Here S AÃ is the von Neumann entropy of ρ AÃ = Tr BB [|Ψ AÃBB Ψ AÃBB |]. Thus the EoP can be understood as a minimal amount of quantum entanglement between AÃ and BB in the extended system.
which does not involve minimization procedures. However, this quantity does not satisfy all properties required for entanglement measures. Also, it does not coincide with the (von Neumann) entanglement entropy when the total system is pure. Moreover, it is natural to expect that this quantity will not have a simple holographic dual in terms of a tractable geometric quantity in generic setups, especially in higher dimensions. This is partly because it coincides not with the von Neumann (n = 1) but the Rényi entropy at n = 1/2 when the system is pure. See also recent discussions in e.g. [29].

JHEP04(2018)132
The general properties of EoP are intensively studied in [20] (See also [13,32]). We briefly review a part of them. First, as we already noted in the introduction, the EoP itself is not just a measure of quantum entanglement between A and B, but is a measure of both classical/quantum correlations between them. In other words, EoP always vanishes for all product states (ρ AB = ρ A ⊗ ρ B ) and is strictly positive for any non-product states. Moreover, E P (ρ AB ) coincides with the entanglement entropy S A (= S B ) when ρ AB is pure (i.e. when there is no classical correlation between A and B). This fact allow us to regard the EoP as a generalization of the entanglement entropy to a measure of correlation for mixed states.
There are several inequalities that the EoP enjoys. For instance, the EoP is always bounded from above by the von Neumann entropies, and from below by a half of the mutual information: Here we simply write E P (A : B) ≡ E P (ρ AB ). Similarly, the EoP satisfies the following inequality for all tripartite states ρ AB 1 B 2 : The mutual informations on the left-hand side are based on the reduced density matrices The EoP on the right-hand side measures the correlation between A and B 1 B 2 .
In particular, if the ρ AB 1 B 2 is pure, the EoP satisfies the polygamy inequality: On the other hand, the reverse of (2.5) is called monogamy and the EoP sometimes satisfies this for mixed states. 2 We will discuss this more in section 6. Furthermore, as expected to be true for any correlation measures, the EoP never increases upon discarding ancilla for any states (sometime called as extensivity): (2.6)

Holographic entanglement of purification
In [14,15] the holographic counterpart of EoP was proposed in the context of the AdS/CFT correspondence in the classical gravity limit. This is the entanglement wedge cross-section denoted by E W (ρ AB ). It represents the minimal cross-section of entanglement wedge [17][18][19] in the bulk AdS spacetime, refer to figure 1. This gives a generalization of the holographic formula of entanglement entropy [8,9]. This E P = E W (or holographic entanglement of purification) conjecture is supported by many facts, including the coincidence of all properties discussed in the previous section, as well as the heuristic derivation based on the tensor network description of the AdS/CFT correspondence. It has also an interesting connection to the bit threads picture [34]. A generalization of this conjecture was also discussed in [16] and the results further support it.
JHEP04(2018)132 * Figure 1. Holographic entanglement of purification. The shaded region is the entanglement wedge of the subsystems A and B in holographic CFTs (we take a constant time slice of global AdS). The dotted lines are the minimal surface whose area gives S AB . The entanglement wedge cross-section E W (A : B) is defined by the minimal area (divided by 4G N ) of codimension-2 surfaces which divide the entanglement wedge into two parts. In this figure this minimal surface is denoted by Σ * AB and A phase transition occurs for the holographic entanglement of purification when we change the distance between A and B in holographic states. For example, in the Poincaré AdS 3 geometry which is dual to a 2d CFT on an infinite space, E W (A : B) can be explicitly written as where c is the central charge of 2d CFT and d is the distance between A and B. We set both the sizes of A and B to be l for simplicity. At the transition d * = ( √ 2 − 1)l the value of the EoP jumps, thereby providing a non-zero gap: ∆E W = c 6 log[3 + 2 We plot a typical behavior near to the transition point in figure 2. Mutual information I(A : B) also exhibits a phase transition [35] at the same point d * as described in the figure 2. However, unlike EoP, the mutual information smoothly goes to zero.
It also allows us to read off the properties of the mutual informations for A, B,Ã,B. Let us consider them assuming a non-trivial situation E W (A : B) > 0. First, we observe that SÃ is the area ofÃ itself 3 divided by 4G N . Then it immediately follows that I(Ã : B) = SÃ + SB − SÃB = 0. On the other hand, I(A :Ã) = S A + SÃ − S AÃ will be UV divergent because S A and SÃ are itself divergent. Note that the entanglement wedge cross section E W (A : B) = S AÃ is always finite (assuming A ∩ B is empty). Thus subtracting this term does not make I(A :Ã) finite. With the simple setup described above, it can be written explicitly by   We regard AÃBB as a new boundary of bulk spacetime defining an extended field theory. The subsystemsÃ andB, lying on the minimal surface used for computing S AB , are identified with the ancilla system. The dashed lines denotes the minimal surfaces whose areas give S A or S B , respectively. Now we have to minimize the S AÃ and that is achieved by minimizing the crosssection of the wedge and that surface is denoted by the thick green line.
where is the UV cutoff. After the phase transition, we get a constant I(A :Ã) = 2S A = The holographic entanglement of purification also satisfies an inequality called the strong superadditivity [14]. This property is not satisfied by the entanglement of purification for generic quantum states. Therefore this property can be regarded as a special property for holographic states. We will discuss this later in section 6.

Computing EoP in free scalar field theory
Here we present a general strategy to calculate the EoP in the ground state of a 1 + 1 dimensional free scalar field theory. We discretize the field theory on a lattice and compute the EoP numerically. Our basic assumption is that since ground state wave functionals of free field theories are Gaussian, the wave functionals which appear after the purifications are also Gaussian. We also choose the minimal size ansatz. Under this assumption, we can calculate the EoP from matrix computations as we will explain below.

Free scalar field theory and discretization
Consider a free massive scalar field theory in 1 + 1 dimension defined by the standard Hamiltonian: We consider its lattice regularization by identifying x = an, where a is the lattice spacing and n = 1, 2, · · ·, N describes the position of each site (see e.g. [42][43][44]). We define the discretized scalar field and its momentum at n-th site: φ n = φ(na) and π n = a · π(na), which satisfy the canonical quantization condition [φ n , π n ] = iδ n,n . We impose the periodic boundary condition φ n+N = φ n and π n+N = π n . Then the rescaled Hamiltonian H = aH 0 reads where the N × N matrix V is given by The ground state wave function Ψ 0 of this lattice scalar theory is computed as where the matrix W is given by √ V or more explicitly: Note that W is a symmetric and real valued matrix. In the present paper, we will set a = 1 by rescaling the definition of the mass parameter m. In our actual numerical computations we will always choose N = 60 and consider five different masses m = 0.0001, 0.001, 0.01, 0.1, 1. A sketch for N = 16 can be found in figure 5.

Calculation of entanglement entropy
We will follow the analysis in [1,[42][43][44] of computation of entanglement entropy in free scalar models. We decompose the Hilbert space H tot as H tot = H A ⊗ H B by choosing the subregion A and its complement B in a lattice system. The numbers of sites in A and B are called |A| and |B|. Consider a gaussian state |Ψ AB in H tot , which is in general written as follow: We define the matrix W and its inverse: where we have the obvious relations Note that for physically acceptable quantum states, the wave function should be normalizable i.e. W should be positive definite.

JHEP04(2018)132
In this setup the entanglement entropy S A = S B = −Tr[ρ A log ρ A ] is computed as follows [1,[42][43][44]. First we compute the eigenvalues {λ i } of the matrix Λ defined by which is positive definite. The entanglement entropy is then computed by the formula (3.11)

Calculations of EoP
Now we are in a position to present how to calculate the EoP: E P (A : B) = E P (ρ AB ) defined by (2.2) for the ground state Ψ 0 in our free scalar lattice model. We divide the total lattice system into subregions A, B and C such that We defined their lattice sizes to be |A|, |B| and |C|. In this setp, we would like to compute the EoP which measures a correlation between A and B. First, we write the ground state wave functionals in the following form: Note that the matrices P, Q, R are all real valued; P and R are symmetric matrices. Then the reduced density matrix ρ AB = TrÃB |Ψ AÃBB Ψ AÃBB | is obtained by integrating out C: Our basic and crucial assumption is that the optimal purified state |Ψ AÃBB in each setup, which minimizes S AÃ , is a gaussian state, described by the gaussian wave functional where J and L are real symmetric matrices and K is a real matrix. For later use, we introduce the matrix S:

JHEP04(2018)132
Since the reduced density matrix ρ AB should agree with (3.13), we find the following two constraints: With these constraints (3.16) imposed, we can calculate the entanglement entropy S AÃ = S BB from the total wave function Ψ AÃBB (3.14) and minimize its value against the parameters in K and L. This is our basic strategy to calculate the EoP.
Here the gaussian ansatz of the purified state (3.14) is just an assumption which we cannot justify with any solid argument. However, it is natural to expect that the class of gaussian wave functionals are closed in themselves and that we may have only to take the minimization of S AÃ within this class. Indeed as we will present below, this ansatz produces many reasonable results, being consistent with the general properties of EoP. Even if our expectation fails, our "minimal gaussian EoP" provides at least a useful upper bound of the actual EoP, which is defined by minimizations over all possible purifications.

Symmetry transformation
In our computation of EoP, we can identify a symmetry transformation of the matrices K and L which do not change the value of S AÃ .
We take P and Q to be two non-degenerate matrices with the sizes |Ã| and |B|, respectively. We also introduce related matricesP (size |A|+|Ã|) andQ (size |B|+|B|) defined bŷ where I |A|,|B| are the identity matrices. The symmetry transformation is given by To see if these transformations indeed do not change the entanglement entropy S AÃ , we can look at the matrix W obtained by rearranging (J, K, L) as follows: where we decompose (J, K, L) based on the indices A,Ã, B andB in an obvious way. The sizes of the matrices A, B and C are (|A| + |Ã|) × (|A| + |Ã|), (|A| + |Ã|) × (|B| + |B|), (|B| + |B|) × (|B| + |B|), respectively. In terms of (A, B, C), the transformations are expressed as -10 -

JHEP04(2018)132
Thus Λ = −E · B T is mapped by the similarity transformation Λ → (P T ) −1 ΛP T and thus S AÃ , computed from the formula (3.10), does not change. By using this symmetry, we can reduce the number of parameters in K and L which we need to minimize to |Ã| 2 + |B| 2 .

Minimal Gaussian ansatz
Even if we assume the Gaussian ansatz, still it looks hopeless to numerically calculate the EoP because the sizes of matrices K and L can be infinite. Therefore we adopt a finite size ansatz, especially the minimal size one given by |Ã| = |A| and |B| = |B|. We call this the minimal Gaussian ansatz. This minimal ansatz is employed to produce our numerical results of EoP, which will be presented in coming sections.
Even though we do not have a full justification of this ansatz, we have numerical supporting evidence that this ansatz can give an exact answer: even if we start with larger sizes of the purification spaces |Ã| > |A| and |B| > |B|, we will get back to the minimal one |Ã| = |A| and |B| = |B| after the minimization, as we will present in the section 4.4.
In this minimal ansatz, we can reduce the matrix K into the following form by taking advantage of the symmetry transformation (3.18): which has 2|A||B| parameters. The matrix L is completely determined from K by the constraint (3.16). Thus the numerical computation of EoP in our setup requires the minimization of S AÃ over the 2|A||B| parameters. In our explicit numerical analysis presented below we will focus on the cases (|A|, |B|) = (1, 1), (1, 2) and (2, 2) with the total number of lattice sites N = 60.

Numerical results of EoP
Now we are prepared to present our numerical results of EoP in our free scalar theory. We choose the total lattice size to be N = 60 and the subsystem sizes to be (|A|, |B|) = (1, 1), (1, 2), (2,2). We perform the numerical computation of EoP E P (A : B) for five different scalar field masses m = 0.0001, 0.001, 0.01, 0.1, 1 (we set a = 1). We are interested in how E P (A : B) depends on the distance d between A and B (refer to figure 5). We employ the minimal Gaussian ansatz (3.21). Thus we have only to minimize S AÃ with respect to the 2|A||B| parameters in K AB and K BÃ as the matrix L is completely determined by K.
In the final subsection, we will present some evidence that supports the minimal ansatz. Let us start with the smallest subsystems |A| = |B| = 1. In this case, we need to minimize with respect to two real parameters K AB = x 1 and K BÃ = x 2 . In our explicit numerical calculations, we always find x 1 = x 2 at any minimum points. This can be understood from the obvious Z 2 symmetry in the original system which replaces A with B and vice-versa. This symmetry leads to the symmetry which exchanges (A,Ã) ↔ (B,B) in the purified system. Our numerical results of EoP and a half of the mutual information are plotted in figure 6. Note that the former should always be larger than the latter as in (2.6) and this is indeed true in our numerical results. As is clear from the graphs, both of EoP and mutual information are monotonically decreasing as the distance N d gets larger. As the mass gets larger, both graphs change from the power law decay to the exponential decay, as naturally expected.
We also plotted the values of x 1 = x 2 which minimize S AÃ in figure 7. It is intriguing to notice that each graph has always a peak at d = 2. We can explain this behavior as follows. When d = 2, A and B are the next to the nearest neighbor and there is a single lattice site, called C, between A and B. It is clear that in the original wave functional, the entanglement between A and C and that between B and C are both equally very strong. Therefore in the purified state |Ψ AÃBB , we can expect that bothÃ andB are closely related to the site C. This means that the correlation between A andB and the one between B andÃ get enhanced for d = 2. On the other hand, when d = 1 and d ≥ 3, a similar consideration does not lead to any clear enhancement. Indeed the parameters K AB = x 1 and K BÃ = x 2 are obviously responsible for these correlations. This explains the peak at d = 2. Next we proceed to the case |A| = 1, |B| = 2. The two sites in B are separately called B 1 and B 2 . In this case we minimize S AÃ with respect to the four parameters (x, y, z, w) in    the matrix K of the form (the indices of K are arranged in the order The results of EoP and a half of mutual information are plotted in figure 8. Qualitative behaviors are very similar to the previous ones for |A| = |B| = 1. As follows from the extensivity of EoP (2.6), the result for |A| = 1, |B| = 2 is larger than that for |A| = |B| = 1 with the same mass m and d.
We also plotted the values of (x, y, z, w) where S AÃ gets minimized in figure 9. As the graphs show, the behaviors of x and z are similar to that of x in the previous case  |A| = |B| = 1. Since x (and z) here are related to the correlation between A andB 2 (and B 2 andÃ), this is enhanced because there is only a single site between A and B 2 (refer to figure 5) in the same way as before. On the other hand, the values y and w are related to the correlations which are rather suppressed by this effect.

Example 3: |A| = 2, |B| = 2
Next we proceed to the case |A| = |B| = 2. The two sites in A and B are separately called A 1 , A 2 and B 1 , B 2 (refer to figure 5 again). Note that, constrained by the resources at our disposal, this is the largest size of subsystems we consider in this paper for our convenience. We expect this example can have some features of field theory limits more the than other examples already discussed. In this case we minimize S AÃ with respect to the matrix K of the form (the indices of K are arranged in the order A 1 A 2 B 1 B 2 ×Ã 1Ã2B1B2 ) given by: (4.2) As we can also confirm numerically, the symmetry which exchanges A and B allows us to set K AB = K BÃ , or equally x = x , y = y , z = z and w = w . Thus to calculate the EoP, we need to minimize S AÃ against four parameters (x, y, z, w). After this minimization, we obtain the results of EoP in figure 10. By comparing them with previous ones, we can confirm the extensivity of EoP (2.6). Also both the EoP and mutual information are again monotonically decreasing. As the mass increases, the power law decay gets changed into an exponential decay. However, we now notice an important difference between the EoP and the mutual information: the values of EoP at d = 1 and d = 2 are almost the same, while those of the mutual information are different. This  plateaux in the EOP qualitatively looks similar to what we observe in the holographic EoP, where there occurs a phase transition (refer to figure 2). Indeed the phase transition point is d * = ( √ 2 − 1)l, where l is the size of the subsystem A and B. In our example, we took l = 2 and thus d c ∼ 1, which is consistent with the above behavior.
It is also intriguing to examine the behavior of parameters (x, y, z, w) at the minimum points. They are plotted in figure 11. First of all, we note that w has a clear peak at d = 2 as in figure 7. The reason for this peak is the same as that of x in |A| = |B| = 1 case: A 2 get strongly correlated withB 2 through the vacant site. This effect highly reduces the correlation between A 2 andB 1 and thus the absolute values of z behave in an opposite way. The behavior of x and y are roughly in the middle between these two. We also find from our numerical data that as d gets larger, the four parameters get closer to each other x y z w. This can be easily understood because when A and B are most separated, all four possible correlations A 1 − B 1 , A 1 − B 2 , A 2 − B 1 and A 2 − B 2 should be strong in the same magnitude.

Numerical evidences for minimal ansatz
In this section, we would like to present numerical evidence that the minimal ansatz (3.21) discussed in sections 3.5, which was employed for all our numerical computations, is sufficient to produce the correct EoP. We discuss this in details for |A| = |B| = 1 case first. For the other two cases, the arguments will follow analogously.
For this case |Ã| = |B| = 1 is the minimal ansatz. Now we try to increase the dimensions of this auxiliary Hilbert space. We consider |Ã| = |B| = 2. The matrix W in (3.19) takes the following form ( elements are in the order AÃBB), We can easily see that the minimal ansatz is contained in this. We set some of the entries to zero such that the W matrix takes the following form, From this it is evident that,Ã 2 andB 2 do not remain entangled with the composite A, B,Ã 1 ,B 1 system. One can then recover the previous results for EoP by setting K A,Ã 1 = K B,B 1 = 1 and K A,B 1 = K B,Ã 1 = x, LÃ 2 ,B 2 = 0 and minimizing over the parameter x regardless of the values of LÃ 2 ,Ã 2 , LÃ 2B2 . ( LÃ 2B2 = 0 makesÃ 2 independent ofB 2 hence the S AÃ 1Ã2 naturally coincides with S AÃ 1 .) Now we want to check that even if we start from (4.3), minimization of S AÃ 1Ã2 will demand thatÃ 2 ,B 2 should decouple from the A, B,Ã 1 ,B 1 . To check this numerically we adopt the following strategy. For our case the K matrix is the following,

JHEP04(2018)132
Then we set, where superscript 0 denotes the minimal ansatz value. We varies all these terms with superscript 1 around zero in some small steps and compute the corresponding values of S AÃ 1Ã2 . Also first of these two constraints in (3.16) fixes all J's. The second one fixes some of the components of L's. For |A| = |B| = 1 the matrix QR −1 Q T is 2 × 2 symmetric matrix. For our case the L matrix is, This is s a symmetric matrix and hence we will have 10 parameters. Using the constraints, KL −1 K T = QR −1 Q T we can determine 3 of them. For our case we determine LÃ 1Ã1 , LÃ 1B1 and LB 1B1 . So we have total of 15 parameters (8 K's and 7 L's) and we vary them around their minimal ansatz values in some smaller steps. 5 From this we find that the value of S AÃ 1Ã2 is always greater than the minimal ansatz value obtained in the previous section for all non trivial values of these extra parameters. So this shows that our minimal ansatz is good enough to produce the correct EoP. We gave a sample plot in figure 12 demonstrating this result. We choose d = 1, N = 60, m = 0.0001. Then we set LÃ 2Ã2 = LB 2B2 = 0.01, LÃ 2B2 = 10 −7 + i where i is the parameter. Also we set, We vary i between −0.004 to 0.004 in the steps of 0.00007 and plot the S AÃ 1Ã2 w.r.t. i. From this plot it is evident that S AÃ 1Ã2 is greater than the corresponding minimal ansatz value which is 2.85393.
Similarly we check this numerically for |A| = 1, |B| = 2 case also and confirm that the minimal ansatz is sufficient to reproduce the correct EoP. 5 We here note that inside this range sometimes it may happen that for some combinations of values of some of the parameters of the matrix W , where the elements are in order A, B,Ã1Ã2,B1B2 doesn't remain positive definite. We exclude such combinations.  Figure 12. The variation of S AÃ1Ã2 w.r.t. the parameter i corresponding to the bigger ansatz and its shows that it is always greater than the corresponding value coming from the minimal ansatz.

Mutual information
In this section we compute various types of mutual information. As discussed in the section 2.2, these different types of the mutual information have also interesting behaviors in holographic computations. The analysis in this section will also serve as a good consistency check for our previous results of EoP. Notice the following useful relations: These suggest that if we want to minimize S AÃ we need to make both I(Ã :B) and I(A :B) small. Our holographic analysis in section 2.2 for the current setup, actually predicts I(Ã :B) hol = I(A :B) hol = 0. Thus in the holographic EoP, the minimization procedure is realized maximally. For non-holographic quantum states, we do not expect such an extreme situation as is so in our results shown below. We also want to mention that we confirmed the inequalities (5.1) and (5.2) against our numerical results.

Analysis of I(Ã :B)
We compute the mutual information between two subsystemsÃ andB in the auxiliary Hilbert space. We plot I(Ã :B) for m = 0.0001, 0.001, 0.01 and 0.1 against the distance d between A and B for |A| = |B| = 1 , |A| = 1, |B| = 2 and |A| = |B| = 2 in the figure 13 using the results obtained in the previous sections. From figure 13, it is evident that irrespective of the size of A and B there is a slight increase in I(Ã :B) around d = 2 and then it decreases monotonically. This is consistent with our previous results for the values of parameters which minimize S AÃ (see e.g. figure 7). For d = 2, there is a lattice point between A and B. As we have argued in the section 4.1, JHEP04(2018)132   this enhances the correlation between A andB and the one between B andÃ gets enhanced and this fact gets reflected in the increase of I(Ã :B) around d = 2. This further shows consistency of our results for EoP. Non-vanishing values of I(Ã :B) deviates from the holographic prediction and the argument of (5.1) implies that the minimization in our free scalar field is not as optimal as the one in holographic CFTs.

Analysis of I(A :B)
Next, we compute the mutual information between A andB and plot them against the distance d between A and B. From the figure 14, it is evident apart from a small hump around d = 2 I(A :B) gradually decreases. We omit the details because its behavior and interpretation are very similar to the previous one I(Ã :B).

Analysis of I(A :Ã)
Lastly, we compute the mutual information between A andÃ. We demonstrate this in the plots below. Again we plot I(A :Ã) against the distance between A and B for different mass. From the figure 15, it is evident that, I(A :Ã) rather increases for all the cases unlike the previous ones. This qualitative agrees with the holographic results (see the right graph of figure 4) if we remember that in our setup of the free scalar model, the sizes of subsystems compete with the lattice spacing and also that we do not expect any phase transitions.

Violations of monogamy and strong superadditivity
In this section we study whether the inequalities of monogamy and strong superadditivity are satisfied or not in our numerical computations. For a (bipartite) correlation measure      E # (including EoP or mutual information), the monogamy inequality [45] for a tripartite state ρ AB 1 B 2 is defined by The strong superadditivity for a 4-partite state ρ A 1 A 2 B 1 B 2 is defined by Note that if a measure E # satisfies the monogamy for all tripartite states, then it also satisfies the strong superadditivity. These inequalities are regarded as desirable features of measures of quantum entanglement for mixed states [33,46]. It is known that the mutual information always satisfies the monogamy 6 (and thus the strong superadditivity) for holographic states [48]. Note that the monogamy of mutual information is not satisfied for generic quantum states [48,49], for example, in free scalar and fermion field theories.
The holographic EoP always satisfies the strong superadditivity [14]. On the other hand, it is known that the EoP does not always satisfy the strong superadditivity for generic states [20,32]. Therefore both monogamy of mutual information and strong superadditivity of EoP will be useful to characterize holographic states of classical gravity duals.

Monogamy/polygamy
We define the ratio of the r.h.s./l.h.s. of (6.1) as .
Note that the extensivity property (2.6) of both EoP and mutual information tells us the ratios R mon is always bounded from below: R mon ≥ 1 2 . As the state gets more quantum correlations, we expect this ratio increases. The lower bound R mon = 1 2 occurs when the state is classical.
We compute this ratio R mon for both the EoP and mutual information (with A = A 1 or A = A 1 A 2 , while each subsystem denotes a single site as in the previous section). We plot this in figure 16 against the distance d between A and B.
It shows that the monogamy of EoP and that of mutual information are always violated except the very massive case m = 1. It can also be seen that heavier mass makes the state more monogamous. Note that, as we discussed in the section 2.2, EoP satisfies the polygamy rather than the monogamy for any tripartite pure states. Hence it is natural to observe the polygamous behavior of EoP for ordinary mixed states.
For the mutual information, the difference between r.h.s. and l.h.s. of (6.1), so called the tripartite information I 3 (= r.h.s. − l.h.s.), was already computed in [49] for a free massive scalar field theory. Their results show that as the mass goes to zero, I 3 gets positively JHEP04(2018)132 Rmon(EoP, A=A1)  Rmon(MI, A=A1)  Figure 16. The ratio R mon of the monogamy inequality of EoP (upper graphs) or mutual information (lower graphs). R mon < 1 means the violation of the monogamy. For all massless cases m = 0.1, 0.01, 0.001, 0.0001 the monogamy is violated. For massive case m = 1 the monogamy can be satisfied. Note that R mon is always bounded from below by 1/2. divergent i.e. the monogamy is maximally violated, which is because the zero mode φ 0 of the massless scalar leads to a classically maximally correlated state ρ AB ∼ dφ 0 |φ 0 φ 0 |. Our numerical results of mutual information for our lattice free scalar model indeed reproduce the same behavior. The EoP shows a similar behavior for large d and this will be explained by the same argument of the scalar field zero mode. The new feature of EoP is that there is a peak at d = 2 in the ratio R mod . This will be again explained by the strong quantum correlation between A 2 and B 2 with the vacant site between A and B, as we already observed the similar peaks in other quantities.

Strong superadditivity
Next we also define the ratio of (6.2) as .
Note that the ratio has the lower bound R SSA ≥ 1 2 . The violation of strong superadditivity is equivalent to R SSA < 1.
We plot this ratio R SSA for the EoP and mutual information in the same manner. The results are essentially the same as that of monogamy. Except for the very massive case, the strong superadditivity is staisfied and there is a peak at d = 2. Naturally this behavior can be interpreted as that of the monogamy violation.   Figure 17. The ratio R SSA of the strong superadditivity of EoP (left) or mutual information (right). R SSA < 1 means the violation of the strong superadditivity. For almost all massless cases m = 0.1, 0.01, 0.001, 0.0001 the strong superadditivity is violated, and for massive case m = 1 this can be satisfied. Note that R SSA is always bounded from below by 1/2.

Conclusions and discussions
In this paper, we calculated the entanglement of purification (EoP) E P (ρ AB ) for the ground state of a 1 + 1 dimensional free scalar field theory. We assumed that the purified state |Ψ AÃBB , which gives the minimum of entanglement entropy S AÃ , is described by a minimal Gaussian wave functional. Thus our numerical results at least give upper bounds of the actual EoP, defined by minimizing against all possible purifications. However, since ground states of free field theories are given by Gaussian wave functionals, there is a chance that our ansatz may provide the correct values of EoP. We presented numerical evidence that justify our minimal ansatz |Ã| = |A| and |B| = |B|. However, we would like to leave for a future problem the final answer to the question whether our ansatz can give exact results for the EoP. In our explicit computations, we focused on the three cases (|A|, |B|) = (1, 1), (1, 2) and (2, 2) with the total lattice size N = 60. The subsystems A and B can be separated by an arbitrary distance d and we studied the behavior of the EoP as a function of d.
Our results show that the EoP is monotonically decreasing when we increase the distance d. As we raise the mass of the scalar field, a power law decay is changed into an exponential decay. These are consistent with the fact that the EoP is a measure of correlation between A and B. We noted that the mutual information I(A : B) is also monotonically decreasing as d increases. However, especially in the case of |A| = |B| = 2, we found an interesting difference between the EoP and mutual information. The EoP has a plateau for 1 ≤ d ≤ 2, while the mutual information does not. We argued that this plateau is qualitatively analogous to the one in the holographic EoP, which is missing for the holographic mutual information. It would be an intriguing future problem to confirm this behavior for a larger subsystem A, which will require more powerful numerical computations with more sophisticated numerical algorithms.
We also studied more details of our computations such as the values of parameters which specify purified states with the minimal S AÃ and the mutual informations for other subsystems I(Ã :B), I(A :B) and I(A :Ã). We noticed that they are correlated in interesting ways. In particular, some of the parameters and I(Ã :B) and I(A :B) take -23 -

JHEP04(2018)132
maximal values at d = 2. We argued that this behavior occurs because when there is an empty site between A and B, the purified siteÃ andB are both strongly correlated with that site. On the other hand, I(A :Ã) is a monotonically increasing function which is similar to holographic calculations.
Moreover we examined whether the inequalities known as monogamy and strong superadditivity are satisfied or not for our numerical EoP results. It is known that in general these inequalities are not satisfied by EoP [20,32]. On the other hand, in the holographic EoP, the latter property turns out to be true [14]. In our analysis of the free scalar field model, we found that either of them are violated for a broad range of masses, including the massless limit. When the mass gets as large as the cut off scale, we found that both monogamy and strong superadditivity are satisfied. This behavior is similar to that for mutual information. Interestingly we observed a clear difference between them: only for EoP, there is a enhancement of monogamy and strong superadditvity at d = 2. We interpreted this as the quantum correlation effect via an empty site which appears in several other quantities studied in this paper.
It would be interesting to perform similar computations based on Gaussian assumptions for other measures, especially those quantifies quantum entanglement, such as entanglement of formation and squashed entanglement. It is also an obviously important future problem to calculate the EoP in conformal field theories directly in the continuum limit as we usually do in the replica method calculation of entanglement entropy.