Entanglement Content of Quantum Particle Excitations I. Free Field Theory

We evaluate the entanglement entropy of a single connected region in excited states of one-dimensional massive free theories with finite numbers of particles, in the limit of large volume and region length. For this purpose, we use finite-volume form factor expansions of branch-point twist field two-point functions. We find that the additive contribution to the entanglement due to the presence of particles has a simple"qubit"interpretation, and is largely independent of momenta: it only depends on the numbers of groups of particles with equal momenta. We conjecture that at large momenta, the same result holds for any volume and region lengths, including at small scales. We provide accurate numerical verifications.


Introduction
Measures of entanglement, such as the entanglement entropy, have attracted much attention in recent years, particularly in the context of one-dimensional many body quantum systems (see e.g. review articles in [1,2,3]). Among such systems, those enjoying conformal invariance in the scaling limit are of particular interest as they provide a theoretical and universal description of critical phenomena. In their seminal work Calabrese and Cardy [4] used principles of conformal field theory (CFT) to study the entanglement entropy (EE) [5] of quantum critical systems. Their results generalised previous work [6], provided theoretical support for numerical observations in critical quantum spin chains [7] and highlighted the fact that the EE encodes universal information about quantum critical points, such as the central change of the corresponding CFT and, in more complex setups, about the full primary operator content of CFT [8,9,10].
The von Neumann and Rényi entanglement entropies are measures of the amount of quantum entanglement, in a pure quantum state, between the degrees of freedom associated to two sets of independent observables whose union is complete on the Hilbert space H = H A ⊗ H B .
In particular, the combination of a geometric description, Riemann uniformization techniques and standard expressions for CFT partition functions is very fruitful. Beyond criticality, EEs are accessible by means of the branch point twist field approach introduced in [24] and also through numerical techniques.
Consider a bipartition where the two sets of observables correspond to the local observables in two finite-size complementary connected regions, A and B (see for instance Fig. 1). Let the system by in a state |Ψ L , then the von Neumann entropy associated to region A is where ρ A = Tr B (|Ψ LL Ψ|) is the reduced density matrix associated to subsystem A and the trace (1.1) is over the degrees of freedom in subsystem A. One may obtain the entropy S Ψ 1 ( , L) as a limiting case of the sequence of nth Rényi entropies defined as One may also consider the so-called single-copy entropy [25,26,27], defined as Much of the work carried out so far deals with the entanglement properties of the ground state (mostly, but not always, in infinite systems). In conformal field theory, universal results 1 Starting from a lattice system with a critical point for some value of a parameter λ = λc, the scaling limit to a critical point described by CFT may be taken by first setting λ = λc so the correlation length ξ → ∞ and then taking the thermodynamic limit L → ∞. The near-critical behaviour of massive QFT is recovered by taking the limit λ → λc and L, → ∞ simultaneously, whilst keeping L/ξ and /ξ fixed. This is the regime we consider in this paper.
1 for certain types of excited states are known: in [28,29], the increment of Rényi entropy in an excited state |Υ with respect to the ground state of a CFT for the configuration of Fig. 1 was found to be S Υ n (r) − S 0 n (r) = (1 + n)(h +h) 3n (πr) 2 + O r 2∆ ψ , (1.5) for small values of r = L . The excited state was defined as where Υ(ξ,ξ) is a CFT field, h,h are its holomorphic and antiholomorphic dimensions, ξ,ξ are coordinates on the cylinder, and ∆ ψ = h ψ +h ψ is the smallest scaling dimension of any field in the theory. Therefore, a measurement of the EE of a low-lying excited state in CFT at finite volume can provide information about the primary field content of the theory. The most extensive numerical study of other kinds of excited states in critical systems was performed in [30]. In this work a very detailed study of the excited states of the XY model in a transverse field and the XXZ Heisenberg spin-chain was carried out. The authors focussed on the case when L 1 and on excited states that are macroscopically different from the ground-state (we will consider instead zero-density states). The EE of excited states with finite energy density in quantum field theory (QFT) or quantum lattice models is very simple by the eigenstate thermalization hypothesis (or its extension to integrable systems): it is dominated by the thermodynamic entropy of the corresponding Gibbs (or generalized Gibbs) state and is known to satisfy a volume law [3]. However, little is known so far about the EE of zero-density excited states in gapped systems. The most extensive numerical study in gapped quantum spin chains was carried out in [31] where some of the results we obtain here (see Section 2) were proposed as describing the "semi-classical" limit of the EE. Indeed the authors observe how the EE of certain excited states approaches such semi-classical limit for large enough volumes and appropriate correlation lengths. In our work [32] we have shown that these bounds, and generalisations, provide, in fact, exact large-volume predictions that are much more widely applicable.
In the present paper, we provide full analytical computations supporting some of the results in [32]. We consider excited states of 1+1-dimensional massive QFT with zero energy density: those formed of finite numbers of asymptotic particles, at various momenta. We consider the situation depicted in Fig. 1, in the limit where both the systems size L and the length of region A are large and in fixed proportion Let |Ψ L be such an excited state. Employing the branch point twist field approach [24], we compute the difference between the Rényi entropy in the excited state and in the ground state, in this limit, lim where T is the branch point twist field andT is its hermitian conjugate [24]. Recall that branch point twist fields are local fields of the n-copy "replica" QFT, the theory constructed as n not mutually-interacting copies of the model under study. In the replica theory, the state |Ψ L has the structure |Ψ L = |Ψ 1 L ⊗ |Ψ 2 L ⊗ · · · ⊗ |Ψ n L , (1.10) where |Ψ i L is an excited state of the i-th single-copy theory in finite volume L. We concentrate on the (uncompactified) massive free real boson and free Majorana fermion models. The techniques that we use -based on form factors of branch point twist fields -have been chosen so that they are (hopefully) generalizable to integrable models, in view of extending our results in a future work.
The results we find are very surprising, for various reasons: • All results are independent of the momenta of the excitations, except for the sole condition of coincidence or not of rapidities, and are independent of the model under consideration.
• The structure of all functions ∆S Ψ n (r) is extremely simple. They in fact admit a combinatorial, or qubit interpretation, related to counting all possible configurations with various numbers of excitations (particles) "located" in the region A and outside of it.
• Our numerical analysis also suggests that the formulae above hold very precisely even for arbitrary systems size L, no matter how small, if the momenta of the excitations are large (even though our calculation methodology employs a large volume expansion).
• Additional numerical analysis presented in [32] has shown they hold also in higher dimensional free theories and at least some states of interacting quantum spin chains.
The paper is organized as follows: In Section 2 we review our results for the increment of EE for states with a finite number of excitations. The formulae presented in Section 2 as well as their "qubit" interpretation appeared first in [32]. Here we present a more general discussion of the "qubit" interpretation. In Section 3 we review the connection between branch point twist fields in replica theories and Rényi entropy. We also highlight the challenges of generalizing such connection to finite volume and excited states. We explain how these challenges may be resolved in the case of the massive free boson theory and introduce the "doubling trick" in this context. In Section 4 we derive the general formulae for the Rényi entropy of a single-particle excited state, a k-particle excited state involving distinct momenta only, and a k-particle excited state consisting of equal momenta. We provide concrete examples of all three cases for the 2nd Rényi entropy of the massive free boson theory. We compare the analytical results to numerical results obtained by employing the wave functional method. In Section 5 we generalize the results of the previous section to the massive free fermion. We find that the expressions for the EEs of states with distinct momenta are identical to those in the free boson theory, even if there are technical differences in the computations involved. In Section 6 we present our conclusions and outlook.
In Appendix A we review the wave functional method and its application to the computation of the Rényi entropies of the harmonic chain. In Appendix B we present a derivation of the selection rules which single out those terms in the form factor expansion that provide the leading large-volume contribution to the Rényi entropies. In Appendix C we prove some properties of the functions g n p (r) in terms of which all EEs can be expressed.

Summary of the Main Results
The computation of the ratio (1.9) for a generic k-particle excited state of a massive free theory in finite volume involves the use of a considerable number of techniques we will be presenting in the next sections: the form factor programme for branch point twist fields [24], the generalization of this programme for finite volume correlators following the ideas of [33,34], the rewriting of the branch point twist field in terms of U (1) fields of the replica free theory by employing the "doubling trick" introduced in [35]. We then use a new numerical technique based on wave functionals in order to test our analytical results. This is therefore a rather technical work. However, the results that we have obtained are surprisingly simple and can be easily summarized. They have been shown to hold more widely in [32].

Main Formulae
Consider a state consisting of a single particle excitation. Let us denote the entropy increments of such a state by ∆S 1 n (r). We find that The increment of von Neumann entropies is given by 2) and the increment of single-copy entropies has the form For excited states consisting of a finite number k of excitations of distinct momenta the results are simply as above, multiplied by k. In the free boson, we may also consider states consisting of k particles of equal momenta. We will denote the entropy increments of such states by ∆S k n (r). We find ∆S k n (r) = The single-copy entropy is a function which is non-differentiable at k points in the interval r ∈ (0, 1) (generalizing (2.3) which has one non-differentiable point). The positions of these singularities are given by the values and the single-copy entropy is given by Therefore, if the rapidities are distinct, the contribution to the entanglement entropy of k particles is exactly k times the contribution of a single particle excitation, while if they are equal, this is not true: the contribution is in fact smaller. For generic states containing a mixture of excitations with equal and distinct rapidities we find formulae which are sums of those above. Denoting by ∆S k 1 ,k 2 ,··· n (r) the Rényi entropies of an excited state consisting of k i particles of momentum p i with p i = p j for i = j we find   .7) for a state of 10 equal momenta and for three "mixed" states with some equal and some distinct momenta. We plot the Rényi entropies for n = 2, 3, 5, 8, 11, 17 and the von Neumann and single-copy entropies. In each figure, the dashed (outer-most) curve is the von Neumann entropy and the dot dashed (inner-most) curve is the single-copy entropy.
It is easy to show that all the differences of von Neumann entropies have their maximum value at r = 1 2 . For states with k distinct rapidities this maximum value is given simply by k log 2 so that a k-particle excited state may at most add k qubits to the entanglement entropy with respect to its ground state value. This fact was discussed in [36] for one-particle excitations and shown to hold beyond free theories, for integrable and non-integrable theories.
For states with some equal rapidities, the maximum is lower. In particular for k coinciding rapidities, it is given by ∆S k

Qubit Interpretation
It turns out that the general formulae (2.1)-(2.7) have interpretations as the entanglement entropies of simple states formed of qubits, and are easily understandable from a quasi-classical particle picture of the actual QFT states considered. This was discussed in [32], and we give here slightly more precision.
In order to explain this, consider a bipartite Hilbert space H = H int ⊗ H ext . Each factor H int H ext is the Hilbert space for N j distinguishable sets each of j indistinguishable qubits, for j = 1, 2, 3, . . .. Making the relation with the entanglement problem described above, we associate H int with the interior of the entanglement region of length and H ext with its exterior, and we identify the qubit state 1 with the presence of a particle and 0 with its absence. With k particles lying on (0, L), we construct the state |Ψ qb ∈ H by the (naive) picture according to which equal-rapidity particles are indistinguishable, and a particle can lie anywhere in (0, L) with flat probability: any given particle has probability r of lying in the entanglement region, and 1 − r of lying outside. We make a linear combination of qubit states following this picture, with coefficients that are (square roots of) the total probability of a given qubit configuration, taking proper care of (in)distinguishability. Then, the Rényi and von Neumann entanglement entropies of |Ψ qb are given exactly by the formulae seen earlier. In general 10) and the statement is that S Ψ qb n (r) = ∆S Ψ n (r) for some excited state |Ψ L corresponding to the probability distribution described above.
More precisely, Here C j+1 is the Hilbert space of j indistinguishable qubits, with basis elements |q , q = 0, 1, . . . , j labelled by the number of qubits that are in their state 1. One can also write H int where N is the total number of groups, N = j≥1 N j , and j i take values j 1 = · · · = j N 1 = 1, j N 1 +1 = · · · = j N 1 +N 2 = 2, etc. We denote the basis of vectors in H int H ext by |q for q = (q i : i = 1, . . . , N ) ∈ j≥1 {0, 1, . . . , j} N j . We use the notationq = (j i − q i : i = 1, . . . , N ) for the state where the qubits are inverted. We then construct where p q is the probability of finding the particle configuration q in the entanglement region according to the naive picture above, given by For instance, if a single particle is present, then the state is as either the particle is in the region, with probability r, or outside of it, with probability 1 − r.
If two particles of coinciding rapidities are present, then we have as either the two particles are in the region, with probability r 2 , or one is in the region and one outside of it (no matter which one), with probability 2r(1 − r), or both are outside the region, with probability (1 − r) 2 . For two particles of different rapidities,

Rényi Entropies and Branch Point Twist Fields
It has been known for some time that several entanglement measures, including the Rényi entropies, can be expressed in terms of correlation functions of a special class of local fields T which have been termed branch point twist fields in [24]. Branch point twist fields are, on the one hand, twist fields in the broader sense, that is, fields associated with an internal symmetry of the theory under consideration, and on the other hand related to branch points of multi-sheeted Riemann surfaces. They are twist fields associated to the cyclic permutation symmetry of a model composed of n copies of the original model, with exchange relations where O i (y) is any local field on copy number i, and with O n+1 (y) = O 1 (y). The idea of quantum fields associated with branch points of Riemann surfaces in the context of entanglement appeared first in [4], where their scaling dimension was evaluated in CFT (see also [37] for an earlier work concerned with similar ideas in the context of orbifold CFT).
Here c is the central charge and n is the number of sheets in the Riemann surface. The general description in terms of branch point twist fields as symmetry fields associated to cyclic permutation symmetry of the n Riemann surface's sheets, as per (3.1), was given in [24], where they were studied in integrable massive QFT. This description is however independent of integrability, and it was first used in massive QFT outside of integrability in [38]. The missing logical link that connects the Riemann surface structure mentioned above with a computation of entanglement measures comes through a result commonly known as the replica trick. Mathematically speaking, the replica trick is simply the statement (1.3) with (1.2). However, the word "replica" originates from the fact that the object Trρ n A which features in (1.2) can be interpreted as the partition function of a replica QFT understood as n non-interacting copies of the original QFT. In the limit L → ∞ (for the configuration in Fig. 1 with L → ∞), this partition function is evaluated precisely on a Riemann surface with n sheets as described earlier, with a branch cut of length across which Riemann sheets are connected cyclically (when the branch cut starts at the origin and L → ∞ this is exactly the structure of the Riemann surface of the function n z z− ). Hence the number of sheets and the number of replicas are both n. For finite volume L, the Riemann sheets are replaced by cylinders of circumference L cyclically connected along a branch cut in the compactified (space) direction. In this picture, the nth Rényi entropy with the partitioning protocol presented in Fig. 1 is given by: where |Ψ L is an excited state of a finite number of excitations in the finite-volume L, replica QFT. The structure of the state is as reported in (1.10),T = T † is the hermitian conjugate of the branch point twist field T , and ε is a non-universal short-distance cut-off. Notice that the dependance on the cut-off ε cancels out when considering the entropy increment (1.9). 7

Challenges Posed by the Treatment of Excited States
For L → ∞ in the ground state the function (3.4) has been extensively investigated, both from the point of view of its universal features [24,38] and for particular models [39,40,41,42,43].
The study of excited states however presents new challenges. First, in the context of integrable models of massive QFT, it is natural to evaluate twopoint functions by inserting a sum over a complete set of states and then computing the matrix elements of local operators that are the building blocks of the resulting sum. Schematically we can represent this process by writing (3.5) The advantage of this decomposition is that for integrable models at least, there exist effective methods to exactly compute the matrix elements L Ψ|T (0)|Φ L . In infinite volume such methods are usually refered to as the form factor programme [44,45] and they provide the most powerful and successful approach to the computation of correlation functions, both analytically and numerically. For branch point twist fields, the form factor programme was generalized in [24]. For finite volume and local fields O -excluding twist fields -matrix elements of the type L Ψ|O(0)|Φ L are also well understood [33,34]. In 1+1 dimensions, the states |Φ L and |Ψ L are characterized by a discrete set of rapidities (or momenta). Should any of the rapidities in one set coincide with some in the other set, the matrix element L Ψ|O(0)|Φ L for L → ∞ will develop, in the usual infinite-volume normalization of the states, δ-function singularities. Considering instead finite volume form factors, provides a natural regularization scheme to deal with such singularities. Indeed, for local operators, a systematic prescription exists to compute the "physical part" of matrix elements such as L Ψ|O(0)|Φ L by subtracting the contributions of any occurring singularities in a way which is controlled by the particular pole structure of the form factors of local fields.
In our case however, we face the challenge that the branch point twist fields are not local in the sense required to apply the techniques of [33,34]. Although they are local with respect to the Lagrangian density of the replica model (as they implement a symmetry) they are non-local with respect to the fundamental fields of the theory (those whose associated modes create and annihilate the physical particles). It is however, still very plausible that the standard general ideas for the computation of finite-volume non-diagonal form factors will be applicable to branch point twist fields. We confirm this below by analytical and numerical results in free theories.
Second, branch point twist fields sit at the origin of branch cuts which, in the standard prescription, originate at the twist field and extend indefinitely in the space direction. For the two-point function, the two branch cuts emerging from the twist field and its hermitian conjugate combine to create a branch cut of finite length which is interpreted as the length of subsystem A. However, once we write down the expansion (3.5) we need to evaluate the matrix elements L Ψ|T (0)|Φ L and L Φ|T ( )|Ψ L . For these matrix elements, an infinitely long branch cut extending in space is incompatible with working in finite volume L. This conflict can be resolved by adopting an approach which is reminiscent of that taken in [46] for the Ising field theory and the matrix elements of its Z 2 twist field σ. We may use the fact that the branch cut can be continuously deformed without changing the value of the correlation function. Therefore we may continuously "stretch" the branch cut along the time direction as indicated in Fig. 3. The result is a product of fields with branch cuts extending in the time direction. In this configuration, the fields are well defined in the quantization on the circle, where they are intertwining operators. The operator ordering of the two-point function in the quantization scheme on the circle, is implemented in the path integral by a time ordering: an infinitesimal shift τ along the cylinder, as in Fig. 3. In parallel to the situation in [46], the Hilbert space of quantization on the circle is divided into sectors characterized by periodicity conditions: if an internal symmetry σ exists, then the Hilbert space H σ is that of field configurations with the quasi-periodicity condition O(x + L) = σ · O(x). For the Ising model, the Z 2 symmetry leads to two sectors, Ramond-Ramond and Neveu-Schwarz. In the case of our replica model, we have in particular n sectors labelled by cyclic elements of the permutation group. The intertwining operators corresponding to the branch point twist fields act as follows: where ω is the elementary cyclic permutation symmetry of the n-copy replica model, taking copy i to copy i + 1 mod n. This is seen as follows: the condition (3.1) imposes continuity between O i below and O i+1 above the branch extending towards the right. After the deformation as in Fig. 3, this becomes continuity between O i on the left and O i+1 on the right of the branch extending towards negative times. This adds a factor of ω on the Hilbert space on which T acts, or equivalently, a factor ω −1 on the image Hilbert space. Therefore, in the matrix elements L Ψ|T (0)|Φ L and L Φ|T ( )|Ψ L , the state |Ψ L is in a different sector than the state |Φ L . In the cylinder picture of Fig. 3, the state |Φ lies between the twist fields, in the time slice of extent τ introduced by the operator ordering. Finally, the question arises as to how the matrix elements of branch point twist fields with states in different sectors can be computed. Answering this question in general integrable QFT is somewhat complicated and will be discussed in a future work. However, for free theories there are additional resources at our disposal. More precisely, for free theories, it is possible to express the branch point twist fields in terms of simpler U (1) twist fields, where the permutation symmetry has been diagonalized. This is achieved by employing the so-called doubling trick introduced in [35] and employed successfully in the branch point twist field context in [24,47], where it allowed for the computation of the vacuum expectation value of the branch point twist field. A similar idea was also used in [48] in the study of the EE of free theories.
The doubling trick is the simple idea that a real free fermion (Majorana) and a real free boson theory can be doubled to construct a complex free fermion (Dirac) and a complex free boson theory. This doubling induces a U (1) symmetry in the new theory to which a U (1) twist field is associated. In a replica theory whose fundamental building blocks are doubled free theories, the U (1) symmetry on each individual copy is extended to a U (n) symmetry, which includes cyclic permutation of the copies. Diagonalizing the cyclic permutation, in the new basis the branch point twist field is then expressed as a product of n individual U (1) twist fields T p for p = 1, . . . , n.
Having summarized the main challenges and techniques involved in the computation of Rényi entropies of excited states in finite volume, we proceed now to present these techniques in some detail for the case at hand.

Doubling Trick and Replica Free Boson Model
In this and the remaining subsections, we concentrate on the free boson model. We then generalize the construction to the free fermion in Section 5.
In [35] Fonseca and Zamolodchikov introduced the "doubling trick". There, it was employed to find differential equations that are satisfied by certain combinations of correlation functions in the Ising model. This technique was later used in order to obtain vacuum expectation values T in infinite volumes in the works [24] (free fermion) and [47] (free boson).
The idea is to "double" the free theory in order to have an additional continuous symmetry. Let φ a and φ b be two independent free massive real bosons. We construct a free massive complex boson as: which has an internal continuous U (1) symmetry. This symmetry can then be exploited in order to obtain information about the original (not doubled) theory. In the context of the branch point twist field, the doubling trick is used as follows. In the doubled replica model, the combination of the U (1) symmetry of the complex field on each replica, and of the permutation symmetry of the replica, implies the existence of a U (n) symmetry of the model. Cyclic permutations form a subgroup of the U (1) symmetry group of rotations amongst the copies, which can be diagonalized. The diagonal basis is a new set of n independent complex free bosons, each of which has its own U (1) symmetry, and the cyclic permutation action is expressed as a product of U (1) actions on each of these bosons. Therefore, the branch-point twist field acts as a product of U (1) twist fields in the diagonal basis.
In the replica theory we have n copies of the complex free boson, Φ j with j = 1, . . . , n. Since the components φ a,j , φ b,j are commuting fields and the permutation symmetry ω acts in a factorized way as ω a × ω b , the branch point twist field also factorises, Therefore, correlators of T in any state that is factorized into the copies a and b, also factorize into those of T a and T b in the real boson theory. The idea is to perform calculations in the replica complex free boson theory and interpret the results in terms of the real free boson using this factorization.
In matrix form, the permutation symmetry ω acts as The eigenvalues of this matrix are exactly the nth roots of unity λ p = e 2πip n for p = 1, . . . , n. The cyclic permutation action is diagonalized by the fields which are themselves canonically normalized complex free boson fields. Since ω acts diagonally on the basisΦ p , it can be factorized into a product of U (1) fields. We denote by T p the U (1)field acting nontrivially on sector p, andT p its hermitian conjugate. The field T p has exchange relations for q, p = 1, . . . , n with q ≡ q + n and p ≡ p + n, andT p has similar exchange relations with complex conjugate phase. Then, where, by definition, the field T n is the identity field. For free bosons, such U (1) fields have been studied and it is known that they have scaling dimensions [49] which coincides with (3.3) for c = 2 (the central charge of the complex free boson). In order to study the entanglement entropy of excited states in finite volume L, we will consider states of the complex boson theory which are k-particle states in copy a times the vacuum in copy b, In this factorized state, we have The second factor is the vacuum expectation value, which is known. We therefore obtain the required real free boson result as In order to describe the many-particle states |k L more precisely, we introduce the creation and annihilation operators (a ± j ) † (θ) and a ± j (θ), respectively, of the complex free boson Φ j ; these create / annihilate a particle of rapidity θ and U (1) charge ± in replica copy j. The creation operator on doubling-trick copy a and replica copy j is expressed as Therefore, the normalized k-particle state (3.15) is, explicitly in the case of distinct rapidities, In the free boson theory, one may consider states with some coinciding rapidities, in which case the normalization of the state needs to be slightly modified. This will be discussed in more detail in Subsection 4.2.2.
On the other hand, in the diagonal basis (3.10), the operatorsã ± p (θ) (and hermitian conjugate) are given by Expressing the operators a ± j (θ) in terms of the tilde operators through (3.20) leads, after some manipulations, to This is a useful expression, because thanks to (3.12), correlation functions of twist fields in the diagonalized basis factorize into correlations on the sectors p = 1, . . . , n. Let us introduce the short-hand notationã Then, the following state factorizes as where we write |0 L = ⊗ n j=1 |0 j;L . Using this, for n = 3 we have for instance where we used the fact that since T 3 is the identity field for n = 3. In this way, any two-point function can be expressed as a sum of factorized correlators involving only particles and U (1) twist fields acting on a particular sector of the theory. A detailed computation for k-particle states of equal and distinct momenta will be presented below. The computation of matrix elements such as (3.24) requires two additional ingredients: first, the introduction of finite volume form factors, and second, the understanding of how particle rapidities are quantized in finite volume. We address these questions in the next two subsections.

Infinite Volume Form Factors of U (1) Fields
As explained in the previous subsection, explicit computations of the Rényi entropy may be obtained by computing matrix elements of U (1) twist fields. Let us review here some of the properties of these form factors in the free boson theory. The form factors have been known in the literature for quite some time [50,47]. We define the two particle form factors of the p-th U (1) field as The last two form factors are vanishing for symmetry reasons (the twist field preserves the total U (1) charge). The form factor programme for quasi-local fields [44,45,51] tells us that the nonvanishing form factors may be computed as the solutions to a set of three equations. First, Watson's equations where γ ± p are the factors of local commutativity associated to the bosons ±. From the exchange relations (3.11) we expect that n . Finally, the kinematic residue equation is is the vacuum expectation value. Based on the equations above it is easy to make a general ansatz: where A and a are constants to be determined. It is then easy to show that the equations are satisfied if a = p n − 1 2 and A = −τ p sin πp n . (3.31) This gives the solution (3.32) Another solution can be obtained by shifting j → j + n but if we assume p ≤ n the solution above is singled out. Since the theory is free, higher particle form factors can be obtained by simply employing Wick's theorem. For the complex free boson they have the structure where we introduced the normalized two-particle form factor and σ are all elements of the permutation group S m of m symbols.
In what follows we will require the form factors (3.33) as well as slightly more general matrix elements. These can be related to form factors as as long as θ i = θ i and β i = β i for all i. That is, any matrix element can be written in terms of form factors as long as there are no repeated rapidities leading to additional singularities [45].

Finite Volume Matrix Elements: A Simple Example
Once the correlation function has been expressed in terms of correlators acting on a particular sector, the latter can be computed by the insertion of a complete set of states. In finite volume both the rapidities of the excited state and intermediate states are quantized. We will use the following simple example to explain what these quantization conditions are in general. Consider a simple matrix element on sector p of the form We will think of this matrix element as a particular building block of a more complicated twopoint function. This means that the external state k i=1 a † p (θ i )|0 p;L depends on rapidities {θ i } which are the same rapidities of the original excited state |k L in (3.19). Here |q p;L are q-particle states of the form 37) and the sum over intermediate states is a sum over q = 0, . . . , ∞ and over β i s. Charge conservation requires that 2s − q = k. (3.38) In finite volume L one must choose a quantization sector in order to determine the set of values the rapidities {θ i } and {β i } may take. Below we choose the state |k L to be in the trivial quantization sector, where the field is periodic, Φ j (x + L) = Φ j (x) for all j. In each copy this generates the Hilbert space H 1 . According to (3.6), the twist fields T andT change quantization sector as follows:T : where H ω is the Hilbert space with quasi-periodicity condition Φ i (x + L) = Φ i+1 (x). Therefore, in the two-point function (3.36), the intermediate states are in the quantization sector H ω . As per (3.11), in the diagonal basis, H ω has quasi-periodicity conditionΦ p (x + L) = e 2πip nΦ p (x). This means that the quantization of momenta (rapidities) is as follows: for the external state (as these are the rapidities of the excited state), and for the intermediate states (3.37). Note that the different signs in (3.41)-(3.42) are associated with particles created by operators a † p (β i ) and b † p (β i ), respectively. These quantization conditions provide the generalization of the Bethe-Yang equations [52,53] (in the free case) in the presence of the branch cut induced by the U (1) twist field T p and can be naturally extended to more general external states.
With this information the finite-volume correlator can be expanded as (the full details of this expansion will be discussed in Section 4) Although (3.43) only shows the form factor expansion of a particular correlator, the above analysis easily extends to any other cases. We note that the expansion (3.43) may alternatively be expressed by replacing the sums {J ± i } by a set of contour integrals such that the sum over residues enclosed by the contours reproduces the original sum. This technique turns out to be rather useful in order to generalize the computation above to any external state. We will make full use of it in Subsection 4.1.2.
The final ingredient needed to evaluate (3.43) are the finite-volume non-diagonal form factors inside the sums. Fortunately, it is known [33,34] that such matrix elements can generically be related to the infinite-volume form factors (3.33) simply as , up to exponentially decaying corrections O(e −µL ). The functions in the denominator are the so-called density functions of the left-and right-states, respectively. In general, these can be computed from the Bethe-Yang equations [52,53]. However, for free theories they are simply products over the particle energies times the volume, that is, with E(θ) = m cosh θ. The form factor in the numerator is exactly the same function as in the infinite volume expression (3.35) up to the quantization conditions on the rapidities discussed earlier.
We now know that finite-volume form factors are proportional to infinite-volume ones up to quantization of the rapitidities. It is worth noting here an important property of the form factor (3.32), namely its leading behaviour near the kinematic singularity. Consider the form factor f n p (β 1 − θ 1 + iπ) and suppose that the rapidites are quantized through Bethe-Yang equations of the form (3.40) for θ 1 and (3.41) for β 1 . Then the leading contribution for θ 1 ≈ β 1 can be expressed as . (3.47) Later computations will often involve the evaluation of the modulus square of f n p (θ) near a kinematic pole, giving rise to sums of the form A proof of the equality (3.48) and a discussion of some other properties of the functions g n p (r) is presented in Appendix C.

Rényi Entropy in the Massive Free Boson
We have now reviewed all the techniques necessary to perform a computation of the Rényi entropy of the configuration in Fig. 1 in a generic zero density excited state of the massive free boson theory. Below, we describe in detail the computation for the cases of a single-particle excitation, a k-particle excitation with distinct rapidities, and a k-particle excitation with equal rapidities. In each case we illustrate our method for n = 2 and, for many-particle excitations, we choose the simplest state with k = 2.

Single-Particle Excited States
We will start by considering the simplest type of excited state, namely a one-particle excited state of rapidity θ, that satisfies the Bethe-Yang equation (3.40) with quantum number I. The excited state (3.19) for k = 1 has the form As explained in the previous section, such a state admits a more intuitive expression after changing to the new basis of creation operators (3.20), as per (3.21). Here we write it as where the C n ({N ± }) coefficients contain all the phase factors from the transformation (3.20), and the summation runs over the integer sets These are the boson occupation numbers of particles/antiparticles in each sector. As seen before, both the branch point twist fields and generic states factorize into sectors so that the two-point function of branch point twist fields in the excited state (4.1) at finite volume may be expressed using (4.2) as where * denotes complex conjugation, and is the finite-volume two-point function in sector p. Note that both here and later, the order of the creation and annihilation operators is irrelevant as they all commute in the free boson case. In sector n, the U (1) twist-fields coincide with the identity, hence the two-point function is only nonzero if N ± n =Ñ ± n , and its value is just the normalization of the finite-volume states For other sectors however, the matrix elements (4.5) are non-trivial. As standard, they can be obtained by inserting a complete set of states between the two fields so that (4.5) becomes a sum over products of the form factors (3.33). Explicitly, without any restriction. We can now insert the complete set of states (4.8) into the two-point function (4.5). Employing the action of the translation operator on energy eigenstates, and the finite-volume form factor formulae (3.44), we arrive to wherex := x + iπ. As seen earlier in (3.33) the form factors above are only non vanishing if Note that the order of rapidities is chosen as in the definition (3.33) (this is just for convenience as for free bosons the order is irrelevant).
We will now take the expression (4.9) and evaluate its leading large-volume behaviour. There are two equivalent ways of doing this which we present below.

Computation by Exact Summation over Quantum Numbers
For large volume, the density factors in the denominator of (4.9) become large. However, if some rapidity of the intermediate states approaches the rapidity of the excited state, the kinematic poles of the form factors will give rise, due to (3.47), to positive powers of the volume. The powers in the numerator and denominator will combine to give an overall power of the volume L. In this section we will show that the largest such power is zero. Therefore, as L → ∞ the two-point function (4.9) tends to a volume-independent value. There are three different cases we should investigate for a given rapidity θ i or β i . Recall that from (3.33) each of the form factors above consists of a large sum of products over two-particle form factors.
The first case of interest occurs when we consider the contribution to (4.9) of those terms where the same rapidity θ i is paired up with the rapidity θ (in the Wick-contraction sense of (3.33)) in a two-particle form factor coming from each of the form factors in (4.9) . If θ i ∼ θ, then the form factor product above will be dominated by the contribution around the corresponding kinematic poles and we can write and, similarly . . . ,θ 1 , . . . ,θ i , . . .θ m + , . . .) , (4.12) wherex means that the variable x is no longer present in the form factor. Above we kept implicit the dependence of the form factors on sets of rapidities not involved in the contraction. The combinatorial factors N + p andÑ + p come from the many pairings of θ i with θ as per the permutation in (3.33).
The leading large-volume term from the summation over the quantum number J + i , pertaining to the rapidity θ i , is where, as before, r = L and we used the Bethe-Yang equations (3.40) and (3.41) to express the rapidites in terms of the associated quantum numbers. Here g n p (r) are the functions (3.48). Note that since the sum is over all integers, the value of the integer I has no effect on the outcome of the sum. In other words, the result is independent of the value of the rapidity θ. Similarly, for the case when some β i is paired up with θ in both the form factors we get As a consequence, if a rapidity is paired up with θ in both the form factors, the summation gives an (mL) 2 factor. The second case of interest occurs when none of the rapidities θ i , β i are paired up in any of the form factors with θ. In this case, the large volume limit is regular, there is no kinematic singularity playing a role, and we can replace the summation over quantum numbers by integration This operation generates additional factors of order mL for each integral. Finally, there is a third case which is a mixture of the previous two, namley when θ i or β i is paired up with θ in one of the form factors but with a different rapidity in the other. Due to the shifts in the Bethe-Yang equations (3.40), (3.41) and (3.42), the summation is not singular at any value of the volume, and it can be rewritten with principal value integral giving once more an mL factor. By successively using the expansion (4.11), (4.12) with the summations (4.13) and (4.14) we can calculate the overall leading large-volume contribution to the two-point function. Indeed, in Appendix B we show that this leading large-volume contribution is of order L 0 and is obtained exactly when N ± p =Ñ ± p with N ± p ≤ m ± , and N + p (N − p ) intermediate rapidities θ i (β i ) are paired up with θ in both form factors. Each pairing of the rapidities gives rise to a sum of the type (3.48) with the remaining, unpaired rapidities giving rise to form factors dependant on a smaller set of variables. Explicitly where q ± = m ± − N ± p is the number of remaining intermediate state rapidities after the contractions. The factorials in the denominator are just m ± !, that came from the complete set of state insertion. Out of m ± original intermediate rapidities, N ± p are paired up with the rapidity θ in the sense described earlier. All particular pairing choices are equivalent to each other under relabelling of the rapidities, as they are all integrated over, that is counted by the binomial factors. The N ± p ! combinatorial factors arise from the pairing of the chosen intermediate rapidities to θ in the form factors, as explained in (4.11) and (4.12). Once all possible contractions with a rapidity θ have been carried out, two form factors will still remain depending on q + + q − rapidities. In addition, we know from (3.33) that only form factors with q + = q − = q are non-vanishing. Simplifying we obtain p , the expression above exactly reproduces the finite-volume vacuum two-point function in the given sector, i.e. p,L 0|T p (0)T p ( )|0 p,L . As a consequence, our end result for the finite-volume two-point function can be expressed as In particular, for p = n, the factor reproduces the norm of the finite-volume state as expected, since g ±n (r) = 1 and n;L 0|T n (0)T n ( )|0 n;L = 1.

Computation by Contour Integration
Another way of calculating the leading large-volume term of the two-point function in a given sector (4.9) is to transform the summation over quantum numbers of the intermediate states into contour integrals. This approach not only leads to the same result (4.19) but seems more amenable to generalization to interacting theories, something we would like to attempt in future work. Consider generic sums of the form and where h ± are functions that are regular at the positions θ i , β i , respectively and C J ± i is a small contour encircling θ i , β i with positive orientation, and the denominators inside the integrals are the exponential form of the Bethe-Yang equations (3.41) and (3.42), that is zero at every solution of the equations. From now onwards we will omit the tilde on the integration variables.

20
Transforming every sum in (4.9) into a contour integral we obtain the expression Our next step is to combine the small contours around the Bethe-Yang solutions into a contour encircling the real axis for each variable. While doing so, the contour will cross the kinematic poles of the form factors, whenever θ i = θ or β i = θ for some i, and we need to account for the residues of these poles. It is easy to see from (3.28), that the contribution from residues at θ coming from a single kinematic singularity is of order L 0 in the volume and therefore they will be strongly suppressed by the power of L in the denominator of (4.22). However, if we consider terms where both form factors have a kinematic pole at the same location θ i = θ or β i = θ, then we have to calculate the residue of a second order pole, and this can change the order in the volume. Let us calculate this residue for a particular rapidity θ i Recall that here, as earlier hatted variables are variables shifted by iπ. From the kinematic residue equation (3.28) it follows that near the kinematic poles the integrand may be approximated as Evaluating the corresponding residue we obtain where the checked variables are absent. Simplifying, the final result is where we also used the Bethe-Yang equation (3.40), and the N + p ,Ñ + p combinatorial factors are the result of the pairing of θ i with the θs. It is important to note, that the result is proportional to the volume, and also to the function g n p (r) introduced in (3.48). An entirely similar computation, for a rapidity β i gives the result As a consequence of these residues, we get the leading large-volume contribution to the twopoint function, if we pick up the largest possible number of residues of second order poles which are enveloped as the contour is deformed. The maximum number of second order poles is min(N ± p ,Ñ ± p ), that implies, that m ± ≥ max(N ± p ,Ñ ± p ). These terms have an e mLR dependence on the volume with (4.28) As argued more generally in Appendix B (the formula above can be seen as an especialization of equation (B.4) in Appendix B) the leading contribution is obtained when N ± j =Ñ ± j , and in 22 that case R = 0. The leading large-volume term of the two-point function then becomes ×F p,n q + +q − (θ 1 , . . . , θ q + ; β 1 , . . . , β q − )F n−p,n q + +q − (β 1 , . . . ,β q − ;θ 1 , . . . ,θ q + ) , where C denotes the contour encircling the real axis, q ± = m ± − N ± p , and the combinatorial factors came from counting the various choices of intermediate rapidities giving rise to double pole residue integrals. Simplifying the combinatorial factors and noticing that q + = q − = q for the form factors above to be non-vanishing, we can easily factor out the vacuum two-point function from the expression above and we obtain once more the result (4.19).

Example: 2nd Rényi Entropy of a Single-Particle Excitation
Let us illustrate the general methods above with the simplest example: we compute the 2nd Rényi Entropy, i.e n = 2, of a single-particle excited state. From (4.1) we can easily write down the state where we used the fact that g 2 2 (r) = g 2 −2 (r) = 1 and g 2 1 (r) = g 2 −1 (r) = 1 − 2r. Therefore the difference of Rényi entropies is which agrees with the expression (2.1) for n = 2. This is also exactly the second Rényi entropy of the two qubit state (2.13).

Multi-Particle Excited States
In this section we adapt the techniques presented for the one-particle excited state case to more general states involving both distinct and equal rapidities. As we will see the essential ideas are the same but the state is more involved which makes the combinatorics of the problem more complicated.

Distinct Rapidities
Let us denote a general k-particle state (3.19) involving only distinct rapidity excitations as | 1, 1, . . . , 1 k L . It can be expressed similarly as (4.2) in the form where the C n ({N q,± }) coefficients are all identical for each value of q (the state is invariant under relabelling of the rapidities). For fixed q they are exactly the same as for the one-particle state. We have the following restrictions for the integers n p=1 =± N q, p = n , where To find the leading contribution in the volume to F p {N q,± p }, {Ñ q,± p } , we follow the same steps as in Section 4.1. As seen in Subsections 4.1.1 and 4.1.2, we need to focus on the contributions arising when some intermediate rapidity approaches one of the rapidities of the excited state in both of the form factors. In other words, we need to pair up the intermediate rapidities with the same rapidity of the excited state from the in-and out-states. This mechanism singles out the leading large-volume contribution as corresponding to N q,± p =Ñ q,± p for every q. Carrying out the calculation, the combinatorial factors simplify to yield the result As a consequence, in the infinite volume limit, the result for a state involving k distinct rapidities factorizes into k single-particle state contributions. That is This in turn leads to the relation which is a special case of the formula (2.8).

Coinciding rapidities
The simple result (4.41) no longer holds if all or some rapidities of the excited state coincide. Let us consider a k-particle excited state where all the rapidities coincide, and are denoted by θ.
In this case the norm of the k-particle state as written in (4.34) is k! n , thus the normalization needs to be appropriately modified. The properly normalized state can then be written as which looks very much like the one-particle state (4.2). This is not too surprising as both states depend on a single rapidity variable. The coefficients D k n ({N ± }) are related to the coefficients C n ({N ± }) of the previous subsections by This relation is of practical use when evaluating our formulae with the help of algebraic manipulation software. The two point function is then where F p is the same function as for the one-particle case (4.5), but now the integers N ± p obey the selection rule

The General Case
The techniques we have just presented for states of distinct and equal rapidities can be easily adapted to deal with more general states: states where some rapidities are equal and other distinct. As expected, the EE difference for a multi-particle mixed state is a sum over the EEs of simpler states associated with groups of coinciding rapidities. This result is expressed by the formula (2.8).
Regarding the results of this section overall, it is worth noting that we do not yet have closed formula for coefficients C n ({N ± }) and D k n ({N ± }) for general n, however it is straightforward to calculate them systematically on the computer and we have done this up to n = 6 for two coinciding rapidities and up to smaller values of n as the number of coinciding rapidities was increased to k = 6. Once the coefficients are known we can easily evaluate formula (4.46) for several values of k, and we observe that the results are always polynomials that have r ↔ 1 − r symmetry as expected. It was by working out such particular examples that we were eventually able to establish the general pattern (2.1)-(2.7).

Example: 2nd Rényi Entropy of a Two-Particle Excitation
In order to make the results above more concrete, we will now consider the EE of two-particle excited states both with distinct and with equal rapidities. The non-trivial part of the computation is in the characterization of the states, namely the computation of the coefficients C n ({N ± }) and D k n ({N ± }) as arising in the formulae (4.34) and (4.46). Once these are know the EEs can be systematically obtained for any state.
Let us consider a two-particle excited state with distinct rapidities which we represent as |1, 1 L . From the general expression (3.21) it is easy to see that The state can be fully characterized by the coefficients C 2 ({N q,± }) with q = 1, 2 and these give two copies of the coefficients (4.31) of the one-particle state (4.30). Substituting these values into the formula we obtain exactly the square of (4.32), that is Consider instead a two-particle excited state of equal rapidities. The state may be written as where the last line follows from noting once more that g 2 2 (r) = g 2 −2 (r) = 1 and g 2 1 (r) = g 2 −1 (r) = 1 − 2r. This then gives the expression

Numerical results: the harmonic chain
The formulae (2.1)-(2.7) are somewhat surprising for their simplicity and their qubit and semiclassical interpretations, especially as that they emerge from an exact, involved QFT computation. It is therefore important to convince ourselves that this is indeed the behaviour of entanglement that emerges when explicitly carrying out the scaling and thermodynamic limit of a discrete quantum mechanical system. In the free boson case the ideal model on which to test our formulae is the harmonic chain. The numerical method that we have employed is a wave functional method and we present the details in Appendix A. This is a method based on the exact inversion of a matrix, and gives machine-precision results for the EE of excited states in the harmonic chain in finite volume. Some of the results are presented in this section, see Figs. 4, 5. In all cases, there is excellent agreement between the numerical computation in the limit of large volume and region length L, and small lattice spacing ∆x (the large-volume scaling limit L, m −1 ∆x) and the analytical large volume results (2.1)-(2.7).
As was explained in [32], the results are in fact expected to be correct in a regime of parameters that goes beyond the universal scaling regime of QFT. The condition, expressed in full generality in [32], is that the minimum of the maximal De Broglie wavelength 2π/P i of all particles, and the correlation length ξ = 1/m, must be much smaller than the minimum of and L − . This include large momenta regions, beyond the low-energy QFT regime, and holds independently of the value of ∆x. Below we present some large-momenta results that confirm this.
The results presented in Fig. 4 are as follows. On the left panel a series of the Rényi entropies (2.1) is presented in the case of a single particle, k = 1. In the cases n = 2, 3, 4, 5, 6, 11, both the analytic (continuous curves) and numerical (dots, squares, triangles etc.) results are presented. All curves have a single maximum at r = 1 2 . The numerical results are in perfect agreement with the analytic results, with relative errors less than 10 −7 . Numerical results are obtained for mL = 5 and with the largest momentum allowed by the chosen lattice spacing (∆x = 0.01), which is in the middle of the Brillouin zone.  Figure 4: Comparison between analytic results (continuous curves) and numerical values (dots) of the Rényi entropies. Left: a single particle, Rényi entropies from n = 2 (red) to n = 11 (orange), with momentum P = 100π. Right: two particle states, n = 2, with distinct momenta given by P 1 ≈ 30, P 2 ≈ 45 (squares, red curve) and with equal momenta P 1 = P 2 ≈ 50 (dots, blue curve). Additional choices of the momenta are explored in Tables 1 and 2. The right panel in Fig. 4 shows the 2nd Rényi entropy for a two-particle excited state. The outer-most curve is twice the function (2.1), that is This is twice the second Rényi entropy of a single excitation. The squares exactly fitting this curve are the numerical values for volume L = 10, m = 1 and a particular choice of (relatively large) distinct momenta. It is interesting to investigate how the chosen values of the momenta affect the accuracy of the fit. Table 1 shows an additional example for distinct (small) momenta P 1 ≈ 0.6 and P 2 ≈ 2.

28
The inner-most curve (with the lowest maximum) is the function (2.4) with k = 2, that is This describes the entanglement of a two-particle excited state with particles of the same momentum. Numerical results are presented with L = 10, m = 1 and P 1 = P 2 = 50. Table 2 shows additional values for the same quantity and momenta P 1 = P 2 = 2 and P 1 = P 2 = 10. High precision is obtained even for relatively small momenta.  Table 1: The difference of 2nd Rényi entropies of a two-particle excited state with distinct momenta. The second row shows the exact values of the function (4.53). The third row shows the numerical values for the given momenta. The other parameters are m = 1, L = 10 and ∆x = 0.01. We see that agreement is not as good as for the data in Fig. 4 (top right), especially for small . This is due to momenta being too small. More precisely min(2π/P 1 , 2π/P 2 , ξ) = 1 which is larger than some of the values of considered, a regime in which we do not expect our formulae to hold. However, even for such small momenta the disagreement with (4.53) is at worse around 10%.  Table 2: The difference of 2nd Rényi entropies of two-particle excited states with equal momenta. The second row shows the exact values of the function (4.54). The third and fourth rows show numerical values for the given momenta. The other parameters are m = 1, L = 10 and ∆x = 0.01. For P 1 = P 2 = 2 agreement is poorer, especially for small due to the momenta being too small. More precisely min(2π/P 1 , ξ) = 1 which is larger than some of the values of considered, a regime where we do not expect our formulae to hold. However the disagreement with (4.54), even for such small momenta is relatively small. For P 1 = P 2 = 10 (as for 50, in Fig. 4) agreement is excellent for all values of r.
which is the second Rényi entropy of a three-particle excited state with equal momenta.
Finally, the right panel of Fig. 5 is the 2nd Rényi entropy of a four-particle excited state. Here four cases are shown: the outer-most curve is the case where all momenta are distinct, corresponding to the function 4∆S 1 2 (r); the curve with the second highest maximum is the case where particles are divided into two distinct-momentum groups of two equal-momentum particles, corresponding to the function 2∆S 2 2 (r); the curve with the third highest maximum is the case where three particles have equal momenta and the fourth particle has a different momentum, corresponding to the function ∆S 1 2 (r) + ∆S 3 2 (r); finally, the inner-most curve is the second Rényi entropy of a four-particle excited state with all rapidities equal. This is given by the function In all cases the volume is again mL = 10 and momenta are chosen high enough. As discussed earlier we observe that as the momentum increases, the agreement between numerics and analytical functions becomes better. At large volume, it is possible to reach very high precision with momenta that are large enough while still within the QFT regime, where the dispersion relation is relativistic. We also studied momenta beyond the QFT regime see Fig. 4 (left), towards the middle of the Brillouin zone (where the energy is maximal), P ≈ π/∆x = 100π. There, we not only observed near machine precision, but also, the condition of the volume L being much larger than the correlation length m −1 is no longer necessary: results keep near machine precision for any values of L, , ∆x with L, ∆x, even with m −1 L, (large correlation lengths). We do not currently have a derivation of this result. Intuitively, this indicates that when the wave function of the excited state presents a large number of oscillations within each subregion, then the entanglement behaves as that of the qubit system explained in Subsection 2.2: the large number of oscillations guarantees that the particle is "evenly distributed" within the subregions.
It is interesting to numerically study the finite-volume corrections to our formulae (2.1)-(2.7) and to compare the results to a QFT computation. We expect to investigate this problem in a future work. Some results were reported in the supplementary material of [32] which, for the harmonic chain, where compatible with integer power law corrections in L.

Excited State Entropies of the Massive Free Fermion
Technically speaking the computations presented in the previous few sections follow through with few but important changes for the free fermion theory. Interestingly however, the results (2.1)-(2.3) hold unchanged for free fermions. For free fermions states involving two identical creation operators have zero norm and therefore the more involved cases (2.4)-(2.7) do not arise in this case. Instead, for a state |1, 1, . . . L of k particles of distinct rapidities the results (2.1)-(2.3) hold as well upon multiplication by k (as for free bosons). As we will see later, in some respects, the free fermion theory is easier to treat by the techniques outlined in this paper simply because states have a simpler structure. In this section we review those technical features that are different for free fermions and present a detailed computation of the case of a one-particle excitation.

Doubling Trick and Replica Free Fermion Model
In this section we develop similar ideas as in Section 3.2. Consider two copies of a real (Majorana) fermion labeled by a and b. This gives us our "doubled theory" which we can now regard as a single complex (Dirac) fermion. The suitably normalized spinor components of this complex fermion are and, if ψ a,b ,ψ a,b are real, then . At the level of creation (annihilation) operators there exists a similar relation: (a (a) ) † (θ) = 1 and, considering now n-copies of such a real fermion in the replica theory, labelled by an index k we similarly have As noted in [24,48] where the ground state entanglement of free fermions was studied by employing similar ideas, it is possible to diagonalize the branch point twist field as well but it is important to make a distinction between n even and n odd. More precisely, the relation (3.9) generalizes to (5.5)

31
and similarly for the fields Ψ L,j . Note that, unlike for the free boson case, the matrix above is different depending on whether n is even or odd, a feature that has been discussed in [24,48]. The eigenvalues of this matrix are λ p = e 2πip n for p = − n−1 2 , · · · , n−1 2 , that is the nth roots of unity for n odd the nth roots of −1 for n even. The cyclic permutation action is diagonalized by the fieldsΨ 6) and the creation operators satisfy the relations and , a j 2 (β)} = 0 for all j 1 , j 2 = 1, . . . , n. The relation can also be inverted to and For free fermions, the U (1) fields associated to these generators have been also studied (see e.g. [54]) and it is known that they have scaling dimensions note that for the massless Dirac fermion c = 1. The form factors of these U (1) fields are also discussed in [54] and they are very similar to those found for free bosons. The two particle form factors have the same structure: and satisfy The two last form factors are vanishing for symmetry reasons. The form factor programme for quasi-local fields [44,45,51] tells us that these form factors are solutions to a set of three equations. First, Watson's equations where γ ± p = e ± 2πip n are the factors of local commutativity associated to the fermions ±. Finally, the kinematic residue equation tells us that where Since the theory is free, higher particle form factors can be obtained by simply employing Wick's theorem. For the Dirac fermion they have the structure where once again f n p (θ) is the normalized two-particle form factor and σ is an element of the permutation group S m of m symbols and sign(σ) is the sign of the permutation σ.
An important property of the form factor (3.32) is its leading behaviour near the kinematic singularity. Consider the form factor f n p (θ 1 −β 1 +iπ) and suppose that the rapidites are quantized through Bethe-Yang equations of the form Then the leading contribution for θ 1 ≈ β 1 can be expressed as .

(5.20)
Note that for free fermions it is common to distinguish between periodic and anti-periodic boundary conditions for the Bethe wave function. These lead to quantization conditions (5.19) which either require I, J ∈ Z or I, J ∈ Z + 1 2 . In our particular computation this choice makes no difference to the final result as we will obtain expressions such as (3.47) which only depend on quantum number differences. In addition, the U (1) twist fields do not change the Z 2 sector (contrary to σ field in the Ising model). For this reason and without loss of generality we consider the quantization condition (5.19) only.

EE of Single-Particle Excitations
Given the relations (5.4) we can represent a replica one-particle excited state in a free fermion theory as In the basis of the generators a j (θ) =ã + j (θ) and b j (θ) =ã − j (θ) this state becomes where ω = e − 2πi n . For instance, for n = 2 it is easy to show that the state takes simply the form where we introduced the notation S(n) to denote the sum over creation operators. For n = 3 we have instead We note that the main difference between n even and n odd is that for n even there is no "trivial" sector with index 0.
These particular examples illustrate the general structure of the states. For both n even and odd, they can be constructed recursively starting from the two simple examples just discussed. The state for a given n can be obtained from the state for n − 2 as follows where α is a phase which can be determined for every n. Its determination is actually a rather non-trivial problem but, as the states (5.22) have norm one by construction, we know it must be a real number. Its value has no effect on subsequent computations as only the norm of e iα will be involved.

Leading Contribution to the Rényi Entropy
The leading contribution to the Rényi entropy can be easily evaluated as all correlators emerging from the states above have a very simple factorized structure. For instance, for n = 2 the leading contribution will come from the matrix elements whereas for n = 3 we have instead By leading contribution we mean here that non-diagonal matrix elements (involving different states on the right and left) have been neglected as the arguments presented in Appendix B show that these, even when non-vanishing, will give sub-leading contributions in the volume. As can be seen from these examples, the building blocks of the correlation function are generally matrix elements of the form Matrix elements of the form p;L 0|a p (θ)b p (θ)T p (0)T p ( )b † p (θ)a † p (θ)|0 p;L have leading large L behaviours which are identical to those of so they involve once more matrix elements of the type (5.28).
The leading large volume contribution to such correlators can be evaluated along the same lines presented for the free boson theory. For instance, let us take one particular example: Recall that the sets {J ± i } are integers corresponding to the quantization of rapidities {θ i }. In finite (large) volume we can write as usual From here, once more the leading contribution will come from terms in the form factor squared such that the rapidity θ + iπ is "contracted" with the same rapidity θ 1 , . . . , θ s in both form factors. Such terms (there are s + 1 such choices) contribute a two-particle form factor squared times the vacuum two-point function, which once more factors out. This gives both for n even and odd. The fact that this gives the same entanglement entropy as the free boson is mathematically very interesting in the sense that in this case it comes from a single product of functions g n p (r) whereas for the free boson it was the result of adding together a constant plus various powers and products of these same functions. It is also not difficult to see that this same structure is recovered when considering multi-particle states of distinct rapidities.

Conclusions and Outlook
In this paper we have studied the n th Rényi entropy increment of a single-interval in one space dimension and its limits n → 1 (von Neumann entropy) and n → ∞ (single-copy entropy). Our work has focussed on a very particular class of QFTs and excited states |Ψ : the former are massive free QFTs in 1+1 dimensions and the latter are zero-density states, populated by finite numbers of particles. We have considered the particular limit , L → ∞ with r := L finite. It is well-known that the EE of finite-density excited states in gapped systems satisfies a volume law [3]. In the current work we have shown that for zero-density excited states in infinite volume the EE of one interval saturates to a value which, upon subtracting the ground state contribution, is a simple non-negative function of the ratio r. More precisely, for 0 < r < 1 the excited state provides a net positive additive contribution to the saturation value of the entanglement entropy. For any zero-density states and entropies, this contribution is maximal for r = 1/2. Moreover, for excited states consisting of k excitations of distinct rapidities, the maximum is k log 2, that is, every excitation "adds" exactly log 2 to the entanglement entropy of the ground state. The simple form of our results makes them amenable to a qubit interpretation in which each k-particle excited state is associated with an entangled qubit state with coefficients that are probabilities of finding q excitations in region A and k − q in region B (see figure above) for q = 0, . . . , k.
Some of our results have previously appeared in the literature (see e.g. [31,36]) and have been described as semi-classical limits of the EE. Our work, together with the companion paper [32], strongly suggest that the results (2.1)-(2.7) apply much more generally, in fact, to any situations were one can reasonably speak of localized quantum excitations. It is also worth emphasizing that our derivation is the only analytic explicit computation we know of, leading to the formulae (2.1)-(2.7).
The domain of applicability of (2.1)-(2.7) may be formally characterized by the condition: where P is the largest momentum of any of the excitations in the state |Ψ L and 2π P can be interpreted as the De Broglie wave length associated to that particular excitation, whereas ξ = m −1 is the system's correlation length. Interestingly, this condition implies that we may have a situation where the correlation length of the system is very large and P is also very large and yet still find the same results. Indeed, we provided numerical evidence of this in Fig. 4 and also for higher dimensions in [32].
This work offers ample scope for generalization and extension. It is reasonable to expect that the same results should also hold for interacting integrable models of QFT. There are three main reasons for this expectation. First, technically, the key mathematical property leading to formulae (2.1)-(2.7) from the form factor calculations of Sections 4 and 5, is the kinematic pole structure of the form factors. However, this form is rather universal and not exclusive to free theories. Second, from [32] and [31] there is evidence that the same results hold in gapped interacting quantum spin chains whose thermodynamic limit should be described by integrable QFT. Finally, the qubit interpretation is quite universal (at least, as long as there is no particle production) so that we see no reason why results should change in more general theories. However, it would be nice to have a rigorous derivation of this result and we hope to provide this in a future work.
Another interesting problem is the investigation of finite volume corrections to (2.1)-(2.7). These can be computed both from the form factor expansion and numerically form the wave functional method presented earlier. Some numerical analysis of such corrections was presented in the supplementary material of [32] but a more detailed analysis of how the corrections depend on the energy of excitations, the value of r and the replica number n would be very interesting. According to our general arguments in Appendix B we expect the next-to-leading order correction the entropy increment to be of order 1/L in the volume so that, for a generic state we should have where f (n, r, {θ i }) is some function of n, the region size, and the rapidities of the excitations, which can be computed from a form factor expansion. It would also be interesting to extend the analysis to higher dimensions. For critical systems it has been shown that the EE contains information about the shape of the regions (e.g. the number of vertices) [55,56,57,58,59,60,61] and we would like to investigate whether or not such information can also be red off from the finite volume corrections. At present we cannot compute these exactly in higher dimensional gapped QFT, but for free theories, we can use the wave functional method to investigate the problem numerically as in [32].
To conclude, our results provide further evidence that measures of entanglement encode universal information about quantum models, be it their universality class [6,4,62], operator content [8,9,28,29], particle spectrum [24,38,63] or, as in this case, the number and nature of their excitations above the ground state. These results come at an exciting time in the understanding of entanglement measures as experimental results for particular Rényi entropies have recently become available [64,65]. It would be extremely interesting to connect our results to experiments and to understand their implications in the wider quantum information context.

A Wave Functional Method
In this appendix, we describe how to evaluate numerically traces of the n th powers of reduced density matrices for few-particle excited states in the quantum free boson model. We use the wave functional method, which is based on completely different principles than methods using form factors explained in the main text, thus offering an independent verification of our results. After discretizing the model to a finite chain of size N , the method reduces the problem to the inversion of a single nN by nN matrix, which can be performed numerically.
Consider the real free boson, with hamiltonian where Φ(x) and Π(x) are hermitian canonically conjugate fields, [Φ(x), Π(x )] = iδ(x − x ). The wave functional of the ground state can be obtained by methods similar to those used for the ordinary harmonic oscillator in quantum mechanics. The annihilation and creation operators are A p and A † p for p ∈ (2π/L)Z with We use the representation of wave functionals Ψ[ϕ] = ϕ|Ψ , with wave functionals taking as arguments fields ϕ : [0, L] → R. In this representation, The vacuum satisfies A p Ψ vac = 0, which gives where N is a normalization factor. Excited states are obtained by acting with the creation operator, giving for instance We denote the reduced density matrix of the vacuum state as ρ B|vac , and that of the excited state as ρ B|{p j } . We are interested in the ratio Tr(ρ n B|{p j } ) Tr(ρ n B|vac ) . (A.10) By using the Gaussian form of the vacuum wave functional (A.4) and the fact that excited states are obtained by multiplying by polynomial functionals of the fields, as in (A.7), we see that (A.10) is the average in a Gaussian measure over the fields ϕ j , of a product of the monomials α p j . In order to evaluate numerically this average, we discretize space. For this purpose, we choose ∆x = L/N , (A. 11) for some N ∈ N, restrict space and momentum variables to We also change the action to its discrete version, which gives rise to the following change in the equations of motion, This induces a change in the dispersion relation, the new energy function being Putting these ingredients together, some calculations show that the final result can be expressed as follows. Define The average · · · is over the Gaussian measure given by the discretized vacuum wave functional, O[ϕ 1 , . . . , ϕ n ] = Dϕ 1 · · · Dϕ n O[ϕ 1 , . . . , ϕ n ] exp − 1 2 M Dϕ 1 · · · Dϕ n exp − The matrix M is an nN by nN matrix, and the inverse matrix M −1 can easily be evaluated numerically. Schematically, the matrix M has the following block structure where the matrices K Q 1 Q 2 have entries (K Q 1 Q 2 ) ij := (L/N ) 2 K(x i − x j ) with x i ∈ Q 1 and x j ∈ Q 2 .

B Selection Rules for Leading Terms in the Form Factor Expansion
In this appendix, we identify the terms in the form factor expansion that contribute in the limit of large system size L. We show that these terms contribute to order L 0 (that is, are finite and nonzero), and that all other terms contribute to orders L −1 or less (that is, vanish as L → ∞).
The leading terms are analyzed in the main text, and give rise to the main results of this paper. For simplicity, we will consider the case where the excited state depends on a single rapidity value: either it is a single particle state, or a many-particle state, where all particles have the same rapidity θ (this is of course only possible in the free boson case). The general case, involving many distinct rapidities, can be understood along similar lines.
Consider a generic term in the form factor expansion (4.9). A generic term is characterized by a number N of particles in the (bra) state on the left, a numberÑ of particles in the (ket) state on the right, the set B = We first establish the leading power of L corresponding to (B.1). Due to (3.44) a finite-volume form factor contributes a factor 1/ √ L for each rapidity: Each particle in the intermediate state that is not contracted with a particle in left or right states (and is, each particle with label in B \ (A ∪Ã)) contributes a factor of L, as for such particles, the sum is evaluated by transforming it into an integral, θ ∼ L dθ: Finally, each element in A contributes a factor L, and each element inÃ also contributes a factor of L. This accounts for two situations. First, a particle may be contracted with one in the left (or right) state but not with any particle in the right (or left) state, j ∈ A and j =∈Ã (or vice versa). In this case, the contraction gives rise to a single pole. The sum over θ j can then be transformed into a converging, principal-value integral L P dθ j , giving a factor of L. Second, a particle may be contracted both with one in the state on the left, and one in the state on the right, j ∈ A and j ∈Ã. In this case, the leading contribution is obtained by "zooming in" onto the second-order pole that develops, and summing the resulting secondorder pole contribution without transforming the sum into an integral. This sum is convergent, and results in a factor L 2 using the fact that momenta are proportional to 1/L. For instance θ j 1/(θ j − θ) 2 ∼ I j ∈Z L 2 /(I j − I − q) 2 for some I ∈ Z and q ∈ (0, 1). The factor of L 2 indicates that we must count a factor of L for the particle both as an element of A and as an element ofÃ. Thus, we have L |A|+|Ã| .
In order to find the leading behaviour, we must therefore maximize Thus R will be maximized whenever the cardinality of A ∩Ã is maximized. This occurs when either A ⊆Ã orÃ ⊆ A, giving Finally, this is maximized by taking N =Ñ . In this case, we have |A| = |Ã| and thus A =Ã, and we find R = 0. This shows the claim at the beginning of this Appendix. Moreover, the argument can be easily generalized to states consisting of various particle types.
C The Functions g n p (r) Throughout this paper we have used the relations The fact that the sum above is a simple polynomial in r can be of course checked numerically. It can also be shown analytically, for instance, by showing that the second derivative with respect to r is zero. We compute