Multi-charged moments of two intervals in conformal field theory

We study the multi-charged moments for two disjoint intervals in the ground state of two $1+1$ dimensional CFTs with central charge $c=1$ and global $U(1)$ symmetry: the massless Dirac field theory and the compact boson (Luttinger liquid). For this purpose, we compute the partition function on the higher genus Riemann surface arising from the replica method in the presence of background magnetic fluxes between the sheets of the surface. We consider the general situation in which the fluxes generate different twisted boundary conditions at each branch point. The obtained multi-charged moments allow us to derive the symmetry resolution of the R\'enyi entanglement entropies and the mutual information for non complementary bipartitions. We check our findings against exact numerical results for the tight-binding model, which is a lattice realisation of the massless Dirac theory.


Introduction
As Schrödinger already recognised one century ago, entanglement is at the core of quantum mechanics.Nowadays it turns out to be the fundamental notion behind many quantum phenomena, from quantum algorithms [1] to gravity [2,3], passing by critical phenomena and topological phases of matter [4][5][6][7], triggering unexpected connections between apparently far branches of physics.At the center of all these ideas, we find the (Rényi) entanglement entropies which are powerful entanglement measures that provide fundamental insights about the investigated system or theory.They are defined as follows.Let us consider an extended quantum system in a pure state |Ψ and a spatial bipartition into A and B. The subsystem A is described by the reduced density matrix ρ A = Tr B |Ψ Ψ| and the associated Rényi entropies are given by the moments of ρ A as where we assume that n is an integer number.After the analytic continuation to complex values of n, the limit n → 1 of Eq. (1.1) yields the von Neumann entanglement entropy For bipartite systems in a pure state, the von Neumann and Rényi entropies can be used as measures of the entanglement shared between the two complementary parts.One of the most interesting properties of entanglement entropies is that they are sensitive to criticality.
In particular, for one-dimensional gapless systems, if A is a single interval, then the ground state entanglement entropy breaks the area law and is proportional to the central charge of the 1+1 dimensional CFT that describes the low-energy spectrum of the system [8][9][10][11].
In the case considered in this work in which A consists of two subsystems A 1 and A 2 , i.e.A = A 1 ∪ A 2 , the ground-state entanglement entropy depends on the full operator content of the CFT, encoding all the conformal data of the model [12][13][14][15].It is important to remark that, in this situation, the entanglement entropies quantify the entanglement between A and B but not between the two parts of A, for which one must resort to other entanglement measures such as negativity [16][17][18][19][20][21][22].Nevertheless, from the entanglement entropies, it is possible to construct the following quantity, dubbed mutual information, which is a measure of the total correlations between A 1 and A 2 .The computation of twointerval Rényi entanglement entropies is a difficult problem, even for minimal CFTs [21,23], as it boils down in general to determine the partition function of the theory on a higher genus n-sheeted Riemann surface [14,15].In fact, exact analytic expressions are only available for the free theories or special limits [13][14][15][24][25][26][27][28][29][30][31][32][33][34][35][36][37].Moreover, the analytic continuation in n to obtain Eq. (1.2) is still a challenging open issue.
The symmetry resolution of entanglement in the two-interval case has not been much explored in CFT.Ref. [78] studies it at large central charge, in the context of holography, while, in Ref. [74], the charged Rényi negativity is analysed for the complex free boson.Here we take a different charge for each part of A, which leads to introduce the multi-charged moments of ρ A .This non-trivial generalisation of the charged moments, first considered in Ref. [82] in the context of quench dynamics, is the main subject of this work.In CFT, they correspond to the partition function on the n-sheeted Riemann surface, but with the insertion of a different magnetic flux across each subset (interval) of A. We compute the multi-charged moments analytically for the ground state of two bidimensional CFTs with central charge c = 1 and global U (1) symmetry-the massless Dirac field theory and the free compact boson-generalising the expressions for the (neutral) Rényi entropies found in Refs.[24] and [14] respectively.From the multi-charged moments, we derive the ground state symmetry-resolved entanglement entropy and mutual information of two disjoint intervals.
The paper is organised as follows: in Sec. 2, we define the symmetry-resolved entanglement and mutual information as well as the multi-charged moments, and we briefly describe the general approach to compute the latter in CFTs.We then move on to calculate the multi-charged moments for the ground state of the massless Dirac field theory in Sec. 3 and of the free compact boson in Sec. 4. In Sec. 5, we apply the previous results to obtain the symmetry-resolution of the mutual information in these theories.When possible, we benchmark the analytic expressions with exact numerical calculations for lattice models in the same universality class.We draw our conclusions in Sec.6 and we include three appendices, with more details about the analytical and numerical computations.

Definitions
In this section, we first give the definition of the symmetry-resolved entanglement entropy and mutual information for a subsystem composed of two disjoint regions.We explain their relation with the multi-charged moments of the reduced density matrix, and we introduce the replica method to calculate them in CFTs.

Symmetry-resolved entanglement entropies and mutual information
As we already pointed out in Sec. 1, we take a spatial bipartition A ∪ B of an extended quantum system in a pure state |Ψ , with A made of two disconnected regions, A = A 1 ∪A 2 .We assume that the system is endowed with a global U (1) symmetry generated by a local charge Q.Given the partition of the system in different subsets, we can consider the charge operator in each of them; for example, in region A, it can be obtained as and, by taking the trace over B, we find that [Q A , ρ A ] = 0.This implies that the reduced density matrix ρ A presents a block diagonal structure, in which each block corresponds to an eigenvalue q ∈ Z of Q A .That is, where Π q is the projector onto the eigenspace associated to the eigenvalue q and p(q) = Tr (Π q ρ A ) is the probability of obtaining q as the outcome of a measurement of Q A .Notice that Eq. (2.1) guarantees the normalisation Tr[ρ A (q)] = 1 for any q.
The amount of entanglement between A and B in each symmetry sector can be quantified by the symmetry-resolved Rényi entropies, defined as (2.2) Taking the limit n → 1 in this expression, we obtain the symmetry-resolved entanglement entropy, According to the decomposition of Eq. (2.1), the total entanglement entropy in Eq. (1.2) can be written as [1] where S c is known as configurational entropy and quantifies the average contribution to the total entanglement of all the charge sectors [41,[50][51][52], while S num is called number entropy and takes into account the entanglement due to the fluctuations of the value of the charge within the subsystem A [41,88,89,92,[106][107][108].
Since the total charge in A is the sum of the charge in A 1 and A 2 , , then the reduced density matrices ρ A 1 , ρ A 2 of A 1 and A 2 can be independently decomposed in charged sectors as we did for ρ A in Eq. (2.1).Therefore, we can define the symmetryresolved entropies S A 1 n (q 1 ), S A 2 n (q 2 ) for the regions A 1 and A 2 analogous to Eq. (2.2) for A, with q = q 1 + q 2 .In Ref. [82], it has been proposed to define the symmetry-resolved mutual information as The quantity p(q 1 , q − q 1 ), normalised as is the probability that a simultaneous measurement of the charges Q A 1 and Q A 2 yields q 1 and q − q 1 , respectively, while the charge of the whole system A is fixed to q.Although Eq. (2.5) is a natural definition, I A 1 :A 2 (q) is not in general a good measure of the total correlations between A 1 and A 2 within each charge sector since, in some cases, it can be negative [82].Nevertheless, we find interesting to investigate this quantity given that it provides a decomposition for the total mutual information (2.5) similar to the one reported in Eq. (2.4) for the entanglement entropy, where num is the number mutual information.

Charged moments and symmetry resolution
The computation of the symmetry-resolved entanglement entropies and mutual information from the definitions (2.2) and (2.5) requires the knowledge of the entanglement spectra of ρ A , ρ A 1 and ρ A 2 and their symmetry resolution.However, this is usually a very difficult task, in particular if one is interested in analytical expressions.Alternatively, one can employ the charged moments of the reduced density matrices.For ρ A , they are defined as Similar quantities can also be introduced for the two subsystems A 1 and A 2 that constitute A by replacing ρ A and Q A by ρ Ap and Q Ap , p = 1, 2. If we take now their Fourier transform, the symmetry-resolved entanglement entropies of the subsystem A are given by [39,40] n . (2.10) In a similar manner, replacing A with A 1 and A 2 , we can obtain the symmetry-resolved entanglement entropies of the two components of A.
Notice that, for computing the symmetry-resolved mutual information of Eq. (2.5), we need to determine p(q 1 , q − q 1 ), i.e. the probability that a measurement of Q A 1 and Q A 2 gives q 1 and q − q 1 respectively, with Q A fixed to q.In order to calculate it, we consider the generalisation of the charged moments in Eq. (2.8) introduced for the first time in Ref. [82], We refer to them as multi-charged moments.When α = β, Eq. (2.11) reduces to the charged moments of A = A 1 ∪ A 2 of Eq. (2.8).If we take the Fourier transform of Eq. (2.11), then Z A 1 :A 2 1 (q 1 , q 2 ) can be interpreted as the probability of having q 1 and q 2 as outcomes of a measurement of Q A 1 and Q A 2 respectively, independently of the value of Q A .Therefore, it satisfies the normalisation and p(q 1 , q − q 1 ) can be calculated as the conditional probability which fulfills Eq. (2.6).

Charged moments in CFT
In the rest of the paper, we will analyse the previous quantities in 1 + 1-dimensional CFTs with a global U (1) symmetry.We will assume that the entire system is in the ground state and that the spatial dimension is an infinite line which we will divide into two parts A and B, with A made up of two disjoint intervals, namely If we denote by 1 and 2 the lengths of the two intervals and d their separation, we have where we have also introduced the cross ratio x of the four end-points, which takes values between 0 and 1.
As explained in detail in Refs.[10,11], using the path integral representation of ρ A , the moments Z A n (0) = Tr[ρ n A ] are equal to the partition function of the CFT on a Riemann surface, which we call Σ n , obtained as follows.We take the complex plane where the CFT is originally defined and we perform two cuts along the intervals A 1 = (u 1 , v 1 ) and A 2 = (u 2 , v 2 ).Then we replicate n times the cut plane and we glue the copies together along the cuts in a cyclical way as we illustrate in Fig. 1.We eventually obtain an n-sheeted Riemann surface of genus n − 1, which is symmetric under the Z n cyclic permutation of the sheets.
Alternatively, instead of replicating the space-time where the CFT is initially defined, one can take n copies of the CFT on the complex plane and quotient it by the Z n symmetry under the cyclic exchange of the copies.We then get the orbifold theory CFT ⊗n /Z n .The moments Z A n (0) are equal to the four-point correlation function on the complex plane [10,109], where τ n and τn are dubbed as twist and anti-twist fields [10,[109][110][111].They implement in the orbifold the multivaluedness of the correlation functions on the surface Σ n when we go around its branch points.In fact, the winding around the point where τ n (τ n ) is inserted maps a field O k living in the copy k of the orbifold into the copy k + 1 (k − 1), that is (2.17) The twist and anti-twist fields are spinless primaries with conformal weight where c is the central charge of the initial CFT.
< l a t e x i t s h a 1 _ b a s e 6 4 = " a l T M P Y s J q z F n W g n e 1 X J 6 p g V g n b g = " o u S I 3 U C S P X 5 I 7 c k w f r 1 n q 0 n q z n j 9 E 5 a 5 b Z I d 9 g v b w D 9 S q V v A = = < / l a t e x i t > ↵ < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 L 6 D s v 9 u g j b 1 C / q 6 Z 0 9 c 9 b j w H d E = " > A A A B / n i c d V D L S s N A F J 3 4 r P V V d e l m s A g u N C R a t N 0 V 3 b i s Y B / Q h n I z v W 2 H T J J h Z i K U U P A X 3 O r e n b j 1 V 9 z 6 J a a 1 g s + z O p x z D / f e 4 0 v B t X G c V 2 t u f m F x a T m 3 k l 9 The charged moments of ρ A can also be computed employing the previous frameworks.As argued in Ref. [39], the operator e iαQ A can be interpreted as a magnetic flux between the sheets of the surface Σ n , such that a charged particle moving along a closed path that crosses all the sheets acquires a phase e iα .For the multi-charged moments introduced in Eq. (2.11), we have to insert two different magnetic fluxes α and β at the cuts A 1 and A 2 respectively, as we pictorially show in Fig. 1.They can be implemented by a local U (1) operator V α (x) that generates a phase shift e iα along the real interval [x, ∞).Then the charged moments are equal to the four-point correlation function on the surface Σ n In the orbifold theory, the magnetic flux can be incorporated by considering the composite twist field τ n,α ≡ τ n • V α .Thus, if we take a field O k in the copy k of the orbifold, then the winding (z − u) → e 2πi (z − u) around the point u where τ n,α is inserted gives rise to a phase e iα/n , (2.20) The same applies to the composite anti-twist field τn,α ≡ τn • V α , which takes a field from the copy k to k − 1 adding a phase e iα/n .Therefore, Eq. (2.19) can be re-expressed as the four-point function on the complex plane In Ref. [39], it is shown that, if V α is a spinless primary operator with conformal weight h V α , then so are the composite twist and anti-twist fields, with conformal weights One can further consider other configurations for the magnetic fluxes between the sheets of the Riemann surface Σ n .In general, if we assume that a particle gets a different phase e iα j when it goes around each branch point, provided they satisfy the neutrality condition α 1 + α 2 + α 3 + α 4 = 0, then the partition function of this theory is given by or, in terms of the composite twist fields, by Then the multi-charged moments Z A 1 :A 2 n (α, β) can be treated as the particular case in which α 1 = −α 2 = α and, due to the neutrality condition, In Sec. 3, we will compute the charged moments Z A 1 :A 2 n (α, β)-and more in general the partition functions Z A 1 :A 2 n ({α j })-for the massless Dirac fermion using the orbifold theory CFT ⊗n /Z n .On the other hand, in Sec. 4, we will adopt a geometric approach to obtain the multi-charged moments of the compact boson from the correlation function on the Riemann surface Σ n of Eq. (2.23).

Free massless Dirac field theory
The massless Dirac field theory is described by the action where ψ = ψ † γ 0 .The γ µ matrices can be represented in terms of the Pauli matrices as γ 0 = σ 1 and γ 1 = σ 2 .The action of Eq. (3.1) exhibits a global U (1) symmetry: it is invariant if the fields are multiplied by a phase, i.e. ψ → e iα ψ and ψ → e −iα ψ.By Noether's theorem, this symmetry is related to the conservation of the charge The ground state entanglement of a subsystem A made up of multiple disjoint intervals in the ground state of this theory was first investigated in Ref. [24].For the case of two disjoint intervals, A = A 1 ∪ A 2 , it was found that the moments of ρ A are where c n is a non-universal constant.
In this section, we will compute the multi-charged moments of Eq. (2.11) in the ground state of the massless Dirac field theory.We will extend the approach introduced in Ref. [24] for the moments of Eq. (3.2).Similar techniques have been exploited in Ref. [60] for studying the charged moments of Eq. (2.8), when A is a single interval, in two dimensional free massless Dirac theories and in Ref. [67] in the context of the charge imbalance resolved negativity.We will benchmark our analytical results with exact numerical calculations in a lattice model.

Charged moments
In Sec. 2, we explained that the partition function Z A 1 :A 2 n ({α j }) can be obtained either by considering the theory on a complicated Riemann surface or by replicating it n-times and working with the orbifold on the complex plane.For the massless Dirac field theory, the latter approach is more convenient.Thus, let us take the n-component field where ψ j is the Dirac field on the j-th copy of the system.Eq. (2.20) describes the effect of the composite twist fields on the components of Ψ when going around the end-points of the subsystem This transformation can be encoded in the matrix In the general case of Eq. (2.24), Ψ transforms according to T α 2p−1 when winding around the point u p and to the transpose matrix T t α 2p when going around the point v p , with p = 1, 2. The matrix T a in Eq. (3.4), sometimes called twist matrix, was introduced for the case a = 0 in Refs.[24,109] and for general a in Ref. [60].Its eigenvalues are of the form By simultaneously diagonalising all the T α j with a unitary transformation (which is independent of α j ), we can recast the replicated theory in n decoupled fields ψ k on the plane, which are multi-valued, Notice that this technique, known as diagonalisation in the replica space, can be applied only to free theories since, otherwise, the k-modes do not decouple.For the free massless Dirac theory, this allows us to write the Lagrangian of the replicated theory as Following this approach, the partition function of Eq. (2.24) factorises into where ({α j }) is the partition function for a Dirac field ψ k with the boundary conditions of Eq. (3.6).
The main difference between the partition functions Z A 1 :A 2 k,n ({α j }) and the standard computation of Ref. [24] for Rényi entropies is that the boundary conditions of the multivalued fields around the branch points now depend on the phases α j .This multivaluedness can be removed, as done in [24] for α j = 0, by introducing an external U (1) gauge field A k µ coupled to single-valued fields ψk .In fact, if we apply the singular gauge transformation then the Lagrangian for the k-th mode can be rewritten as with the advantage of absorbing the phase around the end-points of A 1 ∪ A 2 into the gauge field.The only requirement that A k µ in Eq. (3.9) must satisfy is that, integrated along any closed curve C that encircles the end-points of A, the boundary conditions of Eq. (3.6) for ψ k must be reproduced.For this purpose, we require where C up and C vp are closed contours around the end-points of the p-th interval.Moreover, we have to impose that, if C does not enclose any end-point, then C dy µ A k µ = 0. Applying the Stoke's theorem, the conditions of Eqs.(3.11) can be expressed in differential form, Once the transformation of Eq. (3.9) is performed, the partition function ({α j }) of the k-th mode is equal to the vacuum expectation value where j µ k ≡ ψk γ µ ψk is the conserved Dirac current for each mode.Eq. (3.13) can be easily computed via bosonisation [24], which allows us to write the current in terms of the dual scalar field φ k such that j µ k = µν ∂ ν φ k / √ π.If we use this result in Eq. (3.13), and we apply Eq. (3.12), then Notice that the neutrality condition α 1 + α 2 + α 3 + α 4 = 0 ensures that the latter correlator does not vanish.The correlation function of vertex operators in the complex plane is wellknown (see, for instance, Ref. [113]) and, therefore, Eq. (3.14) can be easily calculated.
Plugging the result into Eq.(3.8) and performing the product over k, we obtain where x is the cross-ratio defined in Eq. (2.15).When we take α 1 = −α 2 = α and α 3 = −α 4 = β in this expression, we get the multi-charged moments in Eq. (2.11) as We assume that all the length scales in this formula have been regularised through a UV cutoff which is included in the multiplicative constant c n;α,β .An interesting case to analyse is when the two intervals A 1 and A 2 become adjacent; that is, when d → 0. In that limit, the cross-ratio x tends to one such that and Eq.(3.16) vanishes.Nevertheless, in this regime, the distance d must be regarded as another UV cutoff, which can be absorbed in the multiplicative constant c n;α,β .Therefore, from Eqs. (3.17) and (3.16), one expects which agrees with the fact that, according to Eq. (2.21), the multi-charged moments In Fig. 2, we check the expressions obtained in Eqs.(3.16) and (3.18) for Z A 1 :A 2 n (α, β) with exact numerical calculations performed in the tight-binding model, which is a chain of non-relativistic free fermions whose scaling limit is described by the massless Dirac field theory.The details of the numerical techniques employed are given in Appendix A. In order to compare Eqs.(3.16) and (3.18) with the numerical data, we need the concrete expression of the non-universal factors c n;α,β and cn;α,β for this particular model.When the two intervals A 1 and A 2 are far away, that is in the limit d → ∞, the multi-charged moments of A 1 ∪ A 2 factorise into those of A 1 and A 2 , Therefore, one expects c n;α,β to be the product of the two non-universal constants associated to A 1 and A 2 as single intervals.The latter were obtained for the tight-binding model in Ref. [45] by exploiting the asymptotic properties of Toeplitz determinants.We can also apply here those results, taking into account that each interval is associated to a different flux, either α or β.Then we have where [45]  Expanding Υ(n, α) up to quadratic order in α, then Eq. (3.20) can be approximated as where ζ n = ln 2 − 2π 2 nγ 2 (n) and . (3.23) In the limit of adjacent intervals, given by Eq. (3.18), the multiplicative constant cn;α,β cannot be determined using the known results for Toeplitz determinants.However, we can conjecture an analytical approximation for it at quadratic order in α and β.When d → 0, we can associate to each end-point of the intervals u 1 , v 1 = u 2 and v 2 the fluxes α, β −α and −β respectively.From the results for one interval of Ref. [45], one can conjecture that each end-point with flux α contributes with a factor e − ζn 4π 2 n α 2 to the constant cn;α,β , if we restrict to terms up to order α 2 .Therefore, the combination of all the fluxes in our case should contribute with a total factor e − ζn 4π 2 n (α 2 +(β−α) 2 +β 2 ) .We then expect that cn;α,β should be well be approximated by When α = β = 0, the expression above simplifies to the multiplicative constant for the moments of a single interval obtained in Ref. [112].In spite of the heuristic reasoning of this result, in Fig. 2, we check its validity by comparing it against exact numerical data.
In the top and middle panels of Fig. 2, we study Z A 1 :A 2 1 (α, β) as function of α, for various values of β, 1 , 2 and d > 0. The points correspond to the exact numerical values obtained as we described in Appendix A, while the solid curves are the analytic prediction of Eq. (3.16), taking as multiplicative constant c n,α,β that in Eq. (3.20).We find an excellent agreement.In the bottom left panel, we analyse the case of adjacent intervals (d = 0).The numerical data for Z A 1 :A 2 1 (α, β) are in very good agreement with the analytical expression (solid curves) of Eq. (3.18), using Eq.(3.24) as multiplicative constant.Finally, in the bottom right panel, we plot Z A 1 :A 2 1 (α, β) as function of the cross-ratio x, for various values of α, β, 1 and 2 .The curves represent the prediction of Eq. (3.16).The continuous ones correspond to take as non-universal constant c n;α,β the full expression of Eq. (3.20), while for the dashed ones we have used the quadratic approximation of Eq. (3.22).The agreement between the analytic prediction and the numerical data is extremely good, even considering the quadratic approximation for the non-universal constant.As expected, this agreement is better for small values of α and β, while, around ±π, we need to take larger subsystem sizes 1 , 2 in order to suppress the finite-size corrections which are well-known and characterised for the charged moments with a single flux insertion [45].

Free compact boson
The second theory we focus on in this manuscript is the free compact boson, which is the CFT of the Luttinger liquid, and whose action reads The target space of the real field ϕ is compactified on a circle of radius R, i.e. ϕ ∼ ϕ+2πmR with m ∈ Z.The compactification radius R is related to the Luttinger parameter K as R = 2/K.The action of Eq. (4.1) is invariant under the transformation ϕ → ϕ + α which, due to the compact nature of ϕ, realises a U (1) global symmetry.The associated conserved charge is The moments of the ground state reduced density matrix of this theory are well-known.An exact analytic expression for the two-interval case was obtained in Ref. [14], which was generalised to an arbitrary number of disjoint intervals in Ref. [30].In particular, for two intervals, it reads [14] where c n is a non-universal constant and We denote by Θ the Riemann-Siegel Theta function with characteristics ε, δ ∈ (Z/2) n−1 , u ∈ C n−1 and Ω a complex (n − 1) × (n − 1) matrix.In Eq. ( 4.3), the characteristics are zero 0 = (0, . . ., 0) and, therefore, we have used the standard shorthand notation Θ(u|Ω) ≡ Θ[ 0 0 ](u|Ω).The matrix Γ(x) in Eq. ( 4.3) has entries given by and with Although the moments of ρ A are known for all the integer n, its analytic continuation to complex n and, consequently, the von Neumann entropy of Eq. (1.2) is still not available for all the values of the Luttinger parameter.
A remarkable contact point between the theories described by the actions of Eqs.(3.1) and (4.1) is the case K = 1.Notice that, when the Luttinger parameter takes this value, the function F n (x) in Eq. ( 4.3) simplifies to F n (x)=1 and the moments of Eq. (4.2) for the massless compact boson present the same universal dependence on 1 , 2 and x as the ones in Eq. (3.2) for the massless Dirac fermion.A detailed discussion on this identity can be found in Ref. [29], where it is explained the reason why, although these two theories are not related by a duality, their partition functions on the Riemann surfaces Σ n arising in the two-interval replica method are actually equal.Here we find that this identity extends to the partition functions Z A 1 :A 2 n ({α j }) on the surface Σ n with different twisted boundary conditions at each branch point.In general, when A is made up of more than two intervals and the Rényi index n is larger than two, the moments of the reduced density matrix in these CFTs (and the corresponding Rényi entropies) are different [29].

Charged moments
We now generalise the result (4.2) to the multi-charged moments in Eq. (2.11).Starting from Eq. (2.19), we will compute them as the four-point function of the field V α on the Riemann surface Σ n .Since the U (1) conserved current is proportional to ∂ x 1 ϕ, V α is identified in this case with the vertex operator [39] V α (z) = e i α 2π ϕ(z) , (4.7) which has conformal dimensions In the following, it will be useful to introduce the rescaled Luttinger paratemeter η = K/(2π 2 ) in order to lighten the expressions.Without loss of generality, let us consider that the end-points of subsystem A are u 1 = 0, v 1 = x, u 2 = 1 and v 2 = ∞.Using Eq. (2.24), and given that the composite twist fields are primaries, we can eventually obtain the expression for an arbitrary set of end-points through a global conformal transformation.Therefore, according to Eq. (2.19), the multi-charged moments can be derived from the four-point correlation function of the vertex operators of Eq (4.7) on the n-sheeted Riemann surface Σ n (x) with branch points at 0, x, 1 and ∞.This surface of genus n − 1 can be described by the complex curve The correlator of vertex operators on a general Riemann surface of arbitrary genus was obtained in Ref. [115].In order to give the explicit expression in our case, we need to introduce some notions about Riemann surfaces [114].
There are different parameterisations of the moduli space of genus n − 1 Riemann surfaces.One possibility is through the matrix of periods, which we denote by Γ.This is a (n − 1) × (n − 1) symmetric matrix with positive definite imaginary part.Notice that, according to Eq. (4.10), the Riemann surface Σ n (x) is parametrised by the cross-ratio x.Therefore, the corresponding matrix of periods only depends on x, i.e.Γ = Γ(x).In order to define it, we need first to specify a particular homology basis for Σ n (x), i.e. a basis of 2(n−1) oriented non-contractible curves on the surface, which we denote by a r and b r , with r = 1, . . ., n − 1.The detailed description of the specific basis that we consider is given in Appendix B. We also have to choose a basis of holomorphic differentials ν r , r = 1, . . ., n−1, normalised with respect to the a r cycles.That is, ar dzν s (z) = δ r,s , r, s = 1, . . ., n − 1. (4.11) Then the matrix of periods is defined as For the surface Σ n (x), the normalised holomorphic differentials read In Appendix B, we thoroughly explain the derivation of the expression for ν r .Inserting it in Eq. ( 4.12), it is then easy to show that the entries of the matrix of periods Γ(x) are precisely those of Eq. (4.5).
If we now consider four vertex operators inserted at generic points in the surface Σ n (x) and with arbitrary dimensions satisfying the neutrality condition α 1 + α 2 + α 3 + α 4 = 0, then its correlation function is of the form [115] In this expression, we denote by E(z, z ) the prime form of the surface Σ n (x), which we will define precisely later, and w(z) = (w 1 (z), . . ., w n−1 (z)) is the Abel-Jacobi map, which relates a point z in the surface Σ n (x) to a point w(z) in the genus n − 1 complex torus This map can be written in terms of the normalised holomorphic differentials of Eq. (4.13) as where we have taken as origin the branch point z = 0.The images under the Abel-Jacobi map of the points z = 0, x, 1, and ∞, where the vertex operators in Eq. (4.9) are inserted, can be easily computed using Eq.(4.13) and applying the identities of Eq. (B.9).Then we find w(x) = q, (4.17) where q = (1/n, . . ., 1/n) and p(x) = (p 1 (x), . . ., p n−1 (x)) with Therefore, for the case z 1 = 0, z 2 = x, z 3 = 1, z 4 = ∞, Eq. (4.14) simplifies to where Let us now focus on the prime form E(z, z ).It can be defined as [114] E(z, z where Θ 1 2 is a shorthand notation for the Theta function of Eq. (4.4) with both characteristics equal to (1/2, 0, . . ., 0) ∈ (Z/2) n−1 and g(z) is Notice that the holomorphic differentials ν r (z) in Eq. (4.13) and, therefore, g(z) are singular at the branch points of the curve that defines the surface Σ n (x).This means that the correlation function (4.14) is in principle not well-defined.In order to solve this issue, the vertex operators inserted at the branch points have to be regularised by redefining them as a proper limit from a non-singular point.We will first extract and remove from g(z) the divergent terms at the branch points.Then, by considering the limit in which the distance between A 1 and A 2 tends to infinity, we will fix the correct definition of the regularised vertex operators at the branch points.
Close to the branch points z = 0, x, 1, and ∞, the holomorphic normalised differentials ν r (z) of Eq. (4.13) behave as with and Observe that, in the four singularities, the divergent term when → 0 is a global factor and, once we take it out, the subleading corrections in vanish.Therefore, these singularities can be removed in the correlation function of Eq (4.14) by defining the vertex operators at the branch points as the limit In this definition, we have included a possible global rescaling factor κ n , which may depend on the genus of the surface, and we will adjust by studying the limit of large separation between the two intervals.If we replace in Eq. (4.21) the vertex operators by the regularised ones introduced in Eq. (4.28), then the resulting correlation function can be written in the form where h T = h α 1 + h α 2 + h α 3 + h α 4 and E ( * ) (z j , z j ) stands for the regularised prime form with and the expressions of ν ( * ) (z) at the branch points are those given in Eq. (4.26).
In Appendix C, we conjecture and numerically check the following identities for the regularised prime forms that appear in Eq. (4.29), and Plugging them into Eq.(4.29), we finally find In Eq. (4.9), once the vertex operators V α are replaced by the regularised ones V ( * ) α , we can exploit Eq. (4.36) to get In particular, when α 1 = −α 2 = α and α 3 = −α 4 = β, we get the multi-charged moments After a global conformal transformation to a subsystem A with arbitrary end-points (u 1 , v 1 , u 2 , v 2 ), we obtain the following result Note that the rescaling factor κ n , which was introduced in the definition of the regularised vertex operators at the branch points, is still undetermined.We can fix it by analysing the limit in which the two intervals A 1 and A 2 are far, i.e. d → ∞, as done for the Dirac theory.In that case, the charged moments Z A 1 :A 2 n (α, β) must verify Eq. (3.19).Since F n (0) = 1 and the constant c n;α,β factorises into those for the intervals A 1 and A 2 , then Eq. (4.38) satisfies the limit d → ∞ of Eq. (3.19) In conclusion, for the massless compact boson, the partition function on the surface Σ n with general twisted boundary conditions is of the form and the multi-charged moments are When the Luttinger parameter is K = 1, then F n (x) = 1, and the partition function of Eq. (4.40) is equal to the one obtained in Eq. (3.15) for the massless Dirac field, as we anticipated at the beginning of this section.Interestingly, the factor in Eq. (4.40) due to the magnetic fluxes is the same, when α = β, as the one derived in Ref. [78] for a large central charge CFT with the Luttinger parameter K replaced by the level of the Kac-Moody algebra of that theory.

Symmetry resolution
In this section, we apply the approach described in Sec.2.1 in order to evaluate the symmetry resolution of the mutual information in the two CFTs analysed in Secs. 3 and 4 from the expressions obtained there for their multi-charged moments Z A 1 :A 2 n (α, β).

Fourier transforms
The first step is to determine the Fourier transform (2.12) of the multi-charged moments.We need to know how the non-universal constant c n;α,β does depend on α and β.In Sec. 3, we have concluded that, for the tight-binding model, it can be well approximated if we only take into account the quadratic terms in α and β.In the following, we will assume that this is in general a good approximation [40].Therefore, we will take In the case of the tight-binding model (K = 1), we obtained in Eq. (3.22) that λ n = e ζn .Therefore, applying Eq. (5.1) in the result of Eq. (4.40), the multi-charged moments can be rewritten as where ˜ p = λ n p .The evaluation of Eq. (2.12) using the expression above yields the following multivariate Gaussian function for the Fourier modes of the multi-charged moments . (5.3) Notice that the Luttinger parameter K enters in the Gaussian factor as an overall rescaling of its variance.
In the limit of large separation between the intervals, i.e. d → ∞ (x → 0), Eq. ( 5.3) tends to namely Z A 1 :A 2 n (q 1 , q 2 ) factorises into the contributions of A 1 and A 2 .This is consistent with the probabilistic interpretation for the case n = 1: the outcomes of the charge measurements in the two intervals are independently distributed when the separation between A 1 and A 2 is large enough.On the other hand, in the limit of two adjacent intervals, i.e. d → 0 (x → 1), the multi-charged moments have the form (see also Eq. (3.18)) whose Fourier transform is (5.6) Setting α = β in Eq. ( 5.2), we obtain the charged moments (2.8) with a single flux In this case, performing the Fourier transform of Eq. (2.9), we end up with Taking n = 1 in this expression, we obtain the probability p(q) of having charge q in the subsystem A, namely p(q) = Z 1 (q).Now we can plug it together with the result for Z A 1 :A 2 n (q 1 , q 2 ) found in Eq. ( 5.3) into Eq.(2.14) to obtain the conditional probability p(q 1 , q 2 ) of having charge q 1 and q 2 = q − q 1 in the intervals A 1 and (5.9)In this expression, ˜ p = λ 1 p ; in particular, for the tight-binding model λ 1 = e ln 2+1+γ E , with γ E the Euler-Mascheroni constant.As a non-trivial consistency check, we have verified that the probability functions we have obtained satisfy the normalisation conditions be derived from the Fourier transform of the charged moments determined in Eq. (5.8) by applying Eq. (2.10).The three symmetry-resolved entropies can eventually be written in the form (5.11) where S X 1 is the total entanglement entropy of subsystem X.In this expression, when X = A p , we have to take Λ = p and σ = 1 while, if X = A 1 ∪ A 2 , then Λ = 1 2 (1 − x) and σ = 2.The auxiliary quantities δ and ξ in Eq. (5.11) are defined in terms of λ n as For the tight-binding model, we know the explicit value of these non-universal constants, From Eqs. (5.9) and (5.11), we can now obtain an explicit expression for the symmetryresolved mutual information.Since the conditional probability p(q 1 , q − q 1 ) satisfies Eq.
(5.10), we have where I A 1 :A 2 is the total mutual information of Eq. (1.3) and we have introduced the rescaled subsystem length ˜ δ p = δ p .We plot this function in Fig. 5.As we anticipated, the symmetry-resolved mutual information is not a good measure of the total correlations between A 1 and A 2 in each symmetry sector since it can assume negative values.Finally, we can also derive the number mutual information defined in Eq. (2.7).Applying in that formula the result for I A 1 :A 2 (q) obtained in Eq. (5.14), we have (5.15) then Eq. (5.15) becomes . (5.17) In the limit 1 , 2 , d → ∞, this expression behaves as This result resembles the one for the number entropy of a single interval (see e.g.[45]), where a double logarithmic correction in the subsystem length also appears, even though, in our case, the dependence on the parameters 1 , 2 , d is more involved.On the other hand, it is a simple function of the Luttinger parameter K since, as we already pointed out, the only effect of K in the Gaussian factor of the Fourier transform of the multi-charged moments is renormalising its variance.

Conclusions
In this manuscript, we have computed the multi-charged moments of two intervals in the ground state of the free massless Dirac field and the massless compact boson, with arbitrary compactification radius.Using the replica approach, the multi-charged moments are given by the partition function of the theory on a higher genus Riemann surface with a different magnetic flux inserted in each interval.We have carried out the analysis of such partition function for the two CFTs under investigation in full generality, allowing the background magnetic flux to generate a different twisted boundary condition at each end-point of the intervals.In the case of the Dirac field, we have adapted the diagonalisation in the replica space method of Ref. [24], to account for the different monodromy of the fields at each endpoint.In the compact boson theory, we have chosen a geometric approach, and we have directly considered the four-point correlator on the Riemann surface of the vertex operators that implement the flux.It turns out that the known formulae for such correlator diverge in our case [115].Once we properly regularised them, we have obtained a cumbersome expression for the multi-charged moments in terms of Riemann-Siegel Theta functions.Nevertheless, we have found several remarkable identities concerning the prime forms of the Riemann surface that allow to dramatically simplify the final result for the multicharged moments.The factor due to the magnetic fluxes is eventually an algebraic function of the lengths of the intervals and their separation-identical to the one obtained in the Dirac theory.In fact, for a certain value of the compactification radius, the multi-charged moments of both theories are equal, generalising the known identity between their twointerval Rényi entropies [29].
Given the simple expression obtained for the multi-charged moments, we can easily calculate their Fourier transform, which has a Gaussian form.From it, we have finally derived formulae for several interesting quantities such as the joint probability distribution of the charge for simultaneous measures in the two intervals, the symmetry-resolved mutual information [82] and the number mutual information.
Let us conclude this manuscript discussing few outlooks.The multi-charged moments analysed here can be used to study the symmetry decomposition of the negativity in imbalance sectors [65].This is a measure of entanglement in mixed states which involves the partial transposition of the reduced density matrix.In the replica approach, this operation can be performed by properly fixing the different fluxes and exchanging the end-points of the transposed interval.A further generalisation is to identify the holographic dual of the multicharged moments, which would be the starting point to compute the symmetry-resolved mutual information in the AdS/CFT correspondence, as done for the entanglement entropy in Ref. [77].Partition functions with twisted boundary conditions, as the ones considered here, have been also proposed as non-local order parameters to distinguish various topological phases of spin chains [104,105].We think that our analysis for the multi-charged moments can be useful to make progresses also in that direction.Finally, even though far beyond the scope of this manuscript, it would be interesting to obtain a rigorous proof of the identities for the prime forms that we have found and numerically checked here.
Applying the rules for the composition of Gaussian operators [121] and introducing W = 2C A − I, we get [82] Z A 1 :A 2 1 (α, β) = (e −iα/2 + e iα/2 ) 1 (e −iβ/2 + e iβ/2 ) 2 det and the notation W pp , p, p = 1, 2, refers to correlations between sites in A p and A p .This result allows to exactly compute the multi-charged moments in the tight-binding model for different values of α and β, as showed in Fig. 2. The Fourier transform of Z A 1 :A 2 1 (α, β) gives the quantities Z A 1 :A 2 1 (q 1 , q 2 ) analysed in Fig. 3 and Fig. 4.

B Derivation of the normalised holomorphic differentials
In this Appendix, we present the detailed derivation of the expression given in Eq. (4.13) for the normalised holomorphic differentials ν r (z) of the Riemann surface Σ n (x).Recall that, according to Eq. (4.11), they are normalised with respect to the homology basis {a r , b r }.In order to define this basis, it is convenient to introduce an auxiliary basis of non-contractible cycles, which we call ãr and br .The cycle ãr encloses anticlockwise the branch cut (u 1 , v 1 ) in the r-sheet of the surface.The dual cycle br connects anticlockwise the r and n sheets through the branch cut (u 1 , v 1 ) and then it goes back to the r sheet through the branch cut (u 2 , v 2 ).In Fig. 6, we draw the auxiliary homology basis in the case n = 3.The cycles a r and b r are defined in terms of the auxiliary ones by We can now obtain ν r (z) by following the usual procedure considered in the literature, see e.g.[117,118].We first take the basis of canonical holomorphic differentials of the surface Σ n (x),

C Prime form identities
The identities for the regularised prime forms E ( * ) (z, z ) of Eqs.The results for n = 2 led us to conjecture the generalisation of Eqs.(4.32)-(4.35)for higher genus.Unfortunately, we have not been able to find a proof for them, although they can be tested numerically with great accuracy, as we show in Fig. 7.The dots are the result of evaluating numerically the definition in Eq. (4.30) of E ( * ) (z, z ), while the solid lines are the functions on the right hand side of the identities in Eqs.

Figure 1 .
Figure 1.Representation of the n-sheeted Riemann surface Σ n for n = 3.The red edge of each cut is identified with the blue edge of the corresponding cut in the lower copy.The calculation of the multi-charged moments in Eq. (2.11) requires to insert different magnetic fluxes between the sheets, which we indicate by the arrows.The operator e iαQ A 1 is implemented by the flux insertions α and −α along the left interval, while e iβQ A 2 corresponds to the fluxes β and −β at the right interval. )

Figure 2 .
Figure 2. Analysis of Z A1:A2 1 (α, β) for the tight-binding model.In the top and middle panels, we show Z A1:A2 1 (α, β) as a function of α at fixed β for different intervals of lengths 1 , 2 , separated by a distance d > 0. In this case, the solid lines are the theoretical predictions of Eq. (3.16), taking Eq.(3.20) for c n;α,β .In the bottom left panel, we repeat the analysis but considering adjacent intervals (d = 0) and the solid curves correspond to Eq. (3.18), with cn;α,β conjectured in Eq. (3.24).In the bottom right panel, we plot Z A1:A2 1 (α, β) as a function of the cross-ratio x.Here the solid lines correspond to Eq. (3.16) using the exact expression for c n;α,β of Eq. (3.20), while for the dashed curves we have considered instead the quadratic approximation for this constant of Eq.(3.22).In all the cases, the points are the exact numerical values for Z A1:A2 1 (α, β) calculated as described in Appendix A.

Figure 5 .
Figure 5. Symmetry-resolved mutual information of Eq. (5.14) in the tight-binding model.We plot our analytical prediction for different combinations of 1 , 2 as a function of the cross-ratio x (left panel) and q (right panel).

Figure 7 .
Figure 7. Numerical check of Eqs.(4.32)-(4.35)involving the regularised prime form E ( * ) (z, z ).We test them for different values of the genus n − 1.The points correspond to the direct numerical evaluation of the definition in Eq. (4.30) of E ( * ) (z, z ).The continuous lines are the plot of the functions on the right hand side of Eqs.(4.32)-(4.35).