{'fraction_non_alphanumeric': 0.03867188300116753, 'fraction_numerical': 0.011921099526842957, 'mean_word_length': 4.336904241364232, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Symmetry Resolved Entanglement of Excited States in Quantum Field Theory II: Numerics, Interacting Theories and Higher Dimensions

In a recent paper we studied the entanglement content of zero-density excited states in complex free quantum field theories, focusing on the symmetry resolved entanglement entropy (SREE). By zero-density states we mean states consisting of a fixed, finite number of excitations above the ground state in an infinite-volume system. The SREE is defined for theories that possess an internal symmetry and provides a measure of the contribution to the total entanglement of each symmetry sector. In our work, we showed that the ratio of Fourier-transforms of the SREEs (i.e. the ratio of charged moments) takes a very simple and universal form for these states, which depends only on the number, statistics and symmetry charge of the excitations as well as the relative size of the entanglement region with respect to the whole system's size. In this paper we provide numerical evidence for our formulae by computing functions of the charged moments in two free lattice theories: a 1D Fermi gas and a complex harmonic chain. We also extend our results in two directions: by showing that they apply also to excited states of interacting theories (i.e. magnon states) and by developing a higher dimensional generalisation of the branch point twist field picture, leading to results in (interacting) higher-dimensional models.


Introduction
The symmetry resolved entanglement entropy (SREE) is a relatively new measure of entanglement for many-body quantum systems. The concept was introduced in [1][2][3] although the name SREE first appeared in [4], where it was defined as a new measure of entanglement associated to theories that posses an internal discrete or continuous symmetry. The role of symmetries in the structure of the entanglement entropy and the contributions of symmetry sectors to the total entropy were also studied in [5], simultaneously and independently of [4]. The approach of [4], which we adopt here, falls within the framework of branch point twist fields [6][7][8][9][10] and their composite generalisations [4,[11][12][13]. These fields have only been formally defined for 1+1D theories and have proven particularly useful in the study of entanglement measures for conformal field theories (CFTs) [8,[14][15][16] and integrable quantum field theories [9,[17][18][19][20]. In this context, typical examples of theories that posses internal symmetries are the Ising field theory, with (discrete) Z 2 -symmetry, and the sine-Gordon model which has (continuous) U (1)-symmetry.
Let us recall the general definition of the SREEs. Consider a 1+1D quantum field theory and a bipartition of space into two complementary regions A andĀ so that the Hilbert space of the theory H also decomposes into a direct product H A ⊗ HĀ. Assume that the theory is in a pure state |Ψ . We define the reduced density matrix of subsystem A as and the standard von Neumann and nth Rényi entropy of subsystem A are defined as S = −Tr A (ρ A log ρ A ) and S n = log(Tr A ρ n A ) 1 − n with S = lim n→1 S n , where Tr A ρ n A := Z n /Z n 1 is the normalized partition function of a theory constructed from n non-interacting copies or replicas of the original model.
Assume now that the theory has an internal symmetry, with symmetry operator Q whose projection onto subsystem A we call Q A . By construction, we have that [Q A , ρ A ] = 0. Let q be the eigenvalue of operator Q A in a particular symmetry sector. Then Z n (q) = Tr A (ρ n A P(q)) , with P(q) the projector on the symmetry sector of charge q, is the symmetry resolved partition function. In terms of this object, the SREEs can be written as S n (q) = 1 1 − n log Z n (q) Z n 1 (q) and S(q) = lim n→1 S n (q) .
In particular, for n = 1 the partition function (3) represents the probability of obtaining the value q upon measuring the charge.
In part I of this series [61], we presented a study of these SREEs by focusing on complex free theories and qubit states. We considered a class of excited states |Ψ consisting of a finite number of excitations and a particular scaling limit, namely one in which parts A andĀ are infinitely large, but the ratio of their lengths by the total length of the system is r and 1 − r respectively, with r ∈ [0, 1]. Since the ratio of particle number to volume is zero, we call these states zero-density excited states.
As discussed in [4] the partition function (3) can best be obtained in terms of its Fourier modes, the charged moments Z n (α) = Tr A (ρ n A e 2πiαQ A ) .
If the symmetry is U (1) the following relationship holds which we employed in [61] in order to compute the SREEs. However, similar formulae hold for other symmetries and, in particular, the integral is replaced by a sum if the symmetry is discrete (see e.g. [39]).
Let us now summarise the main results of [61]. We computed the ratio of charge moments between the excited and ground states. As anticipated, the results are functions of the state |Ψ , the parameter α and the ratio r. Formally, we define the ratios generalising the notation in the definition (5) to include the dependence in r.
We found that for a state consisting of a single particle excitation of charge = ±1 above the ground state |Ψ = |1 , the ratio is M 1 n (r; α) = e 2πi α r n + (1 − r) n , and for a state of two excitations with opposite charges ± we have M 1 1 − n (r; α) = (e 2πi α r n + (1 − r) n )(e −2πi α r n + (1 − r) n ) .
For a state of k identical excitations, that is, with the same charges and momenta, we obtain where f k j (r) := k C j r j (1 − r) k−j and k C j = k! j!(k−j)! is the binomial coefficient. This formula acts as the fundamental building block for all other cases. A generic state comprising s groups of k i i identical particles of charge i will have In the present follow up paper we test these results numerically as well as extend them to interacting and higher dimensional theories.
Note that the possibility of having identical excitations is excluded for fermionic theories. For α = 0 these formulae reduce to those found in [62,63], later generalised to entanglement measures of multiple disconnected regions [64] and to higher dimensions for free bosons in [65]. These results in turn have been extended in a series of works [66][67][68][69][70][71] to deal with finite volume corrections and non-localised excitations. More recently, some of the α = 0 results were recovered as a semiclassical limit in the presence of an interaction potential [72]. This semiclassical picture had already been invoked much earlier, see for instance [73]. However, it is worth emphasising that our formulae are not merely semiclassical limits but hold for genuine quantum theories. The quantum nature of the model is encapsulated by the symmetry resolved entanglement entropy of the ground state (and its associated moments), which are highly non-trivial for quantum models. In other words, it is only the ratios (7) that are simple, not the individual charged moments. This paper is organized as follows: In Section 2 we present numerical results for two (free) discrete systems: a 1D lattice Fermi gas and a 1D (complex) harmonic chain. We find that in both cases the formulae above are reproduced with great precision, even if the scaling limit of the Fermi theory is a massless free fermion whereas for the harmonic chain it is a complex massive boson. For the complex free boson, we show how the SREE can be computed by employing a wave-functional method partly presented in this section with technical details in Appendix A. In Section 3 we show that our formulae also hold for some magnon states in interacting theories, including those of two interacting magnons. In Section 4 we present an extension of our results to higher-dimensional quantum field theories, which draws on the connection between entanglement measures and twist fields and is close in spirit to the wave-functional approach employed for the complex free boson. We conclude in Section 5.

Numerical Results
In this section we present numerical results for two very different discrete models. First we consider a 1D lattice Fermi gas, which has critical features but also possesses highly excited states whose entanglement is well described by our formulae, and then we look at the harmonic chain, whose scaling limit is a massive free boson. Whereas for the first model we can only consider distinct excitations, for the second we consider also states of identical excitations. The good agreement found confirms the more general picture put forward in [62] that these kinds of formulae hold under the broad assumption of localised excitations. These are present both in gapped systems due to finite mass scale/correlation length, and in critical systems, when the De Broglie wave length of the excitations (which is inversely proportional to their momentum) is sufficiently small compared to subsystem size.

1D Lattice Fermi Gas
In this section we analyse a particle-hole excited state of a 1D lattice Fermi gas, comparing our analytical predictions with the numerical data. Even though the model is critical, it was realised in [68] that certain highly energetic quasiparticle excitations still have a universal entanglement content. More precisely, if one assumes that a set of quasiparticles with small enough De Broglie wavelengths (compared to the typical geometric lengths) is present and their momenta are sufficiently separated, then the quasiparticles will be essentially uncorrelated with each other and with respect to zero-point fluctuations. We refer the interested reader to [66][67][68][69][70][71] for further details about the universal entanglement content of quasiparticles in critical systems.
Our goal here is to briefly review the numerical techniques involved in the characterisation of fermionic Gaussian states [74] and their application to the computation of symmetry resolved measures. We start by considering the Hamiltonian of free spinless fermions on a circle of length where µ is the chemical potential and {f j } j=1,...,L , {f † j } j=1,...,L are the ladder operators obeying the standard anticommutation relations When |µ| < 1 the theory is gapless, and the ground state is a Fermi sea with Fermi momentum k F = arccos(µ). The two-point function evaluated in the ground-state at Fermi momentum k F takes the following form Here, we analyse the quasiparticle excited state described by the following two-point function It corresponds to the insertion of a fermion of momentum k = k F + π 4 above the ground state and the removal of another fermion (or equivalently, the insertion of a hole) at k = k F − π 4 . The choice of the momentum shift |k − k F | = π 4 is not important in the continuum limit, where the only necessary condition is that |k − k F | remains finite when L → ∞ 1 . We now have to specify the symmetry of the model. The Hamiltonian (12) is invariant under an internal U (1) symmetry associated to the number of fermions generated by and it clearly satisfies the locality condition Q = Q A + QĀ, with As subsystem A we consider the segment of length , that is the sites j = 1, . . . , and investigate its entanglement properties with the complementary regionĀ containing sites j = + 1, . . . , L.
We denote by C A 0 and C A the × matrices resulting from projection of the matrices C 0 and C (defined by Eqs. (14) and (15) respectively) onto subsystem A, keeping only j = 1, . . . , as spacial indices. Following [27] we express the charged moments of the particle-hole state and the ground state by means of the determinants with ρ A and ρ A,0 the respective reduced density matrices. According to our analytical predictions, we expect that the ratio of the charged moments takes the following universal form that is, the expression for two distinct excitations with charges ±1, which contribute to the charged moment with an Aharonov-Bohm phase e ±2πiα . We write to indicate that equality is only expected in the scaling limit of the lattice model.
To test the validity of Eq. (20) we consider two entanglement measures, namely the excess of (total) Rényi entropy and the so-called (following the terminology of [27]) "excess of variance". The excess of entropy is recovered from our formulae for α = 0 and for two distinct excitations takes the simple form We define the variance 2 associated to ρ A as Similarly, we denote by ∆Q 2 A n,0 the variance of the ground state ρ A,0 . From (20) it then follows that the excess of variance is given by A way to physically interpret the result of Eq. (23) is to regard this excess of variance as twice the contribution associated to a single quasiparticle, since particles and antiparticles contribute in the same way. The latter is just the variance of a Bernoulli random variable with success probability given by namely the probability one associates to the presence of a quasiparticle in A computed with the density matrix ρ n A . Since the variance of a Bernoulli variable with probability p is just p(1 − p), we get Eq. (23).
In Fig. 2.1 we report the numerical values of ∆S n and ∆Q 2 A n − ∆Q 2 A n,0 , computed from (18) and (19) using exact diagonalisation of the correlation matrices C A , C A,0 , and our analytical predictions. We keep L fixed, analysing different values of r = /L. Our choice is motivated by the expectation that these plots should be "universal"at large L, meaning different data obtained with different L should collapse to the same universal prediction (independent of lattice details as k F ) when L → ∞. As we see from the plots in Fig. , the match between numerics and analytics is really good.

Complex Free Bosons and the Wave-Functional Method
In this section we consider a complex massive free boson. Unlike the 1D Fermi gas, this model and its lattice version allow us to test formulae for states containing two or more identical excitations. Furthermore, thanks to the computation technique that we outline below, we are able to directly access the ratios M Ψ n (r; α) defined in (7) for different values of α. Our numerical computation is based on the wave-functional method introduced in [63] (see Appendix A of that paper). Here we need to extend the technique to a complex theory and to the symmetry resolved moments. These extensions are not entirely trivial and for that reason we review the wave-functional method in detail.
Let us consider a 1D complex massive boson on the line [0, L] with Hamiltonian: where Alternatively, we can introduce a pair of real bosons Φ 1 , Φ 2 , and express Φ and Π as so that the Hamiltonian becomes that of two real bosons. The only non-vanishing equal-time commutators are: Since space is compact, there are discrete energy levels with dispersion relation: and the Hamiltonian is diagonalised via the introduction of two sets of creation/annihilation operators. The annihilation operators are: and as usual A † p (B † p resp.) creates a positively (negatively) U (1)-charged particle with momentum p from the vacuum. These operators satisfy: Φ(x), Π(x) then admit the usual Fourier decomposition Finally, the charge operator corresponding to the U (1) symmetry of the theory is: where the normal ordering means that we subtract every non-operatorial term resulting from the commutation relations (i.e. terms involving δ-functions).

The Wave-Functional Method
A very useful way to represent the states in the theory is provided by the wave-functional formalism. In this approach, we associate to every state a functional Ψ acting on the space of classical field configurations φ , φ † : [0, L] → C and formally defined by: The action of the operators Φ, Π on the (rather abstract) state |φ, φ † mimics the action of the position and momentum operators on the state |x in non-relativistic quantum mechanics, so that in the wave-functional representation: and analogously Notice that, for consistency, when iΠ, iΠ † act on a ket |φ, φ † there is a minus sign in front of the functional derivative. The functional of the vacuum state is defined by: and the only solution to these functional differential equations up to normalisation is the Gaussian functional: Notice that K(x) is a real, even function of x. The functionals of the positively and negatively charged one-particle states are obtained through the action of A † p and B † p : and the functional for a state with k + positive excitations and k − negative excitations (all with different momenta) is: A correct choice of the normalisation in (40) ensures that the functional above has unit norm with respect to the bra-ket product. However, the normalisation must be modified when some of the particles' momenta are equal, according to (32). If there are k + i (k − i ) positively (negatively) charged particles with momentum p i (q i ), for i = 1, . . . , m + (m − ) we define: The action of the charge operator on a wave-functional immediately follows from that of the fields and their conjugates: and in particular one finds that the vacuum and the Fock state functionals are charge eigenstates: However, because |φ, φ † is not associated to any charged state in the Fock space, it is not a charge eigenstate. The exponential of the charge operator acts on |φ, φ † by introducing phases (notice the minus sign in front of the integral): (47) This can be argued, for instance, by looking at the corresponding infinitesimal translation: Employing these results it is possible to show (see Appendix A for the derivation) that where G α is a known gaussian functional of the fields φ(x), φ † (x) given in (128). Results for the harmonic chain can then be obtained by discretisation, as we see in the next subsection.

The Harmonic Chain
Since the Hamiltonian (25) reduces to the sum of two Hamiltonians for the real bosons Φ 1 , Φ 2 with prefactors 1 2 , the discretisation proceeds exactly as for the real boson [63]. We divide the interval [0, L] in N parts introducing a spacing: and we define x = L Nx ,x ∈ {0, 1, . . . , N − 1}, so that we can replace every integral with a sum: If we discretise the Laplace operator as: and impose periodic boundary conditions Φ(0) = Φ(L), Φ † (0) = Φ † (L), the Hamiltonian (25) reduces (upon integration by parts) to two independent harmonic chains for real fields. The set of momenta is now restricted to the first Brillouin zone, p = 2π Lp ,p ∈ {0, 1, . . . , N − 1}, and the dispersion relation becomes: from which the relativistic relation E 2 p = m 2 + p 2 is obtained when pL 2N 1. Notice that since we restrict the set of momenta, the function K(x) defined in (40) becomes a finite sum: thus it is no longer an even function of x, though the property K * (x) = K(−x) still holds.
For the sake of simplicity we will take Φ and Φ † to be the fundamental degrees of freedom in the following, while keeping in mind that the real degrees of freedom are recovered using (27). In the formula (131) the functions U ± i , V ± i are modified by simply replacing the integrals with sums following the prescription (51). On the other hand, discretisation of the measure G α leads to a finite-dimensional (nN ) × (nN ) matrix G which couples the fields φ † i (x) and φ j (y): where we explicitly wrote the complex conjugate in (128) before discretising the integrals. Wick's theorem ensures that the Gaussian average in (131) can be computed from the contractions of pairs of fields, which are in turn obtained via the inversion of the matrix G: The matrix G has a block structure, consisting of n 2 blocks G i·,j· , each of which is an N × N matrix. From the above expression we see that the only non-vanishing blocks are either in the diagonal G i·,i· or just off the diagonal, G i·,(i±1)· . Each block G admits a sub-block structure in terms of the matrices K Q 1 Q 2 , Q 1 , Q 2 ∈ {A,Ā}, whose elements are: Notice that K AA is a square matrix with dimensions L N × L N , K AĀ has dimensions L N × L− L N and so on. From (55) we obtain the following basic structures: Note that the block structure is different from that in [63] because the roles of regions A and A are now exchanged. Although this exchange has no effect on the form of the entanglement entropy, for the SREE it makes a difference as the symmetry between A,Ā is broken when we choose to place the charge in subsystem A. In terms of the N × N blocks above, we can schematically write the matrix G as follows: where we introduced λ α = e 2πiα .

Numerical Results
In Appendix A we have explicitly derived the ratio of charged moments for excited states. We thus have all the ingredients needed to obtain numerical results. Let us take a bi-partition where A is a segment made of N A ≤ N consecutive sites, with N A /N = r and analyse the behaviour of M Ψ n (r; α). In Fig. 2 we compare results for two kinds of two-particle excited states: those of particles with identical charges and either distinct or equal momenta p 1 and p 2 . Our analytical predictions for M n (r; α) are In our numerics we have chosen L = N = 30, so that the lattice spacing L/N = 1. We also fix the mass scale to m = 0.1, which corresponds to a typical correlation length of ξ = m −1 = 10 sites. Finally we choose either p 1 = p 2 = π or p 1 = π and p 2 = 2π 5 , both in units of the lattice spacing.
In all our figures we chose non-negative values of α. Given the formulae above, taking α < 0 is equivalent to complex conjugation with α positive, so the figures for negative α are identical except for a change of sign in the imaginary part of all functions. We have also considered the value α = 0 (in green) which is the limit where there is no flux. As expected, in this case our formulae recover those for the excess Renyi Entropies in [62,63], which are symmetric in r and The numerics for the top row figures employ momenta p 1 = p 2 = p 3 = π, for the central row p 1 = p 2 = π, p 3 = π/3, whereas for the bottom row we took p 1 = π, p 2 = π/3, p 3 = π/5. have vanishing imaginary part. Despite the fact that the correlation length is not particularly small with respect to the system's size L (ξ 0.33L), we took highly energetic states (momenta being fixed in the large-volume limit) and we thus expect the validity of our predictions.
In both Fig. 2 and 3, we plot the numerical data (triangles) against analytical predictions ( (58) and (60)) as functions of r fixing n = 2 for several values of α between 0 and 1 2 , which correspond to flux ±1, respectively. At these two points, the ratio becomes purely real. The figures show excellent agreement between numerical data and analytical predictions.

Magnon States
The numerical results of the previous section provide convincing evidence that the results of [61] are correct. We now consider how our results might be applicable in a broader context. A natural starting point are magnon states. Such states describe the eigenstates of a variety of spin chain Hamiltonians, with or without interactions. They admit a simple explicit form in the spin basis so that entanglement computations are easy to perform. We also know from [63,71] that the total entanglement entropy of magnon states is described by our formulae with α = 0. As we see below, even in the presence of non-trivial scattering, the agreement extends to α = 0.
The main idea behind this construction is somewhat similar in spirit to the qubit picture [61,63], namely, that the entanglement content of quasiparticles can be easily understood if one factors out the zero-point fluctuations. In other words, instead of considering the full quantum theory where the quasiparticles are constructed on top of a nontrivial ground state, which in general has its own entanglement content, we consider a simpler theory in which particles are constructed above a trivial ground state. It turns out that the entanglement of this simpler model keeps track of the exact entanglement of the quasiparticle and discards explicitly the entanglement of the true ground state. Our magnon states belong to a Fock space generated by multiparticle configurations. In this case, the symmetry we construct will be "internal", meaning that it acts just as a phase on each multiparticle state; this phase is directly related to the quantum numbers of each particle.

One-Magnon States
We firstly focus on a single magnon state on the lattice, belonging to the one-particle sector of a quantum spin-1 2 chain of length L. A one-magnon state can be written as where |j is the state of a localised magnon of momentum p in position j. If one imposes boundary conditions on the chain, the momentum p is quantised as follows We introduce the action of the symmetry operator e 2πiαQ , where Q is associated with an internal symmetry. For our purposes we just need to specify its action on the vacuum state |0 and on the one-particle sector. In addition we assume that the magnon is charged with respect to Q, and it has charge +1.
We are interested in the entanglement between spins 1, 2, · · · , in region A and the rest of the system. Associated to this region there is a restricted symmetry generator e 2πiαQ A , which acts as where δ j∈A gives 1 if j is in A and 0 otherwise. The reduced density matrix of the region A is These two terms appearing in the formula above are interpreted as the contributions associated to the presence/absence of magnon in subsystem A, respectively. It is easy to show that a relation which shows explicitly the presence of the symmetry for this state ([ρ A , e 2πiαQ A ] = 0). After a straightforward calculation one gets Putting the previous two pieces together, we arrive at the expected final result Tr ρ n A e 2πiαQ A = r n e 2πiα + (1 − r) n , with r = /L, which provides the exact charged moments of a single magnon state.

Two-Magnon States
In the following we consider a state of two magnons with the same symmetry charge. This example is more interesting because it allows us to test whether the presence of non-trivial interaction changes our results. Given a pair of quasimomenta p and p , we parametrise this state in the following way where S is a scattering matrix and |jj is the state with two localised magnons in sites j and j . The choice of the S-matrix is not really relevant for our purpose, but for the sake of concreteness we set using the same conventions as in [62]. The action of the restricted symmetry operator e 2πiαQ A on the two-particle sector of the Hilbert space is e 2πiαQ A |jj = e 2πiα(δ j∈A +δ j ∈A ) |jj .
It is possible to decompose ρ A ≡ TrĀ |Ψ 2 Ψ 2 | as follows where ρ A is the two-particle contribution (both particles in A), ρ A is the vacuum contribution (no particles in A) and ρ

(3)
A is the one-particle contribution (one particle in A and one inĀ). The introduction of the flux gives rise to the following relation No approximation was made up to this point, but the explicit expressions of ρ

(j)
A , given in [62], are cumbersome and not particularly enlightening for our purpose. However, one can show that in the limit L → ∞ and /L, p = p kept fixed, Tr A (ρ (j) A ) n simplifies drastically: Putting all the pieces together one finally gets This computation shows that in this particular scaling limit the interaction between particles has no effect on the final result, and the total charged moment is just a product of two singleparticle charged moments. A different result is obtained if p = p and fixed. In that case, the magnons are indistinguishable and one can prove that so that Tr A ρ n A e 2πiαQ A r 2n e 4πiα + 2 n r n (1 − r) n e 2πiα + (1 − r) 2n , which no longer factorises into one-magnon contributions. Both results are special cases of (8) and (10). The results of this section generalise previous work for the excess entanglement entropy of excited states [63] and are also related to the results of [71] where the entanglement of magnon states was considered more generally. In particular, it was shown that for states consisting of several magnons, entanglement will factorise into the contributions of groups of magnons which are well-separated from each other in momentum space (that is their momentum difference is of order 1 rather than 1/L, where L is the length of the system). Such results also apply to the present case up to the introduction of the appropriate phases.

Further Generalisations: Higher Dimensions
So far we have derived the behaviour of the charged moments of the SREE of quasiparticle excited states making use of two different formalisms: the form factor expansion in 1+1 integrable QFTs [61] and the analysis of qubit/magnon states on the lattice. Unfortunately these techniques are specially suited for 1+1D theories. In this section we want to consider instead a generic QFT in higher dimensions. To this aim we introduce a slightly different approach. Indeed, on the one hand, the description of the Rényi entropy as a correlation functions of branch point twist fields inserted at some points of the space is special of 1+1 dimensional QFT [8,9]. On the other hand, the description of excited states as magnon/qubit states, even if generalisable to higher dimensions, has the disadvantage that it does not take into account the zero-point fluctuations.
Despite these technical limitations, we expect that in the particular scaling limit we are considering the universal entanglement content of the Rényi entropy, together with its symmetry resolved version, should not depend on dimensionality, presence or not of interactions or even integrability. To support these claims, Section 3 on magnonic states provided evidence that the same results tested numerically for free theories in Section 2 can be recovered for interacting systems. As for higher dimensions, at least one precedent for this generalisation already exists. In a previous work [65] the excitations of the free massive boson in D := d+1 dimensions have been extensively analysed and their Rényi entropy was computed in terms of graph partition functions, and found once more to fit the same formulae, with r replaced by the ratio of generalised volumes.
In this section we slightly generalise the formalism of [65], to take into account possible interactions, and provide, as a proof of concept, a simple calculation of symmetry resolved entanglement of a single-particle excited state. The key ingredients we need are the description of the excited states as local operators acting on a vacuum state and a twist field operator, which generalizes the composite branch point twist field to higher dimensional settings (an introduction to composite twist fields can be found in [61]). The only strong assumption we make in our derivation is the presence of a finite mass gap m, whose inverse m −1 is much smaller than the typical lengths of the system. 3 We anticipate here that our formulae (8)- (11) are unchanged in higher dimensional theories, up to the identification which is the ratio of (generalised) volumes between subsystem A and the total system. It may seem surprising that results should only depend on r, as defined in (77), and not on other features of the entanglement region, such as the connectivity and smoothness of its boundary. Indeed, the charged moments and symmetry resolved entropies of both the ground state and of excited states will depend on such properties, as would finite volume corrections to our results. However, our computations deliver results for the ratio of charged moments between the excited and ground states in the infinite-volume limit, and it is this ratio in this limit which is universal and independent of boundary features, not the charged moments themselves. This independence of boundary features has been analytically shown for one simple example, namely the case of one-dimensional disconnected regions, where the same formulae as for one connected region were found to apply, with r the sum of the lengths of all disconnected parts [63].

Excited States and Operator Algebra
Let us consider the vacuum state |0 of a Hilbert space H, together with an algebra A of observables 4 acting on H which has |0 as a cylic vector (see [75] for a modern review of this algebraic viewpoint in QFT). This allows us to represent any state |Ψ of the Hilbert space as We would like to assume further that the vacuum state is translation invariant, namely that it is invariant under a certain faithful representation of the translation group in d dimensions. Strictly speaking, since we consider a finite-size system, we have to slightly modify this requirement. Specifically, we put our system on a d-dimensional torus M of volume V and we require that |0 is invariant under the isometries of the torus. Other boundary conditions can be considered too, but they do not change the picture in the scaling limit we are interested in. We also require locality of the observables, asking that A is generated by a set of fields {O(x)}, which can be applied at any point of M.

For any field O(x) one can construct its Fourier transform O(p) as
and these transformed fields are building blocks for the following set of translation invariant states The state above corresponds physically to k particles distributed on M with momenta {p j } j=1,...,k , and the choice of the fields {O j } may depend on the particle species and quantum numbers. At first sight, this construction is similar to the usual way of generating particle states in free theories acting with creation operators on the vacuum on a Fock space. However, the real advantage of our formulation is that it is directly related to local observables, a property which is fundamental to correctly define entanglement measures. Let us take a set of orthogonal fields {O j }, so that the correlation function O † i (x)O j (x ) vanishes for i = j. In other words the fusion is present only if i = j and their operator product expansion (OPE) can be expressed formally as where we neglected explicitly the contributions coming from less relevant operators (with respect to the identity). The exact evaluation of the correlation function above can be hard, but the assumption of a finite gap m ensures that it vanishes exponentially for |x − x | m −1 . This is the only property we really need in our subsequent discussion.
Consider a smeared version of the modes O(p), with support in a subsystem only, that is a region of space. To each spacial region A ⊆ M and field O(x), we associate Given any two regions A and A , we , making use of some approximations. First, we consider only the most relevant term in the fields OPE, namely Second, since we are working in the limit of small correlation length (compared to the geometry), the leading contribution comes from the insertion of the fields at small distances, which is present if x, x ∈ A∩A ; this observation motivates the change of variable x = x −x, and the subsequent approximation The second integral may be difficult to compute and in principle it could require a UV regularisation for |x | < m −1 . However it does not depend on the regions A, A and in our computation appears only as a multiplicative constant. In conclusion, we end up with where V A∩A is the volume of A ∩ A , which is the main result of this subsection. Since the volume in (86) emerges from the integrals (85), which involve a Fourier transform, we require that subsystem A ∩ A consists of a finite number of disconnected regions, whose boundaries are piecewise smooth.
It is natural to ask how our discussion above would be modified for a vanishing gap m = 0. The main change is in the scaling of correlations functions: exponential localization of the correlation function in a region of typical length m −1 does not hold any longer, due to the long algebraic tails of the correlation functions. We conjecture that, as long as the momenta are fixed in the infinite-volume limit, the main conclusion (86) is unchanged. A qualitative argument is that in this case the inverse momentum, say the De Broglie length, plays the role of typical length scale. In order to make this consideration more precise, let us analyse Eq. (85) for a 1+1D CFT, where O is a field of conformal dimension ∆ O . We focus on the following integral 5 One should note that hermitian conjugation and Fourier transform do not commute. Indeed, we have that which we regulate both in the UV, with a cutoff , and in the IR, with a cutoff L, as followŝ This integral can be explicitly computed. However, the important feature is that for ∆ O > 0, p > 0 and > 0 all fixed, the integral converges to a finite value when L → +∞. This is no longer the case if p ∼ 1/L in the infinite-volume limit. In practice, this means that for small momentum and scaling dimension 0 < ∆ O ≤ 1 the considerations we made so far regarding the scaling at large sizes cannot be applied. As a matter of fact, for free CFTs the scaling dimensions of the fundamental fields are smaller than 1: the fermionic field Ψ has dimension 1/2 while the derivative of a compact boson ∂ x Φ has dimension 1.
While these considerations are not mathematically rigorous in establishing convergence of the OPE expansion in the large volume limit, they are sufficient to explain why low-energy states of gapless theories, or multiparticle states with small momenta difference, are not well captured by our predictions. Indeed, for such states the excess entanglement was computed in [76,77] and is clearly different from the formulae in [62,63].

Replica Construction for Symmetry Resolved Entanglement
Consider now a replica version of the theory, consisting of n copies. For any state |Ψ we consider its replicated version |Ψ n . Our goal is to define a composite twist operator, which generalises the composite branch point twist field as defined in [4], generalising [9], to higher dimensional theories. It is well known from a large body of literature e.g. [4,[36][37][38][39] that the charged moments of the symmetry resolved Rényi entropies are given by correlation functions of composite twist fields (CTF). A short review of their main properties can be found in [61].
In this section, we want to extend the notion of CTFs to higher dimensions. To this aim we will introduce a composite twist operator whose expectation value over |Ψ n gives exactly the charged moments Z n (α). This type of operator was already considered in the literature, especially in the absence of the flux insertion (see for example [65,[78][79][80]), but here we are mostly interested in its relationship with the algebra of local operators.
The first point we have to clarify regards the symmetry, and its action on the space of fields. Starting from e 2πiαQ , the global generator of U (1) symmetry in the non-replicated theory, we Since one can decompose the space of fields as irreducible representations of U (1), we restrict our analysis to charged fields. Going back to the replicated theory, we consider the algebra of replicated observables A n as the algebra generated by the tensor product of n observables in A. Thus, to any field O(x) ∈ A, we associate O j (x) ∈ A n defined as where O(x) lives only in the jth replica. Consider a spacial region A, and its complementĀ. We define a composite twist operator T α A which implements the structure of the n-sheeted, cyclically connected, Riemann surface where the replica theory is defined. That is T α A "implements" the gluing of the replicas along A with an additional flux insertion due to the action of the U (1) symmetry. As for a standard CTF, it does so via its commutation relations with any charged field O j (x), which we require to be Using this choice, the flux is inserted only between the nth and the first replica. We would like to emphasize that a similar definition has already appeared in the context of 1+1D integrable QFTs (see [36][37][38][39]). In particular for A = [0, ], one can identify and the commutation relations for T α A can be expressed as commutation relations for the CTF T α n . The only novelty here is that for higher-dimensional theories it is not clear what the replacement for CTFs is, thus how to represent T α n . Regarding the definition of the twist operator T α A , it is worth mentioning also that in QFT it is known [81,82] that the twist operators are not local observables of the algebra A n . They are rather observables of the orbifolded algebra usually written as A n /Z n , which extends A n and is obtained by taking the coset over the cyclic symmetry of replicas.
We are now ready to relate the twist operator to the symmetry resolved entanglement. The charged moments of |Ψ are given by The above definition, together with the commutation relations (91) and the OPE of Eq. (86), should be enough to prove the explicit analytical expression of the ratio of charged moments between |Ψ and the ground state defined in (7). We now show how these ideas come together for a simple example.

Single-Particle State
In this section we analyse an excited state |Ψ made of a single quasiparticle with momentum p generated by a charged field O. Its explicit expression is given by and the replicated version is just For the sake of convenience, we split so that its commutation relations with T α A become more transparent. Indeed, using just (91) one can express Up to now everything looks exact, while the approximate evaluation of n Ψ|T α A |Ψ n makes use of the OPE contraction in (86). Note also how the phase e 2πiα is only present for nth factor above, similar to the Gaussian measure found by employing wave-functionals in Section 2. Among all the terms which are generated, all but two are vanishing and they give Similarly, we can evaluate the norm n Ψ|Ψ n which does not require the splitting of O j (p) In the evaluation of the ratio n Ψ|T α A |Ψ n n Ψ|Ψ n (100) the proportionality constant (which is non-universal and could be absorbed in a redefinition of the field) cancels out, and one can write with r = V A V . In the expression above, the first piece is universal while the second is not, it is just the nth charged moment of the ground-state. Taking the ratio with the ground state contribution, we finally arrive to the desired result Results for multiparticle states can be obtained in a similar fashion.
To conclude, the striking simplicity of these results, and the possible generalisations to multiparticle states, rely especially on the truncation of the OPE in (82), which is expected to become exact in the limit m d V 1. We further expect that for finite m d V further contributions in the OPE can be recast as a (possibly non-integer) power series in (m d V ) −1 , which generalises the explicit (mL) −1 power expansion that is obtained for 1+1D free theories using form factor techniques [61]. In massless theories, we instead expect corrections as a power series in (|p| d V ) −1 .
The explicit evaluation of these corrections, which are expected to be non-universal, that is, momentum and QFT-dependent, and any possible issues regarding the convergence of these power series are all beyond the scope of this work.

Conclusions
In this paper we have extended the results of [61] by first, providing numerical evidence for their validity and second, showing that they apply much more broadly than the form factor computation of [61] would suggest, to interacting and higher dimensional theories.
Regarding the numerics, we have considered a 1D Fermi gas and a complex harmonic chain. Although both these models are discrete and amenable to numerical computations, their microscopic features are extremely different. Whereas the 1D Fermi gas is a fermionic theory whose continuous limit is a massless complex free fermion, the complex harmonic chain's continuous limit is a complex massive free boson. It is therefore quite remarkable that the same set of formulae for the ratios of charged moments should apply for both theories. This is nonetheless the case. Whereas for the 1D Fermi gas our formulae hold for highly excited states containing excitations of large momenta, thus small De Bloglie wavelengths, for the complex harmonic charge, the formulae hold as long as the correlation length is small compare to subsystem size and the excited states are usually of low energy. The common feature of both types of states is that they are characterised by localised excitations.
A simple class of states in one space dimension are magnon states. These are typical excited states of quantum spin chains, with and without interaction. They have a simple representation in terms of spin degrees of freedom and their entanglement entropies can be computed analytically. We have shown that, irrespective of interaction, the formulae found in [61] still apply, in line with observations made in [63,71].
Finally, we have proposed a method to compute the charged moments of zero-density excited states in higher-dimensions with or without interactions. The generalisation relies on generalising the action of branch point twist fields to higher dimensions (without generalising the fields themselves) and on standard assumptions about the asymptotics of correlations functions of local fields.
There are several extensions of this work that are still outstanding: the study of finite-volume corrections to these results in a QFT setting, the extension of our formulae to other symmetry resolved measures of entanglement, such as the logarithmic negativity [22,53,64], and the study of theories where quasiparticle excitations are not localised (such as free fermion chains where the local degrees of freedom are the spins, which are non-local with respect to fermions). We hope to return to these problems in future work.

A Trace Calculations via Wave Functional Method
In this Appendix we present explicit calculations of the ratio of charged moments for excited states in a complex free boson theory using the wave functional method introduced in Section 2.2. The discretisation of these results leads to the formulae for the complex harmonic chain presented in Subsection 2.3.

A.1 Zero flux
We now wish to compute the ratio of charged moments where ρ A and ρ 0,A are the reduced density matrices of the excited and ground states, respectively. We define restricted wave functionals which take as arguments (complex) fields that either have support on region A ≡ [0, ) or onĀ ≡ [ , L]. Leaving the dependence on the conjugate fields implicit, we write: and similarly for the action of Φ † (x), iΠ † (x). The action of the charge operator Q A on the wave functional is simply: while from (47) it follows that The reduced density matrix then admits the functional representation where Dφ is shorthand for DφDφ † . We start with the simple case where there is no symmetry resolution, that is, we compute (103) for α = 0, following [63]. For n integer, we can insert the resolution of the identity n times, with the identification n + 1 ≡ 1, and we get: The product of the diagonal terms, i.e. those in which the fields act on the same copy for both subsystems is with K(x) defined in (40). For the non-diagonal terms we notice that since K(x) * = K(−x) we have Ψ vac [φ i+1,A , φ i,Ā ] * = Ψ vac [φ i+1,A , φ i,Ā ] and: Putting all the terms together, we end up with the Gaussian measure: where The numerator of (103) instead is for a state of k ± particles of charges ±1 and momenta p j ± (assuming p j ± = p i ± for j ± = i ± ), and using (43): Putting everything together we have that where the correlation function is defined with respect to the Gaussian measure: and the operators are If the excited state is of the form (44), the result (120) is minimally modified. The terms inside the correlator are exactly the same, except that now the range of the indices j ± is that specified after (44) and the prefactor is modified to:

A.2 non-trivial flux insertion
We now come to the quantity (103) with α = 0, assuming the excited state to be of the form (43). Let us consider the denominator first. Because of the flux insertion to the right of the nth operator ρ 0,A , everything is the same as in (109) except for the last resolution of the identity, which produces a term: Ψ|φ nĀ e 2πiαQ A |φ 1A = Ψ|e 2πiα φ 1A , φ nĀ = Ψ[φ 1A e 2πiα , φ nĀ ] * .
Thus, the Gaussian measure in the presence of the charge is modified as follows: with As for the numerator of (103), we have similarly The product in the first line gives what we already had found in the α = 0 case, with the replica index running up to n − 1 only. The product of the last two functionals gives: and from the definitions of the α and β operators we obtain our final result: where . . . α is defined as in (121) with respect to the Gaussian measure (128). U ± i (p) are defined exactly as in (122), while V ± i (p) are now modified by the presence of the charge as follows: If the excited state is given by (44), the right-hand side of (131) is changed exactly as in the case α = 0.