Highly entangled spin chains and 2D quantum gravity

Motzkin and Fredkin spin chains exhibit the extraordinary amount of entanglement scaling as a square-root of the volume, which is beyond logarithmic scaling in the ordinary critical systems. Intensive study of such spin systems is urged to reveal novel features of quantum entanglement. As a study of the systems from a different viewpoint, we introduce large-N matrix models with so-called ABAB interactions, in which correlation functions reproduce the entanglement scaling in tree and planar Feynman diagrams. Including loop diagrams naturally defines an extension of the Motzkin and Fredkin spin chains. Contribution from the whole loop effects at large N gives the growth of the power of 3/2 (with logarithmic correction), further beyond the square-root scaling. The loop contribution provides fluctuating two-dimensional bulk geometry, and the enhancement of the entanglement is understood as an effect of quantum gravity.


Introduction
Entanglement is one of the most characteristic features of quantum mechanics, which provides correlations between objects that are unable to be explained in classical mechanics. In case that a given system S is divided into two subsystems A and B, the reduced density matrix of A is defined by tracing out the degrees of freedom of B in the density matrix of the total system ρ S : ρ A = Tr B ρ S , where Tr B means the trace over the Hilbert space belonging to B. Even if ρ S is a pure state, i.e., can be expressed as the form ρ S = |ψ ψ| for some state |ψ , ρ A is no longer so in general and takes a form like ρ A = c 1 |ψ 1 ψ 1 | + c 2 |ψ 2 ψ 2 | + · · · (c i 's are positive numbers summed to 1) that is called a mixed state. Entanglement is normally measured by the entanglement entropy (EE): which vanishes for the pure states but not for the mixed states. ρ A carries information of interactions between A and B, some of which can be read off through (1.1). We can say that difference of the behavior of (1.1) reflects difference of dynamical property of the system. Let us consider ground states of quantum many-body systems with local interactions. Normally, their EEs are proportional to the area of the boundaries of A and B (called as area law [1]). This can be naturally understood in gapped systems because the correlation length is finite and relevant interactions to the EE are localized along the boundaries. However, gapless systems are exceptional. For example, in (1 + 1)-dimensional conformal field theory, the EE violates the area law by a logarithmic factor, namely grows as the logarithm of the volume of the subsystem [2][3][4]. Recently, Movassagh and Shor discovered a quantum spin chain (called as Motzkin spin chain), whose EE grows as a square-root of the volume and greatly violates the area law in spite of local interactions [5]. A different spin chain with smaller degrees of freedom but exhibiting the same scaling of the EE, called as Fredkin spin chain, was constructed by Salberger and Korepin [6,7].
In this paper, we introduce large-N matrix models whose correlation functions at the tree and planar level reproduce the square-root scaling of the EEs of the Motzkin and Fredkin spin chains. By including loop contribution, such matrix models naturally give an extension of the spin chains. By analyzing the exact solution of one of the matrix models, we find that analogous quantity to the EE including loop effects scales as the power of 3/2 (with logarithmic correction) beyond the square-root. Whereas the tree diagrams are called as rainbow diagrams and look like skeletons [8], loop effects generate diagrams like fishnets that dominate around a critical point and can be regarded as a random surface. This gives intuitive understanding of the enhancement of the correlation and the entanglement between the subsystems. Since the emerging random surface picture defines quantum gravity on two-dimensional bulk, it would be intriguing to discuss the models from the holographic point of view.
This paper is organized as follows. In sections 2 and 3 the Fredkin and Motzkin spin chains and their EEs of ground states are briefly reviewed. In section 4, we introduce large-N matrix models and their connection to the Fredkin and Motzkin spin chains is discussed.
In section 5, from the exact solution of one of the matrix models, we compute analog of the EE that includes effects of fluctuating bulk geometry, and find the enhancement of the square-root scaling to the power of 3/2. Section 6 is devoted to summarize the result and discuss some future directions. The matrix models have so-called ABAB interactions, which are not soluble in the standard manner. In appendix A, we briefly explain the exact solution obtained by character expansion in [9]. Based on the solution, we compute more nontrivial one-point functions from Schwinger-Dyson (SD) equations in appendix B, which are used in section 5.

Fredkin spin chain
We start with a spin chain of length 2n, where up and down spin degrees of freedom with multiplicity (called as color) s are assigned at each of the lattice sites {1, 2, · · · , 2n}. The up-and down-spin states with color k ∈ {1, 2, · · · , s} at the site i is expressed as u k i and d k i , respectively. The Hamiltonian of the Fredkin spin chain [6,7] is given by the sum of projection operators: The Hamiltonian consists of local interactions ranging up to next-to-nearest neighbors. For colorless case (s = 1), we represent the up-and down-spin states as arrows in the (x, y)-plane pointing to (1, 1) (up-step) and (1, −1) (down step), respectively. Then, a spin configuration of the chain corresponds to a length-2n path consisting of the up-and down-steps. The Hamiltonian (2.1) has a unique ground state at zero energy, which is superposition of spin configurations with equal weight. Each spin configuration appearing in the superposition is identified with each path of length-2n Dyck walks that are random walks starting at the origin, ending at (2n, 0), and restricted to the region y ≥ 0.
For s-color case, the above identification is still valid with additional color degrees of freedom. Each spin configuration of the chain corresponds to a length-2n path consisting the up-and down-steps with color. The ground state is unique, and corresponds to length-2n colored Dyck walks, in which the color of each up-step should be matched with that of the subsequent down-step at the same height. The other is the same as the colorless case.
The ground state is given by where P F, 2n, s denotes the formal sum of length-2n colored Dyck walks, w runs over monomials appearing in P F, 2n, s , and N F, 2n, s stands for the number of the length-2n colored Dyck walks: (2.6) The two states of the summand are drawn as colored Dyck walks in Fig. 1.

EE of the ground state
We divide the total system into two subsystems (called as A and B), and compute the EE by tracing out spins in B. Here, let us take a block of the first (n + r) spins as A and the remaining (n − r) spins as B, and consider the case n ± r = O(n) → ∞. Spin configurations in A correspond to a part of colored Dyck paths from the origin to (n + r, h) in the (x, y)-plane, denoted by P F, n+r, s . The height h takes non-negative integers. Similarly, spin configurations in B correspond to the paths from (n + r, h) to (2n, 0), denoted by P (h→0) F, n−r, s . Note that for any colored Dyck path, the part P    Combinatorial arguments give the numbers of the paths P F, n−r, s . The ground state is decomposed as a linear combination of tensor products of two states belonging to A and B (Schmidt decomposition): Here, From the density matrix of the ground state ρ S = |P F, 2n, s P F, 2n, s | with (2.9), the reduced density matrix is obtained as where we used the orthonormal property: (2.14) Since ρ A is a diagonal form, the EE (1.1) is recast as F, n+r,n−r, s does not depend on κ 1 , · · · , κ h and the sums s κ 1 =1 · · · s κ h =1 yield the factor s h .
We plug (2.4) and (2.8) to (2.12) and evaluate its asymptotic behavior as [6] (see also [10,11]) p (h) F, n+r,n−r, s ∼ s −h 1 + (−1) n+r+h 2 8 √ π n (n + r)(n − r) The leading term scales as a square-root of n and significantly violates the area law in spite of local interactions. Note that this originates from the following part of (2.15): F, n+r,n−r, s is crucial to get the √ n-scaling.

Motzkin spin chain
The Motzkin spin chain [5] has additional spin degrees of freedom (we call zero-spin) at each site compared with the Fredkin spin chain. In total, there are up-and down-spin states with color k ∈ {1, 2, · · · , s} and the zero-spin state at the site i, denoted by u k i , d k i and |0 i , respectively. The Hamiltonian of the Motzkin spin chain of length 2n is given in the form of the sum of projection operators: , (3.4) and the interactions are among nearest neighbors. The Hamiltonian has a unique ground state at zero-energy. We can repeat the same identification of the spins and 2D steps as before with the additional zero-spin corresponding to the arrow (1, 0) (flat-step). For colorless case (s = 1), the ground state is expressed by the equal-weight superposition of length-2n Motzkin walks, which are random walks consisting of up-, down-and flat-steps, starting at the origin, ending at (2n, 0) and not allowing paths to enter y < 0 region. For s-color case (s > 1), the color assigned to each up-step should be matched with that of the subsequent down-step at the same height, which is the same as in the Fredkin spin chain.
The ground state is expressed as where N M, 2n, s in the normalization factor is the number of the length-2n colored Motzkin walks given by where 2ρ stands for the number of the flat-steps. For example, 2n = 4 case reads which corresponds to the colored Motzkin walks in Fig. 3.

EE of the ground state
Computing in the same manner as in the previous section, we obtain (3.10) The asymptotic form of p (h) M, n+r,n−r, s is evaluated as [5] (see also [10,11]) Again, the leading term grows as a square-root of n that is beyond the logarithmic violation of the area law usually seen in critical systems, although interactions of the Hamiltonian (3.1) are local. This behavior originates from the following part of (3.8):

Large-N matrix models
In this section, we consider large-N matrix models which reproduce the √ n-scaling of the EEs in (2.17) and (3.12).

Case of Fredkin spin chain
Let us start with one-to-one correspondence between colored Dyck walks and rainbow diagrams in the Gaussian matrix model of N ×N hermitian matrices M f (f = 1, 2, · · · , s): The one-point function is expressed as the sum of rainbow diagrams in the large-N limit. The operator 1 N tr (M 2n ) makes a length-2n loop with a marked point. Feynman graphs for the one-point function are drawn by all possible pairwise contractions of 2n Ms by the propagator In the large-N limit, there remain only planar diagrams among these that are called as rainbow diagrams [8]. For example, Fig. 4 shows rainbow diagrams for 2n = 4 case. We can see that each semi-circle of the rainbow diagrams has one-to-one correspondence to a color-matched up-and down-spin pairs. Thus, The first and second terms corresponds to the first and second terms in Fig. 1, respectively. A semi-circle corresponds to a color-matched upand down-spins.
Next, for the remaining part of (2.18), we consider the operator by introducing another N × N hermitian matrix X. M n+r and M n−r represent spins in the subsystems A and B respectively, and X borders of the subsystems. We can see that the connected correlation function as N → ∞, as far as we ignore diagrams including any contraction by M f -propagators In Fig. 5, we show the diagram corresponding to the divided path in This observation naturally leads to a matrix model action with so-called ABAB-type interactions: under which the one-point function of (4.5) evaluated by tree and planar diagrams gives Here, tree diagrams do not include any loop composed solely by internal lines (propagators) and the four-point vertices tr (M f XM f X) (f = 1, 2, · · · , s). Note that the vertices should be located along the line connecting two Xs in (4.5). Otherwise, the diagram will not be a tree or planar graph. Finally we find F, n−r, s .  Fig. 6. That can give a generalization of the usual EE that includes fluctuating bulk geometry.

Simpler matrix model
We note that the same conclusion as in the above is reached by a simpler matrix model:

Case of Motzkin spin chain
For the Motzkin spin chain, we can use the same matrix models as in the above with simple modification to operators. Let us consider the operator 1 N tr (M + F ) 2n , where F is some N × N matrix. In the binomial expansion of (M + F ) 2n , we can see that M represents the up-and down-spins with color as before, and F corresponds to the zero-spin. Note that terms including the odd numbers of M vanish in its one-point function due to the Z 2 -symmetry under M f → −M f for each f . As far as concerning the EE, namely the number of the Motzkin paths, we may set F to the identity matrix in the computation. Thus, for the Gaussian matrix model

Matrix model solution and extended EE
The matrix models (4.7) and (4.11) have so-called ABAB interactions, to which we cannot apply the standard method to reduce eigenvalue problems [12,13]. However, exact solutions have been obtained by applying a technique to solve O(n) model on a random surface [14] or by using character expansion [9]. 1 In what follows, we focus on the simpler matrix model (4.11) or equivalently the model defined by the action: In what follows, · ′ s denotes an expectation value evaluated by S ′ s . As discussed in appendix A, this model becomes critical at g = g c ≡ − 2 9 . The coupling constant g counts the number of the vertices in Feynman diagrams. As g approaches to g c , diagrams which consist of large number of the vertices dominantly contribute. Since each vertex appearing in the diagram represented by a plaquette in its dual graph, Feynman diagrams tare interpreted as randomly quadrangulated surfaces by the plaquettes. Namely, this defines a lattice model for two-dimensional quantum gravity [15]. In this section, we always consider the large-N limit, which corresponds to extracting surfaces of planar topology, and suppress the symbol "lim N →∞ " for notational simplicity. We introduce a lattice spacing a (length of the edges of the plaquette), and take the continuum limit of the lattice model as a → 0, g → g c with physical quantities (for example, area) fixed finite. As shown in the end of appendix A, the model (5.1) can describe pure quantum gravity in the two-dimensional bulk with the string susceptibility exponent γ str = −1/2 [16][17][18]. However, it exhibits a different scaling relation between boundary length and bulk area (boundary length) ∼ (area) 3/4 (5. 2) from what we have seen in the pure gravity ( (boundary length) ∼ (area) 1/2 ).

Case of Fredkin spin chain
We compute an analog of the leading term of the EE for the Fredkin spin chain defined by S grav.
t and ζ ′ stand for bulk and boundary cosmological constants in the continuum theory that control area and boundary length, respectively. In order to consider contribution from large area and large boundary length, we ignore terms analytic in t or ζ ′ . Here and below, the ellipsis means such irrelevant terms. Let us see large-order behavior of the expansion of we obtain This should be compared with behavior of (2.4) or (4.12): (5.10) The power of n is common. Effects of fluctuating bulk surface are found in t-dependence in the overall factor and 9 2 s = 4s · 9 8 .

g
N tr (M n+r XM n−r X) s Next, we obtain asymptotic behavior of g ∂ 1, 2). After taking t-derivative, we have We will see that these three terms provide different scalings. Use of (5.6), 1 and min{k, ℓ} leads to the large-order series with z i = √ s z ′ i (i = 1, 2). Since the boundary length scales as a −3/2 from (5.5), we define the length in the continuum: Then, the sum of L in (5.14) is evaluated for k = n+r 2 , ℓ = n−r 2 as for u / ∼0. Reading off the coefficient of 1/(z n+r 1 z n−r 2 ) we find

Extended EE with fluctuating bulk geometry
Now we take the ratio of (5.17) and (5.9), and obtain the extended EE (5.3) as From the derivation, this expression is valid unless u is around the origin or ±b. For |u| being of the same order as b (typically |u| = 1 2 b, 1 3 b, etc.) or smaller, the third term becomes dominant scaling as b 3/2 or b 3/2 ln b, which implies that the bulk quantum gravity greatly enhances the square-root scaling (2.17). It would be intuitively understandable because it seems that fishnet diagrams describing a random surface provide stronger correlations than rainbow diagrams at the tree level. Although r-or u-dependence is different, the first term of (5.18) reproduces the square-root scaling. The first and third terms come from the first and third terms in (5.12), respectively. Difference between the two, namely the factor 1/(ζ ′ 1 + ζ ′ 2 ), is crucial to the enhancement. It seems natural because the factor cannot be factorized as a product of a function of ζ ′ 1 and a function of ζ ′ 2 , representing the entanglement of the subsystems A and B. Whereas (2.17) has the maximum at r = 0 (separation at the middle), (5.18) grows as |u| increases when |u| ≪ b. This suggests that fluctuating geometry provides quite different behavior from (2.17) on the dependence of the separation.

Case of Motzkin spin chain
For case of the Motzkin spin chain, in view of (3.13), (4.17) and (4.18), we compute the extended EE around the critical point. First, a generating function of 1 N tr [(M + 1 N ) 2n ] s is related to ω(z ′ ) as In (5.4), we consider large-order behavior in the expansion of and obtain As comparison, (3.10) or (4.17) behaves as We recognize similar effects to the case of the Fredkin model. Next, for g ∂ ∂g 1 N tr (M + 1 N ) n+r X (M + 1 N ) n−r X s , we do similar computation to the Fredkin case. In ω (2) (z ′ 1 , z ′ 2 ) given in (B.11), we use (5.22), 1 and (5.13) to obtain Converting variables as (5.15), we find g ∂ ∂g for u / ∼0. Combining (5.23) and (5.26), we end up with the extended EE including fluctuating bulk geometry as which has essentially the same structure as in the Fredkin case (5.18).

Discussions
In this paper, we introduce large-N matrix models, which reproduce the leading terms of the EE in highly entangled spin chains in their Feynman diagrams at the tree and planar level. By using the exact solution in one of such models, we compute analogous quantity to the EE including the full loop effects in the large-N limit. Although diagrams look like skeletons at the tree level, a two-dimensional random surface emerges in the bulk by including the loop effects. The effects greatly increase the entanglement from the square-root scaling to the scaling of the power 3/2 (with logarithmic correction), and make change the dependence of the separation: the entanglement grows as the difference of the length of the subsystems A and B increases, as far as the difference is small. An intuitive explanation to the former is that fishnet diagrams describing a random surface provide much more correlation than the skeletons. It will be interesting to understand a physical meaning of the latter property.
Since the (s + 1)-matrix model of M f (f = 1, · · · , s) and X can express more details of spin configurations compared to the two-matrix model of M and X, it will be important to analyze the (s + 1)-matrix model and gain deeper insights into the system. For s ≤ 2 the exact solution is found in [14]. In that case, it will be nice if any technique is developed to obtain the relevant one-point functions 1 N tr (M 2n ) and 1 N tr (M n+r XM n−r X) with M = s f =1 M f . It will also be worth doing analogous investigation for Rényi entanglement entropy [10,11] and extending to deformed Motzkin/Fredkin spin chains in which the EEs grow linearly [19][20][21].
In the matrix models, the operators corresponding to spin configurations are regarded as one-dimensional objects, whereas their Feynman diagrams naturally generate twodimensional surfaces. It seems intriguing to investigate the matrix models from the viewpoint of holographic (random) tensor networks [22][23][24].
integers satisfying In terms of {h i }, χ R (A) = det j, k a h j k /△(a) with a i being eigenvalues of A and △(a) the Van der Monde determinant △(a) ≡ det j, k a N −j The expansion coefficient ζ R (Ng) is given by integrals over U(N) matrices: , ǫ j 1 ···j N is an N-th rank totally antisymmetric tensor normalized by ǫ 1···N = 1, and ⌊x⌋ denotes the greatest integer not exceeding x.
Now by using the formula the partition function reduces to eigenvalue integrals as with C N being a constant depending only on N,

A.3 Large-N analysis for highest weights
In order to consider saddle point equations for highest weights {h i }, we define we assume that h(x) = h ′ i becomes a continuous function in large-N limit. A typical distribution of the highest weights is that m k = m k+1 = · · · = m N = 0 for some k and the rest nonzero. Setting α = h ′ 1 and β = h ′ k , we then see that h(x) decreases along the slope of −x for k N < x < 1, and more rapidly decreases for 0 < x < k N . This implies that the highest weight density ρ(h) = − ∂ ∂h x(h) saturates as ρ(h) = 1 for 0 < h < β and 0 < ρ(h) < 1 for β < h < α. (A.18) can be written as Next, let us introduce a function L(h) such that L(h) has the same cut [0, α] as H(h) and Note that λ(h) ≡ exp 1 2 L(h) and h(λ) are functional inverses of each other as mutlivalued functions as discussed in [26]. The inversion of (A.17) is iteratively done as The large-N saddle point equation of (A.7) with respect to h ′ i reads 3 In terms of the analytic function which asymptotically behaves as Its solution is given by where the integration contour C encloses only the square-root cut [β, α] but not the other singularities. The integral is evaluated by inflating the contour as the suitable solution of which reads behaving as X = 1 + 3g 2 + O(g 4 ) for g ∼ 0. The solution determines α and β as The critical point of g is given by a singular point of (A.30) nearest from the origin: g c = − 2 9 . Near the critical point, expansion in ∆ ≡ 2 3 g−gc −gc leads to The critical point h = h * of (A.36) satisfying λ ′ (h * ) = 0 is The functional inversion of λ(h) around the critical point gives Since h † (λ) is given by (A.38) with ± replaced by ∓. From (A. 16), Choices of the branches are fixed by ρ λ (x) ≥ 0. We should take the "−" branch in (A.40) and (A.38). ω(ζ) has a universal meaning in the critical behavior, which defines the quantity in the continuum theory. Here we find the unusual scaling of bulk and boundary cosmological constants implying (boundary length) ∼ (area) 3/4 [9]. The one-point function 1 N tr (M ′ XM ′ X) is computed as The third term is the leading non-analytic term at t = 0, which is relevant to the critical behavior. The fractional power t 3/2 indicates that the string susceptibility exponent is γ str = −1/2 and the bulk surface is described by the same universality class as the c = 0 pure gravity [16][17][18].

B Schwinger-Dyson equations
In this appendix, we derive several SD equations in the matrix model (A.1), by solving which we obtain one-point functions to be needed to compute a generalized EE including bulk gravity effects. Let T p (p = 1, · · · , N 2 ) a basis of N ×N hermitian matrices satisfying tr (T p T q ) = δ pq , Matrices M ′ and X are expanded by the basis as where M ′ p and X p are expansion coefficients. The SD equation in the large-N limit. By combining two SD equations from the identities in the limit N → ∞. Use of (A.42) with z ′ = 3 The identity leads to the SD equation The second term is relevant to the critical behavior and provides the quantity in the continuum theory.