Entanglement of magnon excitations in spin chains

We calculate exactly the entanglement content of magnon excited states in the integrable spin-1/2 XXX and XXZ chains in the scaling limit. In particular, we show that as far as the number of excited magnons with respect to the size of the system is small one can decompose the entanglement content, in the scaling limit, to the sum of the entanglement of particular excited states of free fermionic or bosonic theories. In addition we conjecture that the entanglement content of the generic translational invariant free fermionic and bosonic Hamiltonians can be also classified, in the scaling limit, with respect to the entanglement content of the fermionic and bosonic chains with the number operator as the Hamiltonian in certain circumstances. Our results effectively classify the entanglement content of wide range of integrable spin chains in the scaling limit.


Introduction
Entanglement between subsystems in extended many-body systems has become one of the key ingredients to understanding of quantum matter [1][2][3][4][5]. A quantitative characterization of quantum entanglement is the entanglement entropy and it is defined as follows: For a quantum system in state |K , one can divide the whole system into two subsystems A and B, integrate out the degrees of freedom of the subsystem B of the total system density matrix ρ K = |K K|, and obtain the reduced density matrix (RDM) ρ A,K = tr B ρ K of the subsystem A. The entanglement entropy is the von Neumann entropy of the RDM It can be calculated as the n → 1 limit of the Rényi entropy S (n) The entanglement entropy and Rényi entropy in various extended quantum systems have been investigated for the ground state , the low-lying excited states  and the typical states in the middle of the spectrum [57][58][59][60] using exact and numerical calculations. See also [61][62][63][64][65][66] for the behavior of the average entanglement entropy of the excited states.
When the system is integrable or in general when there is no particle production one can consider quasiparticle excitations whose entanglement shows a remarkable universal property which has a natural qubit interpretation [41,42], see also [32,34,39] for earlier similar conclusions. In [41,42] it was argued that the entanglement entropy for the quasiparticle excited states composed of finite numbers of quasiparticles with finite De Broglie wavelengths or finite intrinsic correlation length is largely independent of the momenta and masses of the excitations, and of the geometry, dimension and connectedness of the entanglement region. Using exact calculations on free fermions and bosons, i.e.
Hamiltonian being equal to the number operator, these results were later confirmed and generalized in [50,52,54,55]. In particular it was shown that for the excited states with more than one mode excited in the most generic circumstances apart from the universal terms there are extra terms. Especially, these additional terms have nontrivial dependence on the momenta of the excited quasiparticles. These terms are non-negligible when the momentum difference of at least a pair of distinct quasiparticles is small. It was argued that these extra terms are also universal though in a weaker sense. In particular, it was shown numerically that similar equations are valid also for the excess entanglement entropy of the excited state in slightly gapped and critical models.
In this paper we consider a circular chain with L sites and choose its subsystem A = [1, ] with consecutive sites. We calculate the quasiparticle excited state entanglement entropy between A and its complement B = [ + 1, L] in the scaling limit L → +∞ and → +∞ with fixed ratio x = L . We first make a conjecture that the results presented in [52,54] are valid not only for non-coupled systems, i.e. the Hamiltonian being just the totally free number operator, but also for generic translational invariant free fermionic and bosonic chains as far as the extra momenta are big. We also made new conjectures for the differences of excited state entanglement entropy. These conjectures are supported with numerical calculations. Then we calculate exactly the entanglement in magnon excited states of interacting chains such as spin-1/2 XXX and XXZ chains. We prove that when finite number of magnons are excited, in the scaling limit, one can decompose the entanglement content into sum of finite number of contributions that are exactly equal to the ones coming from non-coupled fermionic and bosonic chains with particular momenta. In other words, one can classify the entanglement content of also interacting models with respect to the entanglement of the excited states of the free fermionic and bosonic theories.
The paper is organized as follows: In the next section, i.e. section 2, we review the entanglement content of the excited states of the number operator in fermionic chains. Then from numerical calculations we conjecture that these results are universal in the scaling limit for generic quadratic free fermionic systems. In section 3 we make similar arguments for bosonic chains. In section 4 and 5 we first calculate exactly the entanglement content of a few magnon states in the spin-1/2 XXX and XXZ chains and then show that in the scaling limit for the general magnon state the entanglement content can be written as the sum of fermionic and bosonic terms introduced in sections 2 and 3. Finally in section 6 we conclude the paper with discussions.
The paper is accompanied with numerous appendices which in that we explore the exact calcula- A r a † j a j+r + with the local fermionic modes a j , a † j and the parameters A * r = A L−r , B r = −B L−r . It could be diagonalized in terms of the translational invariant quasiparticle modes c k , c † k in the form where ε k ≥ 0 is the energy of the mode k.
The simplest possible case one can consider for the above Hamiltonian is the number operator: The above Hamiltonian, for example, can be considered as the extremely gapped (infinite transverse field) limit of the spin-1/2 XY chain. For the above Hamiltonian we write the quasiparticle modes as Then the Hamiltonian becomes Here p k is the actual momentum and k is the total number of waves, which is an integer and a halfinteger depending on the boundary conditions. Note that p k ∼ = p k + 2π and k ∼ = k + L. In this paper we will just call p k actual momentum and call k momentum. We only consider the case that L is an even integer and the states in the Ramond sector, i.e. periodic boundary conditions for the spinless fermions a L+1 = a 1 , a † L+1 = a † 1 , and so we have the integer momenta The ground state |G is annihilated by all the lowering operators c k and the general excited state |K is generated by applying the raising operators c † k on the ground state The number of excited quasiparticles is |K| = r.
The entanglement entropy of the above states were studied comprehensively in [52,54]  In the scaling limit, i.e. L → +∞ and → +∞ with fixed ratio x = L , the equations (2.9) and (2.10), or equivalently (B.27) and (B.28), can be calculated easily and they show a remarkable universal property [52,54]. The numerical results in the scaling limit of XY and XX chains suggest that the excess Rényi and entanglement entropies take the forms as long as the momenta in the set K satisfy 1 We use 1 ε k to represent the size of the quasiparticle with momentum k. The above condition is just that all the sizes of excited quasiparticles are much smaller than the sizes of the subsystem A and its complement B. Note that with the new condition (2.13), we no longer require that the gap of the model is large or each excited quasiparticle has a large momentum k. Instead, we only require that each excited quasiparticle has a large energy ε k . 1 In [52,54], the proposed condition is [41,42] where ∆ is the gap of the model and the size of excited quasiparticle with momentum k is represented as L |k| . From the example of the XX chain, now we think it is better to represent the size of excited quasiparticle with momentum k as 1 ε k . The condition is modified as Supposing that ε k is a continuous function of k, we have in the scaling limit We obtain min max which is just the new condition (2.13). The physical meaning of the condition used in [41,42,52,54] is that either the correlation length of the model or the maximal size of the excited quasiparticles is much smaller than the sizes of the subsystem and its complement. Using 1 ε k , instead of L |k| , to represent the size of the quasiparticle with momentum k, we see that the correlation length of the model is always larger than or equal to the maximal size of the excited quasiparticles, and we only need to impose the condition for the maximal size of the excited quasiparticles.
The results in the previous paragraph could be further generalized. Suppose we have two states |K = |k 1 , · · · , k r and |K ∪K = |k 1 , · · · , k r , k 1 , · · · , k r in a general translational invariant fermionic chain of the Hamiltonian (2.1). If the momenta in the set K ∪ K satisfies (2.13), from (2.11) and (2.12) we get S (n),fer If the sets K and K have large momentum differences in the scaling limit L → +∞, by which we mean |k − k| → +∞, for all k ∈ K and all k ∈ K, (2.16) the quasiparticles in K and K would decoherent and contribute independently to the entanglement Then we get the differences of the Rényi and entanglement entropies This is just a corollary of (2.11) and (2.12). We emphasize that the above formulas (2.18) and (2.19) are derived under the condition: • The momenta in the set K ∪ K satisfy (2.13); • The sets K and K have large momentum differences in the scaling limit L → +∞, i.e. (2.16).
We note here that the excited state |K of the Hamiltonian (2.1) with specific (A r , B r ) can be actually the ground state |G of another Hamiltonian with a new (A r , B r ), for some examples see [46].
The prime example is the ground state of the XX chain which is an excited state of the number operator N (2.3) as the Hamiltonian. Then one is tempted to guess that in the scaling limit as far as there is a gap between the two states then the equations (2.9) and (2.10) should work. We conjecture that the formulas (2.18) and (2.19) are still valid under a weaker condition: • The momenta in the set K satisfy (2.13); • The sets K and K have large momentum differences in the scaling limit L → +∞, i.e. (2.16).
In particular, the set K may not satisfy (2.13), and we do not have (2.14) and (2.15), but we still have (2.18) and (2.19). Moreover, the set K could possibly has an infinite number of elements in the scaling limit. This conjecture is beyond the results in [52,54]. To support the conjecture, in figure 1, we show examples of the differences of the Rényi and entanglement entropies in the fermionic chain with Hamiltonian    which is in fact the Jordan-Wigner transformation of the spin-1/2 XY chain.
With the above results, we further conjecture that the formula is valid as long as the sets K and K have large momentum differences in the scaling limit L → +∞, i.e. in figure 2. Note that in panels (b) and (c) of figure 2 the large energy condition (2.13) for K is broken and in panels (a) and (b) the large energy condition for K is also broken, while in all the panels the large momentum difference condition (2.16) for K and K is preserved.

Quadratic bosonic chain
This section is parallel to the previous section in many aspects. We consider the Hamiltonian with the bosonic modes a j , a † j and the parameters A * r = A L−r , B r = B L−r . As before we first study the uncoupled Hamiltonian and then make conjectures regarding more general Hamiltonians. In this case a general excited state takes the form with the quasiparticle modes c k , c † k defined the same as (2.4), the ground state |G defined the same as (2.7), and the normalization factor The subsystem mode method can be also used in this case, see appendix C, and this leads to the same equations of the Rényi and entanglement entropies as (2.9) and (2.10) where the exact form of the matrix R A,K can be found in appendix C.1. Since in this case the excited states are not Gaussian states, instead of the correlation matrix method, one has the wave function method [41,42], which leads to a permanent formula for the Rényi entropy in the uncoupled bosonic with the n|K| × n|K| matrix Ω

(n)
A,K is written in appendix C.2. The formula of the Rényi entropy (3.7) is applicable only for integer n ≥ 2. In appendix C.3 we provide some examples regarding the RDM and the entanglement entropy of the states in which a few modes are excited. Some of the results are beyond the cases that were discussed in [52,54].
Numerical calculations done on the discrete Klein-Gordon field theory presented in [52,54] suggest that, similar to the free fermion case, here too the above equations (3.5) and (3.6) give a good approximation also for the excess entropy in slightly gapped and even critical chains as far as the momenta k's are large. This makes it plausible to guess that the conjecture might be also valid for free bosonic chains given the conditions • The momenta in the set K satisfy (2.13); • The sets K and K have large momentum differences in the scaling limit L → +∞, i.e. (2.16).
We make new numerical calculations on discrete Klein-Gordon field theory with mass m, i.e. the Hamil- as shown in figure 3, supporting this conjecture. Note that in the figure there are no entanglement entropies from the wave function method.
The same as that in the fermionic chain, we further conjecture that there is with the condition (2.16). We check the conjecture in the discrete Klein-Gordon field theory (3.10) in figure 4. In all the panels of figure 4 the large energy condition (2.13) for K is broken, while in all the panels the large momentum difference condition (2.16) is preserved.

XXX chain
In this section we calculate exactly the Rényi and entanglement entropies in the magnon excited states in the XXX chain. We only consider the states with finite numbers of magnons and mainly focus on the scattering states. For the bound states see appendix D. In the scaling limit, the magnon excited state Rényi and entanglement entropies turn out to be simply related to the Rényi and entanglement entropies in the fermionic and bosonic chains with the number operator as the Hamiltonian. This shows remarkably that in the scaling limit the equations (2.9) and (2.10) for free fermions and the equations (3.5) and (3.6) for free bosons are still valid even beyond non-interacting free Hamiltonians.

Magnon excited states
We consider the spin-1/2 XXX chain in transverse field with total number of sites L being four times of an integer and periodic boundary conditions for the Pauli matrices σ x,y,z L+1 = σ x,y,z 1 . We only consider the ferromagnetic phase with h > 0, so that the unique ground state is and low energy excitations are magnons.
The XXX chain can be solved with the Bethe ansatz, for review see for example [67,68]. A general with the normalization factor The state |j 1 · · · j R is the configuration that the spins on the sites j 1 , · · · , j R are flipped. The ansatz for the wave function is where the sum P is over all the R! elements of the permutation group S R . The phase shift θ ii is determined by the actual momenta p i , p i through There is θ ii = −θ i i . The actual momenta p i are related to the Bethe numbers I i ∈ [0, L − 1] as In the XXX chain we use the actual momentum p = 2πk L and the momentum k (which is actually the total wave number) with p ∼ = p + 2π and k ∼ = k + L. We have Note that the Bethe numbers I i are integers. For general L, the momenta k i are not necessarily integers or half-integers, and k i can be complex numbers for bound states. In this paper, we mainly focus on the scattering magnon excited states, for which all the momenta k i are real numbers. We use Bethe quantum numbers of the excited magnons I = {I 1 , I 2 , · · · , I R } to represent the magnon excited states in the XXX chain.

Local mode method
We For each i = 0, · · · , R, we define the new indices characterize the configurations of the subsystems A and B respectively, and write the tensor U as a We write the state (4.9) as Note that the basis is orthonormal Then we get the RDM in an orthonormal basis (4.14) We get the Rényi entropy and entanglement entropy Note that V i is well-defined only for i in the range

Ground state
The ground state (4.2) is a direct product state The Rényi and entanglement entropies in the ground state are just vanishing

Single-magnon states
When only one particle is excited, there is the state These expressions have been obtained in [32,34,39,41,42]. As S

Double-magnon states
The double-magnon states can be scattering states or bound states. We focus on the scattering states in this subsection, and we discuss the bound states in appendix D. The calculations in this subsection have overlaps with [39,41].
We consider the states with Bethe numbers I 1 , I 2 , 24) and the normalization factor Generally there are 0 ≤ I 1 ≤ I 2 ≤ L − 1. The two magnons have physical momenta p 1 , p 2 and momenta Note that k 1 , k 2 may not necessarily be integers and may be possibly complex numbers. The total actual momentum, total momentum, and total Bethe number of the state are The total Bethe number I is an integer in the range [0, 2L − 2]. The angle θ is determined by the equation To the equation (4.29), there are three classes of solutions [68]. The first two classes are scattering states, and the third class are bound states. The class I solutions have I 1 = 0, I 2 = 0, 1, · · · , L − 1, and θ = 0. The class I solutions are just excitations of hardcore bosons. The class II solutions have Most of the class II solutions have I 1 ≤ I 2 − 2, and there are also cases with The class III solutions have I 1 = I 2 or I 1 = I 2 − 1. We rename all the two-magnon solutions. We call the special class I solutions with I 1 = I 2 = 0 case I solutions. We call both the class I solutions with I 2 = 0 and all class II solutions case II solutions. We call all class III solutions case III solutions. We discuss the Rényi and entanglement entropies in case I and case II solutions in the following two subsections, and we discuss the Rényi and entanglement entropies in case III solutions in appendix D.

Case I solution
For the case I solution, there are The state is In the orthonormal basis of the states in the subsystem A we have the diagonal RDM

Case II solutions
For case II solutions, there are Note that The angle θ ∈ [−π, π] is a real solution to the equation (4.29). We use the nonorthonormal basis with the definition of U j 1 j 2 (4.24). For the state |I 1 I 2 there are the 4 × 4 matrices with the definitions α 12 ≡ α k 1 −k 2 , β 12 ≡ β k 1 −k 2 , where α k and β k are defined as (B.4) and (B.5), and as well as p 12 ≡ p 1 − p 2 . The matrix R A,I 1 I 2 (4.42) has four eigenvalues ν i with i = 1, 2, 3, 4 which are The Rényi and entanglement entropies are We are interested in the Rényi and entanglement entropies (4.45) and (4.46) in the scaling limit L → +∞ and → +∞ with fixed x = L . The shift angle becomes with the scaled Bethe numbers 48) and the function There is 0 ≤ ι 1 ≤ ι 2 ≤ 1. As the Bethe numbers satisfy I ∼ = I + L, we have the scaled Bethe numbers where S (n),bos A,k 1 k 2 and S bos A,k 1 k 2 are just the double-particle Rényi and entanglement entropies in the extremely gapped bosonic chain (C.32) and (C.33).
For the cases with ι 1 = ι 2 = ι ∈ (0, 1), there is shift angle θ = π and the states are excitations of two fermions, we have the Rényi and entanglement entropies where S (n),fer A,k 1 k 2 and S fer A,k 1 k 2 are just the double-particle Rényi and entanglement entropies in the extremely gapped bosonic chain (B.47) and (B.48).
For the remaining cases with 0 ≤ ι 1 < ι 2 ≤ 1 excluding the case with ι 1 = 0, ι 2 = 1, there is the shift angle θ ∈ [0, π), whose explicit value is not important to us. The two magnons have a large momentum difference and their contributions to the Rényi and entanglement entropies are independent. We get the Rényi and entanglement entropies S (n),XXX A, They are just the universal Rényi and entanglement entropies obtained in [41,42].
The exact results of the Rényi and entanglement entropies in the case II scattering double-magnon states of the XXX chain could be calculated either numerically from the local mode method (4.15) and (4.16) or from the analytical expressions (4.45) and (4.46). We compare the results with the corresponding analytical ones in the scaling limit in figure 6.

General magnon state
We consider the most general scattering magnon state |I (4.3) with Bethe numbers I = {I 1 , I 2 , · · · , I R }, and we focus on the states with finite numbers of magnons. More precisely, we require . For the exact results in the XXX chain, we have set L = 128.
In the scaling limit, there are fixed scaled Bethe numbers As there are finite number of magnons, the shift angles approach to the same values as those in the double-magnon state with the function f (ι i 1 , ι i 2 ) being the same as (4.49), i.e. that in the scaling limit the shift angles between two magnons only depend on the two relevant Bethe numbers not on the Bethe numbers of other magnons. For example, for the state with seven magnons (I 1 , I 2 , I 3 , I 4 , I 5 , I 6 , I 7 ) = {1, 3, L 4 , L 4 + 2, L 4 + 4, L 2 , L 2 + 2} in the scaling limit, we have the shift angles θ i 1 i 2 , i 1 , i 2 = 1, 2, · · · , 6 forming a 7 × 7 matrix The momenta K a are determined by the Bethe numbers I a = {I ab |b = 1, 2, · · · , β a } from (4.8). We sort the cluster a Bethe numbers I ab in order from the smallest to the largest, and we get K a = {K ab |b = 1, 2, · · · , β a } with Especially, for the case with The formulas (4.60), is more general than (4.60). The result (4.60) is for the case with finite number of magnons, while we conjecture that (4.65) is valid even for the case with infinite number of magnons in the scaling limit.

XXZ chain
In this section we show that the conclusions derived in the spin-1/2 XXX chain the previous section can be extended to more general interacting quantum Bethe integrable chains such as the spin-1/2 XXZ chain. The spin-1/2 XXZ chain in transverse field has the Hamiltonian We only consider the ferromagnetic phase with h > max(0, 1−∆ 2 ), and so that the unique ground state is (4.2) and low energy excitations are magnons.
The XXZ chain can be solved with the Bethe ansatz [67]. The most general magnon excited state |I is still in the form (4.3) with (4.4) and (4.5). We require that the state |I has finite number of magnons. The shift angles are determined by the new equation In the scaling limit L → +∞, the shift angle between two magnons with Bethe numbers I 1 < I 2 is where the scaled Bethe numbers ι 1 ≤ ι 2 are defined as (4.48).
For the Rényi and entanglement entropies in the scaling limit, the only relevant shift angles are the ones with equal scaled Bethe numbers We plot the shift angle θ = f (ι, ι) with equal scaled Bethe numbers ι 1 = ι 2 = ι in the scaling limit of the XXZ chain in figure 7. For ∆ = 1, the XXZ chain becomes the XXX chain, and we have discussed the case in details in the previous section. The same as that in the XXX chain, we group the magnons in the state |I with Bethe numbers I = {I 1 , I 2 , · · · , I R } into α clusters I a , a = 1, 2, · · · , α, according to the values of the scaled Bethe numbers. As before, the magnons with ι = 0 and the magnons with ι = 1 are grouped in the same cluster. In the scaling limit, the contributions of different clusters to the Rényi and entanglement (5.5)      where the momenta K a are determined by (4.8).
In figure 8, we show the Rényi and entanglement entropies S

Conclusion and discussion
In this paper we proposed three conjectures for the quasiparticle excited state Rényi and entanglement entropies of the subsystem A = [1, ] in the scaling limit of the general translational invariant quadratic fermionic and bosonic chains of L sites. The first conjecture is The second conjecture is The RHS (6.1), (6.2), (6.4) and (6.5) are calculated in the non-interacting model with the number operator as the Hamiltonian, and we conjecture that these equations are valid as long as the corresponding condition (6.3) or [(6.6) plus (6.7)] is satisfied. The third conjecture is A,G , (6.8) with the condition |k − k | → +∞, for all k ∈ K and all k ∈ K . (6.10) The second conjecture could be derived from the first conjecture plus the third conjecture.
We also investigated the magnon excited state Rényi and entanglement entropies in the scaling limit of the XXX and XXZ chains. In particular, we showed that as far as the number of excited magnons is small one can decompose the Rényi and entanglement entropies to the sum of the Rényi and entanglement entropies of particular excited states in the fermionic and bosonic chains with the number operator as the Hamiltonian. It would be interesting to see how this picture changes when one looks to the spinon excited states in which the number of excited modes are proportional to the size of the system. The techniques used in the papers [69,70] might be useful in this direction.
In previous works [52,54] and the present paper, we have focused on one-dimensional quantum systems. The entanglement in excited states of higher dimensional free and interacting quantum systems has been investigated in for example [45,51], suggesting that there are simple formulas as in one dimensions. It would be interesting to adapt the subsystem mode method to higher dimensions and check the results in various higher dimensional quantum systems.

A Calculations for states in nonorthonormal basis
In this appendix, we give an efficient procedure to calculate the Rényi and von Neumann entropies for density matrices in a general nonorthonormal basis.
We consider the general density matrix in a general nonorthonormal basis |φ i with the positive matrix It is convenient to define another matrix Here ρ P could be the density matrix of the total system or the RDM of a subsystem. The Rényi entropy and von Neumann entropy could be easily written as The positive matrix Q could be written in terms of the unitary matrix U and the positive diagonal with U and Λ being constructed respectively by the eigenvectors and eigenvalues of Q. Then the nonorthonormal basis |φ i could be written as with the orthonormal basis |ψ i satisfying ψ i |ψ j = δ ij . We further write the density matrix as The Rényi entropy and von Neumann entropy could be easily written as When the matrices P, Q, R, S are block diagonal with one could further write the Rényi and von Neumann entropies as In summary, for every density matrix with matrix P in nonorthonormal basis with matrix Q, we define the matrices R, S. We could calculate the Rényi and von Neumann entropies using either R or S. Actually, anything defined in terms of the density matrices could be calculated from S.

B Free fermions: RDM and entanglement
In this appendix we will review and refine the subsystem mode method for the fermionic chain with the number operator (2.3) as the Hamiltonian, or equivalently the extremely gapped limit λ → +∞ of the Hamiltonian (2.20), and then presents some new forms of the correlation matrix method and finally present some examples.

B.1 Subsystem mode method
The subsystem mode method was used in [52,54] to calculate analytically the Rényi entropy with integer index n ≥ 2 in the quasiparticle excited states of the extremely gapped fermionic and bosonic chain. In this subsection, we will refine the subsystem mode method and show how to use the method to calculate the Rényi and entanglement entropies in the quasiparticle excited states of the extremely gapped fermionic chain. The method is efficient for both analytical and numerical calculations. In fact, the procedure we prescribe is applicable to any quantity that is defined for the RDMs.
There are subsystem A = [1, ] and its complement is B = [ + 1, L]. We will mainly focus on the results in the scaling limit L → +∞, → +∞ with fixed ratio x = L . With the number operator (2.3) as the Hamiltonian, the ground state |G is annihilated by all the local lowering operators a j |G = 0, j = 1, 2, · · · , L, (B.1) and so it could be written as a direct product form |G = |G A ⊗ |G B with a j |G A = 0 for all j ∈ A and a j |G B = 0 for all j ∈ B. We divide the quasiparticle modes as with the subsystem modes There are anti-commutation relations In the scaling limit L → +∞, → +∞ with fixed x = L and fixed k, we have We consider an ordered set of momenta K = {k 1 , · · · , k r } with k 1 < · · · < k r . There is the number of quasiparticles |K| = r. We define the products of subsystem modes We emphasize that all the sets used in this subsection are ordered. The excited state (2.8) is with K\K being the complement set of K contained in K. We have defined the factor which is vanishing or just a sign In the nonorthonormal basis c † A,K |G A with K ⊆ K we write the RDM in the form of (A.1). Note that the set K has 2 |K| different subsets. We have the RDM with the entries of the 2 |K| × 2 |K| matrix P A,K There is also the 2 |K| × 2 |K| matrix Q A,K with entries We need to evaluate the expectation values c A,K 1 c † A,K 2 G and c B,K 1 c † B,K 2 G in the ground state, which are just the determinants where the |K 1 | × |K 2 | matrices A K 1 K 2 and B K 1 K 2 have the entries With the matrices P A,K and Q A,K , we follow the procedure in appendix A, calculate R A,K = P A,K Q A,K and obtain the Rényi and entanglement entropies Note that the matrices P A,K and Q A,K are block diagonal. The method is efficient for both analytical and numerical calculations. When only a few quasiparticles are excited, we could calculate the results analytically, and when more quasiparticles are excited we have to calculate the results numerically.

B.2 Correlation matrix method
To verify the results from the subsystem mode method, we calculate numerically the results using the correlation matrix method [13,[15][16][17][29][30][31]. In the generally gapped fermionic chain, one could define the Majorana modes In the general excited state |K = |k 1 k 2 · · · k r , one defines the 2 × 2 correlation matrix Γ K as The Rényi and entanglement entropies in state |K are With the number operator (2.3) as the Hamiltonian, one defines the × correlation matrix C K as with the function The Rényi and entanglement entropies in state |K are Using the position-momentum duality [71,72], we could also write the Rényi and entanglement entropies as where the |K| × |K| matrix A K has entries with i 1 , i 2 = 1, 2, · · · , |K|, k i 1 , k i 2 ∈ K and the definition of α k in (B.4). We may also define the |K| × |K| matrix B K which has entries with i 1 , i 2 = 1, 2, · · · , |K|, k i 1 , k i 2 ∈ K and the definition of β k in (B.5). In fact, there are A K = A KK , As there is always k i 1 − k i 2 ∈ Z in the fermionic chain, we have Using the determinant of the block matrix 32) and the fact that the matrices A K and B K commute, we could write the Rényi entropy (B.27) as 2 with the n|K| × n|K| matrix Ψ (n) A,K written in n × n blocks as The

B.3.2 Single-particle states
In the single-particle state |k , there is the RDM In the nonorthonormal basis there are matrices The Rényi and entanglement entropies are These results have been obtained in [32,34,39,41,42].

B.3.3 Double-particle states
In the single-particle state |k 1 k 2 , there is the RDM with α 12 ≡ α k 1 −k 2 . We have used β k = −α k for nonvanishing integer k. In the nonorthonormal basis there are matrices Alternatively, in the orthonormal basis we have the RDM In either the nonorthonormal basis (B.43) or the orthonormal basis (B.45), we get the Rényi and entanglement entropies in the double particle state |k 1 k 2 S (n),fer It is remarkable that the Rényi entropy is written as an explicit analytical function of the general

C Free bosons: RDM and entanglement
In this appendix we will review and refine the subsystem mode method for the bosonic chain with the number operator (3.2) as the Hamiltonian, or equivalently the extremely gapped limit m → +∞ of the Hamiltonian (3.10), and then present some new forms of the wave function method and finally present some examples.

C.1 Subsystem mode method
The analytical calculation using the subsystem mode method is similar to that in the fermionic chain in the previous appendix, and we will keep it brief in this subsection. We divide the quasiparticle modes c k , c † k into the subsystem modes There are commutation relations with α k and β k defined the same as those in (B.4) and (B.5).
For an arbitrary set K = {k r 1 1 , · · · , k rs s }, we have the number of momenta We define the products of subsystem modes (C.1) Then there is the excited state with K\K being the complement of K in K and the factor s K,K being defined as Then we get the RDM Note the possible momenta repetitions of the sets and subsets used in this paper, and one need to be careful to calculate the subsets and their complements. The general sets K = {k r 1 1 , · · · , k rs s } has d(K) = s i=1 (r i + 1) different subsets. For example, the set K = {1 2 , 2} has R = 3 quasiparticles and d(K) = 6 different subsets. Its subsets K and their corresponding complements K\K are Then we write the RDM in the form of (A.1). We have with the entries of the d( There is also the d( We need to evaluate the expectation values c A,K 1 c † A,K 2 G and c B,K 1 c † B,K 2 G , which are just the permanents where |K 1 | × |K 2 | matrices A K 1 K 2 and B K 1 K 2 have the same entries as (B.16).
With the matrices P A,K and Q A,K , we get R A,K = P A,K Q A,K and the Rényi and entanglement entropies S (n),bos (C.14)

C.2 Wave function method
We could also calculate the Schatten distance from the wave function method [41,42]. One could also see [54,55]. With the number operator (3.2) as the Hamiltonian, there sets up a permanent formula for the Rényi entropy [54] S (n),bos with the n|K| × n|K| matrix Ω (n) A,K written in n × n blocks as C.3.3 Double-particle states with equal momenta For the double-particle state with the same momenta, there is the RDM In the nonorthonormal basis we have the matrices We get the Rényi and entanglement entropies 26) which are just the universal Rényi and entanglement entropies found in [41,42].

C.3.4 Double-particle states with different momenta
In bosonic chain with the number operator (3.2) as the Hamiltonian, there is the RDM [54] ρ A, In the nonorthonormal basis we have the matrices Alternatively, in the orthonormal basis |ψ 0 = |G A , we have the RDM In either the nonorthonormal basis (C.28) or the orthonormal basis (C.30), we get the Rényi entropy with general index n S (n),bos and the entanglement entropy which are beyond the results in [54].

D Double-magnon bound state in XXX chain
In this appendix, we discuss the Rényi and entanglement entropies in the double-magnon bound state, i.e. the class III solutions in [68] and the case III solutions in this paper, in the XXX chain. This appendix has overlaps with [39]. The new things in this paper are that we give the explicit analytical expressions of the Rényi and entanglement entropies and we also discuss the behaviors of the results in various limits. We show explicitly that the double-magnon bound state Rényi and entanglement entropies approach the single-particle Rényi and entanglement entropies S

(n)
A,k in the strongly bound limit and approach the double-particle Rényi and entanglement entropies in the bosonic chain S (n),bos A,k 2 in the loosely bound limit.

D.1 Case IIIa solutions
There are two subcases in the case III solutions. For case IIIa solutions take the form |I 1 I 2 (4.23) with The parameter 1/v is roughly the size of the two-magnon bound state. When v → 0, the two magnons are loosely bound, and when v → +∞, the two magnons are tightly bound.
Formally, we still work in the nonorthonormal basis (4.39) with (D.1). We get the matrices Then we obtain the matrix R A,I 1 I 2 = P A,I 1 I 2 Q A,I 1 I 2 with the eigenvalues   We still have the scaled total Bethe number defined as (D.4) and behave in the scaling limit as (D.5).
The parameter v is in the range 2π 2 L 2 ≤ v ≤ +∞. (D.20) The parameter 1/v is still the size of the two-magnon bound state.
In the nonorthonormal basis (4.39) with (D.17), we get the matrices Then we obtain the matrix R A,I 1 I 2 = P A,I 1 I 2 Q A,I 1 I 2 with the eigenvalues In the scaling limit L → +∞, → +∞ with fixed x = L and finite v or infinite v, we have (D. 13) and the Rényi and entanglement entropies approach (D.14) and (D.15).
The smallest v behaves as v = w L 2 with fixed w ≥ 2π 2 in the scaling limit L → +∞, and the eigenvalues (D.23) become (D.24) The Rényi and entanglement entropies are If v behaves as v = u L with fixed u ≥ 2π 2 in the scaling limit L → +∞, the eigenvalues (D.23)