U(1) symmetry resolved entanglement in free 1+1 dimensional field theories via form factor bootstrap

We generalise the form factor bootstrap approach to integrable field theories with U(1) symmetry to derive matrix elements of composite branch-point twist fields associated with symmetry resolved entanglement entropies. The bootstrap equations are solved for the free massive Dirac and complex boson theories, which are the simplest theories with U(1) symmetry. We present the exact and complete solution for the bootstrap, including vacuum expectation values and form factors involving any type and arbitrarily number of particles. The non-trivial solutions are carefully cross-checked by performing various limits and by the application of the Delta-theorem. An alternative and compact determination of the novel form factors is also presented. Based on the form factors of the U(1) composite branch-point twist fields, we re-derive earlier results showing entanglement equipartition for an interval in the ground state of the two models.


Introduction
Symmetries play an inevitable role in modern physics and deeply interlace the whole discipline in versatile ways. Restricting to merely a practical point of view, it is generally a very fruitful idea to exploit the additional structures imposed by symmetry for almost any type of physical object.
Interestingly, such an idea has only recently been investigated in the study of entanglement of extended quantum systems. As widely known, when a system is in a pure state, the bipartite entanglement of a subsystem A may be quantified by the Rényi entanglement entropies [1][2][3][4] S n = 1 1 − n ln Trρ n A , defined in terms of the reduced density matrix (RDM) ρ A of the subsystem A. From those, in the replica limit n → 1 the von Neumann entropy is obtained. However, the knowledge of Rényi entropies for different n contains more information than the simple S, as, e.g., it provides the entire spectrum of the reduced density matrix ρ A [5,6].
The recent and explicit idea of considering generally the internal structure of entanglement associated with symmetry was put forward in Refs. [7][8][9][10]. In a symmetric state, the conserved charge corresponding to the symmetryQ commutes with density matrix; under general circumstances also the restriction ofQ to the subsystemQ A , commutes with the RDM as Such commutation implies that ρ A is block-diagonal, each block corresponding to an eigenvalue ofQ A . This fact has the important consequence that the Rényi and von Neumann entropies can be decomposed according to the symmetry sectors ofQ A . The symmetry resolved Rényi and von Neumann entropies are defined as S n (q A ) = 1 1 − n ln Z n (q A ) Z n 1 (q A ) , and S(q A ) = − ∂ ∂n in terms of the symmetry resolved partition sums Z n (q A ) = Tr (ρ n A P(q A )) , where P(q A ) is the projector onto the sector corresponding to the eigenvalue q A .
The calculation of all these symmetry resolved quantities requires in general the diagonalisation of ρ A and the resolution of the spectrum in the conserved charge. An ingenious way to circumvent this difficult path passes through the charged moments [8] Z n (α) = Tr ρ n A e iαQ A , that are nothing but the Fourier transform in of the desired partition sums, i.e. [8] Z n (q A ) = Tr (ρ n A P(q A )) =ˆπ −π dα 2π Z n (α)e −iαq A .
A very standard and powerful trick, used a lot in quantum field theory (QFT) context, is to replace the partition function on the n-sheeted surfaces with an n-copy version of the original model with specific boundary conditions among the replicated fields [45]. Crucially, in 1+1 dimension, these boundary conditions can be implemented via local fields in the n-copy theory that are known as branch-point twist fields [12,46]. The moments of ρ A are eventually equivalent to an appropriate multi-point function of the branch-point twist fields in the n-copy theory. In 2D CFT, the scaling dimensions of these fields are exactly known [12,47,48] and provide their two-point function, i.e.
In 2D CFT, the symmetry resolved entropies are obtained as multipoint correlations of composite branch-point twist fields which fuse the action of the replicas and of the flux of charge [8]. These composite twist fields have been recently identified also in some massive theories: free massive Dirac and complex boson QFT [20], the off-critical Ising and sinh-Gordon theories [21]. As shown in our previous work [21], form factors of the composite twist fields can be determined with the bootstrap program, similarly to the usual branch-point twist fields [46,58,59]. Using these form factors, one can obtain a systematic expansion for the correlation functions of the composite twist fields, which eventually relate to the symmetry resolved entropies.
Our previous work [21] focused on the symmetry resolution of the discrete Z 2 symmetry in the Ising and sinh-Gordon models [21]. Here, instead, we initiate the study of continuous symmetries: we use appropriate bootstrap equations for the composite branch-point twist fields related to the U (1) symmetry and solve these equations. For simplicity, we consider the free massive Dirac theory and the free massive complex boson theory. Although the symmetry resolution of these theories has already been performed [20], it is extremely instructive to study them from the point of view of form factor bootstrap, being the simplest examples of theories with more than one particle species. Indeed, multi-particle theories usually introduce quite a few technicalities, which can be kept at a minimum level in free theories. Nevertheless, as we show here, the form factors of the composite twist fields are non-trivial objects even in free theories. Hence, our main result is the determination and solution of the bootstrap equations for these form factors. We then use them to infer the behaviour of the charged moments when the subsystem is a long interval. Finally, for completeness, we also compute the symmetry resolved partition functions and entropies within the two-particle approximation.
The paper is structured as follows.  [76]. In section 5 a compact, alternative derivation of the 2-particle FFs are presented, which is based on an explicit diagonalisation in the replica space. Section 6 reports general and specific results for U (1) charged moments 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. We finally conclude in section 8, which is followed by the appendices containing some details on the analytic continuation of the replica index n → 1 (appendix A), and the determination of the vacuum expectation value (VEV) of the composite branch-point twist field (appendix B and C) .

Form factors of the branch-point twist fields in integrable models
A primary focus of this work is on the composite branch-point twist field and their FFs in free models, nevertheless, before discussing these objects, it is natural to review the standard branchpoint twist fields and the corresponding FFs first. These FFs are an important reference point for our later investigations. Moreover this brief overview allows us to conveniently introduce some basic ingredients of IQFTs, which are essential building blocks in our derivations. Such central objects in our investigations are the FF bootstrap equations and the FFs themselves. Following closely the logic of Ref. [46], we introduce the bootstrap equations for branch-point twist field and comment on their solution. Since it takes little effort and helps keep connection with interacting IQFTs, we keep the following discussion more general, which thus describes the case of interacting theories with diagonal but non-trivial scattering as well.
In massive field theories, the asymptotic states are spanned by multi-particle excitations whose dispersion relation can be parametrised as (E, p) = (m β i cosh ϑ, m β i sinh ϑ), where β i indicates the particle species and ϑ is the rapidity of the particle. In such models, any multi-particle state can be constructed from vacuum state |0 as where A † s are particle creation operators; in particular the operator A † β i (ϑ i ) creates a particle of species β i with rapidity ϑ i . In an IQFT with factorised scattering, the creation and annihilation operators A † β i (ϑ) and A β i (ϑ) satisfy the Zamolodchikov-Faddeev (ZF) algebra, which, for diagonal scattering, reads where S β i ,β j (ϑ i − ϑ j ) denotes the two-particle S-matrix of the theory that describes the scattering of particles β i and β j .
Turning to an n-copy IQFT one describes the scattering between particles in different and in the same copies as S (β i µ i ),(β j µ j ) (ϑ) = 1, i, j = 1, ..., n and µ i = µ j , where we introduced the replica index µ i , which takes values from 1 to n. As known, in the n-copy QFT the branch-point twist fields are related to the symmetry σΨ i = Ψ i+1 , with n + i ≡ i. When a twist field T (or T n ) is inserted in a correlation function, its action can be summarised as In a similar way, one can also defineT with action The bootstrap equations for branch-point twist fields are natural modifications of the form factor bootstrap equations for conventional fields [77][78][79]. Exploiting the relation of the twist fields to the symmetry σΨ i = Ψ i+1 , the corresponding bootstrap equations can be written as [46] F T |(βµ) k −i Res −i Res where ϑ ij = ϑ i − ϑ j , ϑ and (βµ) are shorthands for ϑ 1 , ϑ 2 , ..., ϑ k and (β 1 µ 1 ), (β 2 µ 2 )...., (β k µ k ) respectively;μ = µ + 1 andβ denotes the anti-particle of species β. In the last equationū δ βγ refers to the position of the pole of the bound state of the particle δ in the S-matrix S β,γ and i Γ δ βγ 2 is the corresponding pole-strength. Furthermore, relativistic invariance implies where s is the Lorentz spin of the operator, which is zero for the branch-point twist fields. The free Dirac and massive complex boson theories have no bound states; consequently Eqs. (14)- (16) and (18) give all the constraints for form factors of the branch-point twist fields.
The above bootstrap equations have many inequivalent solutions. Some additional information is generally needed to sort out a given physical solution. In this respect a very useful tool is the ∆-theorem [76]. Such a theorem states that when at some length scale R the theory is described by a CFT, the difference between the conformal weight of an operator O and its conformal weight in the infrared (IR) limit is where Θ is the trace of the stress-energy tensor. Using the form factors to write the correlation function in the integral above, we have where m is a mass scale r = Rm and mE n are the n-particle energies with E n = n k=1 m β k cosh ϑ k /m. For massive theories, the conformal weights in the IR are zero. Hence taking r = 0 in (20) gives the conformal dimension of the operator O as Consequently, given a set of FFs satisfying the axioms (14)- (17) and (18), the ∆-theorem allows us to identify to which physical operator it corresponds (as long as its conformal dimension is known).
Although in principle the ∆-theorem sum rule requires the knowledge of infinitely many terms and FFs, it is usually a very rapidly converging sum; therefore the first few-particle FFs are enough to identify the UV conformal dimension. Finally, we notice that the UV dimension imposes also a constraint on the growth of FFs for large rapidities. When one rapidity of the FF goes to ±∞, its The FFs of various fields are usually found first solving Eqs.
As a consequence we have and hence so that the only independent quantity is F To solve the last equation, we introduce f 11 (ϑ) satisfying so that Eq. (26) is just the standard equation for minimal form factors of conventional local operators, but with an S-matrix S(nϑ) rather then S(ϑ). Exploiting now the standard parametrisation of S(ϑ) with some function g(t), the minimal FF is where the normalisation N guarantees f 11 (±∞) = 1. Finally we have The minimal form factors are the building blocks to obtain all form factors with particle number k ≥ 2, simplifying the solution of the bootstrap equations. We recall that the zero and the oneparticle FFs must be obtained by other means. In the absence of bound states, the two-particle form factors (which are the most relevant FFs) for the branch-point twist field, satisfying also the kinematic poles axioms, have been determined as [46] F T |(βk),(βj) 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 |(βk),(γj) 2 (ϑ 1 , ϑ 2 , n) depends only on the rapidity difference ϑ 1 − ϑ 2 , justifying writing

Branch-point twist field form factors in the Dirac field theory
The free Dirac theory has a massive fermion particle (denoted by +) and anti-particle (denoted by −) and the simple S-matrix with elements Consequently we have For this model, the FFs of the branch-point twist fields are only non-vanishing for even particle number [46,59]. Moreover, the FFs for any even n can be written as a Pfaffain of the two-particle FF [60].
The FFs of the branch-point twist field is only non-vanishing for neutral states containing an equal number of particles and anti-particles. In the Dirac theory, in particular, we have [46,67] that is F . T D refers to the branch-point twist field of the Dirac theory and accordingly T D,n denotes the VEV of this field in the n-sheeted theory.

Branch-point twist field form factors in the free complex boson theory
Also, the free complex boson theory has a particle (+) and an antiparticle (−). The S-matrix is identically 1 between and within each species, so where ± refers to particles and anti-particles. The minimal form factor is also unity for any n.
The FFs of the branch-point twist field is only non-vanishing for neutral states containing an equal number of particles and anti-particles. In particular, we have the simple expressions [67] In this case we find again that F . T B now refers to the branch-point twist field of the free massive complex boson theory and accordingly T B,n denotes the VEV of this field in the n-sheeted theory.

The form factors of the stress energy tensor Θ in free theories
Later on we will need also the FFs of the Θ field and so we report here some basic details about them. In both the free massive Dirac and the free massive complex boson theory this field has only a non-vanishing 2-particle FF, which reads for the 1-copy Dirac and for the 1-copy complex boson theory, where m is mass of the fermion and boson. In the n-copy theories, the FF of this field behaves in an additive way, that is [46] F Θ|(±j),(∓k) and

Form factors of the composite U(1) branch-point twist field in the Dirac theory
The bootstrap equations for the FFs of the branch-point twist field can be naturally modified to obtain the corresponding quantities of the composite branch-point twist fields [21]. We define the semi-local (or mutual locality) index e iα of an operator O with respect to another field φ via the or, when using the radial quantisation picture, Mutually local operators correspond to e iα = 1, while fields with e iα = 1 are called mutually semilocal. In Ref. [21] it was argued that it is natural to assume that the flux phase e iα can be related with the mutual locality index appearing in the bootstrap equation. This assumption can be based on the intuitive picture associated with the insertion of the Aharonov-Bohm flux on one of the Riemann sheets. In this picture, the flux is carried by the particles of the theory, but Eq. (44) is just an equivalent rephrasing of this idea when φ is an interpolating field associated with creation or annihilation of particles. These ideas can be equivalently rephrased by writing the exchange relations for the composite branch-point twist fields of the n-copy theory as and similarly forT where Ψ p,i (x) is a generic quantum field living on the ith replica and possessing a definite U (1) charge p ∈ Z. Our choice for α/n is dictateted by the requirement, that the total phase picked up by the particle when turning around the entire Riemann surface has to be e iα . Specialising now the bootstrap equations of a semi-local U (1) composite branch-point twist field T α to the Dirac theory, −i Res −i Res where β and β i = ±1 (sometimes simply shortened in ±). T α D denotes the composite U (1) branchpoint twist field in the free massive Dirac QFT. To avoid confusion, the standard U (1) twist field in the Dirac theory is denoted by V α D . This field is also referred to as the vertex operator, but based on its monodromy properties it can be equivalently called a U (1)-twist field. In both cases, the superscript α corresponds to the inserted flux.
The FFs of the composite U (1) branch-point twist fields are non-zero for even number of particles and charge-neutral configurations. This is also true for V α D which can be obtained from the above equations for n = 1. It is instructive to write down its 2-particle FFs, known from earlier investigations [80][81][82], Having obtained the defining equations for the form factors, following the logic of section 2, we can for the minimal form factor F T α D min of the composite branch-point twist field. From this we find and finally we get Akin to the previous case, the only independent quantity is F The solution of F that satisfies Luckily, f T α D (±1),(∓1) can be easily obtained from f 11 by multiplying the latter by an appropriately chosen CDD factor, f . Such a factor must obey (ϑ)f 11 (ϑ) satisfies Eq. (56). The correct choice for and hence the minimal FF is It is easy to check that the ansatz (58) satisfies Eq. (57), and any further CDD ambiguity is excluded by exploiting the ∆-theorem checks performed in the next sections.
With the minimal form factor at hand, we can write the 2-particle form factor as a product where the function P  47) and (48), The only independent quantity is P (+1),(−1) satisfying together with the important property Combining Eq. (49) with the minimal form factor, we have We can easily check that for α = 0, i.e., the case of the conventional branch-point twist field when solves the above equations. When α = 0, we have to account for the complex phase factors but still respecting the cyclic property (63) of P D . This can be achieved multiplying P (ϑ, n, 0) by 2πn periodic functions along the imaginary axis, whose values at iπ and −iπ are e ∓iα/(2n) and e ±iα/(2n) , respectively. The simplest choice for such a function is (ϑ, n, α) = p ± D (ϑ, n, ϑ)P (ϑ, n, 0) satisfies the parity and cyclic property Eqs. (62) and (63) and the equivalent of (49), i.e., (64). Furthermore, p ± D (ϑ, n, α) does not introduce any additional poles on the physical sheet.
Using thus (53) and the product formula (60) for the 2-particle FF, all the defining axioms (47), (48), (49) are satisfied. Unfortunately this solution is not the correct one corresponding to the composite U (1) branch-point twist field. Although for α = 0 the FF of the conventional twist field is recovered, for n → 1 the solution is also expected to simplify to the conventional U (1) twist field form factor, which is not the case. Based on this quantity the right choice for p ± can be determined.
The correct solution for P in which the correct choice for p ± D (ϑ, n, α) does not introduce additional poles on the physical sheet (similarly to our previous naive choice) and for large rapidities it tends to a constant. Using the formula (67) to construct the 2-particle form factors via (60) and (61) all the axioms (47), (48), (49) are satisfied and the α → 0 and n → 1 limits yield the desired results. i.e., the conventional branch-point (35) and the standard U (1) twist field FFs (50), which one can easily check. Performing algebraic simplifications, the full two-particle FF can be therefore written as and, using Eq. (53), we arrive at In all the above treatment α is within the range [−π, π], otherwise descendant fields are obtained.

∆-theorem test and higher-particle FFs
The validity of the solution can be checked by the ∆-theorem, which for the Dirac field can be written as where we (i) used that in the Dirac theory the field Θ has only 2-particle FFs (ii) we performed a changed of variables (iiI) we did a final integration in the k = 2 term of Eq. (21). The chiral conformal dimensions obtained from the appropriate integrals of the 2-particle FFs gives the exact It is also easy to obtain the FFs of higher particle numbers in terms of a Pfaffian. Assuming the ordering of the replica index µ 1 ≥ µ 2 ≥ ... ≥ µ 2k , one has where W is a 2k × 2k anti-symmetric matrix with entries where F Finally, we note that the FFs of the other fieldT α D can be easily obtained from that of T α D . As an example, for the 2-particle case we have 4 Form factors of the composite U(1) branch-point twist field in the complex boson theory In this section we treat the case of the free complex boson for which the bootstrap equations of a semi-local U(1) composite branch-point twist field are −i Res −i Res for the minimal form factor F T α B min of the composite branch-point twist field. From this we find and finally we get Akin to the previous case, the only independent quantity is F The solution of F that satisfies To obtain f T α B (±1),(∓1) it is instructive to study first the conventional U (1) twist field V α B . The FFs of this field were computed in [83], and the two-particle FFs of this field read as where the range for α is now [0, 2π]. From Eq. (85), the minimal FF of the standard U (1) field as from which It is straightforward to check that the ansatz (84) satisfies Eq. (76) and Eq. (77).
Similarly to the case of the Dirac theory, we can again write the 2-particle form factor as a where P satisfies Combining Eq. (78) and the minimal form factor, we have different expressions for the poles compared to the Dirac case, namely We can further factorise P (ϑ, n, α) = p ± B (ϑ, n, α)P (ϑ, n, 0) .
P (ϑ, n, 0) is chosen to be the usual one but we must be more careful with the function p ± B (ϑ, n, ϑ). In fact, to be compatible with (78), it must be 2πn periodic functions along the imaginary axis, whose values at iπ and −iπ are e ∓i(α−π)/(2n) and e ±i(α−π)/(2n) respectively. The most natural choice for such a function is With this choice, we re-obtain the standard U (1) field, when n → 1. We also have to recover the conventional branch-point twist field for α = 0, but unfortunately it is not the case. The latter requirement can, nevertheless, be only satisfied if instead of sin α−π 2n sinh ϑ n sin π n in Eq. (95), we use sin α−π 2n sin π n sinh 2 ϑ 2n /(sinh ϑ n sin 2 π 2n ). Note that, interestingly the function that multiplies sin α−π 2n is exactly the inverse of the corresponding function for the Dirac theory, which multiplies sin α 2n in p ± . Eventually, we arrive at our final expression for P where the range of α is [0, 2π]. We can conform our results to the Dirac case and shift the range of α to [−π, π] by redefining our expression as Finally, using Eq. (53), we get
The FFs of higher particle numbers are obtained by using Wick's theorem. Assuming the ordering where {l, m} runs over the possible pairs and F Similarly to the Dirac theory, the FFs of the fieldT α B can be again easily obtained from that of T α B . As long as the range of α is set to [−π, π], the 2-particle FF can be written as

Form factors of the composite U(1) branch-point twist field via diagonalisation in replica space
In this section, we provide an alternative derivation of the two-particle FF of the modified twist fields, based on the diagonalisation in the space of replicas [45]. This technique has been already employed in [71] for the computation of the FF of the standard branch-point twist fields. It works well for free theories when the S-matrix does not depend explicitly on the rapidities and it is more closely related to the approach in Ref. [20]. We will briefly summarise this formalism, discussing the case in which a non-vanishing flux is inserted.
Let us consider the creation operator A † j,± (ϑ) of a particle/antiparticle in the j-th replica. The cyclic symmetry of replicas can be diagonalised moving to a "replica Fourier space" where we identified j ∼ j + n and k ∼ k + n. In order to account for the correct (anti)commutation relations (and hence S-matrix) of bosons and fermions, the index k should be either integer or half-integer according to the following rules • for bosons (S(ϑ) = 1) k = 0, . . . , n − 1; • for fermions (S(ϑ) = −1) k = − n−1 2 , . . . , n−1 2 .
Actually, different prescriptions can be found in the literature for the diagonalisation in the replica space [20,45,46,84,85]. These prescriptions are all related to each other by unitary transformations and in the end the diagonal modes simply get a phase. Since here we are only interested in the absolute value squared of the form-factors, this issue is irrelevant and any prescription would give the same result. The advantage of this approach is that the modified twist field T α n does not mix different Fourier modes and, consequently, its non-vanishing 2-particle form factors are just These FF have to satisfy the following bootstrap equations which obviously differ from the equations given in [71] for the conventional twist field. The only difference is the additional phase shift e ±iα/n , as a consequence of the fact that T α n introduces a phase e iα/n between each couple of consecutive replicas. (In principle, other choices for the phase of the single field are possible with the constraint that the total phase is e iα ; however, these different choices lead to different Fourier modes (104) for T α n .) For simplicity we will focus on F α k,+ for α > 0, since the other cases are obtained by symmetry. Given the form factors (105) in Fourier space, the physical ones are obtained as sum over the Fourier modes. In particular one has

Complex bosons
A solution of the bootstrap equations for the complex boson is The sum over the modes (109) is easily done using to get which corresponds to the result (97) by standard methods. This is the physical solution only for α > 0, while the one for α < 0 can be obtained by symmetry. One can derive the same expression starting from the ansatz compatible with the monodromy equations, and choosing C 0 , C 1 such that the poles are This result is valid also when analytically continued to n → 1, and provides the FF of V α (nonzero if α = 0). Indeed, in that limit, the poles at ϑ = iπ, 2inπ − iπ collapse together with an additional zero so that a single pole is left at ϑ = iπ, as a consequence of (116) This leftover pole is exactly the one for the single replica model, namely

Dirac fermions
Similar calculations can be done for free fermions, which lead to F α k,+ (ϑ) = i T α n sin α 2n + πk n e αϑ/2πn+kϑ/n cosh ϑ 2 . (118) The sum over the modes is performed using to obtain and this provides an alternative derivation of the FF in Eq. (69).
As for the case of complex boson, this solution can be also obtained making the ansatz and fixing C 0 , C 1 compatibly with the kinematical residues.

U(1) charged moments in free massive QFTs
In this section, we first present some basic and elementary facts about the U (1) charged moments Z n (α) for free theories with U (1) symmetry. Exploiting the QFT scaling form, some of our results are valid for arbitrary massive QFTs with U (1) symmetry with free bosonic/fermionic type UV limiting CFT as well. We restrict our analysis to a subsystem composed of a single interval and the full quantum system is prepared in its ground state. In this setting, the charged moments and later on entropies as well can be calculated from the two-point functions of the U (1) composite twist fields, which we explicitly calculate in the following. Specifying the subsystem as an interval , the charged moments are written as where ε is the UV regulator, ζ α n is a normalisation constant for the charged moments and d α n is the scaling dimension of the composite twist field d α n = 2∆ T α n = 1 12 n − n −1 + α 2 (2π) 2 n , Dirac, To keep track of the leading -dependence as well, we rewrite the charged moments as Using the standard two-particle approximation, the scaling function of the two-point correlation, H α,n (m ), for generic n can be written as due to the expansion of the 2-point function In the above formula f α D/B (ϑ, n) is implicitly defined as To ease our notations, we introduced the index κ which is either D or B referring to the Dirac fermion and complex boson respectively. The calculation of the von-Neumann entropy (and of the corresponding charged moments) requires the correct analytic continuationf α κ (ϑ, n) of the function f α κ (ϑ, n), in such a way to justify the treatment of n as a continuous variable. While for any integer n ≥ 2, f α κ (ϑ, n) =f α κ (ϑ, n), this equality is no longer true at n = 1:f α κ (ϑ, 1) equals f α κ (ϑ, 1) for κ = D, B, any α and ϑ = 0. In other wordsf α κ (ϑ, 1) is not a continuous function in ϑ. An important consequence of this fact is that the derivative of f α κ (ϑ, n) wrt. n contains a δ(ϑ) contribution. See appendix A for more details.
In fact for the leading spatial dependence, we can use the small ϑ behaviour of the f α κ (ϑ, n) derived also in appendix A. Plugging the expansions in Eqs. (162) and (165) into the integral in Eq. (126) and using [21] we immediately have for n > 1 and for both κ = D, B. Notably the leading order of H 2pt α,n (m ) is α-independent. As already mentioned for the limit n → 1 and for the derivative of f α κ (ϑ, n) wrt. n, we have to use the proper the analytic continuationf α κ (ϑ, 1) of f α κ (ϑ, n). Using the Taylor expansion in ϑ of the analytic parts off α κ (ϑ, 1) (that can be found in the appendix as Eqs. (163) and (166)), the leading behaviour is determined by the constant 2 sin 2 α 2 . Therefore we have the κ-independent expression for H 2pt manifesting the discontinuous behaviour compared to Eq. (129). For the derivative of H 2pt α,n , we have again for both statistics κ. Here the first term comes from the δ(ϑ) in lim n→1 ∂ ∂nf α κ (ϑ, n) and in the second line we also used that K 0 (z) ≈ e −z π 2z . Summarising, the leading spatial dependence of charged moments for n ≥ 2 can be written as where κ = D, B, ζ α κ,n is a normalisation constant, ε is a UV regulator and the VEV of the composite twist fields is calculated in appendices B and C for fermions and bosons, respectively. For n = 1, and for the derivative These formulas are the final results of this paper, which are derived based solely on the form factor bootstrap. They perfectly match with the results of Ref. [20] derived by completely different means.

U(1) symmetry resolved partition functions and entropies
For free theories, the symmetry resolved partition functions and entropies have been already calculated in Ref. [20] from the charged moments that were equivalent to those in the previous section. It is however useful to recall some steps of this derivation within the notations of this paper, because the computation of the SR entropies from charged moments is non-trivial.
Let us start recalling the definition of the symmetry resolved partition functions (7) in terms the charged moments (6): To perform the Fourier transform of Eqs. (132)-(134), the knowledge of both T α κ,n and ζ α κ,n is required as well. Although we computed the VEV of the composite twist-fields (see appendices B and C), which is a universal quantity (once the UV normalization is fixed), there is no general recipe to obtain ζ α κ,n since it is non-universal and its knowledge does not rely on QFT techniques. Following Ref. [20], we can write the logarithm of the charged moments (132)-(134) for both fermions/bosons as where ln Z The -independent contributions of order O(1) (in the limit → 0) are neglected explicitly and correspond to non-universal quantities.
From Eq. (136) the symmetry resolved von Neumann (n = 1) and Rényi entropies can be straightforwardly computed following [20]. We only report here the leading and sub-leading contributions to the entropies as for n ≥ 1. The total Rényi and von Neumann entropy has the following large behaviour S n c 6 n − n −1 log(mε) −1 − n n − 1 e −2m 4πm + . . . n > 1, where the non-universal constants are again neglected explicitly and c = 1, 2 for complex fermions and bosons, respectively. We recall that the result for S 1 is already known from Ref. [46]; in particular we get − 1 4 K 0 (2m ) which is the double of what was found in the Ising field theory because two particles are present in the theory. In contrast, for n > 1 the n-th entropy depends strongly on the theory under consideration beyond the leading order reported above. At leading order in the limit → 0, one observes the equipartition of entanglement, namely that S n (q A ) does not depend explicitly on q A . Nevertheless, the equipartition is broken explicitly whenever ε = 0. A careful analysis has been already performed in [20], where it has been shown that the terms which break explicitly equipartition maybe be written as a power series in 1 log(mε) −1 . Finally we want to mention that the total von Neumann entropy can be written as [86] where p(q A ) = Z 1 (q A ) equals the probability of finding q A as the outcome of a measurement of Q A . The contribution S c is called the configurational entanglement entropy and measures the total entropy due to each charge sector (weighted with the corresponding probability) [9,87] 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 [9,33,88,89]. In Eq.
(139) the log-log term is necessary in the SR quantity to cancel the same contribution to the total entropy coming from S f .

Conclusions
In this paper we applied the 1+1D bootstrap approach to compute the form factors of the composite branch-point twist fields which are directly related to symmetry resolved entropies. The technique was initiated in Ref. [21] for discrete symmetries and we generalised here to a U (1) conserved charge. For simplicity, we focused on free theories, namely the free massive Dirac theory and the free massive complex boson theory, both of which admits a U (1) symmetry. The motivation for the study of these free theories is twofold. In fact, despite the absence of interactions, the form factors of the corresponding composite branch-point twist fields are highly non-trivial. As a consequence, these calculations serve both as a warm-up toward interacting theories (without many technical complications like non-diagonal scattering etc.) and as a reference point to test future results for integrable models in the non-interacting limit.
We determined all the form factors of the composite field using two different methods: standard form factor bootstrap approach and diagonalisation in replica space. For the Dirac theory, form factors with higher particle numbers turned out to be Pfaffians of the 2-particle ones. For the complex boson, form factors with many particles can be obtained using Wick's theorem and the 2-particle form factor. The form factors for odd particle numbers and for non-neutral (with respect to the U (1) charge) particle configurations are zero. Our solutions have been tested using the ∆theorem sum rule. Furthermore, we recovered known particular limits: when the flux corresponding to the U (1) charge is zero, the conventional branch-point twist field form factors are re-obtained, whereas for a single replica we got the conventional U (1) twist field from factors. In addition, we also determined the exact vacuum expectation values for the composite fields.
Although the main goal of this work was the complete determination of the non-trivial form factors, we also presented the asymptotic results for the charged moments, the charged partition functions, and, eventually, the symmetry resolved entropies. We used the two-particle approximation for the 2-point correlation function of the composite fields, which gives accurate results for the charged moments when the considered subsystem is a long interval embedded in the ground state of an infinite system. For n = 1, a non-trivial analytic continuation is also necessary. We showed that the charged moments obtained from form factor sums equal those in Ref. [20] obtained by completely different means. Following Ref. [20], we finally re-derived the charged partition functions and the symmetry resolved entropies, including sub-leading corrections that breaks the equipartition of entanglement.
Our results are the starting point to systematically compute higher and higher corrections to the symmetry resolved entanglement and to consider subsystems consisting of disjoint intervals as well.
These directions can be an interesting object of future studies. Finally, the application of these form factor techniques to interacting integrable field theories with U (1) symmetry is on the way and will be presented in forthcoming publications.

Acknowledgments
We are grateful to Giuseppe Di Giulio, Cecilia De Fazio, and Olalla Castro-Alvaredo for useful discussions. All authors acknowledge support from ERC under Consolidator grant number 771536 (NEMO).
A Analytic continuation for f α D/B (ϑ, n) The analytic continuation of the quantity f α κ (ϑ, n) with κ = D, B is a subtle issue. The corresponding quantities f (ϑ, n) for the standard branch-point twist field of the Ising and sinh-Gordon models were carefully analysed in Ref. [46] where it was shown that the analytic continuationf (ϑ, n) with domain n ∈ [1, ∞) can be defined from f (ϑ, n) for n = 2, 3, .... It turned out thatf (ϑ, n) = f (ϑ, n) for integer n ≥ 2, but, for n → 1, f (ϑ, 1) = 0 everywhere except in the origin, where it converges to 1 2 . (To be more precise,f (0, n) as well as f (0, n) are understood as lim ϑ→0f (ϑ, n) and lim ϑ→0 f (ϑ, n). The non-uniform convergence is present in lim n→1 lim ϑ →ϑf (ϑ , n) as a function of ϑ.) The convergence in n is therefore non-uniform, which results in a δ-function in the derivative lim This derivative is an important quantity as it governs the leading long-distance behaviour of the entanglement entropy. However, not all the analysis of Ref. [46] is necessary here to obtainf α κ (ϑ, n): only some of the basic ideas from [46] are needed and presented here accordingly.
We first claim, that for the conventional branch-point twist field, the functionsf κ (ϑ, n) can be simply written as This can be justified following the logic of [46]: recalling the definition off D/B (ϑ, n) one can consider the following contour integral and treat j as a continuous parameter where the contour is a rectangle with vertices (− − iL, n − − iL, n − + iL, − + iL). This contour integral is zero as when L → ∞, the contributions of the horizontal lines vanish and, at least in free field theories, the vertical contributions cancel each other due to the periodicity of s κ (ϑ, z +n) = S 2 κ s κ (ϑ, z) and S 2 κ = 1 for κ = D, B. We can evaluate the integral (145) by the residue theorem; the poles are at the positions z = 1, 2, . . . , n − 1, at z = 1 2 ± ϑ 2πi , and z = n − 1 2 ± ϑ 2πi . Using the explicit values of the residues, we end up with from which Eq. (143) is inferred.
In a similar spirit, we can easily derive the corresponding derivatives forf α κ (ϑ, n). For both the complex boson and the Dirac theory, we can writef α κ (ϑ, n) as and consider the contour integral with the same contour as before. This contour integral is again zero for the same reasons as above: the contributions of the horizontal lines vanish and, in free field theories, the vertical contributions cancel each other because S 2 κ = 1. We can evaluate the integral (151) by the residue theorem. The poles are at the positions z = 1, 2, . . . , n − 1, at z = 1 2 ± ϑ 2πi , and z = n − 1 2 ± ϑ 2πi . Using the explicit values of the residues, we end up with for the Dirac theory and and lim n→1 ∂ ∂n for the complex boson theory.
These expressions lead to the main result of this appendix namely and for the Dirac theory and and for the complex boson.
The computation of further, sub-leading -dependent corrections is not addressed in this work, but can be straightforwardly achieved. To use Eq. (128), one eventually needs to expandf α κ (ϑ, n), The determination of the vacuum expectation value (VEV), i.e., the zero particle FF is generally a difficult task and for the branch-point twist fields, the VEV has been derived only for the Ising model, both for the conventional field [46] and for the composite one [21] and for the conventional twist field in the complex boson theory [83]. Strictly speaking the VEV is not unique, of course, but prescribing the standard conformal normalisation, i.e., requiring that it can be usually unambiguously defined. In this appendix, we derive the VEV T α D,n for the Dirac field using and modifying ideas taken from Refs. [21,45,46,91]. Following these ideas, we eventually proceed in a very similar way to the logic of section 5, that is, we essentially perform a diagonalisation in the space of replicas and exploit the factorisation of the field into Fourier components in the replica space.
Let us consider Dirac fermion fields and denote the one living on the jth replica by Ψ j . Similarly to section 5, we can consider the replica space spanned by vectors of the form (Ψ 1 , ..., Ψ n ) T and search for a matrix τ whose action in the space replica corresponds to the U (1) composite twist field. Since the total phase picked up by the twist field when turning around the entire Riemann surface is e iα , the transformation matrix τ can be obtained by multiplying the transformation matrix of the conventional fields [45] by e iα/n , as done in Ref. [20]. Accordingly the following representation is obtained The transformation matrix τ has to be diagonalised for the determination of the VEV [46]. The eigenvalues of τ can be written as e i2πk/n e iα/n with k ranging from −(n − 1)/2 to (n − 1)/2 for both even and odd n. Using the eigenvectors of τ new fermion fields Ψ k can be defined as well. These new fermion fields satisfy the canonical anticommutation relations The eigenvalues of the transformation τ obtained above are compatible with the four-point function of the composite U (1) branch-point twist field [20] Ψ † k (z)Ψ k (z )T α n (w)T α n (w ) T α n (w)T α n (w ) at the UV critical point. Eq. (172) means that performing a clock-wise turn on Ψ k (z ) around the twist field T α inserted at w, the correct factor of e iπ(2k+α)/n is accumulated. From Eq. (172) the factorisation of the composite U (1) branch-point twist field naturally follows and this fact makes it possible to compute the UV dimensions of the factorised components and eventually to determine the VEV in the massive theory. The factorisation of the composite U (1) branch-point twist field in our case can be written as and the action of T α k,n (w) is non trivial only on the Ψ k and Ψ † k fields. The (chiral) conformal dimension of T α k,n can be can be obtained from the relation [12,47,48] T k (z)T α k,n (w)T α k,n (w ) T α k,n (w)T α k,n (w ) where T k is the stress-energy tensor of the k components, as the stress-energy tensor can be shown to factorise into k sectors as well. The analysis of computing the conformal dimensions was performed in [20] leading us to the result It can now be checked that the total conformal dimension of the composite twist field agrees with the known value 1 24 n − n −1 + α 2 2(2π) 2 n because for even n we have where G(x) is the Barnes G-function. Hence, for the n-copy Dirac theory we have . (180) Using the integral representation of the Barnes G-function, we can rewrite the VEV as To determine the VEV of the composite U (1) branch-point twist field in the free complex boson theory, we can proceed in a similar fashion, and our eventual computation boils down again to writing the VEV of the composite branch-point twist field as a product of VEVs for conventional U (1) twist fields. Contrary to the case of the free Dirac theory, however, we now face an important subtlety when defining the VEVs V ϕ B in the complex boson theory. This theory is not compact and as a consequence the short-distance behaviour of the theory and the proper definition of the VEV are non-trivial [83], because of the presence of a zero mode. Consequently, the explicit value of the VEV is non-universal, in the sense that depends on the employed normalisation. This problem was carefully discussed in Ref. [83], where an expression for the VEV was actually proposed based on a natural regularisation of the fields and on the angular quantisation scheme. In the following derivation, we adopt this convention, which was already used in [83] to derive the VEV of the conventional branch-point twist field (in a similar calculation to what follows).
To proceed with our eventual derivation in the standard way, we first of all need again the transformation matrix τ , whose action in the space replica space (i.e. on the vector (Φ 1 , ..., Φ n ) T ) corresponds to the composite twist field can now be written as [20] τ = e i α The eigenvalues of τ are the roots of unity times e iα/n that is e i2πk/n e iα/n with k ranging from 0 to n − 1. Similarly to the Dirac case, one can introduce new Bose fields, which satisfy the canonical commutation relations [Φ k (x), Φ † k (x)] = δ k,k δ(x − x ), [Φ k (x), Φ k (x )] = 0 and [Φ † k (x), Φ † k (x )] = 0. In complete analogy with the Dirac case, the factorisation of composite twist field is inferred as [20] T α n (w) = n−1 k=0 T α k,n (w) , where the conformal dimension of each component is The total conformal dimension of the composite twist field agrees with the known value for any n, that is Based on the above analysis, the decomposition of the composite branch-point twist fields is assumed again and is rephrased as in terms of conventional bosonic U (1) twist field V ϕ k B , with ϕ k = k n + α 2πn and whose left and right conformal dimensions are exactly 1 2 k n + α 2πn 1 − k n − α 2πn . Assuming that this type of factorisation of the composite branch-point twist field also holds in the off-critical theory we can obtain its vacuum expectation value exploiting the results in Ref. [83] V ϕ B =N (me γ E ) ϕ(1−ϕ)/(2π) 2 exp − 2 πˆ∞ 0 dt sinh t ln(cosh t) (4t 2 + π 2 ) cosh 2 t π cos ϕ 2 sinh t 1 − ϕ π + +2t sin ϕ 2 cosh t 1 − ϕ π , with N = exp − 1 3 (γ E + ln 2) + 1 2 (1 − ln(2π)) −ˆ∞ 0 dt t e −t sinh t − t (e −t − 1)(cosh t − 1) where γ E is Euler's constant.