Symmetry resolved entanglement in integrable field theories via form factor bootstrap

We consider the form factor bootstrap approach of integrable field theories to derive matrix elements of composite branch-point twist fields associated with symmetry resolved entanglement entropies. The bootstrap equations are determined in an intuitive way and their solution is presented for the massive Ising field theory and for the genuinely interacting sinh-Gordon model, both possessing a $\mathbb{Z}_{2}$ symmetry. The solutions are carefully cross-checked by performing various limits and by the application of the $\Delta$-theorem. The issue of symmetry resolution for discrete symmetries is also discussed. We show that entanglement equipartition is generically expected and we identify the first subleading term (in the UV cutoff) breaking it. We also present the complete computation of the symmetry resolved von Neumann entropy for an interval in the ground state of the paramagnetic phase of the Ising model. In particular, we compute the universal functions entering in the charged and symmetry resolved entanglement.


Introduction
Symmetries play a central role in physics and in our understanding of nature. They are important guiding principle when formulating theories, their presence or absence or their breaking have profound consequences on the physical properties of models and real-world systems; last but not least symmetries often provide a larger view in the description of the systems of interest. From a practical perspective, the presence of a symmetry usually leads to some kind of simplifications. In particular, for a quantum system the operator corresponding to the symmetry commutes with the Hamiltonian and hence the two operators have common eigenvectors or, in other words, the eigenstates of the system can be characterised by quantum numbers associated with the symmetry operation. The idea of exploiting the additional structures imposed by symmetry for various physical objects is very fruitful and has been recently extended to the study of entanglement too.
When a system is in a pure state, the bipartite entanglement of a subsystem A may be quantified by the von Neumann entanglement entropy [1][2][3][4]. Denoting the reduced density matrix (RDM) of the subsystem by ρ A , the entanglement entropy is defined as (1.1) Alternatively the Rényi entanglement entropies also provide bipartite entanglement measures in pure states and are related to the von Neumann one by taking the limit n → 1.
The explicit idea of considering generally the internal structure if entanglement associated with symmetry is rather recent [5][6][7][8]. In a symmetric state, the system's density matrix ρ commutes with the conserved chargeQ corresponding to the symmetry; if in additionQ A , the restriction ofQ to this subsystem, satisfies then the RDM ρ A is block-diagonal with respect to the eigenspaces ofQ A and, consequently, the Rényi and von Neumann entropies can be decomposed according to the symmetry sectors. Let us denote with P(q A ) the projectors onto the eigenspace with eigenvalue q A . The symmetry resolved partition functions can be defined as Z n (q A ) = Tr (ρ n A P(q A )) , (1.4) from which the symmetry resolved Rényi entropies S n (q A ) and the symmetry resolved von Neumann entropy S(q A ) can be naturally obtained as , and S(q A ) = − ∂ ∂n respectively. This way the total von Neumann entropy can be written as [9] where p(q A ) = Z 1 (q A ) is the probability of finding q A as the outcome of a measurement ofQ A . The contribution S c denotes the configurational entanglement entropy, which measures the total entropy due to each charge sector (weighted with their probability) [7,10] and S f denotes the fluctuation (or number) entanglement entropy, which instead takes into account the entropy due to the fluctuations of the value of the charge in the subsystem A [7,11,12].
The calculation of the symmetry resolved partition functions and entropies is generally a difficult task; the usual way one proceeds includes the replica method and the computation of the charged moments [6] Z n (α) = Tr ρ n A e iαQ A . (1.7) Considering quantum field theories (QFTs) a natural way of computing the Rényi entropies for integer n is provided by the path-integral formalism: Trρ n A corresponds to the partition function on an n-sheeted Riemann surface R n , which is obtained by joining cyclically the n sheets along the region A [13][14][15]. It was recognised in [6] that the charged moments (1.7) correspond, in the path integral language, to introducing an Aharonov-Bohm flux on one of the sheets of R n . An intuitive picture is given by imagining particles with a specific charge eigenvalue moving from one level of R n to the other until they return to their original sheet [6]; if the charge within the subsystem is q A , the total acquired phase of a given particle is then e iαq A as given by the term e iαQ A in Eq. (1.7).
In a path integral approach to quantum field theories (QFTs), the computation of either Trρ n A or Trρ n A e iαQ A can equivalently proceed for an n-copy QFT, where specific boundary conditions are prescribed for the fields φ 1 , ..., φ n corresponding to the different copies. Crucially, in 1+1 dimensional relativistic QFTs, there exist local fields in the n-copy theory that correspond to the boundary conditions imposed on the fundamental fields in the path integral. These fields have been dubbed branch-point twist fields [14,37]. The nth Rényi entropy of an arbitrary spatial subsystem (i.e. consisting also of disjoint intervals) is equivalent to a multi-point function of the branch-point twist fields in an n-copy theory. Direct access to these fields is established in 2D CFT, where the scaling dimensions of these fields are exactly known [14,38,39]. These dimensions directly provide the scaling of two-points function, corresponding to a single interval for a generic CFT [14]. The behaviour of four-point [40][41][42][43][44][45] and also higher functions [46] of these twist fields are known for special CFTs.
The main subject of this manuscript is however integrable quantum field theories (IQFTs). In these theories, the form factor (FF) bootstrap allows for the calculation of the matrix elements of the twist field [37,47,48]. Via the bootstrap, in principle, all matrix elements can be computed. However, the correlation functions of the fields at large distances are usually well described by the first few members of the form factor series. Such form factor bootstrap program has been used in IQFTs for the calculation of the entanglement entropy in many different situations [49][50][51][52][53][54][55][56][57][58][59].
The symmetry resolved entropies in CFT can be obtained by composite branch-point twist fields in essentially the same way as the conventional entropies [6]. The only price to pay is the introduction of composite twist fields fusing the action of the replicas and of the flux of charge (see below for the precise definition). These new composite twist fields have been identified for Luttinger liquids [6], for the SU (2) k Wess-Zumino-Witten models [6], and for the Ising and Z N parafermion CFT [21]. Furthermore, the existence and applicability of such composite twist fields have been recently demonstrated for the free massive Dirac and complex boson QFT too [19]. These findings suggest that in perturbed QFTs (corresponding to a relevant perturbation of a given CFT), the off critical version of the composite twist field exists. We expect that in IQFTs their form factors can be determined with the bootstrap program, similarly to the usual twist fields [37,47,48].
This paper aims to initiate such a program for interacting IQFTs. In particular, we introduce and discuss appropriate bootstrap equations for the composite branch-point twist fields, find their first few solutions and compute the long-distance leading behaviour of the symmetry resolved entropies (similar twist fields have been introduced for non-unitary QFT [53], but in a completely different context and with different aims). For the sake of simplicity, here we consider the simplest integrable models, namely the Ising field theory, which is equivalent to a free Majorana fermion QFT, and the sinh-Gordon (ShG) model, which is a truly interacting QFT. Both models possess the discrete Z 2 symmetry. While from the point of view of IQFT techniques these models are indeed the simplest possible ones, the resolution of their entanglement in terms of the Z 2 symmetry requires a careful treatment because of the lack of a conserved density (1.3). Integrable QFTs with continuous symmetry present many more technicalities because of their richer particle content and for the presence of non-diagonal scattering. Their analysis is still on the way and will be eventually the subject of subsequent works.
The structure of this paper is as follows. In section 2 the FF approach for conventional branchpoint twist fields is briefly reviewed, focusing on the bootstrap equations and their solution for the Ising and ShG models. In section 3, we show how the bootstrap equations can be modified to obtain solutions for the modified twist fields corresponding to a given symmetry resolution. For the Ising and ShG models, the two-particle FFs of the Z 2 twist fields are determined as well. Sections 4 and 5 are explicitly focused on Ising and ShG models respectively, reporting also ∆-theorem [60] checks of the obtained form factors; for the Ising model the even particle-number FFs are expressed in terms of a Pfaffian involving the two-particle matrix elements. Section 6 reports general results for Z 2 symmetry resolved entropies that can be deduced from the IQFT structure. The leading and sub-leading contributions of the symmetry resolved entanglement are explicitly calculated in section 7 for the paramagnetic ground state of the Ising model. We conclude in section 8, which is followed by the appendices containing the determination of the vacuum expectation value (VEV) of the Ising Z 2 branch-point twist field (appendix A) and some auxiliary calculations.

Form factors of the branch-point twist fields in integrable models
Before presenting our results and discussing the determination of the form factors of modified branchpoint twist fields, it is instructive to give a brief overview of some basic ingredients of IQFTs and in particular on form factors of the conventional branch-point twist fields. Here we mostly follow the logic of Ref. [37] and present some of its results with an emphasis on the bootstrap equation.
Form factors (FF) are matrix elements of (semi-)local operators O(x, t) between the vacuum and asymptotic states, i.e., In massive field theories, the asymptotic states correspond to multi-particle excitations, with dis- where α i indicates the particle species. In such models, any multi-particle state can be constructed from vacuum state by means of the particle creation operators A † α i (ϑ) by where the operator A † α i (ϑ) creates a particle of species α i with rapidity ϑ and |0 is the vacuum state of the theory. In an IQFT with factorized scattering, the creation and annihilation operators where S α i ,α j (ϑ i − ϑ j ) are the two-particle S-matrices of the theory.
Our primary interest now is an n-copy IQFT and the corresponding branch-point twist fields.
For simplicity we assume that there is only one particle in the original theory. Then the scattering between the particles of different and of the same copies is described by S i,j (ϑ) = 1, i, j = 1, ..., n and i = j, S i,i (ϑ) = S(ϑ), i = 1, .., n, (2.4) and the branch-point twist fields are related to the symmetry σΨ i = Ψ i+1 , where n + i ≡ i. The insertion of a twist field T (or T n ) in a correlation function can be summarised as and we can also defineT , whose action is The form factors of the branch-point twist fields satisfy the following relations, which are simple modifications of the form factor bootstrap equations [61][62][63] where µ refers to the replica index of the particle, ϑ ij = ϑ i −ϑ j andμ = µ+1. In addition relativistic invariance implies F T |µ 1 ,µ 2 ,...,µ k k (ϑ 1 + Λ, . . . , ϑ k + Λ) = e sΛ F T |µ 1 ,µ 2 ,...,µ k k (ϑ 1 , . . . , ϑ k ), (2.10) where s is the Lorentz spin of the operator, which is zero for the branch-point twist fields. As the theories we consider in this paper have no bound states, Eqs. (2.7)-(2.9) and (2.10) give all the constraints for form factors of the twist fields.
As usual in this context, the so-called minimal form factor F T |j,k min (ϑ, n) is defined as the solution of the first two equations, Eqs. (2.7) and (2.8). That is, the minimal form factor satisfies It is then easy to show that from which it follows that and hence the only independent quantity is F T |1,1 min (ϑ, n). We can use Eq. (2.12) to determine it, writing (2.14) The solution of the last equation is easily obtained by noticing that if it exists a function f 11 (ϑ) operators, but with an S-matrix S(nϑ) instead of S(ϑ). When S(ϑ) can be parametrised as with some function g(t), the minimal FF is where the normalisation N ensures that f 11 (±∞) = 1 and thus The minimal form factors are very useful to obtain all form factors with particle number k ≥ 2 as they can be used as building blocks, hence simplifying the solution of the bootstrap equations.
The zero and one-particle form factors have to be determined by other means. The most important quantities are usually two-particle form factors. It can be verified that the two-particle form factors for the branch-point twist field, satisfying also the kinematic poles axioms, read [37] F T |j,k 2 (ϑ, n) = T n sin π n 2n sinh iπ(2(j−k)−1)+ϑ where T n = F T 0 is the vacuum expectation value (VEV) of T . Furthermore, relativistic invariance implies that F T |j,k 2 (ϑ 1 , ϑ 2 , n) depends only on the rapidity difference ϑ 1 − ϑ 2 , justifying writing F T |j,k 2 (ϑ 1 − ϑ 2 , n) or merely F T |j,k 2 (ϑ, n). It straightforward to show that forT we have (2.21)

Branch-point twist field form factors in the Ising model
The Ising field theory is surely the easiest integrable field theory. It has one massive particle (a free Majorana fermion) and the simple S-matrix and consequently For this model, it has been shown that the FFs of the branch-point twist fields are only non-vanishing for even particle number [37,48]. Moreover, the FFs for any even n can be written as a Pfaffain of the two-particle FF [49].

Branch-point twist field form factors in the sinh-Gordon model
The sinh-Gordon model, with Euclidean action is arguably the simplest interacting integrable relativistic QFT and for this reason it is often taken as a reference point and has been the subject of an intense research activity since many decades, see, e.g., [64][65][66][67][68][69][70][71][72]. Furthermore, it recently became also experimentally relevant because its nonrelativistic limit is the Lieb-Liniger Bose gas [73], a paradigmatic model for 1D ultracold gases [74].
The spectrum of the model consists of multi-particle states of a single massive bosonic particle.
The function g(t) entering in the parametrisation of the S-matrix (2.17) can be identified with (2.28) It is possible to write down an alternative representation of F T |1,1 min,ShG (ϑ, n) in terms of infinite products [37]. For and efficient numerical computation the following mixed representation is more useful .
Similarly to the Ising model, the FFs of the ShG branch-point twist fields are only non-vanishing for even particle number [37,48].
A very important relation between the ShG and Ising models is that the S-matrix and certain form factors of the ShG theory collapse to that of the Ising model, when the limit B = 1 + i 2 π Θ 0 with Θ 0 → ∞ is taken [68]. It can be checked that both F T |1,1 min,ShG (ϑ, n) and F T |j,k 2,ShG (ϑ, n) in this limit collapse to the corresponding quantities in Ising model. This limit will be an important guide for the case of the composite twist fields discussed below.
3 Form factors of the composite branch-point twist fields for Z 2 symmetry in integrable models After the introduction of the bootstrap equations for the FFs of the branch-point twist field, we now show how these equations can be naturally modified to obtain the corresponding quantities of the composite twist fields. At this point, of course, the existence of such fields is not strictly justified, therefore the formal solutions of the modified bootstrap equations will be subject to subsequent cross-checks.
To achieve our goal, first of all, we define the semi-local (or mutual locality) index e 2πiγ of an operator O with respect to the interpolating field φ via the condition  To be more precise about the connection between e i2πγ and e iα , let us consider briefly a U (1) symmetry for which α is a continuous parameter. From the point of view of the bootstrap equations, it is more convenient not to favour any of the Riemann sheets by adding the flux to it, but rather to divide the flux and introducing it on all sheets. This procedure corresponds to add a phase e iα/n on each sheet and therefore the locality factor e i2πγ and e iα/n must be equal. The further elaboration of the U (1) symmetry will be the subject of a subsequent work because, in this case, the particle content of the IQFT is richer and allows also for non-diagonal scattering leading to more complicated form factors. Here, we focus on the simpler, yet not trivial, analysis of the Z 2 symmetry in models with only one particle species.
However, for the Z 2 symmetry (and more generally for discrete symmetries) there are two subtleties that we cannot avoid mentioning. The first one is rather fundamental: for discrete symmetries Noether's theorem does not guarantee the existence of a conserved density, hence it is not a priori obvious if the reduced density matrix commutes with the symmetry operator. This problem will be discussed in the following sections for the specific cases of the Ising and ShG QFT. The other issue is that the phase is e iπ = −1 cannot be divided as e iπ/n among the various sheets, because e iπ/n no longer corresponds to the Z 2 symmetry of interest. This latter difficulty can be easily overcome by introducing the flux corresponding to the phase e iπ = −1 on all sheets. This step is legitimate if the number of sheets n is odd, as the overall phase acquired by a hypothetical particle winded through all sheets is still (−1) n = −1. Our argument implies that the composite branch-point twist fields associated with the Z 2 symmetry in the Ising and ShG models is a semi-local operator with respect to the fundamental field, with locality index e 2πiγ = −1. Specialising the bootstrap equations of a generic semi-local twist field to the Z 2 case, we have

7)
−i Res where T D denotes the composite branch-point twist field associated with the Z 2 symmetry. Having obtained the defining equations, following the logic of section 2, we can write for the minimal form factor F T D min of the composite twist field T D . From this we find and finally we get (3.10) Akin to the previous case, the only independent quantity is F For even n the above equation is equal to that of F T |1,1 min (ϑ, n), but our analysis is valid only for odd n. The solution of F Luckily, f D 11 can be easily obtained from f 11 by multiplying the latter by an appropriately chosen CDD factor, f CDD . Such a factor must obey The correct choice for f CDD turns out to be It is easy to check that the ansatz (3.15) satisfies Eq. (3.14), but it is not entirely trivial that there is no further ambiguity for the CDD factor and that Eq. Putting the various pieces together, the minimal form factor of the composite twist field is Given this minimal form factor, it is easy to show that Eq. (2.20) for two-particle form factors is It is easy to verify that Eq. (3.17) satisfies the axioms (3.5), (3.6) and (3.7).
Analogously to Eq.
where σ x/z i are the Pauli matrices. The conserved charge corresponding to the Z 2 symmetry is the fermion number parityP Q . HereQ =Q A +QĀ is the fermion number operator, which is clearly additive, andĀ denotes the complement of the region A. Crucially, the parity operator has eigenvalues 0 or 1 and the spacial restriction of this operator is also additive in a mod 2 sense, i.e., where we introduced the shorthandP Q A asP A .
An important quantity directly related toP is (−1)Q. This quantity can be expressed in several ways allowing for the computation of the symmetry resolved entropies in the critical point of the Ising model [6] and in its off-critical, lattice version [21], serving as valuable benchmark for our approach. WritingP as and introducing the disorder operators µ z i = i≤j σ x j and µ x i = σ z i σ z i+1 (satisfying the same algebra of the Pauli matrices), we have when the region A is a single interval from site 1 to . We recall that the disorder operator exists in the continuum limit as well. From Eq. (4.4) it is easy to deduce that the Z 2 branch-point twist field must be related to fusion of the usual branch-point twist field and the disorder operator. This picture is confirmed explicitly at the critical point of the Ising field theory [6], which corresponds to a conformal theory with central charge c = 1 2 . The scaling dimension of µ is ∆ µ =∆ µ = 1 16 and the symmetry resolved Rényi entropies for and interval of length read [6] S n (P A ) = −(n−1/n)/12 1 2 where P A is either 0 or 1.
i.e. when the flux can be inserted on each of the n copies.
The solutions for the bootstrap equations (3.5), (3.6) and (3.7) with locality index e i2πγ = −1 for the Z 2 branch-point twist field in the Ising model are easy to obtain. For the minimal form factor we have is obtained by (3.17). As anticipated, and also confirmed later on in this section, the Z 2 branch-point twist field can be regarded as a fusion of the conventional twist field and the Ising disorder operator (on the same lines of the composite fields for non-unitary theories [53]). In the off-critical theory, the FFs of both fields are non-vanishing only for even particle numbers. It is therefore natural to expect that F T D k is also vanishing for odd k. Nevertheless, even with the presence of FFs for odd particle numbers, their knowledge would be not necessary for any of the considerations of this paper [48] and, in fact, the VEV and the two-particle FFs encode all the physics we are currently interested in.
The FFs for even particle number F T D 2k with 2k ≥ 4 can be written as a Pfaffian of the two-particle FF, similarly to the case of the conventional branch-point twist field. For example, considering the bootstrap equations for particle numbers 2k = 4 and 6, it can be directly verified that F T D k indeed admits a Pfaffian representation. In particular, for j 1 ≥ j 2 ≥ ... ≥ j 2k , one has where W is a 2k × 2k anti-symmetric matrix with entries For general k, the Pfaffian structure (4.8) can be shown by induction, following exactly the same lines of the proof for conventional twist-fields [49]. If the ordering of the indices j i is not the canonical one, using the exchange axiom (3.5) one can reshuffle the particles and their rapidities to have j 1 ≥ j 2 ≥ ... ≥ j 2k so to apply (4.8). When the order of particles with the same replica index is left unchanged, the reshuffling does not introduce any ±1 factors.
Non-trivial checks of the solutions are provided by the limit for n → 1 and the ∆-theorem [60].
For n → 1, one expects to recover the form factors of the disorder operator; in particular for the two-particle case we expect with µ Ising denoting the vacuum expectation value of µ Ising . The limit of the Z 2 branch-point twist field in the Ising model is Since also the FFs of the Ising disorder operator can be cast in a Pfaffian form relying on the two-particle FF, the match between the two-particle FFs implies that The second test for the validity of the solution is given by the ∆-theorem sum rule [60]. The ∆-theorem states that if at some length scale R the theory can be described by a CFT, then the difference of the conformal weight of an operator O and its conformal weight in the infrared (IR) limit can be calculated as (if the integral converges) where Θ is the trace of the stress-energy tensor. Writing the spectral representation of (4.13) in terms of form factors, we have where m is a mass scale r = Rm and mE n are the n-particle energies. For the case of the massive Ising model, the conformal weights in the IR limit are zero. Hence taking r = 0 in (4.14) gives the UV conformal dimension of the operator O as In the Ising field theory, as well as in its n-copy version, the field Θ has non-vanishing form factors only in the two-particle sector, so the sum is terminated by the k = 2 contribution. After easy manipulations, the same as in Ref. [37] for the conventional twist fields, Eq. (4.15) for the Z 2 branch-point twist field can be written as We evaluated the integral in (4.16) numerically for many integer odd n using the FF (3.17). We found that the numerical calculated integrals match perfectly the prediction c 24 n − n −1 + ∆ n [6] with c = 1 2 and ∆ = 1 16 for all the considered n. Such perfect agreement is a strong evidence for the correcteness of the FF F from the FFs of the conventional twist fields by an additional CDD factor (3.15) and a different sign prescription in (3.10). As seen in the previous section, the corresponding solution for the Ising model can be associated with the Z 2 symmetry resolution of entropies. Nevertheless, the question of whether the symmetry resolution is possible, i.e., some/any reduced density matrices commute with the operator corresponding to the Z 2 symmetry is a rather difficult one for the ShG model. In the following, we present a series of arguments to claim that such a symmetry resolution is plausible at least for a single interval in the ground state of the model.
The first argument is based on the application of the Bisognano-Wichmann theorem [82] to the ShG model. This theorem states that for the ground state of a spatially infinite relativistic QFT, the reduced density matrix of a half-infinite line can be written as with the modular (or entanglement) Hamiltonian K where H is the hamiltonian density. For the ShG model, the hamiltonian density H ShG is invariant under the Z 2 transformation ϕ → −ϕ, hence K and ρ commute with the Z 2 symmetry operation.
The ShG model is a massive theory, and hence it is plausible that the RDM of an interval still commutes with the symmetry operation, at least for long enough distance, which is the case for which we eventually apply the novel form factors.
A second argument is given by the conformal limit of the ShG model, which is a free massless conformal boson. For the ground state of CFTs, the modular Hamiltonian is also known for a single interval of length 2R [83][84][85] and reads The Hamiltonian density of the free massless boson is again invariant under the Z 2 transformation ϕ → −ϕ, and, repeating the previous reasoning, the possibility of the symmetry resolution is justified in the UV regime.
Finally, we consider another limit of the ShG theory, namely when B = 1 + i 2 π Θ 0 with Θ 0 → ∞. As already noted, in this limit the form factors of the ShG model reduce to those of the Ising model.
As shown below, F T D |j,k 2,ShG (ϑ,n) is no exception to this rule, because the CDD factor f CDD (ϑ) is the same for the Ising and ShG models and holds: this link between the two models provides another evidence for the plausibility of a Z 2 symmetry resolution of the ShG model.
It is now worth studying some features of these FFs and in particular the two-particle one, F T D |j,k 2,ShG (ϑ, n). First of all, similarly to the Ising model, it is expected that F T D k,ShG vanishes for odd k. The reason is always the same: the Z 2 branch-point twist field can be regarded as a fusion of the conventional ShG twist field and the ShG disorder operator or twist field (which should not be mistaken for the branch-point twist field). In the off-critical theory, the FFs of both fields are non-vanishing only for even particle numbers. Considering now the two-particle FF solution, an interesting insight is given by the n → 1 limit of F T D |j,k 2,ShG (ϑ, n). The first few form factors of the ShG twist field are known and were constructed in [86]. This field can be identified with the off-critical version of the twist field of the massless free boson theory, where a unique field exists which changes the boundary condition of the boson field from periodic to anti-periodic and vice versa. This field has conformal weight ∆ = 1/16 = 0.0625 [87] and can be regarded as bosonic analogue of the fermionic disorder operator.  Table 5.1: The two-particle contributions of the ∆-theorem sum rule compared with the expected conformal dimension of Z 2 and conventional branch-point twist fields in ShG model.
We now show that in the limit n → 1, F T D |j,k 2,ShG (ϑ, n) coincides with F D 2,ShG (ϑ), where F D 2,ShG (ϑ) is the two-particle form factor of ShG twist field (again, the disorder operator, not the branch-point one). According to Ref. [86],   1+cosh(ϑ 1 −ϑ 2 ) to prove our claim. Based on this finding, it is natural to expect that the UV scaling dimension of the ShG Z 2 twist field is c 12 n − n −1 + ∆ n with c = 1 and ∆ = 1/16. We close this section showing that the ∆-theorem [60] is consistent with this assumption. Unlike for the Ising model, the form factors of the stress energy tensor in the ShG model are non-vanishing for the k = 4, 6, ...-particle sectors. In the integral formula of the ∆-theorem only the two-particle contribution is included and so it is not expected to be exact, but still to be a very good approximation. We calculated numerically such total 2-particle contribution for several B confirming such expectation. In the table 5.1 we show such comparison for B = 0.4 and 0.6. Notice that the two-particle contribution is always slightly larger than the expected total value and the difference is larger for larger B (up to B = 1), which is a general feature of the ShG model. This is very similar to what observed for the conventional twist field in Ref. [37] and also the difference is of the same order of magnitude. We stress that the fact that the offset is positive is an error (as the non-ideal name 'sum rule' would suggest): in Eq.
In any 2D QFT, the two (charged and neutral) moments entering in the Rényi entropies of an are written as Z n (0) = Trρ n A = ζ n ε 2dn T n (u, 0)T n (v, 0) , (6.5) where ε is the UV regulator, ζ D n and ζ D the normalisation constants of the composite and conventional twist fields, respectively, and d n and d D n their dimensions, given as where ∆ is the dimension of the field that fuses with the conventional twist-field to give the Z 2 composite one (e.g. the disorder operator in the Ising model or ShG with dimension ∆ = 1/16).
It is then clear that in the two symmetry resolved entropies (6.4), in the QFT regime ε 1, we have Z n (1) Z n (0) because ∆ is positive. Hence we find the 'trivial', yet general, result where S n is the total Rényi entropy. For general n the total Rényi entropy is known for some models, see e.g. [37,48], but its form is rather cumbersome. Instead, in the von Neumann limit, the result considerably simplifies in a generic massive model to [37] where U is a model dependent constant (e.g. calculated for the Ising model in [37]) and m the mass of the lightest particle of the field theory. We anticipate that for n = 1, the corrections in (6.8) gets multiplied by ln ε, as we shall see later in this section.
In spite of its triviality, Eq. However, this is not the end of the story. Eq. (6.8) with (6.4) shows that there are corrections to entanglement equipartition that are calculable within the integrable QFT framework of this paper.
In fact, expanding Eq. (6.4) for Z n (1) Z n (0) we have Notice that for generic n > 1, the ratio Zn(1) Zn(0) is proportional to ε 4∆/n while Z 1 (1) ∝ ε 4∆ and so the former is the leading correction. The two corrections become of the same order in the physically relevant limit n → 1. Notice that these corrections are very much reminiscent of the unusual corrections to the scaling [88,89] as calculated in massive theories [90]. This is not a coincidence since also unusual corrections in field theory come from the fusion of the twist field with a relevant operator [89].
Exploiting Eqs. (6.5) and (6.6), we have . (6.12) This expression provides the leading term breaking equipartition of entanglement for n > 1. With the exception of the normalisation amplitudes ζ n and ζ D n which depend on the precise UV regularisation of the theory (lattice in the following), all the quantities entering in the above ratio are in principle accessible to the bootstrap approach and calculable once the FFs are known.
In the von Neumann limit, n → 1, it is convenient to write down some general formula before taking the limit Z n (1) Z n (0). In general we have where, once again, S is the total entropy, and we defined (6.14) We now take the limit Z n (1) Z n (0) (implying Z 1 (1) 1 and s(1) S), obtaining Here the terms SZ 1 (1) and s(1) behave as ε 4∆ ln ε, while Z 1 (1) is proportional to ε 4∆ . Hence the breaking of equipartition of the von Neumann entanglement entropy at leading order is fully encoded in the quantities Z 1 (1) and s(1) defined above. These are obtainable in the FF approach and we will show with an explicit calculation for the Ising field theory in the next section. Although these terms breaking equipartition are vanishing in the field theory limit, they can be straightforwardly evaluated in any numerical computation (e.g. taking the difference S(+) − S(−) which cancels the leading term and isolate the correction). Such numerical computations can be verified against the predictions after having identified (as e.g. done in the next section for the Ising model) or fitted the non-universal UV cutoff ε. The remaining difference is a universal scaling function of m which is calculable within the FF approach, as again shown for the Ising model in the forthcoming section.

Entropies from two-point functions of the Z 2 branch-point twist field in the Ising model
In this section we show how the calculation of the symmetry resolved von Neumann entropies can be carried out based on the knowledge of the Z 2 branch-point twist field. We restrict our analysis to an interval in the ground state of Ising model in the paramagnetic phase, where the entropies can be calculated from the two-point functions of the conventional and composite twist fields. Our findings will be checked against the continuum limit of the existing results for the lattice model [21].
The calculation follows the logic of Ref. [37] including also steps like the determination of the vacuum expectation value of the Z 2 branch-point twist-field, the analytic continuation of the charged moments, and some further technical, but relatively straightforward, algebraic manipulations. The interested reader is encouraged the consult to corresponding appendices, where we report all the steps not strictly necessary to follow the main ideas.
The symmetry resolved entropies for one interval can be calculated in terms of two-point function of the composite and conventional twist fields just plugging (6.6) and (6.5) into (6.4) and (6.13) (or even to (6.11) and (6.15)). The partition sum Z n (0), i.e., Eq. (6.5), determines the total entropy and all the required quantities for its calculation S n were derived in Ref. [37] (including the analytic continuation). Concerning Z n (1) in Eq. (6.5), the two-point function of the Z 2 twist field and its vacuum expectation value can be determined using purely QFT techniques, whereas the proportionality constant can be fixed by comparing the lattice and QFT results. Explicitly, we Focusing now on the von Neumann entropy, we only need to know Eqs. (6.5) and (6.6) in the vicinity of n = 1. Hence, on top of Z 1 (1) given by Eq. (7.1), we also need its derivative in 1 which we rewrite as We stress that the entire dependence, which is the main focus of this approach, is fully encoded in the universal function H n (m ). The easiest part of the above expressions is dd D n dn , i.e.
In the two following subsections we explicitly calculate all amplitudes and two-point functions of composite twist fields.

Computation of the amplitudes
In Eqs. (7.1) and (7.2), a first ingredient yet to be calculated is the amplitude ζ D n . For n = 1 there is a straightforward way to get it, exploiting the fact that T D 1 equals the standard disorder operator. We can then write where m is the field theoretical mass and v the velocity of light, that in our notation is 1. The where A=1.282427129... is Glaisher's constant. Using now that T D 1 (x, 0) = µ(x, 0), we have The only missing ingredient to find ζ D 1 is the relation between the lattice spacing a and the UV regulator ε that was established in [37] and reads ε = χa, with χ = 0.0566227 . . . . with , (7.13) i.e., the complete elliptic integral. Obviously k 1 = k and k 1 = k . Hence, for n = 1, Eq. (7.11) is (1)| = √ k , that close to the critical point is (2(h − 1)) 1/4 = (2ma) 1/4 . On the other hand, directly in the continuum limit we have Eq. (7.1), which in the limit of large separation and for n = 1 is that provides for ζ D 1 exactly the same result as in Eq. (7.10). The other amplitude to be calculated is ∂ ln ζ D n ∂n n=1 in Eq. (7.2). We can use the last procedure to get this amplitude using (1) derived from Eq. (7.11) in [21], obtaining, for Recalling that, by definition, lim Rearranging the previous expression, one can extract ζ D n and its derivative with respect to n to get Notice that the term in ln(am) cancels, as it should. We also used ε = aχ, cf. Eq (7.9).

The two-point function of composite twist fields
Now we change focus and consider the two-point function entering in Eqs. (7.1) and (7.2). For n = 1, the two-point function of the composite fields in Z 1 (1) is just to the two-point function of the disorder operators, which can be also expressed in terms of a solution of a Painlevé III type differential equation [96]. However, for our purposes, the two-particle approximation of the twopoint functions is more useful because it provides not only the two-point function at n = 1, but also its derivative with respect to n. In this two-particle approximation, the correlation function for generic n can be written as (cf. Eq. (3.17) with (4.7)) We have already argued that the k-particle form factors of the Z 2 twist field vanish for odd k in both the Ising and ShG models. It has been also shown that the possible presence of a one-particle FF is irrelevant for the leading behaviour of the total entropy [48]. Overall, Eq. (7.21) allows us to identify the universal function H n (m ) in Eq. (7.1) in the two-particle approximation as dϑf D (ϑ, n)K 0 (2m cosh(ϑ/2)) , (7.23) an expression that is valid for a generic Z 2 symmetric theory with only the precise form of f D (ϑ, n) depending on the model. Eq. (7.23) with (7.22) provides an explicit final result for the Rényi entropies for any odd integer n ≥ 3 (we recall our FFs are derived for odd n). The calculation of the von Neumann limit n → 1 is more involved because it requires the analytic continuation of Eq.
(7.22) which is not an obvious matter, as we will see soon. However, before embarking in this more difficult calculation, let us consider the explicit form of Z 1 (1). In this case, the form factors of the composite twist field become those of the disorder operator, cf. Eq. (4.10), getting F µ 2 ∝ tanh ϑ/2, cf. Eq. (4.11). Hence we immediately have where the leading term in the m expansion is obtained below, but it can also be extracted using the fact that the integral in (7.24) can be rewritten in terms of the Meijer's G-function (although its form is not illuminating and we do not report it here).
Looking at Eq. (7.2) for s(1), we still need the derivative of both the VEV and of the universal function H 2pt n (m ). The former is rather cumbersome, but does not require any particular care and it is then reported in appendix A, see Eq. (A.32) for the final result. Conversely, the analytic continuation of H 2pt n (m ) is more thoughtful and we report its details in the following. In the two-particle approximation, the required derivative reads where we introducedf D (ϑ, n) which is the analytic continuation of f D (ϑ, n). The evaluation of f D (ϑ, 1) and of its the derivative, nevertheless, involves some subtleties related to the proper analytic continuation in n of the FFs, which is non-trivial as carefully discussed in Ref. [37] for the conventional twist field. For any integer odd n ≥ 3,f D (ϑ, n) = f D (ϑ, n). This is no longer true for −π 2 1 2 δ(ϑ) = 4ϑ sinh 2 (ϑ/2) sinh 3 ϑ −tanh 2 (ϑ/2)−π 2 1 2 δ(ϑ) , (7.27) It follows that the final result for Eq. (7.25) is This term, together with (7.24) includes the entire dependence of the symmetry resolved von Neumann entropies and it represents our final full result.
However, putting the various pieces together is not illuminating without expanding for large m as we are going to do now. The leading term in (7.28) clearly comes from the K 0 (m ) factor, but it is worth discussing a simple method to obtain a systematic large expansion. To obtain the subleading terms by evaluating the integrals in Eqs. (7.28) and (7.24), one first recognises that for large , the integral is dominated by the contribution of the region close to ϑ = 0. One can then expand as a function of ϑ = 0 the function which multiply K 0 (m ) in the integrand, and evaluate the asymptotic behaviour of Expanding arcosh(x) around x = 1, exploiting the asymptotic behaviour of the Bessel function K 0 (z) ≈ e −z π 2z , and keeping the leading x − 1 type terms, we and up with which gives the leading -dependent term for (7.29). In this way, one readily derive the expansion in the rhs of Eq. (7.24) and

Putting the pieces together
In this subsection we put together the different pieces of the symmetry resolved entropies. We first of all write down the expressions for Z 1 (1) and s(1) including the leading corrections and then comment on the symmetry resolved entropy. Z 1 (1) is obtained by plugging Eqs. (7.24) and (7.14) into Eq. (7.1), getting where we introduced the combination of amplitudes 1 + 1 8π (7.33) and (7.35) we also kept the leading and subleading terms accounting for the -dependence.
The analogous term incorporating -dependence has not been derived for the lattice model and represent one of our main achievements. With (7.32) for Z 1 (1) and (7.33) for s(1), we can finally use (6.15) to write down the symmetry resolved entropies including corrections too. Keeping the ε 1/4 ln ε and ε 1/4 terms, we end up with

Conclusions
In this paper, we introduced an approach suited to the computation of symmetry resolved entropies in generic massive (free and interacting) integrable quantum field theories. The essence of the approach is the existence of appropriate modified or composite branch-point twist fields whose twopoint function gives the corresponding charged entropies for a single interval. Then the form factor bootstrap program provides the matrix elements of such fields. In particular, here we discussed the Z 2 symmetry resolution for Ising model in the paramagnetic phase and for the sinh-Gordon quantum field theory.
We wrote down the bootstrap equations for the composite twist fields and provided an intuitive picture behind the choice of the locality factors entering these equations. The two-particle form factors for Z 2 branch-point twist fields were calculated for the Ising both models considered here.
For the Ising model, we were also able to compute the vacuum expectation value, alias the zero particle form factor, we argued that form factors with odd particle number vanish, and finally showed that the form factors for any even particle numbers can are Pfaffian of the two-particle form factors. The obtained form factor solution was cross-checked verifying that for n → 1 the form factors of the disorder operator are recovered and applying the ∆-theorem [60] to reproduce exactly the critical dimensions of the composite fields.
Also the sinh-Gordon form factors have been tested in several ways. First, we considered the limit for the interaction parameter B as B = 1 + i 2 π Θ 0 with Θ 0 → ∞, in which the Z 2 branch-point twist fields for the Ising model are recovered. Then for n → 1, we reproduced the disorder operator of the sinh-Gordon model. Applying the ∆-theorem for the form factors, we recovered the expected UV dimensions with satisfactory precision. The error is ascribed to the fact that, unlike for the Ising model, the ∆-theorem sum rule requires an infinite summation and hence the knowledge of all form factors for the Z 2 branch-point twist field.
The general approach to extract the ground-state symmetry resolved entropies for an interval of length from the two-point function of composite twist fields is discussed in Sec. 6. In particular, we showed that entanglement equipartition follows generically from the property that the UV dimension of the composite twist field is larger than the one for the conventional twist field. The subleading term breaking such equipartition is model dependent. The obtained form factors allow for the complete calculation of the charged and symmetry resolved entropies in the paramagnetic phase of the Ising model which is presented in great detail, with emphasis on the physically relevant von Neumann limit n → 1 (that requires a non-trivial analytic continuation). The final result for the charged partition sum and entropy are reported in Eqs. (7.1) and (7.2) with the various amplitudes computed in Sec. 7.1 and the universal functions of m given in Eqs. (7.24) and (7.28). We stress that these universal functions are the main new physical results of this paper since all other terms could be equivalently calculated by taking the continuum limit of the known results for the Ising chain in Ref. [21]. From Eq. (7.37) we can see that the leading term breaking equipartition scales like ε 1 4 ln ε, as expected. However, Eq. (7.37) also provides the m dependence of this equipartition breaking term. It would be highly desirable to test all these predictions with exact numerical calculations based on the continuum limit of the spin chain.
There are various possible ways this work can be extended. The most natural one is the treatment of models with non-diagonal scattering and continuous symmetries, to which the authors plan to devote another communication. The obtained form factors also allow for the calculation of entropies in excited states, as long as reduced density matrix commutes with the symmetry operator.
Finally, the crossover from critical to massive regime at fixed is a very interesting yet challenging problem, which may require an infinite summation higher particle form factors or the development of alternative techniques.

Acknowledgments
DXH is grateful to Sara Murciano for many useful discussions. The authors are also grateful to Olalla Castro-Alavaredo for fruitful feedbacks on a first version of the manuscript. PC and DXH acknowledge support from ERC under Consolidator grant number 771536 (NEMO).

Ising
Finding the solutions to the FF bootstrap equations is relatively easy. Often it is also not difficult to identify these solutions with the corresponding physical fields. Conversely, the determination of the vacuum expectation value (VEV), i.e., the zero particle FF and the one-particle FF (if nonvanishing) is generally a difficult task. So far, exact expressions are known for all fields in the Ising model and for some in ShG, sine-Gordon, Bullogh-Dodd models, as well as for some of their restrictions, see e.g. [64,[91][92][93]. For the conventional branch-point twist fields, an exact expression for the VEV has been provided only for the Ising model in [37]. In this appendix, we show that for the same model the VEV for T D n can also be exactly determined, under some plausible assumptions. We use and modify ideas borrowed from Refs. [37,94,95]. In this appendix, we work in the fermionic basis and denote the j-th copy of the Majorana fermion as ψ j . We explicitly exploit the property that fermionic and spin entanglement are the same for one interval.
As a first step we search for a matrix τ whose action in the space replica space (i.e. on the vector (ψ 1 , ..., ψ n ) T ) corresponds to the the composite twist field. Given that the total phase accumulated by the field in turning around the entire Riemann surface is −1, the main requirement is τ n ψ j = −ψ j , i.e., τ n = −I, where I is the n × n identity matrix. An easy way to proceed is to modify the transformation matrix for the conventional twist-fields [95], as done in Ref. [19] for the resolution of the U (1) symmetry (both papers consider Dirac fermions, but there is no difference for Majorana except that the phase is fixed). Hence, a first representation of the matrix τ is where it is clear that τ n 1 = −I for odd n. However, it was pointed out in [37] that one has to be careful in the FF approach because fermions of the same copy anticommute, as conventional fermions do, but the fermions of different copies commute (S ij = 1). Conversely, in Refs. [19,95] fermions of different copies anticommute. The anticommutation of fermions on different copies can be achieved in the FF approach by a change of basis as [37] As argued in [37], the action of a permutation on the fields ψ ac j in the new basis is no longer σψ ac j = ψ ac j+1 mod n , but instead When this permutation is applied n times we have σ n ψ ac j = −ψ ac j . Moreover, the eigenvalues of the corresponding matrix equal those of (A.1) for odd n, which the case we are interested in. We can then identify both τ 2 and τ 1 with the transformation matrix that has to be diagonalised for the determination of the VEV [37]. The eigenvectors corresponding to the eigenvalues e i2πk/n are complex conjugate pairs for ±k, except k = n/2 with eigenvalue (−1) and real eigenvector equal to 1 √ n (1, −1, 1, ..., 1). Hence, we can build n−1 2 complex fermions by ψ k and ψ −k as ψ † k = ψ −k for k = 1, . . . , (n − 2) and we are left with one Majorana fermion for k = n/2, which is still a Majorana fermion as ψ † n/2 = ψ n/2 . The anticommutation relations {ψ k , ψ k } = δ k,−k , {ψ k , ψ n/2 } = 0 for k = n/2, and {ψ n/2 , ψ n/2 } = 1 are ensured by our choice for the basis (A.2).
The structure of the eigenvalues of the transformation τ is compatible with the four-point function of the Z 2 twist field at the UV critical point: turning clock-wise ψ k (z ) around the twist field T D at w, the correct factor of e i2πk/n is recovered. Eq. (A.8) is an important formula, which is also proved in Appendix B. It leads to the factorisation of the Z 2 branch-point twist field, it allows for the computation of the UV dimensions of the factorised components, and eventually it leads to the determination of the VEV in the massive theory. The factorisation of the Z 2 twist field can also be inferred from the results of [94], which in our case become where action of T D k,n (w) is non trivial only on the ψ −k and ψ k fields. The scaling dimension of T D k,n can be can be obtained from the relation [14,38,39] where T k is the stress-energy tensor of the ±k components. In fact, using the Ward identity [97] T one can deduce that the coefficient h k in (A.10) equals the conformal dimension of the chiral component of both T D n andT D n . To calculate (A.10), we first show, that the stress-energy tensor can also be factorised into different k-components. We recall that the 2D free massless Dirac theory can be written in terms of the two component Dirac spinor Ψ(z,z) = χ(z) χ(z) , where χ andχ are complex fermion fields. The analytic part of the stress energy tensor is whereas for the neutral Majorana field it reads One Dirac field can be constructed from two Majorana fields as but in our case, as argued before, it is more convenient to use with our Fourier transformed fields ψ k . In this way, the stress-energy tensor of the original n-copy model is decomposed into k sectors each involving complex fermion fields. Using Eq. (A.12), the stress-energy tensor of the ±k components is 16) for k = 1 2 , . . . , n−2 2 and, similarly for k = n The total stress-energy tensor is then Now we explicitly compute the lhs. of Eq. (A.10) to determine h k . We first notice that the action of 1 2πi˛d to the lhs of Eq. This dimension can be also rigorously obtained by applying 1 2πi˛d (A.23) The factor 1 4 in (A.22) compared to 1 2 in (A. 19) is important to obtain the desired − 1 2 ψ n 2 (z)∂ z ψ n 2 (z) with the correct normalisation. The application of (A.22) to (A.8) results in T n 2 (z)T D confirming h n 2 = 1 16 .
Finally, the total dimension of the composite twist field is which is the correct dimension in the Ising CFT as h +h correctly reproduces 1 2 1 12 n − n −1 + 1 8n . We have also seen that, winding the complex fermion field χ k (z) = ψ k (z) around the branchpoint twist field, a phase e iπk/n is accumulated for k = n 2 , which can be attributed to the action of a U (1) composite twist field. A plausible assumption is that the decomposition of branch-point twist fields can be rephrased as Assuming that this type of factorisation of the Z 2 branch-point twist field also holds in the off-critical theory we can obtain its vacuum expectation value exploiting the results in Ref. [91] O α = m 2 where G(x) is the Barnes G-function. Hence, for the n-copy Ising theory we have . (A.28) Using the exact result for µ Ising [96], we can write it as µ Ising = m

B Conformal dimensions
In this appendix we show that Eq. (A.8) holds for Z 2 branch-point twist field in the c = 1 2 CFT. Let us recall what we want to prove here: (B.1) The way we proceed is very similar to Refs. [19,94]. We apply the conformal transformation which maps the R n Riemann surface with branch-points w and w to the complex plane ξ ∈ C.
After this uniformising mapping, the twist fields in Eq. (B.1) do not disappear, but they become the disorder operator of the Ising CFT. This is a manifestation of the fact that T D is the fusion of T and the disorder field µ. To check the validity of this idea, we first compute the scaling dimension of T D along these lines.
We assume that |g D (ϑ, n)| < Ce q|n| for Re(n) > 0 and with q < π 2 ; this assumption is motivated by the fact that both Tr (ρ n A ) and Tr ρ n A (−1) nQ A behave so for finite systems, see again Ref. [37] for a detailed discussion. Then Carlson's theorem can be applied tof D (ϑ, n) −g D (ϑ, n) and implies that the difference is identically zero, i.e. the continuation is unique. To be more precise, we use Carlson theorem in its standard form [99] by applying it tof D (ϑ, 2n + 1) −g D (ϑ, 2n + 1), with n = 1, 2, 3, 4, .... The only price to pay is that the growth on the imaginary axis must be bounded by Ce π 2 |n| rather than the usual restriction Ce π|n| . Anyhow, this is compatible with both the limiting behaviour of f D (ϑ, n) and our motivating assumptions forg D (ϑ, n).