Entanglement and symmetry resolution in two dimensional free quantum field theories

We present a thorough analysis of the entanglement entropies related to different symmetry sectors of free quantum field theories (QFT) with an internal U(1) symmetry. We provide explicit analytic computations for the charged moments of Dirac and complex scalar fields in two spacetime dimensions, both in the massive and massless cases, using two different approaches. The first one is based on the replica trick, the computation of the partition function on Riemann surfaces with the insertion of a flux $\alpha$, and the introduction of properly modified twist fields, whose two-point function directly gives the scaling limit of the charged moments. With the second method, the diagonalisation in replica space maps the problem to the computation of a partition function on a cut plane, that can be written exactly in terms of the solutions of non-linear differential equations of the Painlev\'e V type. Within this approach, we also derive an asymptotic expansion for the short and long distance behaviour of the charged moments. Finally, the Fourier transform provides the desired symmetry resolved entropies: at the leading order, they satisfy entanglement equipartition and we identify the subleading terms that break it. Our analytical findings are tested against exact numerical calculations in lattice models.


Introduction
Symmetries are a pillar of modern physics. They can concern spacetime, as rotational or relativistic invariance, or can be internal symmetries, which do not touch the spacetime coordinates. Their exploration turned out to be a central theme in several fields ranging from elementary particles to the theory of phase transitions, from string theory to solid-state physics. One century ago, Emmy Noether proved that every symmetry of a physical system leads to a corresponding conservation law. For example, the conserved electric charge is the generator of the U (1) gauge symmetry of electromagnetism. Consequently, an evergreen research topic is the characterisation of how the presence of a symmetry influences the properties of a physical system. In particular, this manuscript addresses the question of how the entanglement splits into the different sectors of an internal symmetry.
The von Neumann entropy is the most successful way to characterise the bipartite entanglement of a subsystem A in a pure quantum state [1][2][3][4]. Given the reduced density matrix (RDM) ρ A of a subsystem A, the entanglement entropy is defined as (1.1) A related family of functions, known as Rényi entropies, is (1. 2) The essence of the replica trick is that the von Neumann entropy may be obtained as the limit n → 1 of Eq. (1.2). The reason of this way of proceeding is that for integer n, in the path-integral formalism, Trρ n A is the partition function on an n-sheeted Riemann surface R n obtained by joining cyclically the n sheets along the region A [5,6]. This approach, when applied to critical systems whose low energy physics is described by a (1+1) dimensional conformal field theory (CFT), leads to the famous scaling results [5][6][7][8][9][10][11] when the subsystem A is an interval of length embedded in an infinite one-dimensional system and is an ultraviolet cutoff. The possibility of measuring in an experiment the internal symmetry structure of the entanglement [12] went together with a new theoretical framework developed to address the problem [13,14]. Indeed, in Refs. [13,14] a simple generalisation of the replica trick has been proposed to relate the symmetry resolved quantities to the moments of ρ A on a modified Riemann surface: we refer to them as charged moments. Such technique allowed for the derivation of interesting results about the different symmetry-resolved contributions not only in CFTs, but also in the context of free gapped and gapless systems of bosons and fermions, integrable spin chains, disordered systems and many more (the interested readers can consult the comprehensive literature on the subject ).
The main goal of this manuscript is to investigate the field theoretical techniques for the computation of the charged moments in relativistic free two-dimensional quantum field theories (QFTs). The paper is organised as follows. In Section 2 we provide all the definitions concerning the measures of symmetry resolved entanglement and we briefly recall two tools for the computation of the entanglement in QFT, i.e. the twist fields [7,[36][37][38][39] and the Green's function technique in the replica space [40][41][42]. Our main findings are reported in Sections 3, 4, and 5: in the first we employ the twist fields to compute the charged moments both in the massless and in the massive context (in the limit m 1). These results are extended in sections 4 and 5, where we write down the explicit form of the charged moments for arbitrary m and provide analytic asymptotic expansions valid for large and small m . These outcomes are the starting point for the computation of the symmetry resolved entanglement entropies. Numerical checks for free fermions and bosons on the lattice are also provided as a benchmark of the analytical results. In Sec. 6 we give a result for the charged moments of a free massive scalar theory across a hyperplane in generic Euclidean dimensions d. We draw our conclusions in Section 7. Three appendices are also included: they provide details about the analytical and numerical computations.

Symmetry resolution and QFT techniques
In this section, we provide an overview of notions about symmetry resolved entropies we use throughout the manuscript. We also recall some standard techniques used to compute the entanglement entropy in relativistic free QFT.

Symmetry resolved entanglement
We consider a system with an internal U (1) symmetry and its bipartition into two subsystems, A and B. The charge operator Q is the generator of the symmetry and we assume it obeys Q A ⊕ Q B = Q, where Q i is the charge in the subsystem i. If the system described by the density matrix ρ is in an eigenstate of Q, then [ρ, Q] = 0. Tracing out the degrees of freedom of B, we obtain the RDM of A, ρ A = Tr B ρ. Hence, taking the trace over B of [ρ, Q] = 0, we find that [ρ A , Q A ] = 0. This means that ρ A has a block-diagonal structure where each block corresponds to an eigenvalue q of Q A . The density matrix ρ A (q) corresponding to an eigenvalue q is obtained by projecting ρ A onto the eigenspace of Q A with fixed q, as induced by the projector Π q . Therefore we can write where p A (q) is the probability of finding q in a measurement of Q A in the RDM ρ A , i.e. p A (q) = TrΠ q ρ A . Within this convention, the density matrices ρ A (q) of different blocks are normalised as trρ A (q) = 1. The amount of entanglement shared by A and B in each symmetry sector can be computed through the symmetry resolved Rényi entropies, defined as S n (q) ≡ 1 1 − n ln Trρ n A (q). (2. 2) The limit n → 1 gives the symmetry resolved entanglement entropy, i.e. 3) The total von Neumann entanglement entropy associated to ρ A in Eq. (2.1) splits into [43] S 1 = q p(q)S 1 (q) − q p(q) ln p(q). (2.4) The two contributions are known as configurational and fluctuation (or number) entanglement entropy, respectively [12]. The configurational entropy is also related to the operationally accessible entanglement entropy of Refs. [27][28][29], while the number entropy is the subject of a substantial recent activity [12,17,[30][31][32]. The calculation of the symmetry resolved entropies by the definition (2.2) requires the knowledge of the entanglement spectrum of ρ A and its resolution in the charge sectors. However, this is a difficult task, especially for an analytic derivation. As first proposed in [13], we can rather focus on the charged moments of ρ A Z n (α) ≡ Trρ n A e iQ A α , (2.5) with Z 1 (α = 0) = 1, being Trρ A = 1. Similar charged moments have been already considered in the context of free field theories [44], in holographic settings [45,46], as well as in the study of entanglement in mixed states [47,48]. In this specific case, the charged moments are not the main goal of our computation, but they represent a fundamental tool, because their Fourier transforms are the moments of the RDM restricted to the sector of fixed charge q [13], i.e.
(2.6) (Here we assumed Q to be the generator of a U (1) symmetry and q ∈ Z.) Finally the symmetry resolved entropies are obtained as

Replica method and QFT
In the following sections we will mainly deal with a free fermionic field theory and with a complex scalar one, whose Euclidean actions are given, respectively, by where we employ complex coordinates (z,z) for the 2D spacetime. In S D the fields ψ R/L are the chiral (right-moving R) and anti-chiral (left-moving L) components of the Dirac fermion. In S S the field ϕ is a complex scalar. The actions in (2.8) exhibit a U (1) symmetry, i.e. a symmetry under phase transformations of the fields given, respectively, by By Noether's theorem, this continuous symmetry transformation leads to a conserved quantity, which is the charge Q we introduced before. The actions (2.8) played the role of the simplest massive quantum field theories to study the properties of the entanglement entropy. In the same spirit, they also represent the natural starting point for the field theoretical investigation of the charged moments and, as a consequence, of the symmetry resolved entropies.
In what follows we describe two powerful methods to calculate the entropy for free fields. Starting from the replica trick described in the introduction, the first approach is based on a particular type of twist fields in quantum field theory that are related to branch points in the Riemann surface R n . We denote them by T n andT n . Their action, in operator formalism, is defined by [7,36,37] where z 1 and z 2 are the endpoints of A, z ∈ A and i = 1, . . . , n with n + 1 ≡ 1. The two-point function of the twist fields directly gives [7] Trρ n A ∝ T n (z 1 )T n (z 2 ) . (2.11) In conformal invariant theories (e.g. when the mass terms in the actions (2.8) vanish) the two-point function of twist fields is fixed by their scaling dimension, leading to Eq. (1.3).
In some instances, a simplification arises by the diagonalisation in the replica space: the n-sheeted problem can be mapped to an equivalent one in which one deals with n decoupled and multivalued free fields, generically referred asφ k . Thus, also the twist fields can be written as products of fields acting only onφ k , denoted as T n,k andT n,k . The total partition function is a product of n partition functions, ζ k , each one given by (up to unimportant multiplicative constant) The second approach is the one used in Refs. [40,41] for a fermionic and a complex scalar theory, respectively: it also relies on mapping the problem from the determination of the partition function on R n , to the computation of n partition functions of a free field on a cut plane. However, the difference with respect to the previous approach is that each ζ k is not computed as a two point-function of twist fields, but using the relation between the free energy and the Green's function of each sector k. Denoting by G D the Green's function for the Dirac field and by G S the one for the scalar (in each sector k of the n copies), they are related to the corresponding partition function ζ k by, respectively, The strategy of Refs. [40,41] was to exploit the rotational and translational symmetry of the Helmholtz equations satisfied by G D and G S and analyse their behaviour at the singular endpoints of the cut A so to determine the right hand sides of the above equations. The final expressions for ζ k can be expressed in terms of the solution of second order non linear differential equations of the Painlevé V type. Here we only report the final results for the Rényi entropies of free fields in the limit m → 0 [40,41] (2.14) These formulas have been obtained in the scaling regime with t = m fixed, in the conformal limit m → 0 and after taking the limit of large . Eq. (2.14) shows the leading mass corrections to Eq. (1.3) for the theories in Eq. (2.8) strongly depend on the statistics of the particles. The leading mass correction vanishes for a Dirac field, while it is singular (like ln(− ln(m))) for a Klein-Gordon field (both real and complex). In the literature, this infrared divergence is ascribed to the zero mode of the massless scalar theory [49,50].

Twist Field Approach
In this section we consider 1D critical and close to critical systems. We obtain a general exact result for the conformal invariant charged moments by exploiting the properties of some local operators known as modified or fluxed twist fields [44,45]. This result includes and generalises the ones in Ref. [13]. The same approach also provides the leading asymptotic behaviour of the charged moments for (free) massive field theories.

Modified Twist Fields
In a generic QFT, the replica trick for computing Z n (α) defined in Eq. (2.5) can be implemented by inserting an Aharonov-Bohm flux through a multi-sheeted Riemann surface R n , such that the total phase accumulated by the field upon going through the entire surface is α [13]. The result is that Z n (α) is the partition function on such modified surface, that, following Ref. [13], we dub R n,α . In QFT language, the insertion of the flux corresponds to a twisted boundary condition. This boundary condition fuses with the twist fields at the endpoints of the subsystem A resulting into two local operators T n,α andT n,α . These are modified versions of the standard twist fields T n andT n which take into account not only the internal permutational symmetry among the replicas but also the presence of the flux. Thus, the partition function on R n,α is determined by their two-point correlation function, that is the main object of interest in this section. As already mentioned in section 2.2, rather than dealing with fields defined on a non trivial manifold R n,α , it is more convenient to work on a single plane with a n-component field where φ j is the field on the j-th copy (here the field φ j generically refers to either a scalar field ϕ j or a chiral Dirac one ψ j ; the same applies toφ k that we are going to introduce soon). Upon crossing the cut A, the vector field Φ transforms according to the transformation matrix T α where f = 1 for free Dirac fermions and f = 0 for free complex scalars. When α = 0 we recover the usual transformations for the fields across the different replicas [42]. The matrix T α has eigenvalues f = 0 : λ k = e i α n e 2πi k n , k = 0, . . . , n − 1, By diagonalising T α with a unitary transformation, the problem is reduced to n decoupled fieldsφ k in a two dimensional spacetime. Thus, the total partition function is a product of the partition functions for each k and the twist fields can be written as products of fields each acting on a differentφ k , i.e.
with T n,k,αφk = δ k,k e iα/n e 2πik/nφ k andT n,k,αφk = δ k,k e −iα/n e −2πik/nφ k . Since the partition function on R n,α can be written as the two-point function of the modified twist fields, from  When dealing with a CFT (e.g. when m = 0 in (2.8)) T n,k,α andT n,k,α are primary operators and their two-point function is fixed by conformal invariance to be where (see the Appendix A) (3.7) Let us stress that, in order to have operators with positive conformal dimension, the phase that bosons pick up going around one of the entangling points must be 0 < k n + α 2πn < 1. This can be achieved, since α ∈ [−π, π], by trading α with |α| when we deal with scalar field theories.
Using Eqs. (3.5), (3.6) and (3.7) the logarithm of the partition function on R n,α reads (3.8) The charged moments for the free massless Dirac field theory (f = 1) have been already worked out in the literature with different techniques [13,16,44,45]. Instead, the charged moments for a free massless complex scalar field (f = 0) represent a new result (actually in Appendix A of [45] a result consistent with (3.8) has been obtained using the heat kernel techniques). Let us stress that the presence of a flux in the Riemann surface changes some features of the twist fields in CFT: they remain primary operators (see Appendix A for details), but they do depend on the theory and are not anymore identified only by the central charge (see also [13]).

Massive field theory and flux insertion
In this section we compute the charged moments Z n (α) of a massive relativistic 2D QFT on the infinite line for a bipartition in two semi-infinite lines. Thus, we follow the same logic as in [5] (i.e. the continuum version of Baxter corner transfer matrix approach [51] for the reduced density matrix [5,[52][53][54]), which in turn parallels the proof of the c-theorem by Zamolodchikov [55]. The results of this section are not limited to free theories but hold for generic massive relativistic QFT. Exploiting the rotational invariance about the origin of the Riemann surface R n,α , the expectation values of the stress tensor of a massive euclidean QFT in complex coordinates, T ≡ T zz ,T ≡ T * zz , and the trace, Θ ≡ 4T zz , have the form where Θ 1,α=0 is a non-vanishing constant measuring the explicit breaking of scale invariance in the non-critical system, while T 1,α=0 and T 1,α=0 both vanish. These quantities are related by the conservation equations of the stress-energy tensor (4∂zT + ∂ z Θ = 0) as The conservation equations as well as the rotational invariance are preserved in the presence of the flux α because the Riemann surface R n,α can be thought simply as a complex plane with the insertion of two modified twist fields, as discussed in the previous subsection. Both F n,α and G n,α approach zero for |z| ξ, while when |z| ξ, they approach the CFT values given by the conformal dimension of the modified twist field, ∆ n (α) [13]. Hence we have G n,α → 0, (3.11) and in particular for a massive Dirac field theory (f = 1) and for a complex massive scalar theory (f = 0), using the conformal weights (3.7), we have (3.12) Changing variable to R 2 = zz, we can rewrite the Eq. (3.10) as The corresponding integrated form using the boundary conditions in Eq. (3.14) Taking into account the normalisation of the stress tensor, the definition of G n,α in Eq. (3.9) and that Rn Θ n,α dR 2 corresponds to the variation of the free energy wrt a scale transformation (the mass m in this case) per each sheet of the whole n-sheeted surface, the left hand side of (3.14) is equal to − 2 n m∂ m ln Z n (α). We can therefore integrate this equation at fixed n and α to obtain ln Z n (α). The additive non-universal integration constant can be absorbed in a UV cutoff n,α that consequently depends both on the Rényi index n and the parameter α (consistently with the lattice results in Ref. [17,20] for massless theories). Finally we get 16) or specialising to free Dirac (f = 1) or complex Klein-Gordon (f = 0) fields We should mention that ln Z n (α) for the Klein-Gordon field matches the continuum limit of a chain of complex oscillators obtained through the corner transfer matrix approach [20]. We can specialise Eq. (3.16) to a Luttinger liquid with parameter K, whose underlying field theory is a c = 1 CFT equivalent to a massless compact boson. In this case, one can use the results found in [13] for the conformal dimension of the modified twist field, obtaining which for K = 1 gives the result found for fermions in Eq. (3.17), as it should.

From charged moments to symmetry resolved entropies
Performing the Fourier transforms of the charged moments above, one obtains symmetry resolved moments and entropies. For the Luttinger liquid, which includes Dirac Fermions at K = 1, the α dependence of the leading term is the same as in the massless cases. Hence the analysis is identical to the one of Refs. [13,17] and so we will just sketch the results here. The charged moments (ignoring for the time being the dependence on n and α of n,α in (3.19)), are and hence symmetry resolved entropies with S n the total entropy. Exploiting the knowledge of Z 1 (q) in (3.20) we also easily get the number or fluctuation entropy that in the sum for the total entropy cancels exactly the double logarithmic term in Eq. (3.21).
Although the massive complex boson has been already investigated in Ref. [20], there another critical limit has been taken. Here we are interested in the Fourier transform of Eq. (3.18). In the saddle point approximation, we can neglect the term ∝ α 2 in Eq. (3.18) and the Fourier transform is and hence with S n the total entropy. Also in this case, one easily derives the number entropy from Z 1 (q) obtaining again, at the leading order, the double logarithmic term in S n (q) in Eq.

The Green's function approach: The Dirac field
In this section we derive the charged moments for a massive Dirac field for arbitrary mass and then consider the limits of small and large mass. In Sec. 3.1 we showed that Z n (α) can be written as product of partition functions ζ a on the plane with proper boundary conditions along the cut A, explicitly given bỹ Hence we have Let us introduce the auxiliary universal quantities that, using (4.2), allow us to write the logarithmic derivative of the partition function in R n,α as For n = 1, the function c n (α) is the analogue of Zamolodchikov's c-function [55] in the presence of the flux α. The cutoff , in analogy to what discussed in Sec. 3.2 depends on both α and n, although we almost always omit such a dependence for conciseness. As already discussed in section 2.2, the key observation of this approach relies on the identity between the partition function ζ a and the Green's function in the same geometry (see Eq. (2.13)). Through this relation, the function w a has been already obtained for generic values of a for the massive Dirac fermion [40,42].
The method that we just reviewed provides exact results for the charged moments of a free Dirac field. Indeed, in Ref. [40] it has been shown that the function w a defined in (4.3) can be written as This equation can be straightforwardly solved numerically with any standard algorithm for ordinary differential equations, once we impose the boundary condition as where κ D (a) = − ln 2 + 2γ E + 1 2 (ψ(a) + ψ(−a)), with ψ(z) ≡ Γ (z)/Γ(z) the digamma function and γ E the Euler-Mascheroni constant.
Plugging the numerical solution of the differential equation (4.6) into Eq. (4.3), we obtain the universal constant c n (α). Then, with the further integration (4.4), the desired ln Z n (α) is found to the price of introducing the non-universal cutoff . As examples we report in Fig. 1 the plots of the resulting c n (α) for few values of α and n as functions of t = m . In the figure we also compare our exact solution with the numerical results obtained from a lattice discretisation of the free Dirac theory (see Appendix C for details). The agreement is excellent. We stress that in Fig. 1 there is no free parameter in matching analytical and numerical data for c n (α) (as a difference compared to Z n (α)).
The method we just outlined provides exact results for the desired charged moments and, by Fourier transform, the symmetry resolved entropies. However, the procedure is completely numerical and we would appreciate an analytic handle on the subject. While in general this is not feasible, the limits of small and large t are analytically treatable, as we are going to show.

The expansion close to the conformal point m = 0
Here we use the methods just introduced to derive an asymptotic expansion of the charged moments close to the conformal point, i.e. for t = m → 0. In this limit, the expansion of the function w a (t) has been worked out in Ref. [40], obtaining where we omitted the dependence on a of κ D . In order to compute c n (α) through Eq. (4.3), we again set a = k n + α 2πn and we compute the following sums where ω n = π 2 nΩ n (α) and λ n = π 2 nΛ n (α) so that the remainder functions ρ (4.11) Eq. (4.11) can be now integrated analytically, getting where we defined and ρ z n (α) is defined so that ρ z n (α) = O(α 4 ). Notice, as we already stressed a few times, that in Eq. (4.12) the cutoff comes as an additive constant of integration and it generically depends on both n and α.
Eq. (4.12) represents our final field theoretical result for the charged entropies. We wish to test this prediction against exact lattice computations obtained with the methods in Appendix C. However, in order to perform a direct comparison with lattice data, we have to take into account the additional non-universal contribution coming from the discretisation of the spatial coordinate, i.e. the explicit expression for the cutoff in Eq. (4.12) that, as already mentioned, does depend on α and cannot be simply read off from the result at α = 0. We assume here (as Eq. (4.12) suggests at leading order) that the cutoff does not depend on the mass; consequently we can use the exact value for m = 0 [17] obtained with the use of Fisher-Hartwig techniques. The final result of Ref. [17] may be written as , (4.16) and in particular we will use . (4.17) In Ref. [17] it has been shown that the cutoff in (4.16) is very well described by the quadratic expansion in α and higher corrections O(α 4 ) are negligible for most practical purposes. In Fig. 2 we report the numerical data for the charged moments with the insertion of a flux α for two values of α and m with n = 1. The data are well described by the theoretical prediction (4.12) with the cutoff (4.16). Finally, in order to have a test of the prediction (4.12) that does not rely on an independent lattice calculation we can consider the difference between the charged entropy at finite t (i.e. finite mass) and the massless one. Specifically we consider in which both the cutoff and dependences cancel and it becomes a universal function solely of t (closely related to c n (α)). The results for δZ(α, t) are reported in Fig. 3. The agreement of the numerics with the prediction (4.12) is perfect for small t. Furthermore, the differences emerging for larger t are correctly captured by the numerical exact solution of the Painlevé equation (4.6). The small discrepancies visible in the figure are just finite size effects that are stronger for larger values of n and α.

From the charged moments to symmetry resolution
We are now ready to study the true symmetry resolution by performing the Fourier transform of Z n (α). In this Fourier transform we ultimately use a saddle-point approximation in which Z n (α) is Gaussian and hence we truncate hereafter Eq. (4.12) at quadratic order in α. Consequently, the charged partition function can be well approximated as where we consistently approximated the cutoff at quadratic level and used the lattice cutoff with γ(n) given in Eq. (4.17). A different cutoff just leads to a different additive constant in b n (i.e., a different definition of h n ), but we will use its precise form only for the comparison with numerics and so all the following formulas are completely general. Now we can compute the Fourier transform (2.6) that reads When → ∞, we can perform the integral by saddle point approximation and the integration domain can be extended to the whole real line. We end up in a simple Gaussian integral, obtaining 2bn( ,t) . We check Eq. (4.22) against numerical computations in Figure 4 focusing on n = 1 and the agreement is perfect. We test both the scaling with t = m for fixed q and at fixed t as a function of q. Now we are ready to compute the symmetry resolved Rényi entropies. From the definition (2.7) we have , (4.23) where S D n is the total n-th Rényi entropy for the Dirac fields (cf. Eq. (2.14) up to O(t 2 ln 2 t)). We can further expand the above equation for → ∞ since b n ( , t) diverges logarithmically, obtaining  Figure 5. Symmetry resolved entanglement entropies for a few different values of q and n as functions of . The field theory prediction is tested against exact lattice computations. The agreement with Eq. (4.24), that includes lattice effects, is remarkable. For large |q|, the approximation at the order q 2 is no longer sufficient and neglected corrections to the scaling become important, as well known for the massless case [17]. and ln κ n = −π 2 (h 1 + nh n ) 2 . (4.26) The above formula is valid also for the symmetry resolved Von Neumann entropy taking properly the limits of the various pieces as n → 1. By construction, the total entropy, S D n , coincides with the one obtained in [40] for the massive fermions in the conformal limit up to O(t 2 ).
Let us critically discuss the result in Eq. (4.24). The leading terms for large (up to O((ln ) −2 )) do not depend on q and they are given by the total entropies S D n in Eq. (2.14). We then conclude that at this order, the presence of the mass does not break entanglement equipartition found in conformal field theory [16]. The first term breaking equipartition is at order O((ln ) −2 ) and its amplitude is governed by the constant h n defined in Eq. (4.20). This constant gets contributions both from the non-universal cutoff and from the mass; the two contributions have the same analytic features. In Fig. 5 we test the accuracy of our total prediction against exact lattice numerical calculations. The agreement is remarkable for small values of |q|, but it worsens already at q = 2; this does not come as a surprise since the same trend was already observed in the massless case [17]. Such discrepancies are entirely due to corrections of order o(q 2 ) and are expected to reduce as gets larger. The drawback of the data reported in Fig. 5 is that universal field theory mass contributions  Figure 6. Subtracted symmetry resolved von Neumann entropy δS 1 (q, t) ≡ S 1 (q, t) − S 1 (q, t = 0) for q = 0 (left) and q = 1 (right) as a function of t (fixing = 600 and varying m). This subtracted quantity highlights the mass dependence of symmetry resolved entropies. The continuous lines are just the difference of the same subtracted entropies as obtained from the field theory expansion (4.24). and the lattice non-universal terms are mixed up and the latter are, by far, the largest one. It is then very difficult to observe the dependence on the mass in these plots. An effective and easy way to highlight the role of the mass is to subtract from the symmetry resolved entropies their value for the massless case, i.e. considering the numerical evaluation of the δS n (q, t) ≡ S n (q, t)−S n (q, 0). Such subtracted entropies for n = 1 and q = 0, 1 are reported in Fig. 6, showing that the entropy is a monotonous decreasing function of t (and hence of m at fixed ).

The long distance expansion.
In this subsection we move to the analysis of the charged and symmetry resolved entropies in the limit of large t. The most effective way to proceed is, following Ref. [40], to employ in Eq. (4.6) a boundary condition for t → ∞, that takes the form [40] v a (t) ∼ 2 π sin(πa)K 2a (t), (4.27) where K 2a (t) is the modified Bessel function of the second kind. This is the starting point for a systematic expansion for large t of the solution v a (t) of the differential equation (4.6).
Plugging the resulting expansion into the integral (4.5) for w a (t), we get w a (t) = −e −2t sin 2 (aπ) Summing over a = k n + α 2πn , we obtain the long distance asymptotic expansion for the universal factor c n (α) 29) and for n = 1 This is consistent with the exact result c 1 (0) = 0 coming from the normalisation of the reduced density matrix. For α = 0, Eq. (4.29) reproduces the known results [40]. In Fig. 7 we report the numerical exact solution of the Painlevé equation (4.6) for c n (α); we focus on n = 1, 2 and plot c n (α) as a function of t. For large t, the solutions perfectly match the asymptotic expansions (4.29) and (4.30) (for completeness we also show the small t expansion in Eq. (4.11)). Let us emphasise the presence of a discontinuity in c n (α) for n → 1 as a function of n: it is due to the non-commutativity of the limits n → 1 and t → ∞, as well known and discussed at length in the literature for α = 0 [37,40]. We show here that the presence of α does not cancel such a discontinuity, although for α = 0 the leading term is of the same order e −2t .
The charged entropy is simply given by the integral At large t, the function c n (α) goes to zero exponentially in t for any n; hence the charged entropies approach asymptotically a finite value for large . This saturation value is determined entirely by the infrared physics, i.e. by the value of c n (α) at small t, indeed  Once again, Z n (α) are not continuous functions of n close to n = 1 (as it was already known for α = 0, see [40]) and, above all, the correction of ln Z n (α) does not depend on α for n = 1 at this order. Subleading corrections to Eq. (4.33) can be straightforwardly and systematically worked out, but they are not illuminating, although they do depend on α also for n = 1.
For n = 1, since the leading correction does not depend on α, the Fourier transform is not affected and the symmetry resolved moments with n = 1 just get a multiplicative correction to Z n (q) in Eq. (3.20) (so additive for the logarithm), given by For n = 1 the net effect of the sin 2 (α/2) term in Eq. (4.33) is to renormalise the variance with an exponential additive correction, i.e. the desired probability is The symmetry resolved Rényi entropies with n = 1 are straightforwardly obtained from Eq. (2.7). Indeed, plugging Eqs. (4.34) and (4.35) in (2.7), we get Such a result shows exact equipartition (at this order) which is a clear consequence of the simple form of (4.34). This is reminiscent of the exact results for integrable models studied in Ref. [20]. The limit n → 1 for the von Neumann entropy should be handled with care. We start rewriting Eq. (2.7) as We use this equation to obtain the entire correction in t due to the Bessel function and not only the leading exponential term (as done in Eq. (4.33)). The crucial computation is where a = k n + α 2πn . We can use the integral representation for the Bessel function .
(4.41) We now study the behaviour of F n,α (z) when n → 1. For z = 1, the limit n → 1 is singular. We can isolate this singularity using the polar variables (n − 1, z − 1) → (ρ cos θ, ρ sin θ) and expanding in the radial coordinate ρ. The result of this procedure is whose derivative with respect to n is Plugging this result in Eq. (4.40) and taking the derivative wrt n, we get which, once integrated in t according to Eq. (4.31), gives the full ultraviolet behaviour of ∂ n ln Z n (α) n→1 , i.e.
Plugging the above derivative into Eq. (4.37) finally yields The first two terms in (4.46) are respectively the leading and the subleading terms in the total entanglement entropy of a massive Dirac field, in agreement with the known results in Refs. [36,37,40]. The double logarithmic term appears only in the symmetry resolved result and, as already discussed in Eq. (3.21), it is related to the number entropy. The above derivation clearly highlight this correspondence. As for the Rényi entropy, at this order in ln m , there is perfect entanglement equipartition that will be broken by higher order terms.

The Green's function approach: The complex scalar field
In this section we present a derivation of the charged moments for a complex massive scalar by generalising to α = 0 the results obtained in [41,42]. In Sec. 3.1 we showed that Z n (α) can be written as product of partition functions on the plane with boundary conditions along the cut Ã ϕ k (e 2πi z, e −2πiz ) = e 2πiaφ k (z,z), a = k n + α 2πn , k = 0, · · · n − 1. that, using (5.2), allow us to write the logarithmic derivative of the partition function in R n,α as Even here, for n = 1, the function c n (α) is the analogue of Zamolodchikov's c-function [55] in the presence of the flux α.
As already discussed in section 2.2, the key observation of this approach relies on the identity between the partition function ζ a and the Green's function (see Eq. (2.13)). Through this relation, the function w a has been obtained for generic values of a also for bosonic free massive field theories [42]. As already found in Sec. 3.1 using twist fields, also this approach requires that 0 < a < 1 for the scalar theory (see [41] for details). Thus, in order to compute Z n (α), we will use the results in [42] setting a = k n + |α| 2πn for the complex Klein Gordon theory.
Here we consider the complex massive non-compact bosonic field theory with action given by (2.8) and mass m. The function w a with 0 < a < 1 defined in (5.3) can be written as [41] w a = − ∞ t yu 2 a (y)dy, (5.5) where t = m and u a is the solution of the Painlevé V equation The solution of Eq. (5.6) is showed in Fig. 8: the function w a for a generic value of t can be obtained solving numerically Eq. (5.6) with the initial condition for t → 0 where [42] In the figure we compare the exact result from field theory with numerical data for a chain of complex oscillators, obtained exploiting the techniques reviewed in Appendix C. We have a fairly good agreement between lattice and field theory, although for small values of α the agreement gets worse and one needs a larger and larger subsystem length on the lattice to match the continuum limit. This is not surprising, already in Ref. [20] it was shown for the massless case that the lattice results approach the CFT ones in a non uniform way. In the following we will further discuss this issue in the limits when we have an analytic handle on the problem.

The expansion close to the conformal point
In the conformal limit t → 0 we have that the solution of the Painlevé equation admits the expansion [41]   which diverges as π/|α|. Hence, since ln t grows very slowly with t, the true asymptotic behaviour is attained only for t e π/|α| . For smaller values of t, the mode k = 0 looks almost constant (∼ |α|/π) and similar to the leading term. Exactly for α = 0, the mode k = 0 diverges, so its inverse is just zero and it does not affect the calculation. It is then clear that the approach to the asymptotic behaviour cannot be uniform in α, as already observed numerically in Fig. 8.
This competition among the three terms makes difficult the analytical treatment of the last sum in Eq. (5.11), and, at the same time, the non trivial dependence on the cutoff (that we recall also depends on α and n) complicates the comparison with the numerics. For this reason, we consider only the leading term in Eq. (5.11), which, strictly speaking, is valid in the massless case and in the third regime above. Such a leading term is the same provided as in the twist field approach (cf. Eq. (3.8)) and coincides with some equivalent ones in literature [20,45]. The main advantage of Eq. (5.11) is that it clearly shows what are the problems one faces when considering only the leading term. The comparisons with the numerics are shown in Figure 9: we report the numerical data for different values of n and α; as expected from our previous discussion, the agreement with the predictions is very good for large α, but it worsens as α gets smaller and n gets larger. Smaller is α, larger is the value of on the lattice necessary to observe the true asymptotic behaviour.

Symmetry resolution
The symmetry resolved moments of the RDM can be computed through the Fourier transform of the leading term of the charged moments in Eq. (5.11) where Erfi(x) is the imaginary error function (the overall result is real and positive for q ∈ Z) In the large limit, using the expansion in Eq. (5.16), the charged moments in Eq. (5.15) can be can be approximated as Z n (q) = Z n (0) n ln / n 2 π 2 q 2 + ln 2 / , (5.17) and hence the symmetry resolved entropies are given by S n (q) = 1 1 − n ln Z n (q) Z 1 (q) n S n − ln ln + ln n 1 − n , S 1 (q) S 1 − ln ln − 1.
The leading behaviour is described by the total Rényi entropies, with the usual correction ln ln that is independent on q, confirming the equipartition of the entanglement entropy for a complex massive scalar field theory, in agreement with the result for massive harmonic chains [20] (although the critical limit considered there is different from the one here). Let us mention that a further expansion of Eq. (5.15) leads to subleading corrections behaving as q 2 /(ln ) 2 which explicitly depend on q, breaking the equipartition of the entanglement. This kind of terms has been already found for bosonic systems in [20]. Let us now discuss the effect of the term that we disregarded in Eq. (5.11), namely the sum over k. The mode with k = 0 would provide double logarithimic corrections encountered also in other contexts, like non unitary CFTs [56,57]. These in principle are calculable and partially under control. We mention that such terms have a non-trivial dependence on n in Z n (q) and hence they are responsible of a further breaking of equipartition. Unfortunately, the determination of this correction is not easy because it is influenced by the precise dependence on α and n of the non-universal cutoff (as it should be clear from Eq. (5.11)). Finally, as discussed for the charged moments, the effect of the mode k = 0 is even more dramatic and too difficult to keep under control.

The long distance expansion.
The boundary condition for Eq. (5.6) in the limit in which t → ∞ is [41] The solution of Eq. (5.6) in the long distance regime together with Eq. (5.5) gives w a (t) = −e −2t sin 2 (aπ) π 1 + 3 − 16a + 16a 2 4t . (5.20) Summing over a = k n + |α| 2πn , we get c n (α) = e −2t 2πnt −n 2 t − 8 + n 2 12 + 2|α| π − α 2 π 2 + 2 csc 2 π n − |α| π cos α n + 2α π cot π n sin α n , (5.21) The long distance leading term in Eq. (5.21) is showed in Fig. 10 for two different values of n: it approximates well the solution of the Painlevé equation (5.6) in the regime t 1. The same feature was observed in Sec. 4.3 for the corresponding equations in fermionic systems as also the discontinuity for n → 1, which can be ascribed to the non-commutativity of the limits n → 1 and t → ∞. Also for a complex scalar field, Eqs. (5.21) show that the functions c n (α) vanish for large t and the charged moments stop growing. Hence,  so that all contributions coming from the long-distance behaviour are negligible at order O(1/(ln m )). The resolved entropies are the ones given in (3.24), where S n also takes into account the term ne −2t 4πt . The limit n → 1 can be solved through a technique similar to the one used in Sec. 4.3.

Charged moments across the hyperplane: massive scalar field
In this last section, we provide a result for the charged moments of a free massive scalar theory across a hyperplane in d Euclidean dimensions using a generalisation of the method reported in [5,58]. The action for a free complex massive scalar field is (6.1) We denote the space coordinates by x i , i = 1, · · · , d − 1, and the Euclidean time by x 0 . Let A andĀ be regions with x 1 > 0 and x 1 ≤ 0, respectively. The entangling surface Σ is chosen to be a (d − 2)-dimensional hyperplane at x 1 = 0: Σ = {(x 0 , x i )|x 0 = x 1 = 0}. Let us introduce the metric of the spacetime where we have used the polar coordinates for the plane parametrised by (x 0 , x 1 ). The metric of the n-fold cover M n,α of the original spacetime pierced by a flux α, that is constructed by gluing n copies of the sheet with a cut along A, is still (6.2) with r ≥ 0 and 0 ≤ τ ≤ 2πn. Thus, M n,α = C n,α × R d−2 , where C n,α is the two-dimensional cone parametrized by (r, τ ).
In order to compute the charged moments, we can consider the theory with action in Eq. (6.1) in terms of two real scalar fields ϕ 1 (x) and ϕ 2 (x) [20]. The logarithm of the normalised charged moments ln Z R n (α) of each real scalar field on M n,α is given by where the parameter 2 1 is introduced as a regulator for the UV divergences. Because of the direct product structure of M n,α , the Laplacian decomposes into the sum of those on C n,α and R d−2 : The rotational symmetry of the cone C n,α allows the Fourier decomposition of the real fields by the modes exp(iτ a), where a = l n + |α| 2πn , with integer l, such that, after a period 2πn, they acquire a phase e iα , in agreement with the Aharonov-Bohm effect. Therefore, the eigenfunctions ϕ k,a (r, τ ) of the Laplacian are parametrized by (k, a) satisfying Cn,α ϕ k,a (r, τ ) = −k 2 ϕ k,a (r, τ ), k ∈ R + , ϕ k,a (r, τ ) = k 2πn e iτ a J |a| (kr), where J a is the Bessel function of the first kind. The eigenfunctions form an orthonormal basis on the cone C n,α , namely The orthonormal basis of the eigenfunctions of the Laplacian on R d−2 is spanned by the plane waves, ϕ k ⊥ (y) = exp(ik ⊥ · y)/(2π) (d−2)/2 , with eigenvalues −k 2 ⊥ . Exploiting these two sets of eigenfunctions, the trace of the kernel in Eq. (6.3) is The IR divergence emerges from the following integrals where I a is the modified Bessel function of the first kind. The second term on the right-hand side of Eq. (6.7) is divergent, but it gives rise to a term proportional to n and independent on α, so it does not contribute to ln Z n (α). On the other hand, the UV divergence arises from the summation over the angular momentum l, which can be regularised, for example, by using the Hurwitz zeta function The definition of the Hurwitz zeta function for complex arguments s with Re(s) > 1 and q with Re(q) > 0 is ζ(s, q) = ∞ j=0 1 (q + j) s , (6.9) whose analytic continuation satisfies the remarkable identity The other kernel tr(e −s ) is still divergent because it is proportional to the volume of M n,α but it is also proportional to n coming from the volume of the cone, so, being independent on α, it does not contribute to ln Z n (α). Thus, using Eqs. (6.3), (6.6), and (6.8) we obtain the logarithm of the charged moments across a hyperplane R d−2 for a free complex massive scalar field [20] ln Z n (α) = ln Z R n (α) Very remarkably, the dependence on α in Eq. (6.11) is the same in any dimension and hence the symmetry resolved entropies are also the same in any dimension.
As a consistency check, let us focus our attention on d = 2: for a semi-infinite line, we obtain where we have used the definition of the exponential integral function (6.13) By expanding around = 0, −Ei(−m 2 2 ) γ E + 2 ln(m ) and considering n = 1 we end up into where γ E is the Euler-Mascheroni constant. This coincides with what found previously in Eq. (3.18) up to O(1). In d = 3, Eq. (6.11) coincides with the leading term obtained in [33] for the critical limit of a two-dimensional harmonic lattice.

Conclusions
In this manuscript we characterised the symmetry resolved entanglement for free massive fields in two dimensions, presenting the results for a Dirac field and a complex scalar theory. We showed that two well known techniques in the framework of the replica trick can be adapted -by modifying the n-sheeted Riemann surface and the corresponding partition function-to the calculation of charged moments. Both computations (via modified twist fields and the Green's function approach of Ref. [42]) mainly rely on the boundary conditions of the fields at the endpoints of the entangling region. In the first framework, the conformal dimensions of the twist fields get modified as in Eq. (3.8). In the second setting the change induced by the flux α lies in the precise form of the Painlevé V equations (4.6) and (5.6) providing the generalised partition function. These Painlevé equations are easily solved numerically for arbitrary values of the mass, but they can be also handled analytically in the limit of small masses, leading to the charged moments (4.12) for the Dirac field and (5.11) for the scalar theory. The opposite limit of mass much larger that the interval length can also be treated analytically. For the free complex scalar, we also obtain general results for the charged moments in arbitrary dimension when the entangling surface is an hyperplane.
From the Fourier transform of these charged moments, we extract the symmetry resolved Rényi entropies, stressing the relevant universal aspects. At leading order for small m, the symmetry resolved entropies for both theories satisfy equipartition of entanglement [16]. We also show that the entanglement equipartition is broken by the mass at order (ln ) −2 , which is the same one found in other circumstances [17,18,20,33].
There are two main aspects that our manuscript leave open for further study. The first one concerns the calculation of charged and symmetry resolved entropies in free scalars and fermions in arbitrary dimension and for entangling surfaces that are more complex than the simple hyperplane of Sec 6. To this aim, we expect that some of the existing techniques in the literature, as e.g. in Refs. [59][60][61][62][63][64][65][66][67][68][69], should be readily adapted to our problem. Furthermore, important insights could also come from the holographic correspondence for the entanglement entropy [58,70,71]. The other point is whether interacting QFTs can be handled in two dimensions, e.g. exploiting integrability techniques as those of Refs. [36,37,50,72,73].

Acknowledgments
We thank Luca Capizzi, Dávid Horváth, Pierluigi Niro and Paola Ruggiero for useful discussions. PC and SM acknowledge support from ERC under Consolidator grant number 771536 (NEMO).

A Conformal dimensions of twist fields
The goal of this section is to find the conformal dimension of the twist field T n,k,α defined in Eq. (3.4). We will call it generically T a , where a = k n + α 2πn , with a ∈ [0, 1]. As already discussed in section 3.1, in the neighbourhood of a twist field the k-th component of φ undergoes a phase rotationφ k (e 2πi z, e −2πiz ) = e 2πiaφ k (z,z).
Let us start from the case of the free complex scalar CFT with fields (ϕ k , ϕ * k ) and, following [39], consider the correlation function in the presence of four Z n twist-fields .
Imposing that for z → w we have g(z, w; z i ) ∼ (z − w) −2 and that for z → z i we have g(z, w; z i ) ∼ (z − z j ) −a for j = 1, 3 and g(z, w; z i ) ∼ (z − z j ) −(1−a) for j = 2, 4, we can write (up to an additional constant independent of z and w, A(z j ,z j )) where This is exactly the expectation value of the insertion of the stress energy tensor of the field ϕ k in the four-point correlation function. From the comparison with the conformal Ward identity, we can understand that T a andT a are primary fields with scaling dimensions In order to obtain the conformal dimensions of the twist fields of the free Dirac field theory, let us apply a similar procedure for the chiral or anti-chiral complex fermionic fields, (ψ k , ψ * k ). The scaling dimension of T a can be extracted from the Green's function in presence of two Z n twist fields .
Using the results in [74], the previous expression can be explicitly written as where In the limit w → z This is the expectation value of the insertion of the stress energy tensor of the field ψ k in the two-point correlation function and, as before, the comparison with the conformal Ward identity gives the dimensions of the primary twist fields T a andT a as Putting together the monodromy conditions (A.1) and the scaling dimension of the twist field in Eq. (A.11), we deduce that the twist field of a fermionic field admits a bosonisation formula. We can write the complex fermionic field as ψ k ∼ e iϕ k and the twist field as T n,k,α (z) = e i( k n + α 2πn )ϕk . By introducing the vertex operators V β (z) = e iβϕ(z) , the twist fields take the form T n,k,α (z) = V k n + α 2πn (z) andT n,k,α (z) = V − k n − α 2πn (z) [45,74]. Let us observe that at first sight this result could be misleading since the outcome for bosons in Eq. (A.6) does not appear to agree with that of fermions in Eq. (A.11) given that they are related by bosonisation in 1+1 dimensions [75][76][77]. However, via bosonisation of U (1) complex fermions, the corresponding bosons transform by translation, and thus should instead satisfy the boundary condition ϕ k (e 2πi z) = ϕ k (z) + a. Therefore, our computation for charged bosons is not related to charged fermions by bosonisation.
Before concluding this appendix, let us emphasise that while CFTs are well understood objects, n-copies of a CFT after modding out the Z N symmetry among the replicas form a more complicated object known as orbifold [38,39]. The operator product expansions of the twist fields with other fields have been extensively explored (e.g., see [78][79][80][81][82][83][84][85][86]), but, unless a bosonisation procedure for free theories can be used, as for the compact boson, they remain elusive in general and require a case-by-case study.

B Details for the analytic continuation for the Dirac field
In this appendix we provide some details about the analytic continuation of the quantities defined in Eq. (4.10). First, we rewrite Ω(n, α) (a similar result also holds for Λ(n, α)) as Ω n (α) ≡ 1 n 2 (n−1)/2 k=−(n−1)/2 k + α 2π where we set k ≡ k + n−1 2 . After this manipulation, the functions Ω(n, α) and Λ(n, α) in (4.10) can be split as We used the relation ψ(x) = ψ(x + 1) − 1 x so that all the terms are in a suitable form for the use of the following integral representation of the digamma function, i.e.

(B.5)
The same strategy also works for C j (n, α), j = 1, 2, 3 in Eqs. (B.3), however this computation is longer and cumbersome and we do not report here. In Figure 11 we plot Λ(n, α) (left top panel) and Ω(n, α) (right top panel) as function of α for some n and compare it with the quadratic approximation given in Eqs. (4.10). The closeness of the two curves shows that the quadratic approximation is enough for most of the applications. The quadratic approximation depends non trivially on n, as evident from the plot of ω n and λ n in the bottom panels of Figure 11.

C The lattice models
For the numerical test of our field theory results we consider the following lattice discretisation of the Dirac fermion [40] where the p i and q i satisfy the canonical commutation relations [q i , q j ] = [p i , p j ] = 0 and [q i , p j ] = iδ ij . In the thermodynamic limit, the correlators can be written in terms of the hypergeometric functions [90] q i q j = ζ i−j+1/2 2 i − j − 1/2 i − j 2 F 1 1/2, i − j + 1/2, i − j + 1, ζ 2 , (C.6) where the parameter ζ is defined by Let us denote as X and P the matrices of the correlators of positions and momenta (i.e. X ij = q i q j and P ij = p i p j ). The moments of the reduced density matrix of A can be written in terms of the eigenvalues of √ XP that we call σ i (with i ∈ [1, ]). To have a continuous symmetry, we consider a complex bosonic theory which on the lattice is a chain of complex oscillators. It is the sum of two real harmonic chains in the variables (p (1) , q (1) ) and (p (2) , q (2) ), i.e.
H CB (p (1) + ip (2) , q (1) + iq (2) ) = H B (p (1) , q (1) ) + H B (p (2) , q (2) ). (C.9) The Hamiltonian (C.9) can be also written in terms of particles and antiparticles mode operators, a k and b k , satisfying [a k , a † j ] = δ j,k , [b k , b † j ] = δ j,k . In terms of these operators, the Hamiltonian (C.9) is , ω k = ω 2 0 + 4 sin 2 πk L , (C.10) while the charge operator reads i.e. the total number of particles minus the number of antiparticles. The conserved charge can be as well written in real space and its value in a given subsystem A is the same sum restricted to A, i.e.
For the charged moments, we need to compute Tr(ρ n A e iQ A α ), but using the form in Eq. (C.12) for Q A , the trace factorises as where N a A = j∈A a † j a j and N b A = j∈A b † j b j . Using the relations between the number operator N A and the eigenvalues of the correlation matrix σ i , one finds [20] Z n (α) = (C.14)