Entanglement spectrum of the Heisenberg XXZ chain near the ferromagnetic point

Entanglement is one of the most fascinating concepts in quantum mechanics and the application of entanglement related ideas to condensed matter in recent years has boosted a huge interdisciplinary effort [1]. In particular, there has been a growing interest in studying the properties of the so-called entanglement spectrum (ES) [2]–[10]. The original motivation for the study of the ES was that the ES lies at the heart of the density matrix renormalization group (DMRG) algorithm [2]. More recently, the ES has attracted a lot of attention because of its relation to low-energy boundary modes in topological phases [3].

Considering a bipartition of a system into two parts A and B, the entanglement spectrum (ES), {ξ_i}, is defined in terms of the Schmidt decomposition

$$\begin{eqnarray} &&\vert \psi \rangle =\mathop{\sum }\limits _{i}{\mathrm{e}}^{-{\xi }_{i}/2}\vert \psi _{i}^{A}\rangle \otimes \vert \psi _{i}^{B}\rangle \end{eqnarray} \tag{ 1 }$$

where |ψ〉 is the state of the system and $\vert {\psi }_{i}^{A}\rangle$ ($\vert {\psi }_{i}^{B}\rangle$) form an orthonormal basis for the subsystem A (B). The ES {ξ_i =− logλ_i} can also be thought of in terms of the eigenvalues {λ_i} of the reduced density matrix ρ_A obtained after tracing out the B part of the system density matrix |ψ〉〈ψ|.

In this work we study the ES of the finite spin-$\frac{1}{2}$ anisotropic Heisenberg (XXZ) chain in the vicinity of the ferromagnetic point Δ =− 1. (Here Δ is the anisotropy parameter.) In particular, we focus on the transformation from the gapless phase at  − 1 < Δ ≤ 1, described by a c = 1 conformal field theory (CFT), to the Δ =− 1 point. In the gapless phase (except for the free fermion point Δ = 0 where many results are available [2]) only the distribution of the ES levels is known from conformal invariance [10]. On the other hand, at Δ =− 1 the system is not conformally invariant and a dramatic change occurs. This is reflected in a peculiar behavior of usual entanglement related quantities such as the Renyi entropies and the concurrence [11]–[13]. In particular, the Renyi entropies in the vicinity of Δ =− 1 show the signature of an essential singularity [11, 12]. The purpose of this work is to understand how this change in the behavior of the system is reflected at the level of the ES.

The Δ =− 1 point is unitarily related to the SU(2)-invariant Heisenberg ferromagnetic chain, and hence the ground state is the degenerate SU(2) multiplet corresponding to the highest total spin S_T. At Δ =− 1 the states of the multiplet are flat superpositions of all the possible spin configurations compatible with total magnetization, ${S}_{\mathrm{T}}^{z}$, i.e., the particles are fully delocalized. We therefore refer to the ground state multiplet as the symmetric multiplet. In other terms, at Δ =− 1 the ground state shows a 'mean field' structure, i.e., the system has no notion of distance. This implies, for example, that there is no difference between open and periodic boundary conditions.

While there are well-known many-body model systems lacking notions of distance, such as the Lipkin–Meshkov–Glick model [14]–[17] and the Richardson model [18], these are generally mean field like at the level of the Hamiltonian. The specialty of the present case is that the XXZ interaction is local, and the symmetric structure appears only at a single value (Δ =− 1) of the anisotropy parameter. Varying Δ is a natural way to tune away from the symmetric structure. In this paper, we show how this change from the symmetric to the more usual situation is manifested in the ES.

Although the usual interest is in the entanglement properties of the ground state, in the vicinity of Δ =− 1 it is natural to study the entire set of states which become degenerate at that point. The ES properties we report are found to be common to the entire symmetric sector. For simplicity of demonstration or explicit calculations, we find it convenient to sometimes focus on simpler large-${S}_{\mathrm{T}}^{z}$ members of this manifold (which in the energy spectrum are near the top of this group of states for Δ >− 1) in addition to the ${S}_{\mathrm{T}}^{z}=0$ ground state.

At Δ =− 1, the symmetric structure fully determines the ES [19]–[22]. However, at Δ >− 1 the ES is not restricted, and in the limit Δ → −1⁺ most of the ES levels diverge (i.e. most of the eigenvalues of ρ_A vanish). The divergent levels correspond to components in the ground state wavefunction at Δ >− 1 that do not reflect the symmetric structure present at Δ =− 1. Moreover, while at Δ =− 1 the ES does not contain any information about the boundary conditions or the geometry of the partitioning between A and B, for Δ >− 1 the ES for periodic and open boundary conditions display striking differences. For open boundary conditions the ES is less dense (some of the ES levels are missing). This is related to there being two boundaries between A and B for periodic boundary conditions (and only one for open). This difference is also present in the behavior of the entanglement entropy in the gapless phase (as 1/3logℓ for periodic and as 1/6logℓ for open boundary conditions with ℓ the length of subsystem A). We also find that the multiplicity count of the ES in the vicinity of Δ =− 1 (number of ES levels diverging at different rates) seems to be described by an elegant combinatorial formula which we present as a conjecture.

We also analyze the entanglement eigenfunctions, i.e. the eigenstates of ρ_A. While the eigenfunctions corresponding to the lower (non-diverging) ES levels mimic the behavior of the ground state wavefunction (i.e. they exhibit a symmetric structure in the limit Δ → −1⁺), the others show richer structures and in some cases their leading behavior (in  ≡ Δ + 1) can be given in terms of known polynomial functions.

Finally we investigate the distribution of the ES levels in the vicinity of the ferromagnetic point. Although in the thermodynamic limit this is expected to be a universal function, for finite chains the ferromagnetic point introduces significant corrections in the whole region Δ < 0.

For our analysis we take advantage of the fact that the XXZ is an integrable model and all its eigenstates can be in principle calculated using the Bethe ansatz. We also use exact numerical diagonalization. Since ES levels are determined often up to values of ∼100, the reduced density matrix eigenvalues λ_i = e^−ξ_i need to be calculated with very high precision. Therefore, as in [4], we use arbitrary precision numerics.

In section 2 we discuss the spin-$\frac{1}{2}$ XXZ chain around Δ =− 1, showing the ground state manifold and how the wavefunctions of this symmetric sector behave at Δ → −1⁺. In section 3 we present the main features of the ES, reviewing the exact structure known at Δ =− 1 and then presenting numerical diagonalization data for Δ =− 1 + . In section 4 we treat the ES perturbatively in the distance Δ + 1 from the ferromagnetic point, using the Bethe ansatz. The combinatorial aspects of the ES are also addressed and characterized. Section 5 is devoted to the study of the entanglement eigenfunctions. In section 6 we investigate how the distribution of the ES levels changes as the system is detuned away from the ferromagnetic point. Details of the Bethe ansatz calculation appear in the appendix.

The spin-$\frac{1}{2}$ XXZ open chain of L sites is defined by the Hamiltonian

$$\begin{eqnarray} &&\mathcal{H}=\frac{1}{2}\sum _{i=1}^{L-1}({S}_{i}^{+}{S}_{i+1}^{-}+{S}_{i}^{i}{S}_{i+1}^{+})~+~\Delta \sum _{i=1}^{L-1}{S}_{i}^{z}{S}_{i+1}^{z} \end{eqnarray} \tag{ 2 }$$

and the periodic chain is obtained by adding terms connecting sites L and 1. At Δ =− 1 the model is related to the ferromagnetic Heisenberg model ($-{\sum }_{i}{\vec{S}}_{i}\cdot {\vec{S}}_{i+1}$) by the local unitary transformation S^x,y → −S^x,y applied on every second site of the chain³. The ground state of the Heisenberg ferromagnet is in the SU(2) multiplet with highest total spin S_T (symmetric multiplet), i.e. the one containing the fully polarized state |F〉 ≡ |↑↑↑⋯〉. The L + 1 states of the degenerate multiplet are obtained from |F〉 by successive applications of the lowering operator ${S}_{\mathrm{T}}^{-}\equiv {\sum }_{i}{S}_{i}^{-}$. It is conventional to refer to |F〉 as the 'vacuum' and the overturned spins as 'particles'. The states $({S}_{\mathrm{T}}^{-})^{M}\vert F\rangle$ are equal weight superpositions of all the possible spin configurations with M down spins (particles). The magnetization is ${S}_{\mathrm{T}}^{z}=L/2-M$ in terms of the particle content M; either (${S}_{\mathrm{T}}^{z}$ or M) is a good quantum number for all Δ and labels the states of the multiplet.

In figure 1(a) we show the energies of the states which are connected to the symmetric multiplet at Δ =− 1, as a function of Δ + 1. The energy levels in the vicinity of Δ =− 1 are

$$\begin{eqnarray} &&E=-\frac{L}{4}+\Biggl [\frac{L}{4}-\frac{1}{L-1}\Biggl (\frac{L}{2}-{S}_{\mathrm{T}}^{z}\Biggr )\Biggl (\frac{L}{2}+{S}_{\mathrm{T}}^{z}\Biggr )\Biggr ](\Delta +1). \end{eqnarray} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

The ground state at Δ >− 1 is the one with ${S}_{\mathrm{T}}^{z}=0$, while the other states of the multiplet are excited states. In figure 1(b) we plot the wavefunction of the eigenstate with two particles (which is part of the symmetric multiplet but not the ground state) showing wavefunction components versus the relative position (x₂ − x₁) of the two particles. While at the point Δ =− 1 the wavefunction is completely flat (the two particles are completely delocalized), for Δ >− 1 this property is progressively lost farther away from Δ =− 1.

This behavior is true for any state in the symmetric multiplet, although more complicated to display. For example, in figure 1(c) we show wavefunction components of the ground state (${S}_{\mathrm{T}}^{z}=0$) for several values of Δ. Now there are multiple components for each value of x₂ − x₁. It is clear, however, that in the limit Δ → −1 the wavefunction becomes flatter and eventually becomes an equal weight superposition.

In this section we describe how the structure of the ground state in the vicinity of Δ =− 1 determines the behavior of the ES. To calculate the ES we consider a bipartition of the chain into two parts A and B (of lengths L_A and L_B respectively). The wavefunction |ψ〉 of the chain, given two orthonormal bases for the subsystems $\vert {\varphi }_{i}^{A}\rangle ,\vert {\varphi }_{j}^{B}\rangle$, can be written as

$$\begin{eqnarray} &&\vert \psi \rangle =\mathop{\sum }\limits _{i,j}{\mathbf{M}}_{i j}\vert \varphi _{i}^{A}\rangle \otimes \vert \varphi _{j}^{B}\rangle . \end{eqnarray} \tag{ 4 }$$

The matrix M_ij is a central object in the study of the ES. The reduced density matrix eigenvalues λ_i = e^−ξ_i (equation (1)) are the squares of the Schmidt (singular) values of M. Since the block magnetization ${S}_{A}^{z}$ is a good quantum number for the matrix M we can label the ES levels with ${S}_{A}^{z}$.

3.1. The point Δ =− 1

At Δ =− 1 the flat structure of the states in the symmetric sector makes it possible to construct the matrix M explicitly and to calculate its singular values ${\sqrt{\lambda }}_{i}$. Any eigenstate in the symmetric multiplet is identified given the length of the chain L and the magnetization ${S}_{\mathrm{T}}^{z}$ (particle content). For each of these eigenstates there is one ES level per block magnetization sector ${S}_{A}^{z}$ which is given by [19]

$$\begin{eqnarray} &&{\xi }_{\{ {S}_{\mathrm{T}}^{z},{S}_{A}^{z}\} }^{\{ L,{L}_{A}\} }=-\log \left [\frac{\left({{L}_{A}\atop {L}_{A}/2-{S}_{A}^{z}}\right)\left({L-{L}_{A}\atop (L-{L}_{A})/2-{S}_{\mathrm{T}}^{z}+{S}_{A}^{z}}\right)}{\left({L\atop L/2-{S}_{\mathrm{T}}^{z}}\right)}\right ]. \end{eqnarray} \tag{ 5 }$$

As a consequence of the delocalized structure of the states in the symmetric multiplet, equation (5) holds for both periodic and open boundary conditions.

After restricting to ${S}_{\mathrm{T}}^{z}=0$, that is the ground state at Δ >− 1, and taking the thermodynamic limit L → ∞ one obtains

$$\begin{eqnarray} &&{\xi }_{\{ {S}_{A}^{z}\} }^{\{ {L}_{A}\} }=-\log \left [\frac{1}{{2}^{{L}_{A}}}\left({{L}_{A}\atop {L}_{A}/2-{S}_{A}^{z}}\right)\right ] \end{eqnarray} \tag{ 6 }$$

which has also been obtained in [21].

3.2. The ES in the limit Δ → −1⁺: exact diagonalization results

In figure 2(a) we show the ES in the vicinity of Δ =− 1 (i.e. Δ + 1 ≡  = 0) for a block with L_A = 6 in a chain with periodic boundary conditions and L = 18 sites. For clarity we show only the sector of the ES with ${S}_{A}^{z}=0$, since the results in other ${S}_{A}^{z}$ sectors are qualitatively similar.

Zoom In Zoom Out Reset image size

The ES levels show singular behavior. In the limit Δ → −1⁺ there is only one level converging to a finite value given by equation (5), while all the others diverge. The diverging behavior in the ES signals that in the limit Δ =− 1⁺ the wavefunction has to be a symmetric superposition of all allowed configurations. The components in the total wavefunction which reflect this symmetry give rise to a non-diverging level, while the others contribute to the diverging ones. In section 5 we will show how this distinction is manifested in the structure of entanglement eigenfunctions corresponding to diverging and non-diverging ES levels.

In figure 2(b) we demonstrate that the behavior of the levels in the vicinity of Δ =− 1 is given by  − αlog(Δ + 1) + β, α being an integer. The slopes α, going from bottom to top in the ES, are given by the sequence 0,2,4,6,... where 0 refers of course to the lowest level. A diverging level in the ES as  − αlog(Δ + 1) corresponds to an eigenvalue of the reduced density matrix which vanishes as (Δ+1)^α. Focusing on the number of levels diverging with the same slope, an intriguing multiplicity sequence appears as higher levels in the ES are considered. In figure 2(b) we use the notation α[m_per], m_per being the number of levels with slope α. We observed that m_per(α) is given by the sequence m_per(α) = 1,1,2,3,....

In figure 2(c) we show the same plot but for a chain with open boundary conditions (L = 17). The ES shows the same qualitative behavior, although now the slopes are α = 0,4,8,12,.... The multiplicity structure apparently is described by the sequence m_open(α) = m_per(α/2).

In a generic sector with ${S}_{A}^{z}\not = 0$, the slope sequence (the α values) remains the same (α = 0,2,4,... for periodic boundary conditions and α = 0,4,8,... for open boundary conditions), but the multiplicity structure is different.

The outlined scenario depends only on the symmetric structure of the states in the limit Δ → −1⁺, thus the ES obtained from any state within the symmetric multiplet is expected to show similar features. This observation is supported in figure 3 where we show the full ES obtained from the state in the symmetric multiplet with only two particles (which is not the ground state for Δ >− 1, see figure 1). We show the ES for both periodic (left) and open boundary conditions (right). In both cases there are only three possible sectors ${S}_{A}^{z}$ determined by the number of ways of distributing the two particles in A and B. Nonetheless, the main features observed in figure 2 are present: (a) for each value of ${S}_{A}^{z}$ there is one ES level converging to a finite value in the limit Δ → −1⁺ and this value is given by equation (5); (b) the slope of the only divergent level is 2 for periodic boundary conditions and 4 for open ones, in agreement with the sequences shown in figure 2.

Zoom In Zoom Out Reset image size

The structure of the ES outlined in the previous sections can be quantitatively understood by expanding the state wavefunction for small  ≡ Δ + 1 and then calculating the ES. For the XXZ we can exploit the fact that the model is integrable and use the exact expression of the eigenstate wavefunctions (see appendix for details).

For this description, we choose the orthonormal basis $\{ \vert {\varphi }_{i}^{A}\rangle \} =\{ \vert {x}_{1},{x}_{2},\ldots ,{x}_{{M}_{A}}\rangle \}$ for the A block, where |x₁,x₂,...,x_MA〉 are the block configurations with particles at positions x₁,x₂,...,x_MA. Similarly for the B subsystem we have $\{ \vert {\varphi }_{i}^{B}\rangle \} =\{ \vert {x}_{{M}_{A}+1},{x}_{{M}_{A}+2},\ldots ,{x}_{M}\rangle \}$.

With this notation the elements of the M matrix (equation (4)) are

$$\begin{eqnarray} &&{\mathbf{M}}_{[{x}_{{M}_{A}+1},\ldots ,{x}_{M}]}^{[{x}_{1},{x}_{2},\ldots ,{x}_{{M}_{A}}]}\equiv a({x}_{1},{x}_{2},\ldots ,{x}_{{M}_{A}},{x}_{{M}_{A}+1},\ldots ,{x}_{M}) \end{eqnarray} \tag{ 7 }$$

where [x₁,x₂,...,x_MA] and [x_MA+1,...,x_M] act as row or column indices. Here a(x₁,x₂,...,x_MA,x_MA+1,...,x_M) is the amplitude of the chain configuration with particles at positions x₁,x₂,...,x_MA,x_MA+1,...,x_M. It turns out that the first order in Δ + 1 for the amplitude a in (7) can be obtained analytically for arbitrary M and system sizes L and the resulting ES calculated.

In the three subsections below, we first treat the simplest nontrivial member of the symmetric multiplet (section 4.1), namely the one with only two particles (M = 2), and provide the first order expansion in  ≡ Δ + 1. In section 4.2, we discuss M > 2 and higher orders in , providing a conjecture for the multiplicity sequence of the ES at different orders. In section 4.3 we comment on the open boundary case.

4.1. The sector with two particles (M = 2): first order in Δ + 1

In the sector with two particles in the full chain wavefunction, using the expansion (A.4) and (7) the M matrix can be written explicitly. The possible values of M_A are M_A = 0,1,2. Moreover the blocks of M corresponding to M_A = 0,2 are 1 × 1 (i.e. just numbers). In figure 3 (PBC) these would give the two ES levels with ${S}_{A}^{z}=1,3$, which are not diverging in the limit Δ → −1⁺. For these two levels the first order contribution to the wavefunction gives a first order renormalization of the result (5). After restricting to the case L_A = L_B = L/2 to simplify the expressions that would otherwise be very cumbersome, the ES level in the sector with M_A = 2 is written as

$$\begin{eqnarray} &&-2\log \left [\frac{1}{2}\sqrt{\frac{L-2}{L-1}}\Biggl (1-\frac{L(L-2)}{2 4(L-1)}(\Delta +1)\Biggr )\right ]. \end{eqnarray} \tag{ 8 }$$

Since we chose L_A = L_B the ES level with M_A = 0 is the same.

More interesting is the sector with M_A = 1 (${S}_{A}^{z}=2$ in figure 3), in which case there is a divergent level. The block in matrix M describing this sector reads

$$\begin{eqnarray} &&\fl \mathbf{M}=\sqrt{\frac{2}{L(L-1)}}\Biggl [(1+{A}_{2}\epsilon ){w}^{(0)}\otimes {w}^{(0)}+\frac{\epsilon }{L-1}[{w}^{(0)}\otimes {w}^{(2)}+0\leftrightarrow 2]\nonumber\\ &&-~\frac{2 \epsilon }{L-1}{w}^{(1)}\otimes {w}^{(1)}\Biggr ] \end{eqnarray} \tag{ 9 }$$

where ⊗ is the Kronecker product and we defined the three vectors ${w}_{k}^{(i)}\equiv {k}^{i}$ and  ≡ Δ + 1. The factor A₂ ensures the normalization of the full chain wavefunction (see appendix). M is a symmetric matrix and quadratic polynomial as a function of the coordinates of the two particles. Since the vectors w⁽ⁱ⁾ are linearly independent, the rank of M is three. However, while two of the singular values of (9) are (1) and (), the other one is (²) and it is not meaningful because the wavefunction is accurate only up to first order. The third singular value is thus 'fake' and the effective rank of (9) is two. (See [4] for a detailed discussion of such spurious eigenvalues in perturbative ES calculations.) One ES level of (9) is given by

$$\begin{eqnarray} &&-2\log \left [\frac{1}{2 4\sqrt{2}}\frac{{L}^{2}(L+2)(L-2)}{(L(L-1))^{3/2}}(\Delta +1)\right ] \end{eqnarray} \tag{ 10 }$$

which diverges in the limit Δ → −1 and corresponds to the divergent ES level in figure 3 (PBC). The non-diverging ES level in this sector up to first order is given by

$$\begin{eqnarray} &&-2\log \left [\sqrt{\frac{{L}^{2}}{2 L(L-1)}}+\frac{(L-2)[L(L-1)]^{3/2}}{1 2\sqrt{2}(L-1)^{3}}(\Delta +1)\right ]. \end{eqnarray} \tag{ 11 }$$

Note that the correction to the singular values, at first order in  ≡ (Δ + 1), is always of type  ∼ L for large L (equations (8), (10) and (11)). Thus L = L(Δ + 1) needs to be small for these expansions to be meaningful.

4.2. Many particles (M > 2) and higher orders in Δ + 1: combinatorial structure of the ES

In this section we discuss how the structure of the ES for any state within the symmetric multiplet builds up as the particle number (M) of the wavefunction increases and higher orders in Δ + 1 are taken into account.

In each M_A sector the ES levels are organized in powers of Δ + 1: the lowest ES level converges to a finite value in the limit Δ → −1, while the higher levels diverge as  − αlog(Δ + 1) + β with α = 2,4,6,... as observed in figure 2.

Within a given ${S}_{A}^{z}$ sector and for each value of the slope α there is more than one level, i.e. a nontrivial multiplicity structure arises as higher levels in the ES are considered.

It is natural to investigate how the multiplicity sequence depends on the particle content (M). To this purpose in figure 4 we provide examples of multiplicity sequences in the ES for M_A = 2 and M = 3,4,5 (we consider a chain with L = 18 and L_A = 8). We note that the sequences can be obtained as particular cases from the expansion of the so-called q-deformed binomial which is defined as

$$\begin{eqnarray} &&{C}_{q}(r,s)={\left({r\atop s}\right)}_{q}\equiv \frac{\prod _{k=r}^{r-s+1}(1-{q}^{k})}{\prod _{k=1}^{r}(1-{q}^{k})}. \end{eqnarray} \tag{ 12 }$$

For example, if we fix M = r = 5,M_A = s = 2 we obtain C_q(5,2) = 1 + q + 2q² + 2q³ + 2q⁴ + q⁵ + q⁶. The sequence 1,1,2,2,2,1,1 given by the coefficients of the different monomials matches the one shown in figure 4(c). It can be checked that C_q(4,2),C_q(3,2) give the other sequences in figure 4. We also checked that the result holds for other values of M,M_A. Thus we conjecture that for all the states within the symmetric multiplet the multiplicity of the ES level diverging as  − αlog(Δ + 1) is given by the coefficient of q^α/2 in the expansion of the q-deformed binomial (12) with r = M,s = M_A. This result holds provided that M < L_A.⁴ We observe that for L_A ≤ M ≤ L_B the multiplicity sequence is given by (12) with r = L_A,s = M_A, while if M > L_B it is r = L_A − (M − L_B) s = M_A − (M − L_B).

Zoom In Zoom Out Reset image size

4.3. Open boundary conditions

The XXZ chain with open boundary conditions can also be treated using the Bethe ansatz (details in appendix). The first order expansion of the matrix M is

$$\begin{eqnarray} &&\fl {\mathbf{M}}_{[{x}_{{M}_{A}+1},\ldots ,{x}_{M}]}^{[{x}_{1},{x}_{2},\ldots ,{x}_{{M}_{A}}]}=\mathbf{1}+\epsilon \bigl [{\mathcal{M}}_{1}({x}_{1},{x}_{2},\ldots ,{x}_{{M}_{A}})+{\mathcal{M}}_{2}({x}_{{M}_{A}+1},\ldots ,{x}_{M})\bigr ] \end{eqnarray} \tag{ 13 }$$

where ₁ and ₂ are polynomial functions that can be obtained from (A.8). We stress that (13) does not contain any cross terms of the form x_ix_j with x_i∈A and x_j∈B (note the difference from the result (9) for periodic boundary conditions). This means that at first order in there is no entanglement between A and B.

Despite the different values of α the multiplicity structure in the ES shows similarities with the case of periodic boundary conditions. In figure 5 we show the ES levels in the sector with M_A = 2 and total number of particles M = 3,4,5. As already observed in figure 2 all the ES levels with α = 2k and k odd are missing. Nonetheless the multiplicity sequence m_open(α) is simply related to the one observed for periodic boundary conditions as m_open(α) = m_per(α/2) with m_per given in terms of (12).

Zoom In Zoom Out Reset image size

In this section we examine the entanglement eigenfunctions (eigenvectors of the reduced density matrix). We show contrasting behavior for the eigenfunctions corresponding to the ES levels which diverge in the Δ → −1⁺ limit and those that do not diverge.

We discuss the case of a chain with periodic boundary conditions. We provide information about the structure of the entanglement eigenfunctions in two ways. First, up to first order in Δ + 1 the Bethe ansatz is tractable for calculating the entanglement eigenfunctions for arbitrary M, M_A and L (section 5.1). Second, in the M_A = 1 sector, the numerical data reveal a polynomial sequence for the forms of the wavefunctions, which we express in terms of the discrete Chebyshev polynomials (section 5.2).

5.1. First order in Δ + 1

We start by discussing the sector with M = 2. In this sector the entanglement eigenfunctions reflect the polynomial structure of the matrix M (see (9)).

In figures 6(a) and (b) we show the ES eigenfunctions for M = 2 and M_A = 1 obtained from the expansion (9). The lowest ES level eigenfunction ((a) in figure 6) has a parabolic dependence on x₁ (the position of the particle inside block A) and as the ferromagnetic point is approached it becomes flat. The non-diverging ES levels thus mimic the structure of the full chain eigenstate. The same behavior is present in the sectors with M_A = 0,2, where the ES has only one (non-diverging) level. In contrast, the eigenfunction of the first diverging ES level shows linear behavior as a function of x₁ (figure 6(b)).

Zoom In Zoom Out Reset image size

In figures 6(c)–(f) we show the entanglement eigenfunctions for the case with two particles in the block A (M_A = 2) and only one in B (M = 3). In (c) and (d) we plot the eigenfunction components versus the 'center of mass' coordinate (x₁ + x₂)/2 of the two particles in A. The entanglement eigenfunction of the lowest ES level (c) exhibits the flattening behavior characteristic of the chain wavefunction, while the eigenfunction in (d) is again linear. The same qualitative behavior is observed in the relative distance x₂ − x₁ (see figures 6(e) and (f)).

This result is general for all M_A, M sectors—the entanglement eigenfunctions for the non-diverging ES levels flatten in the limit Δ → −1⁺, while the eigenfunctions for the ES diverging as  − log(Δ + 1) are linear functions of the coordinates x₁,x₂,...,x_MA.

5.2. Entanglement eigenfunctions in the sector with M_A = 1 and arbitrary M

In the M_A = 1 sector the number of ES levels is given by the total number of particles M (equation (12)). These levels are shown in figure 7 (left) for M_A = 1 and M = 5. For each of the ES levels we report (right) the components of the corresponding entanglement eigenfunction versus the position of the particle in block A.

Zoom In Zoom Out Reset image size

We observe that in the limit Δ =− 1⁺ the entanglement eigenfunction of the level appearing at order n (i.e. diverging as  − αlog(Δ + 1) with α = 2n) is a polynomial of order n (as a function of the position of the particle in block A). This is the only information needed in order to know the exact functional form of the entanglement eigenfunctions. The polynomial can be fixed by imposing that the eigenfunctions are mutually orthogonal and are normalized.

For instance, the entanglement eigenfunction of the lowest ES level is the flat superposition |φ₀〉 = ∑_xc|x〉, with |x〉 the configuration with the particle at position x in block A. The coefficient c is determined by imposing the normalization of the eigenfunction. For the first diverging ES level in figure 7 (the red dotted line) one has |φ₁〉 = ∑_x(c⁰ + c¹x)|x〉. The coefficients c⁰,c¹ are fixed by imposing that 〈φ₀|φ₁〉 = 0 and |φ₁〉 is normalized. This gives

$$\begin{eqnarray} &&\vert {\varphi }_{1}\rangle =\sqrt{\frac{1 2}{{L}_{A}({L}_{A}^{2}-1)}}\mathop{\sum }\limits _{x}\Biggl (x-\frac{{L}_{A}-1}{2}\Biggr )\vert x\rangle . \end{eqnarray} \tag{ 14 }$$

The next diverging ES level (the black dashed line) is instead

$$\begin{eqnarray} &&\fl \vert {\varphi }_{2}\rangle =\sqrt{\frac{5}{{L}_{A}({L}_{A}^{4}-5{L}_{A}^{2}+4)}}\mathop{\sum }\limits _{x}(2+L_{A}^{2}+3{L}_{A}+6{x}^{2}-6({L}_{A}-1)x)\vert x\rangle . \end{eqnarray} \tag{ 15 }$$

This procedure allows us to obtain the leading orders of all the entanglement eigenfunctions in the sector M_A = 1. The set of polynomials obtained are known as the discrete Chebyshev polynomials T_n [23] and are defined as

$$\begin{eqnarray} &&{T}_{n}(x)=n !{\mathcal{D}}^{n}\left({x\atop n}\right)\left({x-{L}_{A}\atop n}\right) \end{eqnarray} \tag{ 16 }$$

where is such that f(x) ≡ f(x + 1) − f(x).

In this section we show how the symmetric structure of the Δ =− 1 point affects the distribution of ES levels for Δ >− 1. We use periodic boundary conditions and consider M = L/2 (i.e. the state in the multiplet that is the ground state at Δ >− 1). It has been shown [10] that for 1D gapless systems described by a conformal field theory (CFT) the asymptotic distribution of the ES levels is given by

$$\begin{eqnarray} &&n(\xi )={I}_{0}(2\sqrt{b(\xi -{\xi }_{\min })}) \end{eqnarray} \tag{ 17 }$$

where n(ξ) is the number of ES levels smaller than ξ, ξ_min is the value of the lowest ES level, b = c/6log[Lsin(π/LL_A)], c is the central charge (c = 1 for the XXZ spin chain) and I₀ is the modified Bessel function of order zero. This ES distribution depends only on the central charge (in this sense it is super-universal), i.e. (17) holds irrespective of the value of Δ in the gapless phase.

However, the ferromagnetic point Δ =− 1 (where the model is not conformally invariant) is the endpoint of the c = 1 conformal line which describes the XXZ at  − 1 < Δ ≤ 1. Although in the thermodynamic limit (when all the length scales are sent to infinity) the CFT result is expected to describe the ES distribution in the whole region  − 1 < Δ ≤ 1, for a finite chain deviations from the CFT prediction are expected especially in the vicinity of Δ =− 1. At the Δ =− 1 point itself, the von Neumann and Renyi entropies show scale invariance but not conformal invariance [21]. The asymptotic behaviors are [21]

$$\begin{eqnarray} &&\fl {S}_{A}^{(\gamma )}=\frac{1}{2}\log \Biggl (\frac{\pi {L}_{A}}{2}\Biggr )+\frac{\log (\gamma )}{2(\gamma -1)}+\mathcal{O}(1/{L}_{A}),\tqs {S}_{A}=\frac{1}{2}\log \Biggl (\frac{\pi {L}_{A}}{2}\Biggr )+\frac{1}{2}+\mathcal{O}(1/{L}_{A}).\nonumber\\ \end{eqnarray} \tag{ 18 }$$

Although the scaling is logarithmic, the prefactors differ from the CFT prediction that would give ${S}_{A}^{(\gamma )}=(c/6)(1+1/\gamma )\log {L}_{A}+\mathcal{O}(1)$ and S_A = c/3logL_A + (1).

It is interesting to check how this is reflected at the level of the ES distribution. In figure 8 we show the ES distribution for a chain of L = 18 and block size L_A = 8. Finite size effects are present, but we still notice a good agreement between the numerical data and the CFT prediction (17) for Δ >− 0.15. The numerical data in the range  − 0.15 < Δ < 1 collapse on the same curve, which is well described by the CFT result. However, in the region Δ <− 0.15 no data collapse is observed and the CFT prediction is quite a poor approximation. Thus for these block sizes and system sizes the ferromagnetic point Δ =− 1 appears to dominate the physics of the whole region Δ < 0.

Zoom In Zoom Out Reset image size

In this paper we studied the ES of the spin-$\frac{1}{2}$ XXZ model in the limit Δ → −1⁺ (ferromagnetic point) for a finite chain with periodic and open boundary conditions. To this purpose we used exact diagonalization (with arbitrary precision numerics) and the Bethe ansatz. We characterized the ES for all the chain eigenstates within the highest spin multiplet (symmetric multiplet) in the limit Δ → −1⁺.

At Δ =− 1 these states are flat superpositions of all the spin configurations compatible with the state magnetization ${S}_{\mathrm{T}}^{z}$. In other terms, they show a 'mean field' structure with no notion of distance. As a consequence the ES at Δ =− 1 for each of these eigenstates can be calculated exactly and shows few levels. At Δ >− 1 the 'mean field' structure is progressively lost as we move away from the ferromagnetic point.

We summarize our results as follows.

• (i)  
  In the limit Δ → −1⁺, most of the ES levels diverge. For each value of the block magnetization ${S}_{A}^{z}$ there is only one non-diverging level which reflects components of the state wavefunction that in the limit Δ → −1⁺ are flat superpositions of all the allowed spin configurations. The divergent behavior is given by ξ =− αlog(Δ + 1) + β with α = 0,2,4,... for periodic boundary conditions and α = 0,4,8,... for open ones.
• (ii)  
  The multiplicity of an ES level diverging with slope α is given by the sequences m_per(α) for periodic and m_open(α) = m_per(α/2) for open boundary conditions. We have provided a combinatorial expression (equation (12)) for these sequences.
• (iii)  
  We have studied the entanglement eigenfunctions, i.e., the eigenfunctions of the reduced density matrix ρ_A. The eigenfunctions of the non-diverging ES levels mimic the behavior of the state wavefunction, i.e., in the limit Δ → −1⁺ are flat superpositions of all the possible spin configurations. In contrast, the entanglement eigenfunctions of the diverging levels do not become flat but display rich structures. In particular, in the sector M_A = 1 their leading order in Δ + 1 can be given analytically in terms of the discrete Chebyshev polynomials.
• (iv)  
  Although in the thermodynamic limit in the whole gapless phase at  − 1 < Δ ≤ 1 the ES distribution is described by (17), for a finite chain significant deviations from the CFT prediction and non-universal behavior are observed in the vicinity of Δ =− 1. Within the range of system sizes considered (L_A = 8) the ferromagnetic point Δ =− 1 appears to dominate the whole region Δ < 0.

In this appendix we outline the Bethe ansatz expansions used in the paper.

A.1. Periodic boundary conditions

In the region Δ >− 1 the Bethe equations can be written in the form

$$\begin{eqnarray} &&\arctan \Biggl (\frac{{v}_{l}}{\tan (\zeta /2)}\Biggr )=-\frac{\pi }{L}{I}_{l}+\frac{1}{L}\sum _{m\not = l}^{M}\arctan \Biggl (\frac{{v}_{l}-{v}_{m}}{(1-{v}_{l}{v}_{m})\tan \zeta }\Biggr ) \end{eqnarray} \tag{ A.1 }$$

where M is the number of particles in the chain, v_l (l = 1,2,...,M) are the so-called rapidities and ζ = arccosΔ. I_l are the Bethe numbers and each choice of I_l identifies an eigenstate of the XXZ Hamiltonian. In each ${S}_{\mathrm{T}}^{z}$ sector the lowest energy eigenstate (which is the state in the symmetric multiplet with magnetization ${S}_{\mathrm{T}}^{z}$) is identified by the following choice of the Bethe numbers:

$$\begin{eqnarray} &&{I}_{l}=-(M+1)/2+l\tqs l=1,2,\ldots ,M. \end{eqnarray} \tag{ A.2 }$$

(Note that $M=L/2-{S}_{\mathrm{T}}^{z}$.) From the solutions of the Bethe equations the state can be constructed as

$$\begin{eqnarray} &&\vert \psi \rangle =\mathop{\sum }\limits _{\mathcal{P}}\mathop{\underbrace{(-1)^{{\epsilon }_{\mathcal{P}}}{\mathrm{e}}^{{\sum }_{l}{k}_{{\mathcal{P}}_{l}}{x}_{l}+(i/2){\sum }_{l\lt m}\Phi ({k}_{{\mathcal{P}}_{l}},{k}_{{\mathcal{P}}_{m}})}{}}}\limits _{a({x}_{1},{x}_{2},\ldots ,{x}_{M})}\vert {x}_{1},{x}_{2},\ldots ,{x}_{M}\rangle \end{eqnarray} \tag{ A.3 }$$

where runs over all permutations of the set {x₁,x₂,...,x_M}, x_l being the particle positions. We denote by |x₁,x₂,...,x_M〉 the spin configuration with particles (down spins) at sites x₁,x₂,.... The parity of the permutation is denoted by . We also define k_l ≡ π + 2arctan(v_l/tan(ζ/2)) (mod 2π) and the scattering phase Φ(k_l,k_m) ≡ 2arctan((v_l − v_m)/((1 − v_lv_m)tanζ)).

The idea of the approach is to solve the Bethe equations (A.1) perturbatively order by order in  ≡ Δ + 1 imposing a solution of the form

$$\begin{eqnarray*} &&{v}_{l}={v}_{l}^{(0)}+{v}_{l}^{(1)}\epsilon +\cdots \end{eqnarray*}$$

to obtain a perturbed expression for the state using (A.3). At zeroth order the Bethe equations are given by the set of equations

$$\begin{eqnarray*} &&{v}_{l}^{(0)}=\frac{2}{L}\mathop{\sum }\limits _{m\not = l}\frac{1-v_{l}^{(0)}{v}_{m}^{(0)}}{v_{l}^{(0)}-{v}_{m}^{(0)}}. \end{eqnarray*}$$

A.1.1. First order expansion of the wavefunctions.

Up to first order in the amplitude (see (A.3)) of the wavefunction with two particles can be written analytically as

$$\begin{eqnarray} &&a({x}_{1},{x}_{2})=\sqrt{\frac{2}{L(L-1)}}\Biggl [1+{A}_{2}\epsilon -\frac{\epsilon }{L-1}\Biggl ({x}_{1}-{x}_{2}+\frac{L}{2}\Biggr )^{2}\Biggr ] \end{eqnarray} \tag{ A.4 }$$

where

$$\begin{eqnarray*} &&{A}_{2}=\frac{2}{L(L-1)^{2}}\mathop{\sum }\limits _{{x}_{1}\lt {x}_{2}}({x}_{1}-{x}_{2}+L/2)^{2} \end{eqnarray*}$$

ensures the correct first order normalization. The last term in (A.4) shows a quadratic behavior in the interparticle distance x₂ − x₁ and linear dependence on of the 'curvature' of the wavefunction. Equation (A.4) also shows that in the asymptotic limit  → 0 the wavefunction becomes an equal weight superposition (no dependence on x_i). These features were shown and discussed in figure 1, where we also compared with exact diagonalization data.

In first order it is possible to go beyond M = 2 and derive the amplitude for any state within the symmetric multiplet as

$$\begin{eqnarray} &&a({x}_{1},{x}_{2},\ldots ,{x}_{M})=\frac{1}{\sqrt{\mathcal{I}}}\left [1+{A}_{M}\epsilon -\frac{\epsilon }{L-1}\sum _{i\lt j}^{M}\left ({X}_{i j}+\frac{L}{2}\right )^{2}\right ] \end{eqnarray} \tag{ A.5 }$$

where $\mathcal{I}\equiv 1/M !{\mathop{\prod }\nolimits }_{k=0}^{M-1}(L-k)$, X_ij ≡ x_i − x_j and the normalization is

$$\begin{eqnarray} &&{A}_{M}=\frac{1}{(L-1)\mathcal{I}}\mathop{\sum }\limits _{{x}_{1}\lt {x}_{2}\lt \cdots \lt {x}_{M}}\mathop{\sum }\limits _{i\lt j}\Biggl ({X}_{i j}+\frac{L}{2}\Biggr )^{2}. \end{eqnarray} \tag{ A.6 }$$

While the M = 2 wavefunction corresponds to an excited state at Δ >− 1, the formula above includes as a particular case (M = L/2) the wavefunction of the ground state in the gapless phase. One important property of (A.5) is that it is quadratic and symmetric in all the interparticle distances x_i − x_j. An example of wavefunction amplitudes obtained from (A.5) is shown in figure 1(c).

A.2. Open boundary conditions

It is possible to treat the case of open boundary conditions within the Bethe ansatz formalism [24]. The Bethe equations now read

$$\begin{eqnarray} &&\fl k({v}_{l})=\frac{\pi }{L+1}{I}_{l}-\frac{1}{L+1}\Phi ({v}_{l},-{v}_{l})-\frac{1}{2}\frac{1}{L+1}\mathop{\sum }\limits _{m\not = l}[\Phi ({v}_{l},-{v}_{m})+\Phi ({v}_{l},{v}_{m})] \end{eqnarray} \tag{ A.7 }$$

with Φ(v_l,v_m) ≡ 2arctan((v_l − v_m)/((1 − v_lv_m)tanζ)) and k(v_l) ≡ π + 2arctan(v_ltan(ζ/2)). The Bethe numbers giving the lowest energy state in each magnetization sector (i.e. the states within the symmetric multiplet) are now I_l = {1,2,...,M}. The main difference with respect to the periodic boundary conditions is the presence of the 'reflected' phase shifts Φ(v_l, − v_m).

This fact has dramatic consequences in the wavefunction and therefore in the ES. The general expression of the chain wavefunction has the form

$$\begin{eqnarray} &&\fl \vert \psi \rangle =\mathop{\sum }\limits _{\mathcal{P}}(-1)^{{\epsilon }_{\mathcal{P}}}A({v}_{{\mathcal{P}}_{1}},{v}_{{\mathcal{P}}_{2}},\ldots ,{v}_{{\mathcal{P}}_{M}}){\mathrm{e}}^{\mathrm{i}({k}_{{\mathcal{P}}_{1}}{x}_{1}+{k}_{{\mathcal{P}}_{2}}{x}_{2}+\cdots +{k}_{{\mathcal{P}}_{M}}{x}_{M})}\vert {x}_{1},{x}_{2},\ldots ,{x}_{M}\rangle \end{eqnarray} \tag{ A.8 }$$

where the sum is over all permutations and (arbitrary number of) reflections of the set of rapidities {v₁,v₂,...,v_M}. The amplitude A(v₁,v₂,...,v_M) is not as simple as for periodic boundary conditions and we do not show it. As for periodic boundary conditions we can expand the Bethe equations in the vicinity of Δ =− 1. Quite interestingly now we find that the lowest order solutions are determined by the equations

$$\begin{eqnarray} &&\frac{1}{{v}_{l}^{(0)}}=\frac{1}{L+1}\frac{1+({v}_{l}^{(0)})^{2}}{{v}_{l}^{(0)}}+\frac{1}{L+1}\mathop{\sum }\limits _{m\not = l}\Biggl [\frac{1-v_{l}^{(0)}{v}_{m}^{(0)}}{v_{l}^{(0)}-{v}_{m}^{(0)}}+\frac{1+v_{l}^{(0)}{v}_{m}^{(0)}}{v_{l}^{(0)}+{v}_{m}^{(0)}}\Biggr ]. \end{eqnarray} \tag{ A.9 }$$

By expanding (A.8) in Δ + 1 we find quite remarkably that the sum over the reflections of the rapidities v_l → −v_l cancels at first order in Δ + 1 the products of the form x_ix_j. This fact leads to the suppression of the first order ES levels for open boundary conditions, as we have seen (figures 2 and 5). We provide as an example of this cancellation the first order amplitude a(x₁,x₂) for the (unnormalized) chain wavefunction in the case with two particles, which reads

$$\begin{eqnarray} &&\fl a({x}_{1},{x}_{2})=1-\frac{\epsilon }{L}[{L}^{2}/2+3-L(1/6+{x}_{1}+3{x}_{2})+2{x}_{1}({x}_{1}+1)+2{x}_{2}({x}_{2}+1)]. \end{eqnarray} \tag{ A.10 }$$

Our numerical results indicate that the same type of cancellation will happen at every odd order, but for higher orders this is too cumbersome to perform explicitly.