Symmetry resolved entanglement in gapped integrable systems: a corner transfer matrix approach

We study the symmetry resolved entanglement entropies in gapped integrable lattice models. We use the corner transfer matrix to investigate two prototypical gapped systems: the harmonic chain and the XXZ spin-chain. While the former is a free bosonic system, the latter is genuinely interacting. We focus on a subsystem being half of an infinitely long chain. In both models, we obtain exact expressions for the charged moments and for the symmetry resolved entropies. While for the spin chain we found equipartition of entanglement (i.e. all the symmetry resolved entropies are the same), this is not the case for the harmonic system where equipartition is effectively recovered only in some limits. Exploiting the gaussianity of the harmonic chain, we also develop an exact correlation matrix approach to the symmetry resolved entanglement that allows us to test numerically our analytic results. ar X iv :1 91 1. 09 58 8v 1 [ co nd -m at .s ta tm ec h] 2 1 N ov 2 01 9 Symmetry resolved entanglement in gapped integrable systems: a CTM approach 2


Introduction
Entanglement is a characteristic treat of quantum mechanics since its early days. However, only in the last two decades it became clear that entanglement is an important concept also for many-body systems with ramifications to many different lines of research, ranging from high energy physics and gravity to quantum information and critical or topological extended quantum systems (see e.g. Refs. [1][2][3][4] as reviews). The most successful and used measures of the bipartite entanglement are surely the Rényi and von Neumann entropies, defined as follows. Let |Ψ〉 be a pure state of an extended quantum mechanical system and ρ = |Ψ〉 〈Ψ| its density matrix. Let us consider a bipartition of the system into A and B and define the reduced density matrix (RDM) of the subsystem A as the partial trace over the degrees of freedom of B, i.e. ρ A = Tr B ρ. A measure of the entanglement between A and B is the Rényi entropy of order n The Von Neumann entanglement entropy S 1 ≡ −Trρ A log ρ A is the limit n → 1 of the Rényi entropy. In a quantum field theory, Trρ n A for integer n can be expressed in the path integral formalism as a partition function over suitable n-sheeted Riemann surfaces. For the ground state of critical one-dimensional systems with an underlying conformal field theory, this led to a remarkable universal scaling depending only on the central charge c [5][6][7][8][9][10].
Such a universal behaviour is not strictly a prerogative of the gapless models, but it also occurs for gapped models in the vicinity of a quantum phase transition in the regime in which the correlation length ξ is large but finite [5]. Indeed, using ideas from the famous proof of the c-theorem by Zamolodchikov [11], it has been shown that for a bipartition of an infinite system into two semi-infinite halves, the leading behaviour of entanglement entropies is generically [5] S n c 12 1 This result can be elegantly recovered for integrable lattice models through the Baxter corner transfer matrix (CTM) [12], as reported (and generalised) in many references [5,[13][14][15][16][17][18][19][20][21][22][23]. We will discuss explicitly this technique in the following sections. The CTM approach provided exact results not only close to the critical point, but gave generalisations also to the regime in which the correlation length is small. When the subsystem A is a finite interval of length , as long as ξ, the Rényi entropies are just twice the value in Eq. (2) as a consequence of cluster decomposition in the ground-state of these theories. However, as becomes of the order of ξ, a complicated crossover takes place that is not captured by CTM and requires more complicated techniques, see e.g. Refs. [24][25][26][27].
Only in very recent times, it became clear that it is also important to understand the relation between entanglement and symmetries and in particular how entanglement is shared between the various symmetry sectors of a theory [28,30]. The physical motivations for shading light on the interplay between symmetry and entanglement are manifold. For example, one motivation comes from a recent experiment studying the time evolution of the symmetry resolved entanglement in systems with many body localisation [29]. It has been shown that entanglement has two different contributions, called configurational and fluctuation entanglement (see below, cf. Eq. (4), for a precise definition). These two contributions account for the entanglement within symmetry sectors and fluctuations thereof, respectively. In the presence of both disorder and interaction, their dynamics occur over different time scales: the fluctuation entanglement quickly saturates to an asymptotic value while the configurational one exhibits a slow logarithmic growth [29], providing a nice physical explanation of an older finding [31,32]. The possibility of measuring these quantities sparked the interest in further studying how the entanglement is related to the internal symmetries of a system, leading to many results concerning critical ones [30,[33][34][35][36][37][38][39]. A surprising finding is that conformal invariance forces the entanglement entropy to be equally distributed among the different sectors of a U(1) symmetric theory [34]. It is an open issue to understand whether and when such equipartition of entanglement survives away from criticality. However, to date there are no results concerning gapped systems (with the exception of Ref. [38] for a discrete symmetry, but here we are interested in continuous ones). The goal of this work is to fill this gap and to study how the total entanglement splits into the contributions coming from disjoint symmetry sectors in gapped integrable models, using CTM techniques. We carry out this analysis for two non-critical quantum lattice models with a U(1) symmetry, namely the double or complex harmonic chain (which is a free model) and the XXZ chain (which is genuinely interacting). To this aim, we first calculate the moments of the RDM in the presence of a charge flux, that we call charged moments, and then obtain the contributions of the sectors by Fourier transform.
The manuscript is organised as follows. In Section 2 we briefly review all the quantities of interest and we give an overview of how the RDM of an off-critical quantum chain is related to Baxter's CTM. For integrable models whose weights satisfy a Yang-Baxter relation, the eigenvalues of the RDM can be determined exactly. In Sections 3 and 4 we exploit these exact results for the computation of the symmetry resolved entanglement entropy, for the complex harmonic chain and XXZ spin-chain respectively. We also benchmark the analytic results in Section 3 against exact numerical computations. We conclude in Section 5 with some remarks and discussions. Many technical details of the calculations can be found in four appendices.

Symmetry resolution, flux insertion, and corner transfer matrix
We consider a system with a U(1) symmetry, generated by the charge operator Q, which obeys Q A ⊕ Q B = Q, where Q i is the charge in the subsystem i. If the system described by the density matrix ρ is in an eigenstate of Q, then [ρ, Q] = 0. We are interested in a bipartition of the total system into two semi-infinite halves, A and B, and we denote by ρ A the reduced density matrix of A. Taking the trace over B of [ρ, Q] = 0, we find that [ρ A , Q A ] = 0. This means that ρ A is block-diagonal and each block corresponds to a different charge sector labelled by the eigenvalue q of Q A . Therefore we can write where p A (q) is the probability of finding q in a measurement of Q A in the RDM ρ A , i.e. p(q) = TrΠ q ρ A , where Π q is the projector on the sector of charge q. Within this convention, the density matrices ρ A (q) of different blocks are normalised as trρ A (q) = 1. Now, to understand how the total entanglement arranges into contributions coming from the disjoint charge sectors, we first define the symmetry resolved entanglement entropy as The total von Neumann entanglement entropy associated to ρ A in Eq. (3) can be then written as Let us describe the physical meaning of the two sums in Eq. (5) [29,34,40]. The first contribution is known as configurational entanglement entropy and it depends on the entropy of each charge sector, weighted with its probability. The second contribution is the fluctuation entanglement entropy which is due, as the name says, to the fluctuations of the charge within the subsystem. The configurational entropy is related to the operationally accessible entanglement entropy of Refs. [40][41][42].
For future use, we also define the symmetry resolved Rényi entropies as In order to compute these quantities, following the approach of Ref. [30] we first define the normalised charged moments of ρ A as In a (1+1)-dimensional quantum field theory, this quantity is the partition function on a Riemann surface with the insertion of an Aharonov-Bohm flux α, such that the field acquires a total phase α when moving on the entire worldsheet. Similar charged moments have been already considered in the context of free field theories [43][44][45], in holographic settings [46,47], as well as in the study of entanglement in mixed states [48,49]. The Fourier transforms of the charged moments are just the moments of the RDM restricted to the sector of fixed charge [30], i.e.
Hence the symmetry resolved entropies can be obtained as Finally, also the probability p(q) is simply related to the moments Z n as

The corner transfer matrix and the entanglement entropy
In dealing with the geometric bipartition considered in this paper (i.e. two semi-infinite half lines) the corner transfer matrix provides an exact form for the reduced density matrix [50] and hence it is a formidable tool for the derivation of the charged moments and symmetry resolved entropies. In order to understand how the CTM works, we give a brief review of the construction of the RDM. Generally, a direct computation of the density matrix of a system is tough. A trick to address this problem is to use the fact that the density matrix of the quantum chain is the partition function of a two-dimensional classical system on a strip [51][52][53]. The latter can be solved by means of the transfer matrix T and we can identify the eigenstate |Ψ〉 of T corresponding to its maximal eigenvalue. Given the Hamiltonian of the quantum chain H and its lattice spacing a, the transfer matrix is T = e −aH up to a prefactor; hence |Ψ〉 is the ground state of H. One then obtains the reduced density matrix of a subsystem A of the chain by tracing over all the coordinates belonging to the complement of A. Therefore ρ A is the partition function of two half-infinite strips, one extending from −∞ to 0 and the other from +∞ to 0.
The CTM plays a crucial role: it connects a horizontal row to a vertical one. Choosing the lattice in a clever way [12], when the model is isotropic, the four possible corner transfer matrices [12] are all equivalent and the partition function is just TrÂ 4 , withÂ the CTM. Going back to our quantum problem, the reduced density matrix is [50] We will deal with integrable massive models satisfying the Yang-Baxter equations; in this case, it is possible to show that Eq. (11) has an exponential form given by [50,51] H CTM is known as entanglement or modular Hamiltonian, that in the cases we are interested in can be diagonalised as [50] where n j are number operators and ε j are the single-particle levels of the entanglement Hamiltonian. The result (13) provides exact eigenvalues and degeneracies of the RDM (i.e. the entanglement spectrum of the system [54,55]), from which one calculates straightforwardly the entanglement entropies [5]. However, Eq. (13) contains no information about the distributions of the eigenvalues ε j into the various symmetry sectors (indeed, it has exactly the same form for models with discrete and continuous symmetries). In order to use it to compute the symmetry resolved entropies in gapped integrable models, we should complement Eq. (13) with some other input providing the symmetry resolution, but this should be done on a case by case basis. The rest of the manuscript is devoted exactly to solve this problem for two specific 1D integrable lattice models: the complex harmonic chain and the non-critical XXZ chain in which we will exploit the results of Refs. [56,57] and [58,59] respectively.

The complex harmonic chain
In this section we use the CTM to derive the symmetry resolved entanglement entropies for a double or complex harmonic chain that is U(1) symmetric and its continuum limit is a noncompact massive complex boson, i.e. a Klein-Gordon field theory. We will find an analytic expression for the charged moments as functions of α and we will discuss its limit close to the conformal invariant critical point, when the correlation length ξ is finite but large. Then we will use this result to compute the symmetry resolved entropies. All the analytical results will be compared against exact numerical computations based on correlation matrix techniques [60][61][62].

Brief recap of the free complex scalar field and its lattice discretisation
The meaning of the symmetry of a double harmonic chain is clearer in the field theory language and so we first consider a free complex scalar field φ(x) described by the Euclidean action This action is invariant under U(1), i.e. the field φ can be rotated of an arbitrary phase φ(x) → e iθ φ(x) leaving the action unchanged. The Hamiltonian of this field theory is with Π(x) being the field conjugated to φ(x). We can as well rewrite the model in terms of two scalar real fields φ x (x) and φ y (x) and the same for Π(x). In these variables the U(1) symmetry is an O(2) rotation in the plane (φ x , φ y ). The Hamiltonian (15) in terms of these variables is where in the second line we stressed that it is a sum of two identical Hamiltonians H for the real fields φ x and φ y . One introduces the modes a † i (p) and a i (p) for each field i = x, y and momentum p. The Hamiltonian and the conserved charge are instead better written in terms of particles and antiparticles modes operators The Hamiltonian is (with ε 2 (p) = m 2 + p 2 ) while the conserved charge is i.e. the total number of particles minus the number of antiparticles. The conserved charge can be as well written in real space and its value in a given subsystem A is the same integral restricted to A, i.e.
For the construction of the RDM for the lattice version of the complex Klein-Gordon field theory, we start from discretising each of the two real Hamiltonians in Eq. (17). The lattice discretisation of each of them is the harmonic chain, i.e. a chain of L harmonic oscillators of mass M = 1 with equal frequency ω 0 , coupled together by springs with elastic constant k (hereafter we set ω 0 = 1 − k), i.e. the lattice discretisation of the Hamiltonian H is where variables p i and q i satisfy standard bosonic commutation relations [q i , q j ] = [p i , p j ] = 0 and [q i , p j ] = iδ i j . Hence, the lattice version of the complex field theory is the sum of two of the above harmonic chains in the variables q x and q y , i.e.
which we call complex or double harmonic chain. The reduced density matrix, ρ A , for half of the real harmonic chain was explicitly constructed by Peschel and Chung in [57] in the large L limit. The trick is to relate ρ A to the partition function of a two-dimensional massive Gaussian model in the geometry of an infinite strip of width L with a cut perpendicular to it [56]. Due to the integrability of the Gaussian model, in the case where L is much larger than the correlation length, the H CTM for the harmonic chain may be written in a diagonal form as in Eq. (13), where now we explicitly have Here I(k) is the complete elliptic integral of the first kind, i.e.
and β j , β † j are bosonic annihilation and creation operators (satisfying [β i , β † j ] = δ i, j ). They are related to the ladder operators a i of the original chain by a generalised Bogoliubov transformation [57] as Notice that the transformation mixes a and a † so it does not conserve the number operator. The RDM for the double chain is clearly factorised in x and y part, i.e. the entanglement Hamiltonian is the sum of two H C T M in Eq. (24) one with β x,i and one with β y,i ladder operators. Now we proceed as follows. First we rewrite these two entanglement Hamiltonians in terms of the local ladder operators a x,i and a y,i using the inverse of the Bogoliubov transformation (26). Then, using the lattice analogue of (18), i.e.
we rewrite the entanglement hamiltonian in terms of local ladder operators for particles and antiparticles. This is clearly quadratic (it is the rewriting of a quadratic operator after two linear transformations and so it is quadratic) and commute with the charge operator. Hence, via another Bogoliubov transformation (see Appendix A) which conserve the charge, the entire entanglement Hamiltonian of half-chain is brought into the form The charge operator restricted to the semi-infinite line is just the discretisation of Eq. (21), i.e.
Once we apply the Bogoliubov transformation in Eq. (28), we have up to an unimportant additive constant that we neglect. Since the α i and γ i operators in Eq. (29) commute, the RDM factorises as where we denoted the RDM for α i and γ i with ρ α A and ρ γ A respectively. For the charged moments, we need to compute Trρ n A e iQ A α , but using also that Q A is the difference of the number of α i 's and γ i 's, see Eq. (31), the trace factorises as where N α A = j∈A α † j α j and N γ A = j∈A γ † j γ j . The two factors are equal, except for the sign of α. It is very instructive to see how this factorisation happens for a chain of two oscillators as we report in Appendix A.
If for a single harmonic chain, we introduce the quantity then we have that the charged moments of the complex boson are given by We stress that F n (α) is not the log of a local charged moment because in the single harmonic chain there is no local U(1) symmetry. In the following we show how to compute F n (α) by CTM methods for a single harmonic chain and after we use (35) to get the charged moments.

Charged moments from CTM
Here we first compute the quantity F n (α) for a real harmonic chain and from this Z n (α) is simply derived from Eq. (35). In the above subsection, N A and ρ A for the single chain have been already written in the same basis and the derivation of F n (α) amounts to compute the trace whose logarithm is given by This formula is exact and can be easily computed numerically, since it converges very quickly.
It is plotted in Figure 1 as a function of α for various values of ω 0 and n, but we will discuss its properties later. The charged moments for the complex harmonic chain, cf. Eq. (35), are , (38) where in the last equality we factor out the total partition sum and use the definition (121) for θ 4 (u|q). Notice that the entire α dependence is encoded in the denominator of Eq. (38) and that Z 1 = 1, but Z 1 (α) = 1. Also the total Rényi entropies of the complex harmonic chains are i.e. the double of a real harmonic chain.

Poisson resummation and critical regime.
A drawback of the form (37) is that it does not directly allow a direct expansion in the critical regime, i. e. for small ε. Moreover, we cannot perform an Euler Mac-Laurin summation (as for α = 0, see [5]) since the function f (x) = log(1 − e −2x ) diverges for x → 0. However, following Ref. [14], we can obtain the asymptotic expansion for small ε by using the (generalised) Poisson resummation formula: wheref In order to use this resummation formula for Eq. (37), we must choose a = 1/2, b = 1 and which allows us to rewrite the sum (37) as The cosine-Fourier transform of f n,α (x) iŝ where Φ is the Lerch transcendent function, defined as If α = 0 and n = 1, Eq. (45) simplifies to the known value [14] Plugging Eq. (45) into the Poisson resummation formula, we rewrite log Z n (α) in such a way to isolate the contribution of the term k = 0 which gives the leading divergence in the limit ε → 0, i.e.
Here we have introduced the polylogarithm of order 2 We are now interested in the critical region of the parameter space in which the correlation length ξ (inverse gap) is large but finite. In the critical regime ξ 1 (or equivalently ε 1), the correlation length of the model behaves like Using the results of Ref. [63], the last sum over k in Eq. (48) in the limit ε → 0 behaves like Eq. (52). As discussed in the text, the convergence to the critical result is not uniform and it is slower for smaller α = 0. The function log Z n (α) for the complex chain is twice the real part of F n (α).
and hence the only non-vanishing terms in the asymptotic expansion close to ε = 0 are whose real part is because The charged moments for the complex harmonic chain are now given by Eq. (35), i.e. log Z n (α) = F n (α) + F n (−α) and, in the limit ε → 0, Notice that while F n (α) is generically complex, log Z n (α) for the complex chain is real and even in α.

Discussions.
We concluded our exact computation of the charged moments and we are now ready to critically discuss our findings. Eq. (55) is very suggestive. It tells us that the leading term in the "charged entropies" diverges logarithmically with ξ but with a non-standard prefactor. Indeed, in the conformal field theory of the compactified boson, it has been found that when α = 0, the additional term in the logarithm is proportional to α 2 [30] , while here we also have a linear contribution in α. Obviously the two results are not in contradiction, because the continuous limit of the harmonic chain is non-compact and the prefactor of α 2 in Ref. [30] diverges when the compactification radius is sent to infinity. These results are very intriguing and it would be interesting to recover them directly in a field theory approach; work in this direction is in progress [64].
Another interesting fact is that the limit α → 0 and the expansion for ε around 0 do not commute, as a difference with other known cases (we believe that the origin of the noncommutativity is the non compact nature of the continuum limit). Indeed, if we consider first the limit α → 0 in Eq. (48), the last sum gives leading to the known formula of the Rényi entropies of a real harmonic chain, that in the critical regime ε → 0 is [5,14] (see Eq. (40)) On the other hand, if we invert the order of these two operations, we obtain the divergent term in Eq. (51). Considering now the charged moments of the complex chain, ln Z n (α) = 2ReF n (α), the divergent term (51) cancels, but the finite part is not the total moment ln Z n in Eq. (39). This fact implies that the approach of ln Z n (α) to the critical limit ε → 0 is non-uniform in α: exactly at α = 0 the charged entropy approaches (40), but for any non-zero α the limit is (52) that as a consequence is reached for smaller and smaller ε (i.e. ω 0 ) as α gets closer to 0. All these aspects are evident in Figure 1 where we show (for α ≥ 0 since F n (−α) = F * n (α)) the exact result Eq. (48) (or equivalently (37)) together with its critical limit, Eq. (52). As we discussed above, the former converges to the latter as ω 0 , therefore ε, decreases, but in a non-uniform way. Indeed, while for large α (i.e. close to π) the two curves are very close also when ω 0 is not so small, for smaller and non-zero values of α, we need much smaller ω 0 to approach the critical limit. For α = 0 the limit is different. It is also clear that for higher values of n, the convergence is slower and starts at smaller values of ω 0 . The last observation is a well known fact for α = 0, cf. Ref. [14], and it is not surprising that the effect is amplified in the presence of a flux.

Symmetry resolved moments and entropies via Fourier trasform
The symmetry resolved moments Z n (q) are obtained as Fourier transform of Z n (α) in Eq. (38), i.e.
The integral in the rhs of the above equation can be found in Ref. [65] (exercise 14 at page 489), obtaining which is our final result for the symmetry resolved moments. It is likely that the sum in Eq. (59) can be rewritten in terms of some special functions, but we did not find any particularly useful expression. We define the sum as which can be written in few different equivalent ways that are useful for investigating diverse properties: Clearly in terms of this function we have where we used the explicit form of Z n in Eq. (39). The symmetry resolved Rényi entropies are now easily deduced from Eq. (9), obtaining Taking the limit n → 1, we get the von Neumann entropy The critical limit ε → 0 is easily understood if one focuses on the variation in q of moments and entropies, rather than on their absolute values. Indeed from Eq. (62), it is easy to see that where the last limit is performed by expanding to the second order in ε each term in the sum (60), making carefully the sum in terms of ζ functions, and finally re-exponentiating the result. We stress that this critical limit is not the Fourier transform of the critical limit for Z n (α) in Eq. (55) because the two limiting procedures do not commute. The critical behaviour of the resolved entropies is then easily worked out as which is valid also for n = 1 without any particular limit. Also in the critical limit, it is worth to mention the behaviour which signals the presence of a subleading term proportional to log ε ∼ log(log ξ). Such a term has not a unique interpretation and origin for the (complex) harmonic chain. Indeed, we know that the total entropy of a massive free non-compact boson has such subleading terms in log(log ξ) [66] in the small mass limit, but even that double logarithmic terms appear generically in the symmetry resolution, also for the critical compact boson [30,37]. Let us now critically discuss our results. First of all, there is a very important difference compared to the conformal gapless case [30], i.e. the absence of equipartition of entanglement [34]: the Rényi entropies (63) depend explicitly on q. This dependence is explicitly reported in Figure 2 (a) where, in order to show its variation, we plot it as a continuous function of q, although only integer values are physical. The lack of exact equipartition is not surprising; also in critical models the leading terms for large show equipartition [34], while some subleading terms depend explicitly on q [30,37]. In panel (b) of Figure 2 we focus on the critical limit of Rényi entropies (66) plotting S n (q) − S n (q = 0). As ε → 0, the result approaches the critical form (66), but clearly the convergence is not uniform: it is faster for smaller q and n. Indeed, since this dependence is all encoded in the function Φ q (e −nε ), the parameter that must be small is not ε, but nε. On the other hand, the higher order terms in ε, that have been neglected in (66), become important for large q.
Another interesting feature of the symmetry resolved entropies for this complex harmonic chain is an effective equipartition in two limits. The first one is the limit of large q. Indeed, in Eq. (63) the entire q-dependence is encoded in the function Φ q (e −nε ). Looking at Eq. (61), it should be clear that all the terms with qnε 1 are exponentially suppressed. Practically, the total sum is more or less the same for all q such that nεq 1 (from Eq. (50) this is equivalent to nqπ 2 log ξ in the critical region). Hence, there is an effective equipartition among all q 1/(nε). Actually, since the only physical values of q are the integers, this fact implies that there is an almost exact equipartition (with the exception of S n (0)) of the entropy if nε 1, which corresponds to ω 0 10 −4 (for n = 1). In panel (c) we report the von Neumann entropies S 1 (q) for several values of ω 0 , showing that, as q becomes large enough, the entropies S n (q) do not depend on q anymore. We also explicitly report the (approximate) crossover values for q ∼ 1/ε (as function of ω 0 is given by Eq. (24)), showing that it correctly captures the change of behaviour. Finally, we have effective equipartition also in the critical regime, but in this case also for small q. In fact, Eq. (66) shows that the q-dependent term is proportional to ε 2 , while the leading term of S n (q) (say for q = 0) diverges as ε −1 . Thus the q-dependence is suppressed as ε 3 and there is an effective equipartition. Even if for large q, the expansion (66) breaks down, we do not expect that S n (q) − S n (0) becomes of the order S n (0) and so there is an effective equipartition for all q: the numerical analysis of Eq. (63) seems to confirm this expectation. The functional form of the leading q-dependent term in Eq. (66) is reminiscent of the one found for free fermions [37].

The total entanglement entropy as a consistency check.
As a non-trivial consistency check of our results, we compute the total von Neumann entanglement entropy starting from the symmetry resolved ones using Eq. (5). The probability p(q) is given by Eq. (62) with n = 1 and Eq. (64) provides the symmetry resolved entropies. Plugging these two results into Eq. (5) leads to The last sum over q above can be written as the following derivative where we have used that reflecting that Z 1 (q) is normalised to 1. Taking now the derivative with respect to ε, we finally obtain which is the entanglement entropy of a complex harmonic chain (i.e. the double of a real one).

Numerical checks
In this subsection we test the validity of the results in the previous ones against exact numerical computations. We work only with an infinite real harmonic chain (22) with finite ω 0 . For the complex case, we just combine the results for two real chains. Let us consider a bipartition where the subsystem A is given by contiguous lattice sites. Let us call X A and P A the × matrices of the correlators restricted to the subsystem A, where X i j = 〈q i q j 〉 and P i j = 〈p i p j 〉.
The explicit forms of these correlators in the ground state of the gapped harmonic chain have been already reported many times in the literature (see e.g. Refs. [3,62,67]) and we are not going to rewrite them here. Let us denote by σ k , with k = 1, . . . , , the eigenvalues of the matrix X A P A . We introduce the vectors |n〉 ≡ k=1 |n k 〉, products of Fock states of the number operator in the subsystem A, namely N A . The reduced density matrix of A can be written as [68,69] ρ A = n k=1 where the non-negative integer n k is the k-th element of the -dimensional vector n. Since N A = j∈A n j is the number operator in the orthonormal basis made of the states |n〉, we can write Summing over the possible occupation numbers n k from 0 to ∞, we get This relation holds also in higher dimensions and for a generic shape of the subsystem A provided that is the number of sites in A. Notice the similarity of Eq. (74) We now consider F n (α) = log Tr[ρ n A e iN A α ] for a real harmonic chan. The numerical data for F n (α) for an interval of length should converge to the double (because of the two end-points) of the CTM prediction for the semi-infinite line (with one-endpoint) as soon as becomes larger than the correlation length ξ. In Figure 3 we report the numerical data for (half of) the real and the imaginary parts of F n (α) for different values of n and α. We have set ω 0 = 0.1, so that after a short crossover in , the data saturate. The CTM prediction (37) is also reported for comparison, showing that the analytical result perfectly describes the saturation values. The charged moments for the complex harmonic chain are just log Z n (α) = 2Re[F n (α)] both for numerics and analytics and so Figure 3 is a direct test also for them. We now take the Fourier transform of the numerical data for Z n (α) to test the validity and the accuracy of the CTM predictions for the symmetry resolved moments and entropies. In Figure 4 we report the (square roots of the) numerically calculated symmetry resolved partition sums Z n (q). We compare the data for n = 1, 2, 3 with the CTM prediction (59). The latter perfectly captures the q-dependence, as shown in the panel (a), and gives the value at which the data saturate when studied as functions of , panel (b). Finally, in Figure 5 we report the symmetry resolved entropies for several values of q, n, ω 0 . For large , the numerical data converge to (twice) the CTM predictions in Eqs. (63) and (64). Notice that for the larger values of ω 0 the saturation values do not depend on q because of the effective equipartition, but for smaller ω 0 they clearly do. As ω 0 becomes much smaller (such that ε ∼ 0.1), we expect again effective equipartition, although we do not report such data here because they require very large .

Gapped XXZ spin-chain
In this section we study the symmetry resolved entanglement in the anisotropic Heisenberg model in the gapped antiferromagnetic regime using the CTM approach. The resolved moments are computed starting from the explicit expressions for the eigenvalues of the RDM and their degeneracies. Then the symmetry resolved entropies are deduced and their critical regime is investigated. The discrete Fourier transform of the resolved moments allows us to compute the charged moments and to discuss their behaviour in the critical regime.

Symmetry resolved moments and entropies
The Hamiltonian of the anisotropic Heisenberg model (also known as XXZ chain) is  (63) and (64), to which they clearly approach. Notice that the convergence is slower for smaller ω 0 . For ω 0 = 0.1 we have an approximate equipartition, but this is not the case for ω 0 = 0.01.
where σ i , i = x, y, z are the Pauli matrices. The model has a conformal quantum critical point for ∆ = 1, it is gapless when |∆| ≤ 1 and gapped when |∆| > 1. We consider this model in the antiferromagnetic gapped regime with ∆ > 1.
The XXZ chain is solvable by Bethe Ansatz techniques; unfortunately this framework is not very effective to study the entanglement properties both in the coordinate [9] and in the algebraic [70][71][72][73][74][75][76][77][78] approach. On the other hand, the CTM solution for the XXZ chain is a powerful tool to compute the entanglement entropies; in this approach, the reduced density matrix is related to the partition function of the six-vertex model on a strip with a cut. In Ref. [50] H CTM has been found to be of the form (13) with and n j being some fermionic number operators. Since in the thermodynamic limit, the groundstate of the gapped XXZ spin-chain is doubly degenerate we should clarify which state we are going to deal with in this section. The entanglement Hamiltonian (13) together with (77) selects by construction the ground state that does not break the inversion symmetry, i.e. the one that in the limit of large ∆ is (|N 1 〉 + |N 2 〉)/ 2 where |N i 〉 are the two possible Néel states. However, we prefer to work with the more physical symmetry breaking state |N i 〉. In CTM approach this can be constructed with an entanglement Hamiltonian of the form (13) where the sum over j starts from 1 rather than 0, i.e.
In the remaining part of this section we always focus on the symmetry breaking ground state with the above H CTM . If one is interested into the other state, analogous results may easily be derived. The entanglement spectrum is obtained by filling in all the possible ways the single particle levels in (78) (i.e. setting all n j equal either to 0 or 1). The resulting levels are equally spaced with spacing 2ε and highly degenerate. The degeneracy of the level 2εs, with s = j j (see (77)) is Q(s), the number of partitions of s into smaller non-repeated integers (including zero).
(Notice we use the non-standard symbol Q(s) instead of q(s) to avoid confusion with q, the charge sector. ) We want to characterise how the entanglement of the semi-infinite line A with respect to its complement splits into the different sectors with fixed magnetisation S z ≡ j σ z j /2. We indicate with q the possible values, in the subsystem A, of the difference of the magnetisation with respect to the antiferromagnetic Néel state chosen as a reference configuration. Such variable q is quantised in terms of integer numbers (each spin flip leads to a change of magnetisation of ±1), i.e. q ∈ . With a slight abuse of language, we will refer to q as the magnetisation, although it is a magnetisation difference. To derive the symmetry resolved entanglement, we first write Z n (q), defined in (8), as where λ s are the eigenvalues of the RDM and the sum is restricted to the levels with fixed value of q. Using Eq. (12) and the explicit expression of the entanglement spectrum from Eq. (78), we can write where F(q, s) is the number of eigenvalues at level s with magnetisation q. The degeneracies F(q, s) have been studied in Ref.
where we have also exploited that P(n) is non vanishing only if n is a positive integer.
The two sums in (81) can be conveniently rewritten in terms of generating functions Setting x = e −4nε and y = e −2ε in (82) and plugging them into (81) we obtain We remark that Z 1 (q) is normalised to one, i.e. q∈ Z 1 (q) = 1, as it should be from the definition (79). This is consistent with the interpretation of Z 1 (q) as a probability, see Section 2. The denominator of Eq. (83) can be expressed in terms of elliptic theta functions (see Appendix C) and then Z n (q) reads where κ and κ are defined in (124). Notice that q = 1/4 is exactly the mean magnetisation of the subsystem in the critical limit ε → 0, as we can check by computingq = dqqZ 1 (q), since we are dealing with the symmetry breaking ground state. Notice that the dependence on q in Eq. (84) is entirely encoded in the Gaussian factor and it is symmetric for q → 1/2 − q. Moreover, exploiting the asymptotic behaviours in (130) and (131) in appendix C, we have that in the critical regime Z n (q) becomes where we keep the Gaussian factor in order to have a meaningful result. Once the resolved moments Z n (q) have been worked out, the symmetry resolved entropies follow straightforwardly and, taking the limit n → 1, Notice that as ∆ 1, S n (q) → 0 (see also Figure 6), since in this limit the selected antiferromagnetic ground state is a product state. If we would have considered the non-symmetry breaking ground state (|N 1 〉 + |N 2 〉)/ 2, ∆ 1 we would have found S n (q) → log 2, as for the total entropy [5,14,23]. We stress that although there is entanglement equipartition, the functions S n (q) are not equal to the total entropies S n because there is a non-vanishing fluctuation term like in Eq. (5) for n = 1.
Remarkably, the expressions (86) and (87) for the symmetry resolved Rényi and von Neumann entanglement entropies do not depend on q for any value of n, i.e. they exactly satisfy the equipartition of entanglement for any value of ∆. In the critical case, only the leading terms satisfy such equipartition [34,37].
The relation between the correlation length of the model and ε, in the critical regime ξ 1, is [12] log ξ π 2 2ε which combined with Eqs. (85) provides the expansions of the symmetry resolved entropies in the critical regime We notice that the term − 1 2 − 1 2 log(log ξ/π) appearing in S 1 (q) in Eq. (89) is canceled exactly by the fluctuation entanglement entropy once we consider the total von Neumann entanglement entropy. Indeed, using that the probability is p(q) = Z 1 (q), we write the fluctuation entropy as − dqZ 1 (q) log Z 1 (q). Using (84), computing the gaussian integral in q and then taking the critical limit, we find In the middle panel, we report again Z n (q) at fixed q and as function of ∆ (full lines). As a comparison, we also report the asymptotic expansion (85) for ∆ close to 1 (dashed lines). In the right panel, we report S n (q) and its critical limit, respectively Eq. (86) and Eq. (89), as function of ∆ for n = 1, 2, 3. We recall that S n (q) does not depend on q because of entanglement equipartition.
which exactly cancels the contribution from the configurational entropy. This is in complete analogy with what has been found for critical systems for the log log term [37].
As for the harmonic chain, another useful check is to recover the total von Neumann entanglement entropy from S 1 (q) in Eq. (87). Using the expression of Z 1 (q) = p(q) in Eq. (83) once we set n = 1, the total von Neumann entropy is Let us introduce the constants ∞ k=1 1 − e −4εk = N 1 and ∞ k=1 1 + e −2εk = N 2 . Because of normalisation of Z 1 (q), the first term in Eq. (91) just gives S 1 (q) (since, as already stressed, it does not depend on q), while the second one leads to Performing explicitly the derivative with respect to ε and summing all contributions in Eq. (91), we obtain which is the known entanglement entropy found in Refs. [5,23] for the symmetry breaking ground state.
In Figure 6 we report symmetry resolved moments and entropies. The possible values of q are just integers, but since Z n (q) becomes quickly small as q increases, we consider arbitrary real values. As anticipated, Z n (q) has a peak at q = 1/4 and shows a clear Gaussian shape for all ∆. The exact result (83) is well approximated by its critical limit (85) for ∆ close to 1, but the approach is not uniform and it is worse for larger q (as well as larger n). Clearly, the maximum of Z n (q) is a decreasing function of n. In the last panel of Figure 6, we report the symmetry resolved entropies as functions of ∆ (as we stressed because of equipartition, they do not depend on q). Notice that the window of ∆ for which the critical limit of S n (q) in Eq. (89) is a good approximation of the exact expression (86) is wider for smaller values of q.

Charged moments via Fourier series
The charged moments are obtained from the resolved ones Z n (q) by inverting the formula (8), i.e.
Plugging in the above equation the result for Z n (q) in Eq. (83) and using the definition of the elliptic function θ 3 (z|u) (see Eq. (121) in appendix C), we obtain Setting α = 0 and exploiting the infinite product representation (129) of θ 3 (z|u), we get as found in [5]. As for Z n (q) in Section 4.1, we can express Z n (α) in terms of elliptic functions obtaining Z n (α) in the critical regime is obtained using the asymptotic expansions reported in appendix C (i.e. Eqs. (134), (130) and (131)), finding Taking the logarithm of Z n (α) and using (88) we have log Z n (α) 1 12 Here, the linear term in α is just the mean magnetisation in A,q = 1/4. The leading term in Eq. (99) is very suggestive. Indeed, for the critical compact boson (aka, Luttinger liquid), in the case of A being an interval of length embedded in an infinite 1D system, log Z n (α) diverges logarithmically with as [30] log Z n (α) 1 6  where K is the Luttinger liquid parameter (related to compactification radius). The prefactor of Eq. (99) is exactly half of the conformal result (100) for K = 1 2 , which is the Luttinger parameter at ∆ = 1. The multiplicative factor 1/2 is simply understood because in our geometry there is a single endpoint instead of two as in the conformal case. It is natural to wonder under what hypotheses this can happen since we have seen that it is not true for the harmonic chain. Moreover, for the symmetry resolved entropies, the CFT result is S n (q) − S n = − 1 2 log((2K/π) log ) + O( 0 ) [30,34], which is the same as in Eq. (89) with the replacement → ξ and with K = 1/2.
In Figure 7 we report the plots of the charged moments as functions of α and ∆. Also in this case, the approach to the critical regime is not uniform and it is faster for α closer to 0 (and n close to 1). This is very different compared to what we have seen in the previous section for the harmonic chain for which the limit α → 0 is singular. This is a further confirmation that the anomalous behaviour of the harmonic chain is due to its non-compact nature of the continuum limit.

Conclusions
In this manuscript we found exact results for the symmetry resolved entanglement entropies of half line in infinite integrable systems in the gapped regime. We considered two models for which the RDM (and therefore the entanglement spectrum) of the subsystem can be obtained through the Baxter CTM.
In Section 3 we considered the massive regime of the complex harmonic chain that has a U(1) symmetry corresponding to the conservation of the electric charge. In order to obtain the symmetry resolved entanglement entropies, we first computed the charged moments of the RDM in Eqs. (35), (37), and (48). Their critical behaviour is also discussed and an interesting discontinuity for α → 0 has been pointed out. Then we computed the Fourier transform of the charged moments and the symmetry resolved entanglement entropies (see Eqs. (63) and (64)); we also discussed their leading behaviour in the critical regime, see (66). Interestingly we found that there is no entanglement equipartition, i.e. the symmetry resolved entropies S n (q) explicitly depend on q. However, entanglement equipartition is effectively recovered in two limits: i) for large q, i.e. as soon as q becomes larger than the logarithm of the correlation length and ii) in the critical region for nε 1.
We also derived an exact expression for the charged moments valid for a generic harmonic system in the correlation matrix approach [62]. The final results are the formulas (74) and (75)) which hold in any dimension and for any shape of the subsystem. Here we limit ourselves to use these relations to check numerically the results derived in the CTM approach. We considered a finite interval of length in an infinite chain and we found that for large the results converge to the CTM predictions.
In Section 4 the symmetry resolved entanglement entropies have been computed for the XXZ chain in the antiferromagnetic gapped regime (Eqs. (86) and (87)). Here, the conserved U(1) symmetry corresponds to the rotations in the plane perpendicular to the anisotropy. Somehow surprisingly, for this model, the symmetry resolved entropies exactly satisfy the equipartition of entanglement for any anisotropy ∆ ≥ 1. We found this result very remarkable, although its physical origin is not clear: it would be very interesting to establish a priori which properties guarantee an exact equipartition of entanglement and how they are related to integrability. The computation has been performed exploiting the explicit expressions of the elements of the entanglement spectrum and the degeneracies of each level in a given magnetisation sector [58]. We also computed the charged moments (Eq. (95)) checking that for α = 0 the result of [5,23] was retrieved. We found that Z n (α) have no discontinuities, as a difference with the complex harmonic chain.
Let us conclude this manuscript with some possible directions for future investigations motivated by the results we have found. A first and natural question is whether some of the results we found here may be also recovered in massive two-dimensional field theories both free and integrable. Work in this direction is in progress [64]. It is also interesting to understand what happens when integrability is absent: while a general treatment seems impossible, the results for the entanglement spectrum in Refs. [58,79] suggests that in some non-trivial regimes general results may be derived. Another natural extension is to study symmetry resolved entanglement in higher dimension for which there are only few works for free fermions [38,39]. Our Eq. (74) paves the way for general numerical studies in arbitrary dimension for bosonic systems as well, also in the presence of a spherical constraint [80]. In some cases, also analytical results can be explicitly worked out [81]. Finally, one expects that the symmetry resolved entanglement should help in reconstructing the entanglement (or modular) Hamiltonian, but it is still unclear how. This issue is very timely given the large current effort devoted to understand the structure of the entanglement Hamiltonians both in field theories [82][83][84][85] and lattice models [86][87][88][89][90].

A.1 A two-site chain with complex oscillators
For a single harmonic chain with two sites, the RDM has been worked out e.g. in [57]. The entanglement Hamiltonian of one site is For this site, the β's are related to the a's as which is the specialisation of Eq. (26) to the case of A being one site. Here e θ = (1+ω 2 0 /4) 1/4 , but its explicit value is unimportant. Hence, in terms of the ladder operators a, a † , H A can be rewritten as Rather then one real harmonic oscillator, we consider a complex one, which is the same as two real harmonic oscillators described by the ladder operators a † i , a i , i = x, y such that the only non-vanishing commutators are [a i , a † i ] = 1, i = x, y. Therefore, the entanglement Hamiltonian of these two real harmonic oscillators is simply the sum of two single ones: Let us rewrite Eq. (104) in terms of the particle and antiparticle ladder operators a and b in Eq. (27), i.e.
or, equivalently (up to an additive constant we can absorb in the normalisation factor of the RDM) One can bring Eq. (107) into a diagonal form through Bogoliubov transformations, i.e.

A.2 The Bogoliubov transformation for a chain of arbitrary length
For a real harmonic chain of arbitrary length 2L, the entanglement Hamiltonian for half system is [57] where the eigenvalues j depend on L and in the thermodynamic limit are given by Eq. (24) while for L = 1 by Eq. (101).
The ladder operator β j as function of the local ladder operators are given by Eq. (26), i.e.
Hence, the entanglement Hamiltonian in terms of local operators is Therefore, the entanglement Hamiltonian of a complex chain is just the sum of two real ones with local ladder operators a a, j with a = x, y as in the case of two oscillators in the previous subsection. Such H A can be rewritten in terms of the particle and antiparticle ladder operators in Eq. (27), obtaining (up to constants) which we can put in the diagonal form by the transformation (28).

Appendix B A generalisation of the binomial theorem
In this Appendix we report a proof (based on Refs. [91,92]) of a generalisation of the binomial theorem that has been used in Eq. (135). We also discuss some corollaries of the theorem used in the main text. The generalisation of the binomial theorem is: where n k t is the generating function (in the variable t) for the number of integer partitions with at most k parts, whose largest part is at most n − k, i.e.
We give a combinatorial proof of Eq. (116). Take the left hand side of Eq. (116) and think of it as a polynomial in x (of degree n) with coefficients being polynomials in t, i.e. rewrite it as n k=0 a k (t)x k . Clearly, a k (t) is the generating function for partitions with exactly k parts not exceeding n. In fact, expanding the product on the left hand side, the term x k comes from taking x t j exactly in k factors. In each of them, x comes together with some power of t, which is different for each factor and does not exceed n; hence they are parts of our partition. These partitions can be thought as Young tableaux with k rows and at most n columns. Choosing a given partition, denote as λ j the length of the row j (starting from the bottom). We then have 1 ≤ λ 1 < λ 2 . . . λ k−1 < λ k ≤ n. From this partition, we can produce another one with k rows and at most n − k columns. Just proceed as follows: remove zero boxes from the first row, one box from the second row and, in general, i −1 boxes from the i-th row. So we obtain a partition of µ's, The generating function for µ's is exactly n k t . On the other hand, the generating function for λ's is obtained from the generating function on µ's by multiplying it by t k(k−1)/2 , which takes into account the total number of removed boxes. Therefore we have Another useful property derived from this theorem is the identity [91] n−1 that we used to derive Eq. (135). A final observation is that, through this binomial theorem, one can prove that e −εnq(q−1) / q−1 j=0 (1 − e −2εn( j+1) ) is the generating function of the partitions of an integer into q distinct parts. This result will be useful in appendix D.

Appendix C Some properties of the Jacobi theta functions
In this Appendix we report and discuss some properties of the Jacobi theta functions that we exploited to get some results in the main text.
The Jacobi theta functions θ r (z|u), r = 2, 3, 4 are defined as [65] and we use the standard shorthand θ r (u) ≡ θ r (0|u), r = 2, 3, 4. The functions θ r (u), r = 2, 3, 4 can be expressed in terms of infinite products [65] θ 2 (u) = 2u These three relations allow us to write some particular infinite products in terms of ratios of Jacobi theta functions. An example of such relations is where we defined that can be obtained properly combining the equations in (122). Other formulas that can be derived in this way are and ∞

C.1 Some asymptotic properties of the Jacobi theta functions
In this subsection we report some asymptotic expressions of θ r (z|u), r = 2, 3, 4 in the limit in which the variable u → 1. These formulas are useful to derive results in the critical regime, namely for ε → 0. Let us consider first the case in which the variable z in the theta functions is 0. At the leading order when u → 1, we can write [93] From the definition (124) we therefore obtain at the leading order Two examples in which these asymptotic formulas have been employed in the main text are, that has been exploited to obtain (141), and involved in the computation of the critical limit of (84). The asymptotic expression for u → 1 of θ 3 (z|u) is [93] θ 3 (z|u) π which reduces to the second identity of (130) when z = 0. Plugging (134), setting z = −inε and u = e −4nε , into Eq. (97) we obtain (98).

Appendix D The CTM symmetry resolution
In the main text of the paper, we derived the symmetry resolved entropies for the most interesting case of the conserved charges Q A being the "electrical" charge of the complex harmonic chain and magnetisation of the XXZ chain (equivalently the number operator in fermion language). Being these models integrable, there are many other conservation laws that can be used in place of these, but usually are very difficult to calculate. However, a quantity we can easily deal with in the CTM approach is Q A = j n j = j β † j β j , although it has not a clear physical meaning, if it has one at all. Indeed, since [ρ A , n j ] = 0 for each j, Q A is conserved and the symmetry resolved entanglement for the sectors with different values of this quantity may be studied. We will refer to Q A = j n j as the CTM charge. Although these results have most likely no physical meaning at all, the details of the calculations are rather interesting and worth being presented.

D.1 The CTM symmetry resolution in the harmonic chain
For a single real harmonic chain, the flux resolved partition sum for the CTM charge is just Z n (α) = e F n (α) . Before performing the Fourier transform to get the symmetry resolved mo-ments, it is useful to rewrite e F n (α) as where in the last equality we have used the generalisation of the binomial theorem reported in Appendix B. In addition, Eq. (123) allows us to rewrite the denominator in Eq. (135) in such a way that the Fourier transform Z n (q) is (136) Since q and k are both integer numbers, Eq. (136) simplifies to We also provide the analytic continuation of Eq. (136) to real q where we expressed the finite product in terms of the infinite products: and we introduced the generalised gamma function Eq. (140) reduces to the ordinary gamma function in the limit ε → 0.
In the critical regime ε → 0, as showed in Appendix C, Eq. (137) becomes The symmetry resolved Rényi entropies are easily deduced from Eq. (9), obtaining Taking the limit n → 1, we get the von Neumann entropy The analytic continuations of S n (q) and S 1 (q) to real q are respectively with the leading behaviour for ε → 0 given by where we introduced the functions a(q) = − log Γ (q + 1) and b(q) = a(q) + q. The symmetry resolved entropies do not satisfy entanglement equipartition, like the one for the true charge of the complex chain. However, the breaking of equipartition is rather different: in this case the leading term for ε → 0 which grows linearly in q and is proportional to log ε, while for the complex chain the first term breaking equipartition is subleading and goes like ε 2 (the sums for the entanglement entropies are finite because the probabilities decay fast with q, cf. Eq. (141) for n = 1). Anyhow, from the expressions as sums over q in Eqs. (142) and (143), it is clear that all the terms with 2ε j 1 are exponentially suppressed. Practically, the total sum is more or less the same for all q such that εq 1 (from Eq. (50) this is equivalent to qπ 2 log ξ in the critical region). Hence, there is an effective equipartition among all q 1/ε. Actually, since the only physical values of q are the integers, this fact implies that there is an almost exact equipartition (with the trivial exception of S n (0) = 0) of the entropy if ε 1, which corresponds to ω 0 10 −4 . Some results for the symmetry resolved moments and entropies are reported in Figure 8 as continuous functions of real q, although only the integer values are physical. It is evident from the figure that, as q becomes large enough, the entropies S n (q) do not depend on q anymore, as from the previous argument about effective equipartition. In panel (c) we explicitly report the (approximate) crossover values for q ∼ 1/ε (as function of ω 0 is given by Eq. (24)), showing that it correctly captures the change of behaviour. In panels (a) and (b) we report S 1 (q) and S n (q) respectively, together with the critical limit (146). As expected, the approach to the critical behaviour is highly non-uniform in q: as q becomes larger we need smaller values of ε.
As a final non-trivial consistency check of our results we compute the total von Neumann entanglement entropy starting from the symmetry resolved ones using Eq. (5). The probability p(q) is given by Eq. (137) with n = 1 while the symmetry resolved entropies are in Eq. (143). Plugging these two results into Eq. (5) leads to The sum over q in (147) can be written as the following derivative .
Using (148) in (147) and exchanging the derivative with respect to ε with the sum over q, we can exploit that reflecting that Z 1 (q) is normalised to 1. Taking now the derivative with respect to ε, we finally obtain which is the known result from the CTM calculation in Ref. [5], i.e. Eq. (37) for α = 0.

D.2 The CTM charge for the XXZ spin chain D.2.1 Charged CTM moments.
We now consider the CTM charge in the XXZ spin chain. As a difference compared to the main text, in this appendix we focus on the state that does not break the symmetry, i.e. with entanglement Hamiltonian given by Eq. (13) with the sum over j starting from 0. As usual, we first compute the charged moments Z n (α): where we used that for this model the n j 's are fermionic number operators. Taking The asymptotic expansion for small ε is obtained applying the Poisson resummation formula (41). Defining f n,α (x) as f n,α (x) = log(1 + e −2nx+iα ), where we used f n,α (0) = log(e iα + 1)/2. The cosine-Fourier transform (42) of (153) iŝ with the function Φ defined in (46). For α = 0 and n = 1, it reduces tô We now apply to (154) the Poisson resummation formula (41)  For ε → 0, the leftover sum over k is vanishing. In particular, the last part behaves as Thus, we get log Z n (α) = − Li 2 (−e iα ) 2εn − nπ 2 24ε + log e iα + 1 − n log 2 It is worth to observe that, for this model, the limit α → 0 can be taken after the expansion close to ε = 0 retrieving the result found in [14] log Z n = 1 n − n π 2 24ε + (1 − n) log 2 2 + O(ε).
In Figure 9, we report the α dependence of the charged moments for different values of ∆ and n. We also shows the comparison between the exact result (157) and its critical limit (159). As expected, the latter gets very close to the former as ∆, therefore ε, is close to its critical value.  .
We can check Eqs. (163) and (164) computing Z n (q) directly from the entanglement spectrum, as done in Sec. 4.1 for the case of Q A being the magnetisation. In the symmetry sector with charge q, the degeneracy of the level 2εs is P q (s), i.e. the number of partitions of an integer s in exactly q parts, not exceeding s. The partition function Z n (q) then is Z n (q) = These symmetry resolved entropies have the same form as the ones for the harmonic chain in Eq. (142) except for the explicit expression of ε. Thus, the analytic continuation to real q, the von Neumann limit, and the behaviour in the critical regime are the same as those obtained in the previous section and we do not report here. Notice that these symmetry resolved entropy do not satisfy entanglement equipartition. However, as for the harmonic chain, equipartition is effectively recovered as q 1/ε. Finally, notice the similarity between these symmetry resolved entropies and the ones for the magnetisation in Eq. (86). Apart from a reparametrisation and an additive term, the main difference is that in the case of the CTM charge the sum is up to q and in the magnetisation case it is up to ∞ (and that is why the former does not satisfy equipartition while the latter does). When the upper limits in the former do not matter, the two become practically equivalent.
In Figure 10 we plot Z n (q), showing also a comparison between the exact result (163) and its critical limit, Eq. (164). It is interesting to observe that the maxima of Z n (q) are increasing or decreasing with n depending on the considered values of ∆. In the last panel we report the exact expression of S n (q) and its critical limit, respectively Eq. (144) and Eq. (146), as function of ∆ for different q. The agreement improves for ∆ close to 1, as it should.  Figure 10: Symmetry resolved moments and entropies for the CTM number in the XXZ spin chain. The top two panels report the exact results for Z n (q) (163)-full lines-for ∆ close to 1 and the comparison with the critical limit, Eq. (164)-dashed lines-, as a function of q, for different values of n = 1, 2, 3. At fixed ∆, the approach is not uniform and the smaller values of q converges faster. In the third panel, we report Z n (q) for ∆ far from the critical point, where a peak at q = 1/2 is developed. In the last panel, the exact and critical limit of S n (q), respectively Eqs. (144) and (146), are shown against ∆ for different q.