Avalanche of entanglement and correlations at quantum phase transitions

We study the ground-state entanglement in the quantum Ising model with nearest neighbor ferromagnetic coupling J and find a sequential increase of entanglement depth d with growing J. This entanglement avalanche starts with two-point entanglement, as measured by the concurrence, and continues via the three-tangle and four-tangle, until finally, deep in the ferromagnetic phase for J = ∞, arriving at a pure L-partite (GHZ type) entanglement of all L spins. Comparison with the two, three, and four-point correlations reveals a similar sequence and shows strong ties to the above entanglement measures for small J. However, we also find a partial inversion of the hierarchy, where the four-point correlation exceeds the three- and two-point correlations, well before the critical point is reached. Qualitatively similar behavior is also found for the Bose-Hubbard model, suggesting that this is a general feature of a quantum phase transition. This should be taken into account in the approximations starting from a mean-field limit.

One of the main aims of the present work is to study the relation between correlations and entanglement. First, we consider pure states of q spins with two internal degrees of freedom s = 0, 1 (qubits): (S1) Definitions of the entanglement measures C 2 , τ 3 , τ 4 for two, three and four spin-1/2 systems (qubits) are given in Refs. [S1, S2, S3]. It is shown (see Fig. 2 of the main paper) that √ τ 3 is almost always bounded from above by ||ρ corr 3 || 1 . In the case of four spins, there are three in general different entanglement measures τ (i) 4 , i = 1, 2, 3. Looking at the results presented in Fig. S1 for 9 × 10 5 randomly chosen states |ψ 4 , we observe that for the majority of the states ||ρ corr 4 || 1 gives an upper bound to the four-tangles similarly to the case of three spins just discussed above. Although the values of τ 4 can differ from each other for a particular state |ψ 4 [S5], this difference is only small here such that the overall picture is basically the same and we display in Fig. S1 only the results for τ (1) 4 . Since the four-tangles as well as the one-norms of the correlated density operators are continuous quantities, there should be neither horizontal nor vertical gaps between individual points in Fig. S1. However, those are visible which is an indication of insufficient statistics. We could not manage to overcome this problem probably due to the large dimension of the corresponding Hilbert space.

RANK-TWO APPROXIMATION
The quantum transverse Ising chain is a special case of the transverse XY-chain For γ = 0 it has a quantum phase transition of the Ising type. The reduced density matrices of two (ρ 2 =ρ ij ), three (ρ 3 =ρ ijk ), and four (ρ 4 =ρ ijkl ) neighboring spins [S6, S7] ||ρ corr 4 ||1 τ4 FIG. S1: The four-tangle τ4 is shown here against ||ρ corr 4 ||1 for pure states. We have chosen the pure states statistically from the extended Schmidt form [S4]. Each dot in the figure corresponds to a single state. We do not have enough statistics, as can be seen from our viewgraph. Although states for arbitrary value of τ4 should be there whose value of ||ρ corr 4 ||1 comes arbitrarily close to the line τ4 = ||ρ corr 4 ||1, we don't see any occurrences at larger values of τ4. essentially possess two dominant eigenvalues p 1 and p 2 while the sum of the remaining sub-dominant eigenvalues stays below 2.5%. The second eigenvector interferes strongly with the entanglement of the first, whereas we checked that this is not the case for the third; the remaining error is maximally p I>3 ≈ 0.5%. Therefore we only consider the case of rank two density matrices and neglect the rest of few percents of weight. What one is left with are the two highest weights p 1 and p 2 of the density matrix. We briefly discuss two ways of taking care of them: 1) take p 1 or 2) take p 1 /(p 1 + p 2 ) as the highest weight of the new rank-two density matrix. Whereas in 1) one assumes that the neglected part is as destructive to the entanglement as the second state is, in 2) the remaining states do not enter the calculation at all. Both treatments lead neither to an upper bound nor to a lower bound to the entanglement.
Details on how the approximations work for the concurrence can be seen in Figs. S2 and S3. The procedure 1) works perfectly for nearest neighbors, where the plots can hardly be distinguished for the transverse Ising chain and also for the transverse XY-chain for J up to the factorising field, where they start to deviate considerably. This is different when considering larger distances, where both curves have similar shapes only around the critical point (see Fig. S3). There, it gives good estimates even for the zeros of the exact concurrence and largely avoids over-estimating the entanglement in the state.
The genuine multipartite entanglement content, i.e., the  (1) is plotted together with the two approximations 1) and 2) (see text) for the transverse Ising model (γ = 1). It can be seen that approximation 1) (blue dashed curve; almost invisible here) basically coincides with the exact concurrence (black curve). The approximation following scheme 2) (green dash-dotted curve) slightly lies above the exact curve. This is demonstrated in the inset, where the differences C (1) Diff. := C − C (1) (blue dashed curve) and C (2) Diff. := C (2) − C (green dash-dotted curve) for approximation schemes 1 and 2 are shown. Bottom: The same plots as for the transverse Ising chain are shown here for the transverse XY model with anisotropy parameter γ = 0.5. Here the concurrence C2(1) is well described by approximation scheme 1) up to the factorising field. Beyond this point, the approximation is still reasonable, but lies above the exact curve. three-tangle √ τ 3 and the four-tangle, for the nearestneighboring sites are shown in Fig. 1  for nearest neighbors are roughly the same. Therefore, we do not distinguish between the three four-tangles and drop the upper index. This also indicates that the entanglement would be due to the three-or four-particle GHZ-states, respectively. Such a behaviour goes conform with the expectations for that particular model [S8].  (1)||1 as largest correlation function is an upper bound to the corresponding concurrence. For d = 2, the concurrence drops to 0.04 at its maximum at the critical point, whereas ||ρ corr 2 (2)||1 is about roughly the same as ||ρ corr 2 (1)||1.

RESULTS FOR THE ISING CHAIN
Now we turn to the ground state of the one-dimensional transverse Ising chain in the thermodynamic limit under periodic boundary conditions and study mixed states of a given number of spins which can be located at different distances from each other. First, we observe that strict equality between the concurrence and largest correlation found for pure states turns into an upper bound for the concurrence of mixed states, i.e., C 2 ( 1 , 2 ) ≤ ||ρ corr 1 2 || 1 (see Fig. S4), due to the convex roof involved. (d1, d2)||1 is shown for (d1, d2) from nearest neighbors (1, 1) (black, 0) to (1, 3) (green, 2) together with (2, 2) (blue, 3). It is not much affected as for the two-site case: if we look at their maximal values, it only shrinks to about one third of that for ||ρ corr

3
(1, 1)||1. The maximum is assumed at roughly J 0.945 and is slighly moving towards the critical point Jc = 1 increasing their distance.

Correlation functions
In what follows, we deal with the ground states of onedimensional translationally invariant systems. In this case, the density matrices depend only on the distances between the lattices sites. One can always order the site indices such that 1 < 2 < · · · < q and we can writeρ corr 1... q ≡ ρ corr q (d 1 , . . . , d q−1 ), where d i = i+1 − i > 0 are the corresponding distances. Intuitively, one would expect that (i) the correlations of a fixed number of sites decrease with the distances between the sites and (ii) the correlations for fixed distances decrease with the number of sites. For nearest neighbors, the latter would lead to inequalities ||ρ corr 2 (1)|| p ||ρ corr 3 (1, 1)|| p . . .
Our analysis of two completely different models presented below shows that the expectation (i) is always satisfied whereas (ii) does not necessarily hold. The concurrences C 2 ( 1 , 2 ) ≡ C 2 (d), where d = | 2 − 1 | are shown together with the 1-norm of ρ corr 2 in Fig. S4. Whereas the 1-norm of the correlations has no substantial changes, the concurrence at distance d = 2 modifies to about 2% of the maximal value for nearest neighbors (see inset).
This qualitatively doesn't change much for the 1-norms considering 3 sites at distances (d 1 , d 2 ) which means that when the first particle is at site 1 the next is at site 1 + d 1 and the third one at 1 + d 1 + d 2 (hence, the next neighbor reduced density matrix would beρ 3 (1, 1)). This is shown in Fig. S5, where different distances have been considered for ||ρ corr 3 (d 1 , d 2 )|| 1 : (d 1 , d 2 ) = (1, 1) to (1, 3) and (2, 2). ||ρ corr 3 (d 1 , d 2 )|| 1 decays to about one third of the nearest neighbor situation ||ρ corr 3 (1, 1)|| 1 with a maximum at J 0.945 which is moving tinily up to 0.98. One could extract the tendency that the maximum moves for (1, 1) to (1, d 2 ) and ||1 is shown for distances corresponding to (1, d2, 1) for d = 1 (red, 1) to 3 (blue, 3). The effect of growing d seems to be that the curve for d be an upper limit to m with m < n. The major change happens around the critical value of Jc = 1. The inset shows the distances (2, 1, 2) and (2, 2, 2) as compared to (1, 1, 1). There is not much difference noted.

Difference in the norm
We have seen that the 1-norm serves as an upper bound to the correlation functions in the model. We nevertheless studied also the 2-norm, known as the Hibert-Schmidt or Frobenius norm. The result is shown in Fig S7. It is seen that the 2-norm close to the crtitical pointρ corr 4 (1, 1, 1) is still considerably larger thanρ corr 3 (1, 1); it is however still smaller than ρ corr 2 (1) showing only a partial reordering. The corresponding eigenvalues of the matricesρ corr q (1), q = 2, 3, 4 is shown in Fig. S8.

RESULTS FOR THE BOSE-HUBBARD MODEL
The ground state of the Hamiltonian (4) is obtained numerically for arbitrary J by exact diagonalization in the subspace of the Hilbert space where the total momentum is zero. This allows to calculate exactly the reduced density matrices. In the basis of the occupation numbers n 1 . . . n q , the entries n 1 . . . n q |ρ q (d 1 , . . . , d q−1 )|n 1 . . . n q do not vanish, provided that q i=1 n i = q i=1 n i = n B = 0, . . . , N . Thus, the reduced density matrices possess a block-diagonal structure and the blocks are labeled by n B . The correlated density matrices have a similar structure.
The eigenvalues of the correlated reduced density operatorŝ ρ corr q (1, . . . , 1) are shown in Fig. S9. With the increase of the (1, 1)||2, and for J 0.9 we have the opposite. number of lattice sites q, the number of nonvanishing eigenvalues grow but their magnitudes decrease. This leads to a qualitatively different behavior of the one-and two-norms that are plotted in Fig. S10. The one-norms ||ρ corr q (1, . . . , 1)|| 1 grow monotonically with the increase of J and tend to finite constant values in the limit J → ∞. For small values of J, we indeed observe the hierarchy (S3) but already at moderate values of J the one-norms for different q become comparable to each other. It is quite surprising that ||ρ corr 4 (1, 1, 1)|| 1 becomes quickly larger than ||ρ corr 3 (1, 1)|| 1 and later also larger than ||ρ corr 2 (1)|| 1 . This happens much before the critical point J c of the superfluid-Mott-insulator transition. Hence we observe the same behavior as for the integrable Ising model.
The two-norms ||ρ corr q (1, . . . , 1)|| 2 display completely different behavior because the contribution of small eigenvalues is suppressed. The inequalities (S3) are satisfied for the two-norms at any value of J, although the difference between q = 2, 3, 4 is not very large near and above J c . The two-norms possess broad maxima and approach their asymptotic values (1, 1, 1) (c) for the ground state of the transverse Ising model. at J → ∞ from above.
If we consider one-and two-norms for fixed q but vary the distances between the sites, we find that both norms decrease with the distance which is demonstrated in Fig. S11 for two sites (q = 2). The same was also observed for the transverse Ising model.
In the limit of the ideal Bose gas (J → ∞), the entries of the reduced density matrices depend only on the number of sites q but not on the distances between those: n 1 . . . n q |ρ q |n 1 . . . n q = N ! (N − n B )!n B ! n B ! n 1 ! . . . n q ! 1/2 × n B ! n 1 ! . . . n q !