Boundary effects on symmetry resolved entanglement
BBoundary effects on symmetry resolved entanglement
Riccarda Bonsignori and Pasquale Calabrese , September 21, 2020 International School for Advanced Studies (SISSA) and INFN, Via Bonomea 265, 34136 Trieste,Italy International Centre for Theoretical Physics (ICTP), Strada Costiera 11, 34151 Trieste, Italy
Abstract
We study the symmetry resolved entanglement entropies in one-dimensional systems with bound-aries. We provide some general results for conformal invariant theories and then move to a semi-infinite chain of free fermions. We consider both an interval starting from the boundary and awayfrom it. We derive exact formulas for the charged and symmetry resolved entropies based ontheorems and conjectures about the spectra of Toeplitz+Hankel matrices. En route to charac-terise the interval away from the boundary, we prove a general relation between the eigenvalues ofToeplitz+Hankel matrices and block Toeplitz ones. An important aspect is that the saddle-pointapproximation from charged to symmetry resolved entropies introduces algebraic corrections to thescaling that are much more severe than in systems without boundaries. a r X i v : . [ c ond - m a t . s t a t - m ec h ] S e p ontents A = [1 , (cid:96) ] The characterisation of the interplay between entanglement and internal symmetries has recently be-come the focus of an intense research activity [1–18] aimed to have a deeper resolution of the structureof the reduced density matrix of many-body systems and quantum field theories. The theoretical workin this area has almost entirely focussed on systems with periodic boundary conditions (PBC). How-ever, there are many fundamental reasons to investigate systems with boundaries, in particular openboundary conditions (OBC). Just to quote a few: experimental solid-state systems typically have OBC;in trapped cold atoms, the vanishing of the density beyond a trapping length induces OBC in the inho-mogeneous gas that can be treated with the methods of field theories in curved space [19]; in some nonequilibrium situations such as a quantum quench, the initial state can be treated as a boundary statein imaginary time formalism [20]. The goal of this paper is to study the effects of boundary conditions,in particular open, on symmetry resolved entanglement in conformal field theories and in free fermionschains. The latter can be treated with simple exact methods.The focus of our work is the entanglement entropy . Given a quantum system described by a purestate density matrix ρ = | ψ (cid:105)(cid:104) ψ | and a bipartition of the Hilbert space H = H A ⊗ H B , the entanglement2igure 1: The two geometries we consider in this manuscript for a semi-infinite system. The subsystem A is always an interval of length (cid:96) that starts either from the boundary (top) or at distance (cid:96) from it(bottom with (cid:96) = 3 ).entropy is the Von Neumann entropy of the reduced density matrix ρ A = Tr B ρ , i.e. S vN = − Tr ρ A log ρ A , (1.1)which is the limit for n → of a larger family of entropies, known as Renyi entropies S n = 11 − n log Tr ( ρ nA ) . (1.2)The latter give more information than the entanglement entropy, since their knowledge for different n provides the full spectrum of the reduced density matrix ρ A [21].One of the remarkable results about the (Rényi) entanglement entropy in extended quantum systemsis its universal scaling behaviour in conformal field theory (CFT) when the subsystem A is an intervalof length (cid:96) embedded in the infinite line [22–25] S n = c (cid:18) n (cid:19) log (cid:96) + c (cid:48) n , (1.3)where c is the central charge and c (cid:48) n is a non-universal additive constant. The presence of a boundarystrongly affects the above scaling behaviour. For example, for a semi-infinite boundary CFT starting(say) at x = 0 , the scaling of the Renyi entropies for the interval A = [0 , (cid:96) ] (see Fig. 1) is [22, 23]: S n = c (cid:18) n (cid:19) log (cid:96) + ˜ c (cid:48) n , (1.4)where the non universal constants ˜ c (cid:48) n are related to the c (cid:48) n in Eq. (1.3) by the universal relation [22, 26] ˜ c (cid:48) n − c (cid:48) n g, (1.5)and log g is the boundary entropy [27, 28]. For the tight-binding chain with OBC that we are going toconsider, we have g = 1 [27]. In a microscopic gapless model (e.g., a spin chain) the leading correctionsto the asymptotic conformal behaviour given by Eqs. (1.3) and (1.4) decay as (cid:96) − K/n [30–34] and (cid:96) − K/n [29, 34, 35] for PBC and OBC, respectively ( K is the scaling dimension of a relevant operator aswe will see later on). These corrections oscillate with (cid:96) , except for PBC at n = 1 .When the subsystem A is placed at distance (cid:96) from the boundary (i.e. A = [ (cid:96) , (cid:96) + (cid:96) ] , see Fig. 1),CFT predicts the general scaling form for the entanglement entropies S n = c (cid:18) n (cid:19) log 4 (cid:96) (cid:96) ( (cid:96) + (cid:96) )(2 (cid:96) + (cid:96) ) + c (cid:48) n + log ˜ F n ( x )1 − n , (1.6)3here x ∈ [0 , is the anharmonic ratio x = (cid:96) (2 (cid:96) + (cid:96) ) , (1.7)and ˜ F n ( x ) is a function depending on the full operator content of the considered boundary CFT andnormalised as ˜ F n (0) = 1 (to recover the bulk result (1.3) for (cid:96) (cid:29) (cid:96) , i.e. x → ). Such a universalfunction should be calculated on a case by case basis, analogously to the one for two disjoint intervalsin periodic systems [36–38]. Some general results for ˜ F n ( x ) in the boundary compact boson CFT(Luttinger liquid) appeared only very recently [39, 40].In this work, we consider the symmetry resolved entanglement of CFT and free fermionic systemswith boundaries. In Sec. 2 we provide all definitions about symmetry resolved entanglement and reviewthe necessary results from the literature. In Sec. 3 we present our results within boundary CFT. In Sec.4 we move to free fermions on the semi-infinite line and consider a block starting from the boundarywhich can be analysed with a generalisation of the Fisher-Hartwig formula. In order to study the caseof a block away from the boundary, we first prove in Sec. 5 a general relation between the spectraof certain Toeplitz+Hankel matrices (related to the case of interest) and block Toeplitz ones (relatedto two disjoint blocks in an infinite chain). Hence in Sec. 6 we calculate the charged entropy for twoblocks in an infinite chain, which are then used in Sec. 7 to infer the results for the interval away fromthe boundary. We conclude in Sec. 8 with a summary of the results and further discussions. We consider an extended quantum system possessing an internal U (1) symmetry, generated by a localoperator Q . The system is taken in a pure state described by a density matrix ρ with a definite valueof this conserved charge and hence [ ρ, Q ] = 0 . We consider a bipartition in A and B and the chargeoperator Q itself splits in the sum Q = Q A ⊗ I B + I A ⊗ Q B of the charge operators Q A , Q B associatedto each subsystem. Consequently, [ ρ A , Q A ] = 0 , implying that the reduced density matrix ρ A has ablock diagonal form, in which each block corresponds to an eigenvalue q of Q A , i.e. ρ A = ⊕ q Π q ρ A Π q = ⊕ q [ p ( q ) ρ A ( q )] , (2.1)where Π q is the projector on the eigenspace of the eigenvalue q and p ( q ) = Tr (Π q ρ A ) is the probabilitythat a measurement of Q A gives the eigenvalue q as outcome. Each block ρ A ( q ) of the reduced densitymatrix is normalised so that Tr ρ A ( q ) = 1 . The block decomposition of the reduced density matrixcan be exploited to quantify the contributions of the different charge sectors to the total entanglemententropy. In fact, Eq. (2.1) allows us to rewrite the total entanglement entropy as [41, 42] S vN = (cid:88) q p ( q ) S vN ( q ) − (cid:88) q p ( q ) log( q ) ≡ S c + S f , (2.2)where we have introduced the symmetry resolved entanglement entropy as the entanglement entropyassociated to the block ρ A ( q ) S vN ( q ) = − Tr [ ρ A ( q ) log ρ A ( q )] . (2.3)4he two terms in Eq. (2.2) are called configurational entanglement entropy ( S c ) [42–44] and fluctuation(or number) entanglement entropy ( S f ) [42, 45] respectively. The former measuring the total entropydue to each charge sector weighted with their probability and the latter the one due to the fluctuationsof the charge within the subsystem A . Similarly, we also define the symmetry resolved Rényi entropies as S n ( q ) = 11 − n log Tr [ ρ A ( q ) n ] . (2.4)The evaluation of the symmetry resolved Rényi and entanglement entropies from the previous def-initions would require the knowledge of the resolution of the spectrum of ρ A in Q A , which is notstraightforward because of the nonlocal nature of the projector Π q . An alternative route [2, 3] is basedon the computation of the the charged moments Z n ( α ) ≡ Tr [ ρ nA e iαQ A ] , (2.5)whose Fourier transform Z n ( q ) = (cid:90) π − π dα π e − iqα Z n ( α ) ≡ Tr [Π q ρ nA ] , (2.6)readily provides the symmetry resolved quantities in Eqs. (2.3) and (2.4) as S n ( q ) = 11 − n log (cid:20) Z n ( q ) Z ( q ) n (cid:21) , S vN ( q ) = − ∂ n (cid:20) Z n ( q ) Z ( q ) n (cid:21) n =1 . (2.7)Notice that the probability p ( q ) is nothing but p ( q ) = Z ( q ) . In the next subsection we show how toevaluate the charge moments using the replica trick. In the replica approach, the moments Tr ρ nA can be evaluated for any (1 + 1) -dimensional quantum fieldtheory (QFT) as partition functions over a suitable n -sheeted Riemann surface R n in which the n sheets(replicas) are cyclically joined along the subsystem A [22,23]. Similarly [2], the charged moments find ageometrical interpretation by inserting an Aharonov-Bohm flux through such surface, so that the totalphase accumulated by the field upon going through the entire surface is α . Then the partition functionon such modified surface is nothing but the charged moments Z n ( α ) .This partition function can be rewritten in terms of correlator of properly defined twist fieldsimplementing twisted boundary conditions. Assuming, without loss of generality, that the Aharonov-Bohm flux is inserted between the j -th and ( j + 1) -th replicas, we can write the action of the twistfields as [2] T n,α ( x, τ ) φ i ( x (cid:48) , τ ) = (cid:40) φ i +1 ( x (cid:48) , τ ) e iαδ ij T n,α ( x, τ ) , if x < x (cid:48) ,φ i ( x (cid:48) , τ ) T n,α ( x, τ ) , otherwise. (2.8)In terms of these composite twist fields, the charged moments for a single interval A = [0 , (cid:96) ] embeddedin the infinite line are Z n ( α ) = (cid:104)T n,α ( (cid:96),
0) ˜ T n,α (0 , (cid:105) , (2.9)5here ˜ T is the anti-twist field. In particular, in the case of a (1 + 1) -dimensional CFT, it has beenshown that these fields behave as primary operators with conformal dimension ∆ n,α = ∆ n + ∆ α n , ∆ n = c (cid:18) n − n (cid:19) , (2.10)so that they can be written as T n,α = T n V α , where T n are the standard twist fields with conformaldimension ∆ n and V α , with conformal dimension ∆ α , is a field implementing the insertion of theAharonov-Bohm flux. It follows that Z n ( α ) scales as Z n ( α ) = c n,α (cid:96) − c ( n − n ) − αn , (2.11)where c n,α is the normalisation constant of the composite twist field (with c n, = c n ). The previousarguments apply to a generic CFT. In the case of Luttinger liquid conformal field theories [2] (that are c = 1 free compact scalar bosonic massless theories describing the universality class of many 1D criticalsystems of interest, such as free and interacting spin chains), the operator V α can be identified with thevertex operator e iαϕ ( z ) , so that its conformal dimension is ∆ α = (cid:16) α π (cid:17) K, (2.12)where K is the Luttinger parameter, related to the compactification radius of the bosonic theory.From the charged moments (2.11), we get the symmetry resolved moments Z n ( q ) via Fourier trans-form, that in the limit of large (cid:96) is [2] Z n ( q ) (cid:39) (cid:96) − c ( n − n ) (cid:114) nπ K log (cid:96) e nπ q − ¯ q )22 K log (cid:96) , (2.13)where ¯ q ≡ (cid:104) Q A (cid:105) represents the average number of particles. The latter, being a non-universal quantityof the system, cannot be determined within CFT.From Eq. (2.13) we straightforwardly read the leading order of the symmetry resolved (Rényi)entropy (2.7) as S n ( q ) = S n −
12 log (cid:18) Kπ log (cid:96) (cid:19) + O ( (cid:96) ) , (2.14)where S n is the total entropy (1.3). We observe that at leading orders the symmetry resolved en-tanglement does not depend on q , that is, it has the same value in all the different charge sectorscorresponding to different eigenvalues of the charge operator. This result is known as equipartition ofentanglement [3]. The simplest lattice model described by a Luttinger liquid CFT (with K = 1 ) is represented by freespinless fermions hopping on a 1D lattice, with Hamiltonian H = − (cid:88) l c † l c l +1 + c † l +1 c l + 2 h (cid:18) c † l c l − (cid:19) , (2.15)where the fermionic ladder operators c l obey canonical anti-commutation relations { c l , c † m } = δ l,m and h is the chemical potential. The Hamiltonian is straightforwardly diagonalised in the momentum space6nd its ground state is a Fermi sea with momentum k F = arccos | h | . The conserved U (1) local chargeof the model is given by Q = (cid:80) i c † i c i and it can be split in the sum Q = Q A + Q B for any spatialbipartition. The Jordan-Wigner transformation c l = (cid:32) (cid:89) m 12 [ σ xl σ xl +1 + σ yl σ yl +1 ] − hσ zl , (2.17)where σ x,y,zl are the Pauli matrices at site l . Depending on the boundary conditions, the sum over l inthe Hamiltonians can run over a finite, semi-infinite, or infinite number of sites.The reduced density matrix for an arbitrary spatial subsystem can be obtained using Wick’s theo-rem, and has the following form [46–48] ρ A = det C A exp (cid:88) j,l ∈ A [log( C − A − jl c † j c l , (2.18)where C A is the correlation matrix, i.e. the matrix formed by the correlations (cid:104) c † i c j (cid:105) with i, j ∈ A .The fermionic reduced density matrix (2.18) is equal to the spin reduced density matrix for the samesubsystem A only when A is one interval (starting from the boundary if PBC are not imposed), becausethe Jordan-Wigner transformation is local within a compact subset. Conversely, for a non-compactbipartition A ∪ ¯ A (such as for disjoint intervals with PBC or one interval away from the boundary inan open chain) the non-local nature of the Jordan-Wigner string makes the spin and fermion reduceddensity matrices different [49–51].Denoting with | A | the total number of sites within A (i.e. the length (cid:96) for a single interval), C A is a | A | × | A | real and symmetric matrix. This matrix can be diagonalised by an orthogonal transformation.We write the eigenvalues of C A as (1 + ν k ) / , k = 1 , . . . , | A | . Exploiting the quadratic form of ρ A inEq. (2.18), the charged moments are [2] Z n ( α ) = | A | (cid:89) i =1 (cid:20)(cid:18) ν i (cid:19) n e iα + (cid:18) − ν i (cid:19) n (cid:21) = exp | A | (cid:88) i =1 f n ( ν i , α ) , (2.19)with f n ( x, α ) = log (cid:20)(cid:18) ν i (cid:19) n e iα + (cid:18) − ν i (cid:19) n (cid:21) . (2.20)This formula can be straightforwardly used to evaluate the charged moments numerically in any Gaus-sian state and for any bipartition. For α = 0 , we get back the standard formula for the (neutral) totalmoments.Eq. (2.19) also is the starting point for the analytic computation of Z n ( α ) using the Fisher-Hartwigformula, as done in Ref. [5]. To this aim, one first introduces the determinant [53] D A ( λ ) = det[( λ + 1) I − C A ] = | A | (cid:89) j =1 ( λ − ν j ) ≡ det( G ) , (2.21)7hich is a polynomial of degree (cid:96) in λ whose zeroes are the eigenvalues { ν j , j = 1 , · · · , (cid:96) } . Then werewrite Eq. (2.19) in integral form log Z n ( α ) = 12 πi (cid:73) dλf n ( λ, α ) d log D A ( λ ) dλ , (2.22)where the contour integral encircles the segment [ − , which is the support of the ν j .So far, everything was completely general and applies to an arbitrary Gaussian state for an arbitraryspatial subsystem A . There is even no reference to the boundary conditions. In the following we reviewsome exact results valid for the ground state of an infinite chain, focusing on those aspects we will needfor the generalisation to open chains. We now specialise to the ground state of the infinite free-fermion chain and for the case of A being aninterval of length (cid:96) . The reduced correlation matrix is ( C A ) ij = sin ( k F ( i − j )) π ( i − j ) , (2.23)with i, j = 1 . . . (cid:96) . The matrix G in Eq. (2.21) has a Toeplitz form, i.e. its elements depend onlyon the difference between row and column indices G jk = g j − k . For this matrix, the integral (2.22)can be evaluated analytically in the asymptotic limit (cid:96) → ∞ using the Fisher-Hartwig formula [52], atechnique that has been used extensively to evaluate entanglement in free lattice models [31, 35, 53–59].We briefly recap this derivation in the following in order to illustrate the procedure and set the notationthat will be used for open systems.The Fisher-Hartwig formula is written in terms of the symbol of the Toeplitz matrix, defined as theFourier transform of g l g l = (cid:90) π − π dθ π e ilθ g ( θ ) , (2.24)that in the case considered here is given by g ( θ ) = (cid:26) λ + 1 , θ ∈ [ − π, − k F ] ∪ [ k F , π ] ,λ − , θ ∈ [ − k F , k F ] . (2.25)In the integration domain (2.24) the symbol has two discontinuities, located at θ = − k F and θ = k F .The Fisher-Hartwig formula relies on the possibility to express the symbol in the following form g ( θ ) = f ( θ ) R (cid:89) r =1 e ib r [ θ − θ r − π sgn( θ − θ r )] (2 − θ − θ r )) a r , (2.26)where R is an integer, a r , b r , θ r are constants and f ( θ ) is a smooth function with winding number zero.For the symbol g ( θ ) in Eq. (2.25), there are two discontinuities so that R = 2 , and the constantsassume the values a , = 0 , b = − b = β λ + m and f ( θ ) = f = ( λ + 1) e − ik F m e − ik F β λ , where β λ = 12 πi log (cid:20) λ + 1 λ − (cid:21) , with dβ λ dλ = 1 πi − λ , (2.27)and the integer m ∈ Z labels the different inequivalent representations of the symbol. Usually one referssimply to the Fisher-Hartwig formula when there is a single representation of the symbol and to the8eneralised one when there are multiple representations, as it is the case for us. For a Toeplitz matrix T with a symbol of the form (2.25) without inequivalent representation, the Fisher-Hartwig formulaprovides the large (cid:96) behaviour det T (cid:39) F [ f ( θ )] (cid:96) R (cid:89) j =1 (cid:96) a j − b j , where F [ f ( θ )] = exp (cid:18) π (cid:90) π dθ log f ( θ ) (cid:19) . (2.28)When the symbol has several inequivalent representations, as for our case, one has to sum over all ofthem [31, 52], obtaining, in our specific case, the asymptotic for large (cid:96)D A ( λ ) (cid:39) ( λ + 1) (cid:96) (cid:18) λ + 1 λ − (cid:19) − kF (cid:96)π (cid:88) m ∈ Z (2 (cid:96) | sin k F | ) − m + β λ ) e − ik F m(cid:96) [ G ( m + 1 + β λ ) G (1 − m − β λ )] , (2.29)where G ( z ) is the Barnes G -function.The charged moments Z n ( α ) are evaluated by inserting the result for D A ( λ ) (2.29) into the integral(2.22). It is easy to see that, for α ∈ [ − π, π ] , the leading behaviour for large (cid:96) of such integral is givenby the term with m = 0 D (0) A ( λ ) ≡ ( λ + 1) (cid:96) (cid:18) λ + 1 λ − (cid:19) − kF (cid:96)π (2 (cid:96) | sin k F | ) − β λ [ G (1 + β λ ) G (1 − β λ )] , (2.30)and so the contour integral (2.22) gives [5] log Z (0) n ( α ) = iα k F (cid:96)π − (cid:20) (cid:18) n − n (cid:19) + 2 n (cid:16) α π (cid:17) (cid:21) log(2 (cid:96) | sin k F | ) + Υ( n, α ) . (2.31)In this expression for the charged moment, we recognise immediately the average number of particles(the linear term in α ) given by ¯ q ≡ (cid:104) Q A (cid:105) = k F (cid:96)/π , and the dimension of the modified twist field fromthe term proportional to log(2 (cid:96) | sin k F | ) . The non-universal constant Υ( n, α ) is given by the integral Υ( n, α ) = ni (cid:90) ∞−∞ dw [tanh( πw ) − tanh( πnw + iα/ + iw )Γ( − iw ) , (2.32)that is real as long as α is real. For later purposes, we rewrite it as Υ( n, α ) = Υ( n ) + γ ( n ) α + (cid:15) ( n, α ) , (cid:15) ( n, α ) = O ( α ) , (2.33)where γ ( n ) = ni (cid:90) ∞−∞ dw [tanh ( πnw ) − tanh( πnw )] log Γ( + iw )Γ( − iw ) . (2.34)For α = 0 , the above results reproduce the well known (total) Rényi entropies [31, 53].The symmetry resolved moments are just the Fourier transform Z n ( q ) of Z n ( α ) . In this Fouriertransform we ultimately use a saddle-point approximation in which Z n ( α ) is Gaussian and hence wetruncate hereafter Z (0) n ( α ) in Eq. (2.31) at quadratic order in α . Consequently, the charged partitionfunction can be well approximated as Z n ( α ) = Z n (0) e iα ¯ q − b n α / , (2.35)9here b n = bπ n ln (cid:96) − h n , (2.36)with b = 1 , h n ≡ − π n ln(2 | sin k F | ) + 2 γ ( n ) . (2.37)(Here we slightly change our notations compared to Ref. [5] and follow more closely those in [9].) Inthis way, the Fourier transform is a simple Gaussian integral, with result Z n ( q ) = Z n (0) √ πb n e − ( q − ¯ q )22 bn . (2.38)The symmetry resolved Rényi entanglement entropies are then obtained using Eq. (2.7) as S n ( q ) = S n − 12 ln(2 π ) + 11 − n ln b ( (cid:96), t ) n/ b n ( (cid:96), t ) / − q − n ) (cid:18) b n ( (cid:96), t ) − nb ( (cid:96), t ) (cid:19) , (2.39)where S n is the total entropy. Expanding for large (cid:96) , we have S n ( q ) = S n − 12 ln (cid:16) bπ ln δ n (cid:96) (cid:17) + ln n − n ) − π n ( h − nh n ) − n ) ( b ln (cid:96) ) ++ ( q − ¯ q ) nπ h − nh n − n )( b ln (cid:96)κ n ) + o (ln (cid:96) − ) , (2.40)where, following Ref. [5] we absorbed some subleading corrections in the amplitudes as ln δ n = − π n ( h n − h ) b (1 − n ) , ln κ n = − π ( h + nh n )2 b . (2.41)The above formula is valid also for the symmetry resolved Von Neumann entropy taking properly thelimits of the various pieces as n → .We wrote these formulas in a rather generic fashion as a function of b n because in the other casesstudied here (and elsewhere [8, 9, 13]) only the specific form of b n (or h n ) matters and the final resultis always given by Eq. (2.40) with the minor redefinition of the amplitudes. In this section we present our boundary CFT results for the charged and symmetry resolved entropies.We start from a 1D system in a semi-infinite line [0; ∞ ) and a subsystem A consisting in a finite interval [0; (cid:96) ) as in Fig. 1. From general CFT scaling, one expects the charged moments Z n ( α ) to be Z n ( α ) ≡ Tr (cid:2) ρ nA e iQ A α (cid:3) = (cid:104)T n,α ( (cid:96) ) (cid:105) HP = ˜ c n,α (2 (cid:96) ) − c ( n − n ) − ∆ αn , (3.1)where the subscript HP stands for the average over the (right) half plane z = x + iτ with x ∈ R + and τ ∈ R . Anyhow, it is worth to obtain such (correct) prediction from first principles by merging theboundary CFT approach to the entanglement entropy [22,23] with the insertion of the Aharonov-Bohmflux [2].The n -sheeted Riemann surface then consists of n copies of the half plane x ≥ sewn together along ≤ x ≤ (cid:96), τ = 0 . Once again, the flux between the j -th and ( j + 1) -th replicas can be implemented10y the definition of local composite twist fields T n,α = T n V α at the end-point of the region A , where V α is responsible for the Aharonov-Bohm flux and T n generates the Riemann geometry. The scalingdimension of the composite twist field is obtained by evaluating the expectation value of the totalstress-energy tensor T ( w ) = (cid:80) nj =1 T j ( w ) in the Riemann geometry R n with inserted flux α . First, thetransformation z = (cid:18) (cid:96) − w(cid:96) + w (cid:19) /n , (3.2)maps the whole Riemann surface into the unit disc D = {| z | < } with the flux α . Thus we have (cid:104) T ( w ) (cid:105) R n,α = (cid:88) j (cid:18) dzdw (cid:19) (cid:104) T j ( z ) (cid:105) D ,α + nc { z, w } , (3.3)where the Schwartzian derivative is given in this case by c { z, w } = c (cid:18) − n (cid:19) (2 (cid:96) ) ( w − (cid:96) ) ( w + (cid:96) ) . (3.4)In Eq. (3.3) we added the subscript α to stress that the expectation values are taken in the presenceof a flux. Notice that (cid:104) T j ( z ) (cid:105) D ,α is non-zero for α (cid:54) = 0 . To calculate it, let us start by noticing that thetransformation (3.2) maps the subsystem A (say on the first sheet) into [0 , and the branch point intothe origin. Hence, a closed path encircling the branch point n times is mapped into a single-windingorbit around the origin so that (cid:104) T j ( z ) (cid:105) D ,α = (cid:104) T j ( z ) V α (0) (cid:105) D (cid:104)V α (0) (cid:105) D , (3.5)where the subscript D (without α ) refers to the expectation values on the unit disk in the absence ofthe flux. When V α is a primary operator, the right hand side of Eq. (3.5) is h α /z [60] and hence (cid:104) T j ( z ) (cid:105) D ,α = h α z . (3.6)Plugging the above expression and the Schwartzian derivative (3.4) into Eq. (3.3) ones get (cid:104) T ( w ) T n,α ( (cid:96) ) (cid:105) HP (cid:104)T n,α ( (cid:96) ) (cid:105) HP = (cid:20) c (cid:18) n − n (cid:19) + h α n (cid:21) (2 (cid:96) ) ( w − (cid:96) ) ( w + (cid:96) ) . (3.7)The comparison of Eq. (3.7) with the conformal Ward identity for boundary CFT [60] confirms thatalso in the presence of boundaries, the scaling dimension of the composite twist fields is (2.10), leadingimmediately to the expected result (3.1). Taking the Fourier transform by saddle point approximation,one obtains the asymptotic symmetry resolved moments Z n ( q ) (cid:39) (2 (cid:96) ) − c ( n − n ) (cid:114) nπK log (cid:96) e nπ q − ¯ q )2 K log (cid:96) . (3.8)The other situation of interest is that of an interval of length (cid:96) placed at distance (cid:96) from theboundary. In this case, global conformal invariance fixes the overall scaling. Indeed, the chargedmoment is a two-point function Z n ( α ) = (cid:104)T n,α ( u ) ˜ T n,α ( v ) (cid:105) HP with u = (cid:96) and v = (cid:96) + (cid:96) . By imagecharges technique, this is related to a four-point function on the plane with images at u = − v and v = − u . The scaling of a general four four-point function of composite twist fields is (cid:104)T n,α ( u ) ˜ T n,α ( v ) T n,α ( u ) ˜ T n,α ( v ) (cid:105) = c n,α (cid:18) | u − u || v − v || u − v || u − v || u − v || u − v | (cid:19) n,α F n,α ( x ) , (3.9)11here F n,α ( x ) is a universal function that depends on the operator content of the CFT and x is theanharmonic ratio of the four points. The desired boundary correlation is related to the square root ofthe bulk four-point function, and hence, specifying to the actual values of u i and v i , we have Tr (cid:2) ρ nA e iQ A α (cid:3) = (cid:104)T n,α ( u ) ˜ T n,α ( v ) (cid:105) HP = c n,α (cid:18) (2 (cid:96) + (cid:96) ) (cid:96) (cid:96) ( (cid:96) + (cid:96) ) (cid:19) ∆ n,α ˜ F n,α ( x ) , (3.10)where ˜ F n,α ( x ) is a universal function of the anharmonic ratio x in (1.7) that depends on the operatorcontent of the boundary CFT. As (cid:96) (cid:29) (cid:96) , using ˜ F n,α (0) = 1 , Eq. (3.10) tends to the results for oneinterval in an infinite system (2.11), as it should. Conversely, for (cid:96) (cid:28) (cid:96) , one recovers the result for theinterval close to the boundary, cf. Eq. (3.1), in which the relation between c n,α and ˜ c n,α (generalisationto α (cid:54) = 0 of Eq. (1.5)) ˜ c n,α c / n,α = g − n , (3.11)is provided by the singular behaviour of F n,α ( x ) close to x = 1 .Finally, let us briefly mention what happens for a finite system of length L with OBC on both sides(which is the most relevant situation for physical applications, as e.g., those in Refs. [19, 20]). Theworldsheet of each replica is an infinite (in the time direction) strip of length L . This can be mappedto the half plane by a conformal transformation (the logarithm). Hence the net effect of having a finitesystem is just to replace all separations in Eqs. (3.1) and (3.10) with the appropriate chord distance,e.g. in Eq. (3.1) we have the replacement (cid:96) → Lπ sin π(cid:96)L , (3.12)and similarly in Eq. (3.10). In the same way, one can consider a semi-infinite system at finite tem-perature, with the worldsheet being a semi-infinite cylinder that can be mapped to the half plane bya conformal transformation. Simple algebra leads to the replacement (cid:96) → βπ sinh π(cid:96)β [22, 23], with β being the inverse temperature. A = [1 , (cid:96) ] The results in the previous section are valid for an arbitrary boundary CFT. In this section, we specialiseto a particular microscopic model whose low energy physics is captured by a CFT. We wish to calculatealso the non-universal factors (not predicted by CFT) that enter in a parameter-free comparison withnumerics. We focus on the tight-binding model given by the Hamiltonian (2.15), which describes freefermions hopping on a 1D lattice. We work with a semi-infinite chain, i.e. l ∈ N , and the first site onthe left is free, i.e. we choose open boundary conditions.The fermion correlation function between two arbitrary sites i, j [46] is ( C A ) ij = sin ( k F ( i − j )) π ( i − j ) − sin ( k F ( i + j )) π ( i + j ) . (4.1)For an arbitrary spatial bipartition A ∪ ¯ A , the elements of the correlation matrix C A are just given bythe above correlation with indices i, j restricted in A . The first term in Eq. (4.1) is the same as in the12nfinite system (cf. Eq. (2.23)) and it depends only on the spatial distance between the two points as aconsequence of translational invariance. The second term breaks the latter symmetry and depends onthe average distance from the boundary (and it vanishes for large i, j at fixed separation | i − j | , whenonly the first term survives).Having constructed the subsystem correlation matrix, the entanglement spectrum and hence theRényi entropies (total and symmetry resolved) follow from the general results in Sec. 2.2. In particular,the charged moments are given by Eq. (2.19) in terms of the eigenvalues ν j of the matrix C A . In thefollowing we specialise to two bipartitions that can be handled analytically, which are the ones depictedin Fig. 1 with one block of (cid:96) sites at distance (cid:96) from the boundary (that can be (cid:96) = 0 , recovering aninterval starting from the boundary). The (cid:96) × (cid:96) correlation matrix may be written as ( C A ) ij = f i − j − f i + j +2 (cid:96) , i, j = 1 , . . . , (cid:96), (cid:96) = 0 , , , . . . (4.2)with f k the elements of the correlation matrix for the infinite system in Eq. (2.23). The first term inthe rhs of Eq. (4.2) is a Toeplitz part, since it depends only on the difference of the indices, and isequal to the matrix for an infinite system; the second term is a Hankel matrix since it depends only onthe sum of the indices and it comes from the presence of the boundary. The correlation matrix C A inEq. (4.2) is then equal to the sum of one Toeplitz matrix and a Hankel one. In the case of a block starting from the boundary, i.e. with (cid:96) = 1 , the Toeplitz+Hankel matrix in Eq.(4.2) has been studied a lot in both the mathematical and physics literature [35,61,62]. The final resultof interest for our paper is the conjectured form for a generalised Fisher-Hartwig expansion [35] D A ( λ ) (cid:39) ( λ + 1) (cid:96) (cid:18) λ + 1 λ − (cid:19) − kF (cid:96)π (cid:88) m ∈ Z e i π ( β + m ) (cid:20) (cid:16) (cid:96) + 12 (cid:17) | sin k F | (cid:21) − ( m + β λ ) × e − ik F ( β + m )( (cid:96) +1 / G ( m + 1 + β λ ) G (1 − m − β λ ) , (4.3)where β λ is defined in Eq. (2.27) and G ( z ) is the Barnes function. This formula was conjectured inRef. [35] in such a way to reproduce the results for m = 0 derived rigorously in Ref. [61] and generalisingheuristically to the existence of inequivalent representations of the symbol. The term / added to (cid:96) has been introduced to absorb some of the subleading corrections to better match numerical evaluationof the determinant and has not a rigorous basis (we will see in the following section another reason forits presence). This asymptotic expansion has been used to compute the moments of the reduced densitymatrix Tr ρ nA , for which the leading term corresponds to m = 0 , while the first corrections are given by m = ± . For the charged moments Z n ( α ) , restricting by periodicity to α ∈ [ − π, π ] , the leading termis still given by m = 0 , see also the discussion in Ref. [6].13 .1.1 The leading term The leading behaviour of the determinant D A ( λ ) is given by D (0) A ( λ ) ∼ e i ( π − k F ) β λ ( λ + 1) (cid:18) λ + 1 λ − (cid:19) − kFπ (cid:96) (4 (cid:96) | sin k F | ) − β λ G (1 − β λ ) G (1 + β λ ) , (4.4)(at the leading order we can drop the / in the generalised form) so that the leading contribution tothe charged moments (2.22) is log Z (0) n ( α ) = 12 πi (cid:73) dλf n,α (1 , λ ) d log D (0) (cid:96) ( λ ) dλ = a + a (cid:96) + a log L k + a , (4.5)where a = 12 π (cid:16) π − k F (cid:17) (cid:73) dλf n,α (1 , x ) ddλ β λ ,a = 12 πi (cid:73) dλf n,α (1 , λ ) (cid:18) − k F /π λ − k F /π − λ (cid:19) ,a = 12 πi (cid:73) dλf n,α (1 , λ ) d ( − β λ ) dλ = 1 π (cid:73) dλf n ( λ, α ) β λ − λ ,a = 12 πi (cid:73) dλf n,α (1 , λ ) d log[ G (1 − β λ ) G (1 + β λ )] dλ . (4.6)Putting everything together, the final result is log Z (0) n ( α ) = iαk F (cid:96)π + (cid:20) iαπ (cid:18) k F π − (cid:19)(cid:21) − (cid:20) (cid:18) n − n (cid:19) + 1 n (cid:16) α π (cid:17) (cid:21) log (cid:2) (cid:96) | sin k F | (cid:3) + Υ( n, α )2 , (4.7)where the real constant Υ( n, α ) is the same as the one in Eq. (2.32). Let us now discuss the variousterms, emphasising the differences with the case of the infinite system in Eq. (2.31). The purelyimaginary term gives the average number of particles in the interval (cid:96) which is ¯ q ≡ (cid:104) Q A (cid:105) = k F π (cid:96) + 1 π (cid:18) k F π − (cid:19) + o (1) . (4.8)The leading piece, proportional to (cid:96) , is the average density also present in the infinite system. The O (1) correction is the more interesting and represents the variation of the number due the inhomogeneityclose to the boundary. Indeed the density at the site j is n j = k F /π + sin(2 k F j ) / (2 πj ) (i.e. thecorrelation (4.1) at coinciding points). Thus, the mean number of particles in [1 , (cid:96) ] is (cid:104) Q A (cid:105) − k F π (cid:96) = (cid:96) (cid:88) j =1 sin(2 k F j )2 πj = ∞ (cid:88) j =1 sin(2 k F j )2 πj + O ( (cid:96) − ) = 1 π (cid:18) k F π − (cid:19) + o (1) . (4.9)At half-filling this term is obviously zero since there is no inhomogeneity.Let us now move to the other pieces in Eq. (4.7). The logarithmic term agrees with the boundaryCFT prediction (3.1) for this specific model and confirms the dimension of the modified twist fields.Interesting enough, also the non-universal constant is related to the one in the absence of boundariesand satisfies the universal relation (3.11) with g = 1 , as well known.14igure 2: Real part of the charged Rényi entropies log Z n ( α ) with the insertion of a flux α as functionsof (cid:96) . The exact numerical results (symbols) are shown for different values of α and n at fillings k F = π/ (left) and k F = π/ (right). The numerical data present a slow approach to the asymptotic behaviour(4.7), with large oscillating corrections to the scaling.In Fig. 2 we present the comparison between the analytical formula (4.7) and the exact numericaldata obtained for the charged Rényi entropies for different values of n and α , at filling k F = π/ and k F = π/ . It is evident that the data are slowly approaching the exact predictions, but there are largeoscillating corrections to the scaling whose amplitude increases with n and α . These corrections donot come as a surprise: they have been intensively studied for the total entanglement entropy both forperiodic [30–34] and open boundary conditions [29,35]. For the charged moments with α (cid:54) = 0 and PBC,they have been described and characterised in Refs. [5, 6]. In the following subsection we will deriveand quantify them for OBC. The corrections to the scaling are encoded in the generalised Fisher-Hartwig formula (4.3). The preciseorder in which they appear depends on the quantity to be calculated. It is straightforward to convinceourselves that for < α < π ( − π < α < ) the leading correction is the one with m = 1 ( m = − ).For α = 0 the terms m = ± are of the same order. If not yet clear, all this will be evident andself-consistent at the end of the calculation in this subsection.We perform the calculation keeping only the two terms with m = ± and approximate the charac-teristic polynomial D A ( λ ) as D A ( λ ) (cid:39) D (0) A ( λ ) (cid:26) ie − ik F (cid:96) e − ik F (cid:96) L − − β λ k Γ(1 + β )Γ( − β λ ) − ie ik F (cid:96) e ik F (cid:96) L − β λ k Γ(1 − β λ )Γ( β λ ) (cid:27) (4.10) ≡ D (0) A ( λ )(1 + Ψ (cid:96) ( λ )) , where the leading term D (0) A ( λ ) is the one in Eq. (4.4) and we introduced the shorthand L k = 4 (cid:16) (cid:96) + 12 (cid:17) | sin k F | . (4.11)15he corrections to the scaling are characterised by the difference d n ( α ) ≡ log Z n ( α ) − log Z (0) n ( α ) , (4.12)that for large L k can be written as d n ( α ) (cid:39) πi (cid:73) dλf n,α (1 , x ) d log[1 + Ψ (cid:96) ( λ )] dλ = 12 πi (cid:73) dλf n,α (1 , λ ) d Ψ (cid:96) ( λ ) dλ + · · · . (4.13)The contour integral can be split as the sum of two contributions above and below the segment [ − , d n ( (cid:96) ) (cid:39) πi (cid:20)(cid:90) i(cid:15) − i(cid:15) − (cid:90) − i(cid:15) − − i(cid:15) (cid:21) dλf n,α (1 , λ ) d Ψ (cid:96) ( λ ) dλ . (4.14)The integral is then given by the discontinuity at the branch cut; the only discontinuous function alongthis cut is β λ , that we rewrite as (for − < x < ) β x ± i(cid:15) = − iw ( x ) ∓ , with w ( x ) = 12 π log 1 + x − x . (4.15)Hence the discontinuities across the branch cut given by the two terms in Ψ (cid:96) are respectively (cid:20) L − i − βk Γ(1 + β ) − β (cid:21) β = − iw − − (cid:20) L − i − βk Γ(1 + β ) − β (cid:21) β = − iw + (cid:39) L iwk γ ( w ) , (4.16) (cid:20) L − i +2 βk Γ(1 + β ) − β (cid:21) β = − iw − − (cid:20) L − i +2 βk Γ(1 + β ) − β (cid:21) β = − iw + (cid:39) L − iwk γ ( − w ) , (4.17)where we dropped terms of order O ( L − k ) compared to the leading ones and defined γ ( w ) = Γ (cid:0) − iw (cid:1) Γ (cid:0) + iw (cid:1) . (4.18)We can now perform the change of variable λ = tanh( πw ) , −∞ < w < ∞ , (4.19)and integrate by parts using ddw f n (tanh( πw ) , α ) = πn [tanh( πnw + iα/ − tanh( πw )] , (4.20)to finally get d n ( α ) (cid:39) in (cid:90) ∞−∞ dw (tanh( πw ) − tanh( πnw + iα/ × (cid:104) ie − ik F e − ik F (cid:96) L iwk γ ( w ) + ie ik F e ik F (cid:96) L − iwk γ ( − w ) (cid:105) . (4.21)This last integral can be performed on the complex plane by residue theorem. For the first piece of theintegral in square brackets, we must close the contour in the upper half plane, while for the second piecein the lower half plane. In principle, we should sum over all residues inside the integration contour,but, since we are interested in the limit of large L k , we can limit ourselves to consider the singularities16igure 3: Leading corrections to the scaling for the charged entropies. The plots show the difference d n ( α ) defined in Eq. (4.12). We focus on α = 1 , n = 1 , and different fillings k F = π/ (left), k F = π/ (center) and k F = π/ (right) as function of (cid:96) . The numerical data (symbols) match wellthe calculated leading correction to the scaling (4.22) from generalised Fisher-Hartwig formula both forreal and imaginary parts.closest to the real axis. For the first integral this is at w = i/ (2 n )(1 − α/π ) while for the second one itis at w = − i/ (2 n )(1 + α/π ) . Their sum gives d n ( α ) (cid:39) ie − ik F (2 (cid:96) +1) L − n ( − απ ) k Γ (cid:0) + n − α πn (cid:1) Γ (cid:0) − n + α πn (cid:1) − ie ik F (2 (cid:96) +1) L − n ( απ ) k Γ (cid:0) + n + α πn (cid:1) Γ (cid:0) − n − α πn (cid:1) . (4.22)This expression represents our final result for the oscillating corrections to the scaling. Let us commentit. For α = 0 the corrections are real and give back the result already found in Ref. [35]. Otherwise,Eq. (4.22) is valid only in the range − π ≤ α ≤ π and must be extended periodically outside ofit. Mathematically, one can found the periodic structure using the entire sum in the generalisedFisher-Hartwig formula (4.3). This periodicity in α is identical to what found for PBC in Ref. [6] towhich we remand for an extensive discussion. Within the principal domain − π ≤ α ≤ π , the twocontributions decay with different power laws (as anticipated above), so that only one of the two isdominating, according to the sign of α . However, they become comparable in magnitude as α approacheszero where one should carefully take into account both of them. Conversely, for α close to ± π , thecalculated corrections with m = ± become comparable with the leading term ( m = 0 ) and hence agood agreement between Eq. (4.22) and the data is achieved only at extremely large (cid:96) , since in thederivation we assumed that the leading term is much larger than the corrections.In Fig. 3 we report the exact numerical data for d n ( α ) at fixed α = 1 and n = 1 , . We observethat, for different fillings, the data always match well the analytic prediction (4.22) both for real andimaginary part. Remarkably, the sampling at commensurate points due to the integer values of (cid:96) , matchwell the oscillating structure due to the partial fillings.17igure 4: Von Neumann (left) and second Rényi (right) symmetry resolved entanglement entropies athalf filling k F = π/ for one interval of length (cid:96) starting from the boundary. The analytical predictions(dashed lines) in Eq. (2.40) with the amplitudes in Eq. (4.23) are compared with the numerical data(symbols) for ∆ q = q − ¯ q = 0 , . The agreement of the saddle point results with the numerical datais rather unsatisfactory. The observed deviations are caused by the finite domain of the integral in α ,and are well captured by Eq. (4.24) reported as full lines which falls in the middle between the datafor even and odd (cid:96) . We are finally ready to evaluate the symmetry resolved moments Z n ( q ) as Fourier transform of Z n ( α ) in the saddle point approximation. In this approximation, the leading term is given by Z (0) n ( α ) in Eqs.(2.35) with the variance b n again given by Eq. (2.36), but with the amplitudes b = 12 , h n = − n π log(4 sin k F ) + γ ( n ) . (4.23)The symmetry resolved moments and entropies (for large (cid:96) ) are then given by Eqs. (2.38) and (2.40),respectively, with the amplitudes given in Eq. (4.23).In Fig. 4, we compare the analytical predictions for the symmetry resolved entanglement entropieswith the numerical data, focusing on half filling k F = π/ . The agreement is overall rather unsatis-factory, the asymptotic curves look rather close to data for even (cid:96) , but too far from those at odd (cid:96) .Although the oscillations are very large, we would have expected the analytic curves to be in the middleof the data, as it happens for the charged moments. Actually, the reason of this disagreement is easilyidentified. Indeed, the asymptotic curve in the saddle point approximation is given by the Fouriertransform of Z n ( α ) where the integral has domain on the entire real axis, rather than the interval [ − π, π ] . The difference between the two is exponentially small in the variance b n , and hence is usuallyneglected. However, the variance only grows as log (cid:96) ; consequently the corrections are algebraic in (cid:96) and so rather visible. The data at finite, even large, (cid:96) should be better described by the integral on thedomain [ − π, π ] , i.e. Z n ( q ) (cid:39) Z n (0) (cid:90) π − π dα π e − iα ( q − ¯ q ) − b n α / = Z n (0) e − ∆ q bn Erf (cid:0) i ∆ q + πb n √ b n (cid:1) + Erf (cid:0) − i ∆ q + πb n √ b n (cid:1) √ πb n , (4.24)18here Erf( z ) is the error function. The symmetry resolved entanglement entropies obtained, via Eq.(2.7) are reported for ∆ q = 0 , in Fig. 4 as continuous lines and, as expected, they fall right in themiddle of the numerical data for even and odd (cid:96) . For extremely large (cid:96) , Eq. (4.24) clearly converges tothe asymptotic result in Eq. (2.38). It is instructive to check analytically this slow approach, also tomotivate the strong visible difference compared to the periodic case [5]. To this aim, let us focus on thecase ∆ q = 0 which has the smallest corrections. The relative difference between the truly asymptoticbehaviour and Eq. (4.24) is Z n (0) − Z asy n (0) Z asy n (0) (cid:39) Erfc (cid:16) π (cid:114) b n (cid:17) (cid:39) e − b n π √ π b n ∼ (cid:96) − b n √ log (cid:96) , (4.25)where Erfc( z ) is the complementary error function and in the last equality we only used the leadingbehaviour b n ∼ b log (cid:96) . It is clear that for OBC these corrections are much more severe than for theinfinite system because b = 1 / (instead of b = 1 with PBC). For example, for n = 1 , which is thebetter behaving case among those in Fig. 4, we have corrections going like (cid:96) − / which are extremelyslow. We have seen that in CFT the (charged and neutral) entanglement entropies of one interval in themiddle of a semi-infinite chain are related to the ones of two intervals in an infinite chain. It is rathernatural to wonder whether this correspondence has some akin for free fermions on the lattice. This isindeed the case and it is a consequence of a general relation between Toeplitz+Hankel matrices andblock Toeplitz ones that we are going to show (and that in a similar form was present in Ref. [64], butfor a slightly different class of matrices).Let us start from the correlation matrix of two disjoint intervals of equal length (cid:96) at distance d inan infinite system. We denote this (cid:96) × (cid:96) matrix by C (cid:96),d that is easily obtained from the correlationmatrix of a single interval of length (cid:96) + d given by Eq. (2.23) by erasing all rows and column between (cid:96) + 1 and (cid:96) + d . It has the block structure C (cid:96),d = (cid:18) F E d + (cid:96) ( E d + (cid:96) ) T F (cid:19) , with F i,j = f i − j , E ki,j = f i − j + k , (5.1)and f m = sin( k F m ) πm are the elements of the correlation matrix of the single interval (2.23), although whatfollows is true for arbitrary f m such that f m = f − m . The four blocks in the matrix C (cid:96),d are Toeplitzmatrices, but the total matrix is not. (Notice that by rearranging rows and columns it is possible towrite C (cid:96),d in such a way that it is a Toeplitz matrix in which the elementary block is a × matrix, astandard form in the literature, see e.g. [64]. However this not relevant for us.)Let us denote by C − the (cid:96) × (cid:96) correlation matrix of a block of length (cid:96) at distance (cid:96) from theboundary of a semi-infinite system in Eq. (4.2), which can be written as C − = F − H (cid:96) , with H (cid:96) i,j = f i + j +2 (cid:96) . (5.2)19et us also introduce the matrix C + = F + H (cid:96) , (5.3)which, incidentally, is the correlation matrix of a semi-infinite system with a Neumann boundarycondition on the first site i = 1 (see e.g. [65]). Let us now assume that (cid:126)v − is a (cid:96) -component vector,which is eigenvector of C − , i.e. C − (cid:126)v − = λ v − (cid:126)v − . Let us introduce the reversed of a vector (cid:126)v as (cid:126)v R = ( v (cid:96) , v (cid:96) − , . . . v ) . Then, the (cid:96) -component vector (cid:126)V − ≡ ( (cid:126)v R − , − (cid:126)v − ) is an eigenvector of C (cid:96), (cid:96) +12 with eigenvalue λ v − , as one straightforwardly shows, indeed (cid:18) F E (cid:96) +1+ (cid:96) ( E (cid:96) +1+ (cid:96) ) T F (cid:19) (cid:18) (cid:126)v R − (cid:126)v (cid:19) = (cid:18) F (cid:126)v R − E (cid:96) +1+ (cid:96) (cid:126)v ( E (cid:96) +1+ (cid:96) ) T (cid:126)v R − F (cid:126)v (cid:19) , (5.4)and the components of the top vector are (cid:96) (cid:88) j =1 ( f i − j v (cid:96) +1 − j − f i − j +2 (cid:96) +1+ (cid:96) v j ) = (cid:96) (cid:88) j =1 ( f i − j − f i + j (cid:48) +2 (cid:96) ) v (cid:96) +1 − j , (5.5)which is what we wanted to prove. The same goes on for the bottom vector in Eq. (5.4). Similarly,given an eigenvector (cid:126)v + of C + of eigenvalue λ v + , we have that the vector (cid:126)V + ≡ ( (cid:126)v R + , (cid:126)v + ) is eigenvectorof C (cid:96), (cid:96) +12 with eigenvalue λ v + . By construction, all vectors (cid:126)V ± are orthogonal and hence, since theyare (cid:96) , they form a basis. The main conclusion here is that the spectrum of C is the union of thespectra of C + and C − , a rather remarkable result that likely is known in the literature, but we failedto find it (see however [64]). Thus we can generically relate the spectra of block Toeplitz matrices withToeplitz+Hankel ones. We stress that these results are only valid in the fermionic basis and not in thespin one, where the structure is much more complicated because of the Jordan-Wigner string [51, 66].Let us now move to the main objects of interest here that are the moments Z n ( α ) , both chargedand neutral. Since they are just products of functionals of the eigenvalues of the respective correlationmatrices, cf. Eq. (2.19), we have that the moment of two intervals of length (cid:96) at distance (cid:96) + 1 in aninfinite chain is the product of the moments built from C + and C − . This is valid for arbitrary (cid:96) and (cid:96) . In the scaling limit (cid:96), (cid:96) → ∞ , with arbitrary ratio (cid:96)/(cid:96) , the moments built with C + and C − areexpected to become equal, because both boundary conditions correspond to the same boundary CFTwith g = 1 , see e.g. [65, 67]. Consequently, the scaling form of the moments of one interval at distance (cid:96) from the boundary (with both free and Neumann conditions) is the square root of the one of twointervals of the same length at distance (cid:96) + 1 in an infinite system. En route to the calculation of the symmetry resolved entropies for a block away from the boundary, wefirst discuss the case of two disjoint intervals in an infinite system, since we have just seen the two arerelated. In this case, the correlation matrix is given by Eq. (5.1). For the matrices with this structure,a generalisation of the Fisher-Hartwig formula has been conjectured in Ref. [58] for the evaluation ofthe total Rényi entropies. Here, we just state the final result of Ref. [58] and refer the reader interestedinto more details to the original reference. We write the final result for a generic piecewise constant20ymbol g ( θ ) with R discontinuities at θ = θ r with r = 1 , . . . N , although we are also interested in thecase R = 2 in Eq. (2.25). We introduce t r ≡ g ( θ ) with θ ∈ [ θ r − , θ r ) . We consider the subsystem A = [ u , v ] ∪ [ u , v ] of total length | A | = v − u + v − u . The final result of Ref. [58] is S n ( X ) = A n | A | + B n log ( v − u )( v − u )( v − u )( u − v )( u − u )( v − v ) + 2 C n + . . . , (6.1)where the coefficients depend only on the symbol g ( θ ) and are given by A n = 12 π (cid:90) π − π dθf n ( g ( θ )) , (6.2)with f n ( x ) ≡ f n ( x, in Eq. (2.20), B n = 2 R (cid:88) r =1 J n ( r, r ) , (6.3) C n = R (cid:88) r =1 I n ( r ) − (cid:88) ≤ r (cid:54) = r (cid:48) ≤ R log [2 − θ r − θ r (cid:48) )] J n ( r, r (cid:48) ) , (6.4)where we defined J n ( r, r (cid:48) ) = 12 π (cid:90) t r t r − dλ df n ( λ ) dλ ω r (cid:48) ( λ ) , (6.5) I n ( r ) = 12 πi (cid:90) t r t r − df n ( λ ) dλ log (cid:34) Γ (cid:0) − iω r ( λ ) (cid:1) Γ (cid:0) + iω r ( λ ) (cid:1) (cid:35) dλ, (6.6)with ω r ( λ ) = 12 π log (cid:12)(cid:12)(cid:12) λ − t r λ − t r − (cid:12)(cid:12)(cid:12) . (6.7)For the specific case of two intervals in the ground state with only two singularities at ± k F , one finallyhas A n = 0 , B n = 1 / n − /n ) , and C n = Υ n . This results agrees with the conformal field theoryprediction F n ( x ) = 1 [63].It is completely clear that for the charged entropy of two disjoint intervals, the final result is alwaysgiven by Eq. (6.1) with the minor replacement f n ( λ ) → f n ( λ, α ) . Putting the various pieces togetherwe arrive at log Z n ( α ) = A n,α | A | + B n,α log ( v − u )( v − u )( v − u )( u − v )( u − u )( v − v ) + 2 C n,α , (6.8)with A n,α = iαk F π , (6.9) B n,α = 16 (cid:18) n − n (cid:19) + 2 n (cid:16) α π (cid:17) , (6.10) C n,α = Υ n ( α ) . (6.11)21igure 5: Scaling behaviour of the charged Rényi entropies log Z n ( α ) for several values of α and n . Wefix (cid:96) /(cid:96) = 2 and plot as function of (cid:96) . The numerical data (symbols) are compared with the prediction(7.1). The agreement is satisfactory, despite the presence of strong oscillating corrections. In this section, we finally move to the symmetry resolved entropies of a block A of length (cid:96) placed atdistance (cid:96) from the boundary (set at x = 0 for simplicity). As proved in Sec. 5, the asymptotic scalingbehaviour of Z n ( α ) is just the square root of Eq. (6.8) (with u and v being the mirror images of u = (cid:96) and v = (cid:96) + 2 (cid:96) ). Hence we have log Z n ( α ) = iαk F π (cid:96) − (cid:18) (cid:18) n − n (cid:19) + 1 n (cid:16) α π (cid:17) (cid:19) log (cid:18) (cid:96) (cid:96) ( (cid:96) + (cid:96) )(2 (cid:96) + (cid:96) ) (2 sin k F ) (cid:19) +Υ( n, α )+ o (1) , (7.1)where Υ( n, α ) is the same of Eq. (2.32). The interpretation of the three terms is the same as in Eq.(4.7) for the block starting from the boundary. This result is valid in the limit (cid:96), (cid:96) (cid:29) with theirratio arbitrary and confirms the CFT scaling (3.10). The only new relevant information is that also for α (cid:54) = 0 , the function F n,α ( x ) in Eq. (3.10) is equal to , as could be likely derived by field theoreticalmeans as those in Ref. [63].In Fig. 5, we show the comparison between the numerical data and the analytic prediction (7.1)for the Rényi entropies in the presence of flux for different values of n and α , with fillings k F = π/ and k F = π/ . The data are plot as function of (cid:96) at fixed ratio (cid:96) /(cid:96) = 2 . As expected, the presence ofoscillating corrections (with amplitude that increases as n and α get larger) strongly affect the data,but the agreement is satisfactory.In order to investigate in a better way the correctness of Eq. (7.1) and to isolate the universalpart, it is standard practice to construct combinations of entropies of different regions that cancel allnon-universal amplitudes (for disjoint intervals the mutual information is one of them [36, 37]) and areexplicitly scale invariant, i.e. they depend only on the ratio (cid:96)/(cid:96) or, equivalently, on the anharmonicratio (1.7). To this aim, let us first define the shorthand ζ ( n,α ) A for Z n ( α ) of the subsystem A , because22igure 6: The scale invariant function log t n ( α ) (7.2) at fixed α = 1 plotted as function of theanharmonic ratio x for n = 1 (left), n = 2 (center) and n = 3 (right). The numerical data are obtainedat fixed anharmonic ratio x (1.7) for increasing values of (cid:96) (and (cid:96) ). The data are extrapolated to thescaling limit (cid:96), (cid:96) → ∞ fitting with the first two corrections to the scaling. The obtained results matchremarkably well our analytic prediction (7.4) (continuous line).now we need to refer explicitly to the subsystem. The combination of interest is t n ( α ) ≡ e iα (cid:104) Q [0 ,(cid:96) (cid:105) ζ ( n,α )[ (cid:96) ,(cid:96) + (cid:96) ] ζ ( n,α )[0 ,(cid:96) ] ζ ( n,α )[0 ,(cid:96) + (cid:96) ] , (7.2)where the prefactor is introduced to cancel the mean number of particles in the interval [0 , (cid:96) ] . Thescaling limit of t n ( α ) can be written in terms of twist field as t n ( α ) ≡ e iα (cid:104) Q [0 ,(cid:96) (cid:105) (cid:104)T n,α ( (cid:96) ) ˜ T n,α ( (cid:96) + (cid:96) ) (cid:105)(cid:104)T n,α ( (cid:96) ) (cid:105)(cid:104)T n,α ( (cid:96) + (cid:96) ) (cid:105) , (7.3)in which it is clear that we are taking the ratio of the twist field correlation with its connected part,making it dimensionless and so scale invariant. Notice in particular that for α = 0 , using also the factthat the entanglement of A is the same as the one of the complement, the above ratio corresponds tothe (exponential) of the Rényi mutual information.Combining Eqs. (7.3) and (7.1), it is straightforward to obtain the leading behaviour of t n ( α )log t n ( α ) = − (cid:20) (cid:18) n − n (cid:19) + 1 n (cid:16) α π (cid:17) (cid:21) log (cid:96) (2 (cid:96) + (cid:96) ) = − (cid:20) (cid:18) n − n (cid:19) + 1 n (cid:16) α π (cid:17) (cid:21) log x , (7.4)where in the rightmost hand side we recognised the appearance of the anharmonic ratio x as definedin Eq. (1.7). This result is tested against numerics in Fig. 6. It is clear that the date are plaguedby scaling corrections. In order to achieve the scaling limit we work as follows. For some values of n and α reported in the figure, we computed the ratio t n ( x ) at fixed x and for increasing values of (cid:96) from 10 to 200. The data present finite (cid:96) corrections, which become larger as n and α increase. Theleading corrections to the scaling behave (cid:96) − (1 − α/π ) /n , cf. Eq. (4.22). There are also corrections goingas (cid:96) − (1+ α/π ) /n , (cid:96) − (2 − α/π ) /n , and (cid:96) − and which one is the first subleading depends on n and α . Wethen perform a fit of the finite (cid:96) data, keeping the first two power-law corrections and extrapolating at (cid:96) → ∞ . The data obtained following this procedure are reported in Fig. 6 and they match extremelywell the analytic prediction (7.4). This is particularly remarkable for n = 2 and when the finite (cid:96) data are still very far from their asymptotic values.23igure 7: Symmetry resolved von Neumann (left) and second Rényi (right) entropies at half filling k F = π/ for one interval of length (cid:96) placed at distance (cid:96) from the boundary of a semi-infinite chain.We focus on (cid:96) /(cid:96) = 2 . The analytical predictions (dashed lines) in Eq. (2.40) with the amplitudes inEq. (7.5) are compared with the numerical data (symbols) for ∆ q = q − ¯ q = 0 , . The agreement of thesaddle point results with the numerical data is better than the one in Fig. 4. However, the observeddeviations are again caused by the finite domain of the integral in α , and are well captured by Eq.(4.24), full line, which falls in the middle between the data for even and odd (cid:96) . For the symmetry resolution of moments and entropies, once again, we exploit the stationary phaseapproximation which implies the Gaussian approximation (2.35) for Z n ( α ) . The parameters definingthe variance b n for large (cid:96) , see Eq. (2.36), in this case read b = 1 , h n = − π n (cid:20) ln(2 | sin k F | ) + 12 log 4 (cid:96) ( (cid:96) + (cid:96) )(2 (cid:96) + (cid:96) ) (cid:21) + 2 γ ( n ) . (7.5)Notice that for (cid:96) (cid:29) (cid:96) , h n crosses over to the result in the bulk (2.37), as it should. The symmetryresolved moments and entropies (for large (cid:96) ) are then given by Eqs. (2.38) and (2.40), respectively,with the amplitudes given in Eq. (7.5).In Fig. 7, we compare the analytical predictions for the symmetry resolved entanglement entropieswith the numerical data, focusing on half filling k F = π/ and on (cid:96) /(cid:96) = 2 . Let us critically comparethese results with those with (cid:96) = 0 in Fig. 4. Especially for n = 2 , we observe that the analyticcurves are not in the middle of the data for even and odd (cid:96) , but it is slightly shifted. The shift is muchsmaller than the one for (cid:96) = 0 and it can be described by performing the integral on [ − π, π ] as in Eq.(4.24). Indeed as (cid:96) → ∞ , the data crossover to those for an infinite system, where such an effect wasnegligible [5]. Indeed for the most relevant case n = 1 , the two curves are extremely close.We conclude this section by studying another universal quantity appearing in the symmetry resolu-tion. This is related to the variance of the generalised probability distribution function p n ( q ) ≡ Z n ( q ) Z n .Indeed, as pointed out in Ref. [8], the difference of the variance of a given state with a reference onecancels all non universal factors. In our case, we can take the difference of the variance (cid:104) ∆ q (cid:105) OBC n for24igure 8: Excess of variance D n ( x ) in Eq. (7.6) as function of x at half filling k F = π/ for n = 1 , , respectively obtained from Eq. (7.6) (blue lines) and from numerical evaluation via correlation matrix(symbols). The numerical data are at fixed ratio x for increasing values of (cid:96) (and (cid:96) ). The data areextrapolated to the scaling limit (cid:96), (cid:96) → ∞ fitting with the first two corrections to the scaling. Theobtained results match very well the analytic prediction (7.6) (continuous line).one interval at distance (cid:96) from the boundary with the one of the infinite system (cid:104) ∆ q (cid:105) PBC n . Since,the variance is just the quadratic term in Z n ( α ) , it is clear that the difference between Eqs. (7.5) and(2.37) cancels all non-universal and divergent factors, resulting finally in D n ( x ) ≡ (cid:104) ∆ q (cid:105) OBC n − (cid:104) ∆ q (cid:105) PBC n = 12 nπ log 4 (cid:96) ( (cid:96) + (cid:96) )( (cid:96) + (cid:96) ) = 12 nπ log(1 − x ) , (7.6)where in the rightmost hand side we again recognise the anharmonic ratio (1.7). We can compare theanalytic expression of D n ( x ) with the numerical results that can be obtained from [8] (cid:104) ∆ q (cid:105) n = (cid:88) k ν − k − n + 1) − ν − k − n + 1) , (7.7)where the ν k are again the eigenvalues of the correlation matrix corresponding to the state and bipar-tition of interest. The excess of variance D n ( x ) is then numerically evaluated as the difference betweenEq. (7.7) evaluated for open and periodic boundary conditions. The results are shown in Fig. 8 for n = 1 , , (respectively from left to right) as function of x at half filling. The row numerical data atfinite (cid:96), (cid:96) are far from the analytic prediction (full line) but tend toward them. To be more quantita-tive, we extrapolate to the scaling limit with a fit that takes into account the first two corrections tothe scaling that are expected to go like (cid:96) − /n and (cid:96) − /n respectively. These extrapolations are shownin Fig. 8 and the agreement is very good, especially when thinking that the row data are very far awayfrom the asymptotic ones for n = 2 , . In the present work, we characterised the symmetry resolved entanglement entropies in the presenceof boundaries. We first reported some general results in CFT and then we specialised to spinless freefermions hopping on a 1D lattice. We focused on the case of a single block which can be either placedat the boundary or away from it. The former case can be treated with rigorous methods based onFisher-Hartwig formulas [35, 61] that allow us to estimate also the corrections to the leading behaviour.Conversely, for the interval away from the boundary, we first derived an exact relation with the case25f two intervals on the infinite line and then exploited it to use a recently proposed conjecture for thelatter case [58]. All our analytic results are tested against exact lattice computations. We found thatthe saddle-point approximation from charged to symmetry resolved entropies introduces algebraicallydecaying corrections in (cid:96) that are much more severe than in the periodic case. In all considered instances,we have entanglement equipartition, as expected; however we also identify the first term breaking suchan equipartition.We now briefly mention some possible generalisations and extensions of our results. First, thanksto the correspondence between the correlation matrix of the 1D lattice and the overlap matrix of a freeFermi gas [65, 68], our results immediately generalise to the latter case, although we did not discussit here. Then, we can exploit our results, to infer predictions for higher dimensional free fermioniclattices and gases with boundaries, using dimensional reduction techniques [13]. Less straightforwardgeneralisations concern instead the effect of the boundaries in the symmetry resolved massive fieldtheory. For integrable field theories, one should join the boundary approach of Ref. [69] with the formfactors of the composite twist fields [17]. Instead for free massive theories, we should generalise thetechniques of Refs. [9, 63] to the presence of boundaries. Acknowledgements We are extremely grateful to Sara Murciano for very useful discussions and for comparison of someof the reported results. We also thank Luca Capizzi and Paola Ruggiero for useful discussions andcollaboration on related topics. Both authors acknowledge support from ERC under Consolidatorgrant number 771536 (NEMO). References [1] N. Laflorencie and S. Rachel, Spin-resolved entanglement spectroscopy of critical spin chains andLuttinger liquids , J. Stat. Mech. P11013 (2014).[2] M. Goldstein and E. Sela, Symmetry Resolved Entanglement in Many-Body Systems , Phys. Rev.Lett. , 200602 (2018).[3] J. C. Xavier, F. C. Alcaraz, and G. Sierra, Equipartition of the entanglement entropy , Phys. Rev.B , 041106 (2018).[4] E. Cornfeld, M. Goldstein, and E. Sela, Imbalance Entanglement: Symmetry Decomposition ofNegativity , Phys. Rev. A , 032302 (2018).[5] R. Bonsignori, P. Ruggiero, and P. Calabrese, Symmetry resolved entanglement in free fermionicsystems , J. Phys. A , 475302 (2019).[6] S. Fraenkel and M. Goldstein, Symmetry resolved entanglement: Exact results in 1d and beyond ,J. Stat. Mech. 033106 (2020). 267] N. Feldman and M. Goldstein, Dynamics of Charge-Resolved Entanglement after a Local Quench ,Phys. Rev. B , 235146 (2019).[8] L. Capizzi, P. Ruggiero, and P. Calabrese, Symmetry resolved entanglement entropy of excitedstates in a CFT , arXiv:2003.04670v2.[9] S. Murciano, G. Di Giulio, and P. Calabrese, Entanglement and symmetry resolution in twodimensional free quantum field theories , JHEP 2008 (2020) 073.[10] S. Murciano, G. Di Giulio, and P. Calabrese, Symmetry resolved entanglement in gapped integrablesystems: a corner transfer matrix approach , SciPost Phys. , 046 (2020).[11] P. Calabrese, M. Collura, G. Di Giulio, and S. Murciano, Full counting statistics in the gappedXXZ spin chain , EPL , 60007 (2020).[12] M. T. Tan and S. Ryu, Particle Number Fluctuations, Rényi and Symmetry-resolved EntanglementEntropy in Two-dimensional Fermi Gas from Multi-dimensional bosonisation , Phys. Rev. B ,235169 (2020).[13] S. Murciano, P. Ruggiero, and P. Calabrese, Symmetry resolved entanglement in two-dimensionalsystems via dimensional reduction , J. Stat. Mech. (2020) 083102.[14] X. Turkeshi, P. Ruggiero, V. Alba, and P. Calabrese, Entanglement equipartition in critical ran-dom spin chains , Phys. Rev. B , 014455 (2020).[15] K. Monkman and J. Sirker, Operational Entanglement of Symmetry-Protected Topological EdgeStates , arXiv:2005.13026.[16] E. Cornfeld, L. A. Landau, K. Shtengel, and E. Sela, Entanglement spectroscopy of non-Abeliananyons: Reading off quantum dimensions of individual anyons , Phys. Rev. B , 115429 (2019).[17] D. X. Horváth and P. Calabrese, Symmetry resolved entanglement in integrable field theories viaform factor bootstrap , arXiv:2008.08553.[18] D. Azses and E. Sela, Symmetry resolved entanglement in symmetry protected topological phases ,arXiv:2008.09332.[19] J. Dubail, J.-M. Stéphan, J. Viti, and P. Calabrese, Conformal field theory for inhomogeneousone-dimensional quantum systems: the example of non-interacting Fermi gases , SciPost Phys. ,002 (2017).[20] P. Calabrese and J. Cardy, Evolution of entanglement entropy in one-dimensional systems , J.Stat. Mech. (2005) P04010.[21] P. Calabrese and A. Lefevre, Entanglement spectrum in one-dimensional systems , Phys. Rev. A , 032329 (2008); 27. Alba, P. Calabrese, and E. Tonni, Entanglement spectrum degeneracy and Cardy formula in1+1 dimensional conformal field theories , J. Phys. A , 024001 (2018).[22] P. Calabrese and J. Cardy, Entanglement entropy and quantum field theory , J. Stat. Mech. P06002(2004).[23] P. Calabrese and J. Cardy, Entanglement entropy and conformal field theory , J. Phys. A ,504005 (2009).[24] C. Holzhey, F. Larsen, and F. Wilczek, Geometric and renormalized entropy in conformal fieldtheory , Nucl. Phys. B , 443 (1994).[25] G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Entanglement in quantum critical phenomena ,Phys. Rev. Lett. , 227902 (2003);J. I. Latorre, E. Rico, and G. Vidal, Ground state entanglement in quantum spin chains , Quant.Inf. Comp. , 048 (2004).[26] H.-Q. Zhou, T. Barthel, J. O. Fjaerestad, and U. Schollwoeck, Entanglement and boundary criticalphenomena , Phys. Rev. A , 050305 (2006).[27] I. Affleck and A.W.W. Ludwig, Universal non- integer ground state degeneracy in critical quantumsystems , Phys. Rev. Lett. , 161 (1991).[28] E. Cornfeld and E. Sela, Entanglement entropy and boundary renormalization group flow: Exactresults in the Ising universality class , Phys. Rev. B , 075153 (2017).[29] N. Laflorencie, E. S. Sorensen, M.-S. Chang, and I. Affleck, Boundary effects in the critical scalingof entanglement entropy in 1D systems , Phys. Rev. Lett. , 100603 (2006).[30] P. Calabrese, M. Campostrini, F. Essler, and B. Nienhuis, Parity effect in the scaling of blockentanglement in gapless spin chains , Phys. Rev. Lett , 095701 (2010).[31] P. Calabrese and F. H. L. Essler, Universal corrections to scaling for block entanglement in spin-1/2 XX chains , J. Stat. Mech. (2010) P08029.[32] K. Ohmori and Y. Tachikawa, Physics at the entangling surface , J. Stat. Mech P04010 (2015).[33] P. Calabrese, J. Cardy, and I. Peschel, Corrections to scaling for block entanglement in massivespin-chains , J. Stat. Mech. (2010) P09003.[34] J. Cardy and P. Calabrese, Unusual Corrections to Scaling in Entanglement Entropy , J. Stat.Mech. (2010) P04023.[35] M. Fagotti and P. Calabrese, Universal parity effects in the entanglement entropy of XX chainswith open boundary conditions , J. Stat. Mech. P01017 (2011).2836] S. Furukawa, V. Pasquier, and J. Shirashi, Mutual information and compactification radius in ac=1 critical phase in one dimension , Phys. Rev. Lett. , 170602 (2009).[37] P. Calabrese, J. Cardy, and E. Tonni, Entanglement entropy of two disjoint intervals in conformalfield theory , J. Stat. Mech. P11001 (2009);P. Calabrese, J. Cardy, and E. Tonni, Entanglement entropy of two disjoint intervals in conformalfield theory II , J. Stat. Mech. (2011) P01021.[38] V. Alba, L. Tagliacozzo, and P. Calabrese, Entanglement entropy of two disjoint intervals in c=1theories , J. Stat. Mech. P06012 (2011).[39] A. Bastianello, Rényi entanglement entropies for the compactified massless boson with open bound-ary conditions , JHEP (2019) 141.[40] A. Bastianello, J. Dubail, J.-M. Stephan, Entanglement entropies of inhomogeneous Luttingerliquids , J. Phys. A Quantum computation and quantum information . CambridgeUniversity Press, Cambridge, UK, 10th anniversary ed. (2010).[42] A. Lukin, M. Rispoli, R. Schittko, M. E. Tai, A. M. Kaufman, S. Choi, V. Khemani, J. Leonard,and M. Greiner, Probing entanglement in a many-body localized system , Science , 6437 (2019).[43] H. M. Wiseman and J. A. Vaccaro, Entanglement of Indistinguishable Particles Shared betweenTwo Parties , Phys. Rev. Lett. , 097902 (2003).[44] H. Barghathi, C. M. Herdman, and A. Del Maestro, Rényi Generalization of the Accessible En-tanglement Entropy , Phys. Rev. Lett. , 150501 (2018); H. Barghathi, E. Casiano-Diaz, and A.Del Maestro, Operationally accessible entanglement of one dimensional spinless fermions , Phys.Rev. A , 022324 (2019).[45] M. Kiefer-Emmanouilidis, R. Unanyan, J. Sirker, and M. Fleischhauer, Bounds on the entangle-ment entropy by the number entropy in non-interacting fermionic systems , SciPost Phys. , 083(2020);M. Kiefer-Emmanouilidis, R. Unanyan, J. Sirker, and M. Fleischhauer, Evidence for unboundedgrowth of the number entropy in many-body localized phases , Phys. Rev. Lett. , 243601 (2020).[46] M. C. Chung and I. Peschel, Density-matrix spectra of solvable fermionic systems , Phys. Rev. B64 , 064412 (2001).[47] I. Peschel, Calculation of reduced density matrices from correlation functions , J. Phys. A , L205(2003).[48] I. Peschel and V. Eisler, Reduced density matrices and entanglement entropy in free lattice models ,J. Phys. A , 504003 (2009);I. Peschel, Entanglement in solvable many-particle models , Braz. J. Phys. , 267 (2012).2949] V. Alba, L. Tagliacozzo, and P. Calabrese, Entanglement entropy of two disjoint blocks in criticalIsing models , Phys. Rev. B , 060411 (2010).[50] F. Igloi and I. Peschel, On reduced density matrices for disjoint subsystems , EPL , 40001 (2010).[51] M. Fagotti and P. Calabrese, Entanglement entropy of two disjoint blocks in XY chains , J. Stat.Mech. (2010) P04016.[52] E. L. Basor and C. A. Tracy, The Fisher-Hartwig conjecture and generalizations , Physica A ,167 (1991);E. L. Basor and K. E. Morrison, The Fisher-Hartwig conjecture and Toeplitz eigenvalues , LinearAlgebra and Its Applications , 129 (1994).[53] B.-Q. Jin and V. E. Korepin, Quantum spin chain, Toeplitz determinants and Fisher-Hartwigconjecture , J. Stat. Phys. , 79 (2004).[54] J. P. Keating and F. Mezzadri, Random Matrix Theory and Entanglement in Quantum SpinChains , Commun. Math. Phys. , 543 (2004);J. P. Keating and F. Mezzadri, Entanglement in Quantum Spin Chains, Symmetry Classes ofRandom Matrices, and Conformal Field Theory , Phys. Rev. Lett. , 050501 (2005).[55] V. Alba, M. Fagotti, and P. Calabrese, Entanglement entropy of excited states , J. Stat. Mech.(2009) P10020.[56] F. Franchini, A. R. Its, and V. E. Korepin, Renyi Entropy of the XY Spin Chain , J. Phys. A ,025302 (2008);A. R. Its, B. -Q. Jin, and V. E. Korepin, Entanglement in XY Spin Chain , J. Phys. A , 2975(2005).[57] F. Ares, J. G. Esteve, F. Falceto, and E. Sanchez-Burillo, Excited state entanglement in homoge-neous fermionic chains , J. Phys. A , 245301 (2014).[58] F. Ares, J. G. Esteve, and F. Falceto, Entanglement of several blocks in fermionic chains , Phys.Rev. A , 062321 (2014).[59] F. Ares, J. G. Esteve, F. Falceto, and A. R. de Queiroz, On the Mobius transformation in theentanglement entropy of fermionic chains , J. Stat. Mech. (2016) 043106;F. Ares, J. G. Esteve, F. Falceto, and A. R. de Queiroz, Entanglement entropy and Möbiustransformations for critical fermionic chains , J. Stat. Mech. (2017) 063104.[60] J. L. Cardy, Conformal Invariance and Surface Critical Behaviour , Nucl. Phys. B , 514 (1984).[61] P. Deift, A. Its, and I. Krasovsky, Asymptotics of Toeplitz, Hankel, and Toeplitz+Hankel deter-minants with Fisher-Hartwig singularities , Ann. Math. , 1243 (2011).3062] E. Basor and T. Ehrhardt, Asymptotic Formulas for Determinants of a Special Class of Toeplitz +Hankel Matrices , In: D. Bini, T. Ehrhardt , A. Karlovich, I. Spitkovsky I. (eds) “Large TruncatedToeplitz Matrices, Toeplitz Operators, and Related Topics. Operator Theory: Advances andApplications”, 259 (2017).[63] H. Casini, C. D. Fosco, and M. Huerta, Entanglement and alpha entropies for a massive Diracfield in two dimensions , J. Stat. Mech. (2005) P07007.[64] M. Fagotti and P. Calabrese, Evolution of entanglement entropy following a quantum quench:Analytic results for the XY chain in a transverse magnetic field , Phys. Rev. A , 010306 (2008).[65] P. Calabrese, M. Mintchev, and E. Vicari, The entanglement entropy of 1D systems in continuousand homogenous space , J. Stat. Mech. P09028 (2011).[66] A. Coser, E. Tonni, and P. Calabrese, Spin structures and entanglement of two disjoint intervalsin conformal field theories , J. Stat. Mech. (2016) 053109.[67] P. Fendley, H. Saleur and N. Warner, Exact solution of a massless scalar field with a relevantboundary interaction , Nucl. Phys. B , 577 (1994).[68] P. Calabrese, M. Mintchev, and E. Vicari, The entanglement entropy of one-dimensional gases ,Phys. Rev. Lett. , 020601 (2011);P. Calabrese, M. Mintchev, and E. Vicari, Entanglement entropy of quantum wire junctions , J.Phys. A , 105206 (2012).[69] O. A. Castro-Alvaredo and B. Doyon, Bi-partite entanglement entropy in massive QFT with aboundary: the Ising model , J. Stat. Phys.134