Eigenstate entanglement scaling for critical interacting spin chains

With increasing subsystem size and energy, bipartite entanglement entropies of energy eigenstates cross over from the groundstate scaling to a volume law. In previous work, we pointed out that, when strong or weak eigenstate thermalization (ETH) applies, the entanglement entropies of all or, respectively, almost all eigenstates follow a single crossover function. The crossover functions are determined by the subsystem entropy of thermal states and assume universal scaling forms in quantum-critical regimes. This was demonstrated by field-theoretical arguments and the analysis of large systems of non-interacting fermions and bosons. Here, we substantiate such scaling properties for integrable and non-integrable interacting spin-1/2 chains at criticality using exact diagonalization. In particular, we analyze XXZ and transverse-field Ising models with and without next-nearest-neighbor interactions. Indeed, the crossover of thermal subsystem entropies can be described by a universal scaling function following from conformal field theory. Furthermore, we analyze the validity of ETH for entanglement in these models. Even for the relatively small system sizes that can be simulated, the distributions of eigenstate entanglement entropies are sharply peaked around the subsystem entropies of the corresponding thermal ensembles.


Introduction
Entanglement is a fundamental feature of quantum matter with far-reaching consequences for its macroscopic properties, complexity, and technological potential. It lies at the heart of modern physics and pervades various research fields [1][2][3]. Here, we focus on bipartite entanglement entropies is the reduced density matrix for an energy eigenstate |E of a bipartite system H A ⊗ H B . In the following, we consider compact subsystems A and investigate their entanglement with the complement B as a function of energy and linear subsystem size . The latter is defined such that vol A =: d , where d is the number of spatial dimensions.
In previous work, we argued that, based on the eigenstate thermalization hypothesis (ETH), the crossover of eigenstate entanglement entropy from the groundstate scaling at low energies and small to a volume law at larger energies or is captured by a single crossover function which assumes a universal scaling form in quantum-critical regimes [40]. A system obeys ETH if all eigenstates (strong ETH) or almost all eigenstates (weak ETH) look locally thermal, i.e., if the expectation values of local observables E|O|E are, in the thermodynamic limit, indistinguishable from expectation values Tr( td O) of corresponding thermodynamical ensembles [41][42][43][44][45][46][47][48][49][50][51]. In particular, td = td (E) could be chosen as the (global) microcanonical ensemble for energy E or, due to the equivalence of ensembles in the thermodynamic limit, as any other thermodynamical ensemble with the same energy density. If this ETH equivalence holds for all local observables in subsystem A, it also holds for subsystem density matrices A obtained from |E E| in the sense that A − Tr B td 1 is small. Thus, ETH also holds for their subsystem entropies. Now, the long-range physics of critical condensed matter systems in equilibrium can usually be described by field theories, and their subsystem entropies are then captured by scaling functions that depend on dimensionless parameters like energy ratios. So, in the presence of ETH, the entanglement entropies of (almost) all energy eigenstates are characterized by such a scaling function.
In Ref. [40], we gave the eigenstate entanglement scaling functions in analytical form for critical onedimensional (1d) systems based on conformal field the-ory and for d-dimensional fermionic systems with a Fermi surface. Numerically, we demonstrated the scaling behavior and weak ETH for the entanglement in systems of non-interacting fermions in d = 1, 2, 3 [40], and for the harmonic lattice model (free scalar field theories) in d = 1, 2 [52]. Simulations for such non-interacting systems are efficient, and the results are basically free of finite-size effects.
In this paper, we probe and confirm the described scaling properties and the applicability of ETH for interacting spin-1/2 chains in integrable and nonintegrable regimes. Although the integrable systems do not obey strong ETH, weak ETH still holds, implying that all eigenstates look locally thermal, except for an exponentially small number of untypical ones [51][52][53]. Weak ETH applies very generally to many-body systems and only requires a sufficiently fast decay of correlations [40,44]. Two classes of interacting chains are investigated here with exact diagonalization: the Heisenberg XXZ model with or without next-nearestneighbor interactions (Sec. 3) and the transverse nextnearest-neighbor Ising model as well as its dual (Sec. 4). For each case, we first assert the scaling properties by showing data collapse for subsystem entropies of thermodynamical ensembles. Then, we assess the applicability of ETH for entanglement entropies by computing entanglement for all eigenstates and comparing to the corresponding scaling function.

Numerical simulations
We use exact diagonalization [54,55] for the models to obtain all 2 L eigenvalues and eigenvectors, where L denotes the total system size. Employing periodic boundary conditions, the systems are translation invariant, reflection symmetric, and invariant under spin inversion. Hence, the Hamiltonians are block diagonal in a basis of semi-momentum states with crystal momentum k, reflection parity p r , and spin-inversion parity p s . The XXZ model also conserves the total z magnetization M , which further reduces computation costs. Using full diagonalization in each block, we compute entanglement entropies S(E, ) of all energy eigenstates [Eq. (1)]. For the test of ETH and the analysis of scaling properties, we compute thermal subsystem entropies S td (β, ) = − Tr td,A ln td,A for the ensemble td = 1 Z e −βH with inverse temperature β and subsystem states td,A = Tr B td .
The largest system size in the simulations for the transverse Ising model is L = 18, corresponding to 40 symmetry sectors. The largest two blocks have dimension 7314 each and quantum numbers |k| = 2π/3, p r = ±1, and p s = +1. The largest system size in the simulations for the XXZ model is L = 20, corresponding to 436 symmetry sectors. The largest 36 blocks have dimension 8398 each and quantum numbers |M | = 1, |k| = π/10, . . . , 9π/10, p r = ±1, and p s = ±1.
To properly assess the scaling functions that describe the crossover from small and E to large or E, one needs a sufficient number of low-energy states. Due to the computational restrictions on L, we hence restrict the analysis to critical points/phases of the spin models. Note, however, that eigenstate entanglement in gapped systems may also be described by scaling functions. As demonstrated in Ref. [40], they feature additional parameters, e.g., the mass-temperature ratio.

Scaling of thermal subsystem entropy
For the long-range physics of critical systems, there is just a single energy scale -the temperature β −1 -and we hence expect a universal scaling function for the subsystem entropy (2) in the form S td (β, ) ∼ Φ( /β) [40]. Independent of integrability, the low-energy excitations of the XXZ model in the spin-liquid phase are lowmomentum spinons with a linear dispersion. Because of the resulting rotation and scale-invariance, the corresponding field theory is a 1+1d CFT [60][61][62]. Within CFT, the thermal subsystem entropy can be computed  with the replica trick and analytic continuation [17,63].
We obtain the scaling function [40] Φ with the central charge c and group velocity v such that Here, C(β) = c 3 ln βv πa + c is a subleading term with ultraviolet cutoff 1/a and a nonuniversal constant c . For small subsystem size or small temperature β −1 , one recovers the celebrated log law S( ) ∼ c 3 ln for groundstate entanglement [6,7,[14][15][16][17][18]. The crossover to a volume law occurs around ∼ c := βv/π. For c , Eq. (4) gives an extensive thermodynamic entropy To confirm this field-theoretic prediction, we numerically compute S td (β, ) in the XXZ systems. The group velocity can be calculated from where E 0 (k) denotes the lowest energy in the symmetry sector with crystal momentum k, k gs is the groundstate momentum, and δk = 2π/L is the momentum spacing. We only consider even L. For L = 4n, the groundstate sector has k gs = 0, and even reflection and spininversion parities. For L = 4n + 2, the groundstate has k gs = π, and odd parities. Figure 1 shows the results for both the integrable and non-integrable cases with various inverse temperatures β. The thermal subsystem entropies, plotted as functions of /β after subtraction of the subleading term C(β), show indeed a collapse onto the crossover scaling function (4a). The small deviations at the ends of the constant-β lines are finite-size effects that become prominent for ≈ L/2. The curves for the highest considered temperatures (β = 2, 3) deviate somewhat towards larger S td . This is due to the fact that CFT, employed to derive Eq. (4a), assumes a perfectly linear dispersion for excitations. A linear dispersion up to arbitrarily high momenta can only occur in field-theoretical descriptions of condensed matter systems and results in ultraviolet divergences that need to be regularized. Lattice systems have a natural ultraviolet cutoff which results in a non-linear dispersion at higher energies. In particular, the excitation energy needs to have an extremum at the Brillouin-zone boundary. As shown in Fig. 1, this non-linearity can be compensated for almost entirely by defining an effective temperature through lim →∞ S td (β, )/ =: cπ/[3vβ eff (β)] and plotting with respect to /β eff (β) instead of /β. As the CFT result is very precise for low temperatures where the nonlinearity is irrelevant, Eq. (5) implies that β eff (β) → β at low temperatures. For higher temperatures, β eff has been defined such that replacing β by β eff in Eq. (5) matches the exact thermodynamic entropy density. In Fig. 1, we employ this procedure only for the highest temperatures because β eff is not needed for the lower temperatures, where its proper determination would also be complicated by the finite-size effects.

Entanglement entropy and ETH
We have found that the thermal subsystem entropies S td (β, ) obey the expected scaling behavior. Let us now probe whether ETH applies in the sense that the entanglement entropies S(E, ) of (almost) all eigenstates converge to S td (β, ) in the thermodynamic limit and are hence described by the same crossover scaling function (4).
The XXZ model (3) conserves the total z magnetization M := i σ z i /2. Due to the spin-inversion symmetry, the considered thermodynamic ensemble td = e −βH /Z in Eq. (2) has a vanishing magnetization, Tr( td M ) = 0. Hence, we compare S td to the entanglement entropies of all zero-magnetization eigenstates. Eigenstates with nonzero magnetization could, e.g., be compared with the canonical ensemble for the same M , or with the grand-canonical ensemble e −β(H−hM ) /Z with the right magnetization density fixed through the field h. In the thermodynamic limit, these ensembles are of course equivalent.
For β −1 → ∞, the energy approaches E = 0 because Tr(σ a i σ a j ) = 0 for any i = j and a = x, y, z. Hence, all S td (β(E), ) curves in Fig. 2 end at E = 0. One would reach positive energies by considering negative temperatures or, equivalently, considering −H instead of H. As this does not yield further insights, we refrain from doing so.
For fixed energy E, the distribution of S(E, ) peaks at its upper edge that coincides well with S td (β(E), ). The coincidence of the peak with the thermal subsystem entropy improves with increasing L and decreasing . The weight in the low-S tail of the distribution decreases with increasing L. Not surprisingly, the low-S tails are more prominent for the integrable systems. This is consistent with earlier results which established that standard deviations of observables decay exponentially in L for non-integrable systems (strong ETH) [41,45,64] and algebraically for integrable systems (weak ETH) [44,47,52,65]. Similarly, large deviation bounds suggest that the ratio of untypical (ather- mal) eigenstates decreases double exponentially in L for non-integrable systems and at least exponentially for integrable systems [51][52][53]. In Refs. [40,52], we simulated integrable non-interacting fermionic and bosonic models for which very large system sizes can be reached and found indeed S(E, ) to be very sharply peaked at S td (β(E), ).

Next-nearest neighbor transverse Ising model and its dual
In this section, we study the 1d next-nearest-neighbor Ising model in a transverse field: It is invariant under a π rotation in the xy plane (σ x → −σ x , σ y → −σ y ) and g → −g. Also, the sign of the nearest-neighbor coupling J is inessential as it can be altered by a π rotation in the yz plane (σ y → −σ y , σ z → −σ z ) on every second site, which does not change the other two terms. Hence, we assume g ≥ 0 and J ≥ 0 without loss of generality.
In the framework of the quantum-classical correspondence [67], the model (8) is closely related to the classical axial next-nearest-neighbor Ising (ANNNI) model in 2d, and hence often referred to as the quantum ANNNI model. For κ < J/2, there exists a second-order Isingtype quantum phase transition between a ferromagnetic phase at small g and a paramagnetic phase at large g as shown in Fig. 3. For κ > J/2 there is the so-called antiphase at small g, followed by the critical floating phase, before transitioning to the paramagnetic phase [66,[68][69][70][71][72][73][74][75][76]. All four phases meet at the multicritical point (κ, g) = (J/2, 0). The model is non-integrable except for the three lines J = 0, g = 0, and κ = 0.
And it is frustration-free on the Peschel-Emery onedimensional line g = J 2 /(4κ) − κ [66,73]. Here, we set J = 1 and choose the two points (κ 1 , g 1 ) = (0.21, 0.62) and (κ 2 , g 2 ) = (−1.695, 3.39) on the transition line between the ferromagnetic and paramagnetic phases. Both are described by an Ising CFT with central charge c = 1/2. Ising systems are also often discussed in terms of dual bond variables [77] Together with τ y i := iτ x i τ z i , they obey the Pauli matrix algebra. In terms of these operators, the model (8) reads This dual Hamiltonian has only nearest-neighbor interactions and comprises the integrable transverse Ising and XY models. In this form, the model has, for example, been discussed in Refs. [66,[78][79][80]. While the structure of the phase diagram remains the same, the physical interpretation of the various phases changes.
In particular, the paramagnetic phase becomes a topological phase with two ground states distinguished by parity. With the Jordan-Wigner transformation from {τ x,y,z i } to fermionic operators [81], H tI becomes the Kitaev-Hubbard chain [80].

Scaling of thermal subsystem entropy
The mapping (9) from H tI to the dual H tI is actually an exact unitary transformation if one applies open boundary conditions in the former and suitable fixed boundary conditions in the latter. An up-spin in the original representation maps to a domain wall. Except for unitary string operators, the mapping is local, and the entanglement properties of both models are hence closely related. For the simulations, we employ periodic boundary conditions. While they do not affect the local physics and the structure of the phase diagram, they affect entanglement and thermal subsystem entropies due to boundary contributions. Figure 4 shows that the differences of subsystem entropies of the model (8) and its dual (10) are in fact almost independent of the subsystem size and have only a small temperature dependence. In the following, we study the point (κ 1 , g 1 ) in the next-nearest-neighbor model H tI and the second point, (κ 2 , g 2 ), in the dual model H tI to make sure that ETH and the general scaling arguments from the introduction and Ref. [40] work in both representations.
The group velocity v for the low-energy excitations is again determined according to Eq. (6) in the groundstate sector which, here, has momentum k gs = 0 and even reflection and spin-inversion parities. Figure 5 shows the thermal subsystem entropies for various inverse temperatures β. When subtracting the subleading constant C(β) and plotting as a function of /β, the data collapses very well onto the crossover scaling function (4). Similar to the observations in the XXZ model, for the point (κ 2 , g 2 ) at the highest temperature (gβ = 2), one sees some deviation from the scaling function due to the nonlinearity of the dispersion at higher energies. This effect can again be compensated for by employing the effective temperature β −1 eff . , and the weight in the low-S tails decreases with increasing L. But the distributions are generally broader than for (κ 2 , g 2 ). This should be due to the fact that the point (κ 1 , g 1 ) is relatively close to the two integrable lines at g = 0 and κ = 0, and to the frustration-free Peschel-Emery line; see Fig. 3. So strong ETH should apply in the thermodynamic limit, but, at small L, deviations from the thermal subsystem entropies can be similarly big as in integrable systems, which only obey weak ETH.

Entanglement entropy and ETH
All four plots display arc-like substructures in the eigenstate entanglement, which are more pronounced for small L. Had we not chosen to only plot zeromagnetization data in Fig. 2, similar structures would occur for the XXZ model with one arc for every magnetization sector. In the thermodynamic limit, the magnetization density of almost all H XXZ eigenstates is in a small interval [− , ] around zero (concentration of measure) and, hence, the arc-like structures for the XXZ model disappear for large L. The transverse Ising model H tI and its dual H tI do not have symmetries except for  Fig. 7 displays arc-like substructures that emerge below the curve for the thermal subsystem entropy. They are more pronounced for smaller system sizes L. In the text we explain these features by approximately conserved quantities. For total system size L = 14 and subsystem size = 5, the plots show eigenstate entanglement entropies S(E, ) in the next-nearest-neighbor transverse Ising model H tI (panel a) and its dual H tI (panel b) at the two critical points indicated in Fig. 3. They are compared to expectation values of the number of σ z and τ z domain walls, quantified by i σ z i σ z i+1 and i τ z i τ z i+1 , respectively.
translation, spatial reflection, and spin inversion, which do not explain the arc-like structures in Fig. 6. However, for (κ, g) = (0.21, 0.62), the transverse Ising model (8) is relatively close to the line g = 0, where the total z magnetization is conserved and all eigenstates can be chosen as tensor products of σ z i eigenstates. To first order in g, one needs to superimpose each such state |ψ with the spin-flipped state i σ x i |ψ to obtain the perturbed eigenstates. While the resulting states have zero magnetization, they have a fixed number of σ z domain walls. The eigenstate expectation values of i σ z i σ z i+1 shown in Fig. 7a establish that the number of σ z domain walls is indeed approximately conserved and corresponds exactly to the entanglement arcs. Similar to the case of the XXZ model, the arc features will vanish with increasing L due to the concentration of measure. The situation of the dual model (10) is analogous. For (κ, g) = (−1.695, 3.39), the term τ z i τ z i+1 dominates and i τ z i τ z i+1 is approximately conserved, explaining the arc-like structures in Fig. 7b.

Conclusion
In summary, we have numerically confirmed eigenstate entanglement scaling for critical interacting 1d systems. For (almost) all energy eigenstates, the crossover from the groundstate scaling at low energies E and small subsystem sizes to a volume law at large E or is captured by a universal scaling function derived from CFT. Both integrable and non-integrable models were investigated using exact diagonalization. Although the accessible system size is limited for these interacting systems, the thermal subsystem entropies match the predicted scaling function very well, and the results establish the applicability of ETH for the entanglement entropies. We also shortly addressed differences in subsystem entropies for the next-nearest-neighbor transverse Ising model and its dual. Due to boundary effects, they basically differ by a constant with small temperature dependence. The results here complement numerical investigations of eigenstate entanglement for large non-interacting fermionic and bosonic systems in Refs. [40,52].