More on symmetry resolved operator entanglement

The ‘operator entanglement’ of a quantum operator O is a useful indicator of its complexity, and, in one-dimension, of its approximability by matrix product operators. Here we focus on spin chains with a global U(1) conservation law, and on operators O with a well-defined U(1) charge, for which it is possible to resolve the operator entanglement of O according to the U(1) symmetry. We employ the notion of symmetry resolved operator entanglement (SROE) introduced in Rath et al (2023 PRX Quantum 4 010318) and extend the results of the latter paper in several directions. Using a combination of conformal field theory and of exact analytical and numerical calculations in critical free fermionic chains, we study the SROE of the thermal density matrix ρβ=e−βH and of charged local operators evolving in Heisenberg picture O=eitHOe−itH . Our main results are: i) the SROE of ρβ obeys the operator area law; ii) for free fermions, local operators in Heisenberg picture can have a SROE that grows logarithmically in time or saturates to a constant value; iii) there is equipartition of the entanglement among all the charge sectors except for a pair of fermionic creation and annihilation operators.


Introduction
Quantum entanglement plays a pivotal role to understand emergent phenomena in quantum many-body physics and numerical methods.In this context, the entanglement entropy has received significant attention and has become the most popular measure of bipartite entanglement in quantum systems.It is relevant in various contexts, ranging from high-energy physics [1][2][3][4] to condensed matter theory [5][6][7][8], when studying extended systems like quantum field theories (QFTs) and lattice models.For example, it can be useful to detect and describe phase transitions, even when a conventional order parameter is unavailable.In fact, its behaviour as a function of the subsystem size allows us to discern if the system is in a gapped or gapless phase and what are the universal features of critical systems.It turns out that in the former case, the entanglement follows the area law [9], i.e. it is proportional to the size of the border of the subsystem, in contrast to the thermal entropy which is characterised by a volume law.However, for the ground state of gapless local Hamiltonians in one dimension (1D), the area law is corrected by a logarithmic term [10][11][12][13][14].
One practical implication is that the presence of entanglement makes it difficult to simulate quantum many-body systems on a classical computer.For example, efficient simulations based on Matrix Product States (MPS) work well in non critical 1D systems [15][16][17][18][19], but they are less efficient when we approach a quantum critical point.The counterpart to MPS for approximating operators are Matrix Product Operators (MPOs), which are tensor network representations of operators, and it is natural to ask whether a quantity similar to the entanglement entropy exists to accurately capture the validity of this approximation [20][21][22][23][24][25][26][27].The answer lies in the concept of operator entanglement (OE), which quantifies the entanglement between quantum operators acting on different parts of a quantum system [20][21][22][23][28][29][30][31][32][33][34][35][36][37][38][39][40].The results about the OE depends on the specific operator and on the framework in which it is employed.In this paper we continue the analysis initiated in Ref. [41] of the symmetry resolved operator entanglement (SROE), which extends the notion of OE to consider the entanglement properties of operators with respect to specific symmetries of the quantum system.The interplay between the entanglement of a state and symmetries has been intensively studied in the last years through the symmetry-resolved entropies [42][43][44], both theoretically [45][46][47][48][49][50][51][52][53] and experimentally [41,[54][55][56], for several entanglement measures [57][58][59][60][61].The only symmetry resolution for operators studied so far is the reduced density matrix after a quantum quench from a product state [41] (see also for related quantities in open systems [27,31]).As in the non-resolved case, it exhibits an entanglement barrier [22,62,63]: it grows linearly in time as entanglement builds up between the local degrees of freedom, it then reaches a maximum, and ultimately decays to a small finite value as the reduced density matrix converges to a simple stationary state.
Before moving to the organisation of the paper, we introduce the main concepts: the OE, and how to symmetry resolve it.

Definition of Operator Entanglement
The entanglement properties between a region A and its complementary B are usually defined starting from a state |ψ⟩.Nevertheless, one can also study the entanglement properties of an operator O by vectorizing it.Namely, if the operator O lives in End(H A ) ⊗ End(H B ), we can view it equivalently as a vector in (H A ⊗ H B ) ⊗ ( HA ⊗ HB ), where H A,B denotes the Hilbert space in part A or B, and HA,B its dual.To work with a properly normalised state after the vectorization, we divide the operator O by a normalisation factor Tr(O † O).Operationally, one can pick an orthonormal basis {|i⟩} for (H A ⊗ H B ) and the corresponding basis {|j⟩} for its dual ( HA ⊗ HB ), and write O = ij O ij |i⟩ ⟨j|, where O ij = ⟨i| O |j⟩.Then the normalised operator-state is obtained as [64,65] Importantly, O admits a Schmidt decomposition, where r is the operator Schmidt rank, and the λ i are real positive coefficients that satisfy A,i , O A,j ] = δ ij .This can be seen by performing the ordinary Schmidt decomposition of the pure state |O⟩ and eventually reverting the vectorization to get back to the space of operators.
The OE is defined as follows.From the vectorization |O⟩, we can build the superreduced-density-matrix Tr B⊗B (|O⟩ ⟨O|) -where 'super' refers to the operators in the space of operators-which is a super-operator acting on operators on H A .Then the OE is the Rényi entropy of that super-reduced-density-matrix, Alternatively, the OE is given in terms of the Schmidt values λ i in Eq. ( 2) as As usual, the limit n → 1 produces the von Neumann OE

U (1) charge and symmetry resolution of OE
We now assume that there is a charge operator Q acting on the full Hilbert space H = H A ⊗H B which generates a U (1) symmetry, and which is a sum of the two charge operators acting on subsystems A and B, i.e.
A natural question is how to define a symmetry resolution of the OE of an operator O that possesses a fixed charge q O , i.e.
[Q, O] = q O Q.The answer is based on the symmetry resolution for |O⟩, as we are going to review.This problem has been already addressed in [26,27,41] and we report here the main definitions we will use in this manuscript.
From the commutation relation between O and Q, the terms in the Schmidt decomposition (2) can be reorganised according to their charge q as [41,61] O where A,j = q O (q) A,j ⊗ O . The charge q that appears in these equations can be introduced as the eigenvalue of a charge 'superoperator' Q living in the Hilbert space End(H) ⊗ End( H) Such superoperator satisfies the commutation relation where we have used the vectorization introduced in Eq. ( 1).Using the local structure of Q in A ∪ B, we can write and the following commutation relation holds We can exploit the result above such that Tr B⊗B (|O⟩ ⟨O|) can be decomposed as where Π q is a projector onto the eigenspace of Q A with fixed q and Tr B⊗B (|O⟩ ⟨O|)(q) denotes the (normalised) reduced density matrix built from the vector |O⟩ and restricted to the charge q.The partition function projected in a given charge sector reads and the SROE is given by According to Eq. ( 12), the total von Neumann OE associated to O splits into As it was shown in detail in [41], using the uniqueness of the Schmidt coefficients, the set of all (non-zero) values {λ (q) j } altogether must be the same as the set of values {λ i } from Eq. (2).Therefore, all quantities defined above can be defined in terms of {λ (q) j } as follows

Organization of the paper
Our goal is to extend our previous study on SROE [41].The manuscript is organised as follows.In Section 2 we review the known analytical techniques for computing the (nonresolved) OE in CFT and in free fermion spin chains, and we extend these to include the symmetry resolution.In Section 3 we apply these techniques to the case of the density matrix of a critical 1D system at thermal equilibrium.In Section 4 we study the SROE of local operators evolving in Heisenberg picture for free fermion chains.We draw our conclusions in Section 5.

Techniques to compute the SROE
Before reviewing the technical tools necessary to evaluate the SROE, both in field theory and in free fermionic lattice models, we summarise how to tackle the problem of symmetry resolution of a U (1) symmetric state following Refs.[42,47].

U (1) symmetry resolution
In this section, we explain how to compute the symmetry resolved entanglement entropy in a given charge sector for a state |ψ⟩, denoted by S (n) q (|ψ⟩).Let us consider a system with an internal U (1) symmetry and its bipartition into two subsystems, A and B. Tracing out the degrees of freedom of B, we obtain the reduced density matrix (RDM) of A, ρ A = Tr B |ψ⟩ ⟨ψ|.A measure of the entanglement between A and its complementary part is provided by the Rényi entropies, defined as and the limit n → 1 gives the von Neumann entropy as usual.When |ψ⟩ is an eigenstate of the hermitian charge operator Q, by taking the partial trace of the commutator This means that the reduced density matrix ρ A has a block-diagonal structure where each block corresponds to an eigenvalue where p A (q ′ ) is the probability of finding q ′ in a measurement of Q A in the RDM ρ A , i.e. p A (q ′ ) = Tr Π q ′ ρ A .Within this convention, the density matrices ρ A (q ′ ) of different blocks are normalised as trρ A (q ′ ) = 1.Thus, from the normalised ρ A (q ′ ), we can define the symmetry resolved Rényi entropies as The total von Neumann entanglement entropy associated to ρ A in Eq. ( 18) admits a decomposition as in Eq. (15), The two terms are known as 'configurational entanglement entropy' and 'fluctuation entanglement entropy' (or 'number entanglement entropy') respectively [54].The configurational entropy is also related to the operationally accessible entanglement entropy of Refs.[66][67][68], while the number entropy is the subject of a substantial recent activity [54,[69][70][71][72].
The calculation of the symmetry resolved entropies by the definition ( 19) is a difficult task, especially for an analytic derivation.As proposed in Ref. [42], it is convenient to focus on the charged moments of ρ A , Tr ρ n A e iQ A α [73][74][75][76].Their Fourier transforms with respect to α are the moments of the RDM restricted to the sector of fixed charge q ′ [42], i.e.
Finally the symmetry resolved entropies are obtained as We can apply the same machinery to compute the symmetry resolution of the OE in the charge sectors of the operator in Eq. ( 10) starting from the charged moments of the superdensity-matrix built from |O⟩.They are defined as and their Fourier transform gives where q is the eigenvalue of Q A .Thus, we have described a procedure to compute the SROE without the explicit knowledge of how the spectrum of the operator is resolved into the different charge sectors.
2.2 SROE from the replica trick in conformal field theory (CFT)

Brief review of the replica trick for the OE in CFT
In field theory, the reason for looking at Rényi entropies -instead of focusing directly on the von Neumann entropy-is that, for integer n, Trρ n A can be expressed in path-integral formalism as a partition function on an n-sheeted Riemann surface R n obtained by joining cyclically the n sheets along the region A [11,13].One approach to compute these Rényi entropies is based on a particular type of twist fields in quantum field theory that are associated with the branch points of the Riemann surface R n .We denote them by T n .Their action, in operator formalism, is defined by [14,77,78] where x 1 and x 2 are the endpoints of the interval A = [x 1 , x 2 ] on the real axis, x ′ is a point in A and i = 1, . . ., n indexes the n sheets modulo n.
denotes any operator acting in a single copy.By definition, the two-point function of the twist fields is the partition function on R n which enters the definition of the Rényi entropies [14].
In conformal invariant theories the two-point function is fixed by the scaling dimension of the fields, which here is We can generalise the action of the twist operators to arbitrary permutations σ ∈ S n of the permutation group of the n copies of the field theory, namely we can define a twist operator T σ as the operator with smallest possible scaling dimension such that [22] T σ (x 1 ) If σ has cycles of lengths n 1 + • • • + n c = n, then the scaling dimension of T σ is where c is the central charge of the theory.The Rényi entropy is given in terms of the correlator of the twist fields located at the two endpoints of A [14,22] S n (|ψ⟩ where Here σ has a single cycle of length n, so it has the scaling dimension (28).The OE of an operator O can be computed in a similar fashion.One has to consider Then the analogous expression of Eq. ( 29) reads [22] S We now generalise this formula to include the charge operators needed for the charged moments.

Method for computing the charged moments
The replica trick can be adapted as follows to obtain the charged moments Z n (α).The trick consists in inserting an Aharonov-Bohm flux through the multi-sheeted Riemann surface R n such that the total phase accumulated by the field upon going through the entire surface is α, similarly to what done for the standard entanglement entropy [42].In terms of twist fields, this amounts to computing where is the charge operator acting in the first replica.Let us stress that the operator e iαQ A,1 appears with opposite signs because of the definition of the charge operator in the doubled Hilbert space (see Eq. ( 10)).
In this paper, we will apply formula (31) in the context of a spinless Luttinger liquid, which is equivalent to a c = 1 compactified boson CFT parameterised by a coupling constant K called 'Luttinger parameter' see e.g.Refs.[79,80].Crucially, in the Luttinger liquid, the U (1) symmetry is generated by , where φ(x) is the boson field.Then when the charge operator is restricted to the interval A = [x 1 , x 2 ], the operators e iQ A α appearing in Eq. ( 31) can be written as which is a product of two vertex operators [42].In other words, when one turns around the operator T −α σ −1 (x 1 ) on the Riemann sheet, one goes from replica j to replica σ(j) and one also picks up a phase α when going through the first replica.The scaling dimension of the composite (or charged) twist fields reads [42] In Section 3 we will apply formulas ( 31)-( 32)- (33) to compute the SROE of the thermal density matrix.

Computing the SROE in free fermion chains
It is well established that, for the eigenstates of quadratic Hamiltonians, the entanglement entropies can be efficiently computed in terms of the eigenvalues of the correlation matrix of the subsystem [81][82][83].As pointed out in Ref. [41], this strategy can be adapted to compute the charged moments of the OE, see Eq. ( 23).Let us briefly review the resulting formalism here.The operator O that we are interested in is the Gaussian density matrix ρ of a free fermionic chain of length L. The density matrix can be diagonalised and be put in the form ρ ∝ e − k λ k c † k c k , where e −λ k = n k 1−n k with n k the occupation number of the orbital k = 1, . . ., L.Here c † k (c k ) creates (annihilates) a fermion in the orbital k; the creation/annihilation operators satisfy {c k , c † k ′ } = δ kk ′ .For our purposes it is convenient to write the density matrix ρ as where |0⟩ k (resp.|1⟩ k ) are states where the orbital k is empty (resp.filled), so that by applying the vectorization trick in Eq. ( 1) for ρ, we get |ρ⟩ where the ck operators are copies of the c k 's introduced in the vectorization process, and |0⟩ is the vacuum annihilated by all the c k 's and ck 's.From this pure state, we can build the superreduced-density-matrix Tr B⊗B (|ρ⟩ ⟨ρ|).The correlation matrix of the state |ρ⟩ reads [41,84] In the basis of c k , ck 's, the supercharge operator takes the form At this point, we can compute the 2L × 2L correlation matrix as and by doing a Fourier transform, we can write C in the spatial basis.To evaluate the charged moments in Eq. ( 23), we have to restrict the supercharge operator to Q A and do the Fourier transform of the correlation matrix in Eq. ( 38) to the subspace corresponding to the subsystem A, of size L A .Diagonalizing the latter matrix, we get 2L A real eigenvalues ξ i between 0 and 1.Therefore, one can compute the charged moments of the reduced density matrix built from |ρ⟩ in terms of the eigenvalues ξ i as [41,83] Using Eq. ( 24), we can compute exactly the SROE for the reduced density matrix of a free fermionic chain.We will use similar techniques to compute the SROE of an operator in Heisenberg picture in Section 4.
< l a t e x i t s h a 1 _ b a s e 6 4 = " A P O P j R 0 8 y U G j x S d s b e G a y S 8 j w e a X s 7 Z W r p 5 X S w e G w n h m y R t b J F v H I P j k g x + S E N A g j 9 + S R P J F n 5 8 F 5 c V 6 d t 8 H q m D P 0 r J I f c D 6 + A C n M p W Q = < / l a t e x i t > T ↵ < l a t e x i t s h a 1 _ b a s e 6 4 = " h P 1 w q e a X s 7 Z W r p 5 X S w e G w n h m y R t b J F v H I P j k g x + S E N A g j 9 + S R P J F n 5 8 F 5 c V 6 d t 8 H q m D P 0 r J I f c D 6 + A C n M p W Q = < / l a t e x i t > T ↵ < l a t e x i t s h a 1 _ b a s e 6 4 = " Q q T n 9 n 9 y V S 9 Z 5 q X J T L t Q v 5 / V k y S E 5 I i f E I h e k T q 5 J g z Q J I z F 5 J E / k 2 X g w X o x X 4 2 0 2 m j H m m X 3 y A 8 b 7 F 9 o Q m 4 s = < / l a t e x i t > i v < l a t e x i t s h a 1 _ b a s e 6 4 = " w P g V b J j n The replicated surface that is used to calculate the charged moments in the setup of [42] on a cylinder of circumference βv.The two twist fields T −α σ −1 and T α σ are located at the endpoints of the interval A. Right panel: For the charged moments of the OE, there are four twist fields acting at the endpoint of the interval in H A ⊗ HA .

SROE of a thermal density matrix
In this section we focus on the thermal density matrix, ρ β = e −βH , where β is the inverse temperature.The density matrix is vectorised, ρ β → |ρ β ⟩ ∈ H ⊗ H.For concreteness, consider for instance the free fermion chain (equivalent to the spin-1/2 XX spin chain via a Jordan-Wigner transformation) Then the charge operator Q can be promoted to a 'charge super-operator' Q = Q⊗1−1⊗Q T .Notice that Q |ρ β ⟩ = 0.In this section we compute the SROE of ρ β , relying on the fact that, at low energy, the Hamiltonian (40) corresponds to a free fermionic CFT with c = 1, or equivalently a Luttinger liquid with Luttinger parameter K = 1.We first evaluate the charged moments for the c = 1 free fermion CFT (or Luttinger liquid with K = 1), then extend the results to the case of interacting fermions (i.e. to Luttinger liquid with K ̸ = 1), and use the results to compute the SROE.We then benchmark our analytical results against numerics.

Charged moments for free fermion CFT
As reviewed in the previous section, in CFT the charged moments can be computed as where we have used the charged twist fields T α σ (x 1 ) that cyclically permute the replicas of the subsystem A, leaving B untouched as in Eq. ( 30), and are simultaneously the end points of the charge operator x 1 dx∂φ(x) as in Eq. ( 32).If we view the ratio (41) as a correlation function of twist operators living on an infinitely long cylinder of circumference 2βv, parametrised by the complex coordinate x + iy with (x, y) ∈ R × [0, 2βv], then the four charged twist operators are located at the points 0, L A , iβv and L A + iβv, where L A is the size of the subsystem A (see Fig. 1).
Therefore, the four-point function we need to calculate reads where σ = (1, 2, . . ., n) is a cyclic permutation.The scaling dimension of the composite twist field is given by Eq. ( 33) with K = 1.For free fermions, this object corresponds to the evaluation of the four-point function of charged twist fields on the cylinder of circumference 2βv.The result is given in Ref. [22] for α = 0, and it can be generalised straightforwardly to α ̸ = 0 by replacing the scaling dimension ( 26) by (33), leading to with c n,α a non-universal constant which depends on the microscopic details of the model, and that we will discuss in more detail below.We implicitly assume that the lattice spacing is set equal to 1.We use this result to compute the SROE in Section 3.3.

Generalization to K ̸ = 1
The result (43) for the charged moments can be generalised to a Luttinger liquid -or equivalently a compactified boson CFT-with a Luttinger parameter K ̸ = 1, corresponding to interacting fermions.The scaling dimension of the composite twist field in this case changes according to Eq. ( 33).We are interested in computing The calculation of that four-point function is much more complicated for K ̸ = 1 than for K = 1, as it involves more data from the underlying CFT, including operator content and OPE coefficients (see e.g.[85][86][87][88]).Here we exploit the result of the calculation for two intervals on the infinite line from Ref. [89], and we use the conformal transformation w → z = βv/π log w that maps each sheet in the w-plane into an infinitely long cylinder of circumference 2βv, to obtain the result that we need.It is convenient to introduce the cross-ratio Then the four-point function of the charged twist fields splits into a product of the four-point function of the twist fields in the plane, times the four-point function of vertex operators in the n-sheeted Riemann surface R n [89], Here c n,α denote again ultra-violet non-universal constant.The factor is the partition function ( 41) at α = 0, while the product of vertex operators is readily evaluated [89], leading finally to The partition function Z n (0) defined in Eq. ( 41) with α = 0 can be computed by combining the expression for the OE in Ref. [22] with the result for the four-point function of twist fields in Ref. [86].The result is where F n (x) is the conformal block of twist fields [86], Here Θ(u|Ω) is the Riemann-Siegel Theta function defined for an (n − 1) × (n − 1) complex matrix Ω and an (n − 1)-dimensional vector u: The matrix Γ(x) in Eq. ( 48) has entries given by [86] Γ and with Our final result for the charged moment is then In the two limiting cases L A ≫ βv and L A ≪ βv, where the four-point function factorises into a product of two two-point functions, we see that the result reduces to the expected one, These results for the charged moments have the same dependence on the parameters βv, L A as the one of the OE found in [22], with the scaling dimension ∆ Tσ replaced by ∆ T α σ .Moreover, when K = 1, F n (x) = 1 for all x and the logarithm of Eq. ( 46) coincides with Eq. (43).

Symmetry resolution
To go from the charged moment (52) to the SROE, we need to take the Fourier transform with respect to α.To do this, it is important to know the non universal constant c n,α .We start from the free fermion chain (40), where c n,α is known analytically, and later turn to the interacting case.57) for free fermions on the lattice with dispersion ε(k) = − cos k, for two different values of n.The Fermi velocity is v = 1.We take a chain of L sites at inverse temperature β and cut an interval of length L A .The black symbols represent the extrapolated data for fixed L A /β = 0.5.

SROE in the free fermion chain
The constant c n,α is known in the free fermion chain (40) from the charged moments of the reduced density matrix computed in Ref. [69].The ratio c n,α /c n,0 behaves quadratically at small α as [69] c n,α c n,0 where (This expression is real, as can be checked easily.)For later convenience, we define 2 hn .The constant c n,0 was computed by Korepin and Jin in Ref. [90]: This leads to the following result for the charged moments in the free fermion chain: We check the result (57) against numerics in Fig. 2. The symbols represent the numerical data obtained through exact lattice computations, using the expression for the charged moments described in Section 2.3.The data present finite size corrections, which become larger as n and α increase.In the right panel, in order to achieve the correct scaling limit, we have computed log Z n (α) at fixed α, n = 2 and for fixed L A /β = 0.5.The leading corrections to the scaling behave as (see also [91] for a similar analysis).We then perform a fit of the finite L A data (fox fixed α, n = 2), keeping the first two power-law corrections and extrapolating at L A → ∞.The data obtained following this procedure are reported as black symbols in the right panel of Fig. 2.
Plugging the expression (57) of the charged moments in the Fourier transform (24) and applying the saddle-point approximation for β, L A ≫ 1, we get where and This formula is also valid for the symmetry resolved von Neumann OE taking properly the limits of the various pieces when n → 1.The expression ( 58) is checked against numerics in the free fermion chain in Fig. 3.The small deviations between data and analytics are a consequence of the finite size corrections already present in the charged moments (Fig. 2).

Interacting case
For K ̸ = 1, we do not know a microscopic model where the non-universal constant c n,α is known.Nevertheless, we can generalise the previous results, assuming that c n,α behaves smoothly as a function of α near α = 0, i.e.
for some constants c n,0 and h n that depend on the microscopic model and are, in general, not equal to the ones in the free fermion chain.By repeating the same steps as before, the SROE reads

Discussion
To get insight into the meaning of Eqs. ( 58) and ( 62), we study the following two asymptotic regimes, which are valid for any value of K:  58) (solid lines) for β = 200, 20, n = 1 (upper panel), n = 2 (lower panel).For n = 2, we also check the result for the total OE (green line).We observe that as the temperature increases, the finite size corrections become more important.
We find that the leading order term corresponds to the total OE, which diverges logarithmically in L A at low temperatures, i.e. when the reduced density matrix ρ A is very close to the one of the ground state, while it is bounded (in L A ) at finite temperature β.This is the main striking difference with respect to the usual entanglement, which at finite temperature has an extensive behaviour.The practical consequence of this result is that ρ β can be efficiently represented by an MPO [22].Then S (n) q (ρ β ) presents a double logarithmic correction in L A at low temperature, while it remains bounded at finite temperature, with a double logarithmic correction depending on the inverse temperature β.Finally, the fact that the leading order term of the SROE coincides with the total OE resembles the entanglement equipartition for the symmetry resolution of a U (1)−invariant theory [43], and it can be traced back to the conformal invariance of the system also in this case.

SROE of local operators in Heisenberg picture in free fermion chains
The problem of computing the OE of local operators evolving in the Heisenberg picture, ϕ(t) = e iHt ϕe −iHt -also called 'local operator entanglement' by Kos, Bertini and Prosen [32]-, has attracted some attention recently [22,[30][31][32][33][34][35][36][37].In this section we study the symmetry resolution of the OE of a local operator in Heisenberg picture in the XX chain with periodic boundary conditions, which is mapped to the free fermion chain (40) by the Jordan-Wigner transformation.Following Refs.[20][21][22], we distinguish initial operators that are local in terms of the fermions (such as c † x or c † x c x ) from operators that are attached to a Jordan-Wigner string.

SROE of creation/annihilation operators and other local operators
The goal of this section is to study the SROE of c † x (t), σ z x (t).In the former case, we take advantage of the form in which Eq. ( 6) can be written, while in the latter we exploit the connection between the charged moments of the super-density-matrix and the correlation matrix, similarly to what has been done in Section 2.3.

SROE of creation operator c †
x (t) We start by studying the SROE of a single creation operator c † x (t) (the analysis for c x (t) can be performed in a similar way).We observe that [Q, c † x (t)] = c † x (t), with Q = j c † j c j , so q O = 1 in Eq. ( 6).We choose x = 0 and we can use the Fourier transform, for a ring of L sites, In the thermodynamic limit, the matrix elements U 0m (t) of the unitary evolution operator can be written in terms of Bessel functions, J m (t).We can directly recast c † 0 (t) into the form of Eq. ( 6), which reads This implies that there are only two possible charge sectors, q = 1, 0 and we can also directly write down the Schmidt values λ m (0) = λ m (1) = J m (t).At large times, the Bessel functions satisfy Combining Eq. ( 16) and ( 66), we get the following result at large time: and q p(q)S q (c † 0 (t) This implies that the SROE vanishes in all the charge sectors and the total OE is given only by the number/fluctuation entanglement (the second term of Eq. ( 15)) i.e. by the probability of finding q = 1, 0 as an outcome of Q A .

SROE of local density operator c
We now consider a less trivial example, which is the symmetry resolution of a pair of creation and annihilation operators.In the spin language (up to additive constants) this corresponds to c Let us focus on σ z 0 (t).From the state-operator correspondence, we want to study the entanglement due to the state c † 0 (t)c † 0 (t) |0⟩.Before reporting the details of the computations, we report the main result here.We observe that there are 3 non trivial charge sectors, q = 0, ±1, where and q p(q)S q (σ z 0 (t)) i.e.Eq. ( 15) is satisfied.Therefore, differently from what we found before for one single creation operator c † 0 (t), there exist one charge sector q = 0 where the SROE is different from 0 and both terms in Eq. ( 15) contribute to the total OE.
Given the simplicity of the result, it might be possible to find it following the logic of the previous subsection about c † (t).However, as an illustration of the techniques used in Section 2.3, we follow an alternative path.First of all, we need to compute Therefore, plugging Eq. ( 64) into Eq.( 72), we get By repeating the same computations for c † n c † m , cn c m , cn c † m , we can write down the 2L A × 2L A correlation matrix for a subsystem (−L A + 1, 0) as where C nm (t) = i m−n J m (t)J n (t), while 0 L A and 1 L A denote the L A × L A null and identity matrix, respectively.In the scaling limit t, L A → ∞, we can compute Tr(2C A (t) − 1) j , observing that This can be understood as follows: first of all, because of the construction of C A (t) with the blocks C(t) and 1 − C(t), Tr(2C A (t) − 1) j vanishes for any odd values of j.For the even powers, we can start by studying the case j = 2.We restrict to the upper L A × L A block and we compute the following trace (the trace of the square of the lower block is the same) Therefore, Eq. ( 76) and (66) give  The generalisation to higher powers of j follows the same logic, since all the powers involving terms like j J 2k j (t) vanish for k > 1 in the scaling limit.We can now come back to our main focus of the section, i.e the evaluation of the charged moments built from σ z 0 (t).The charged moments in Eq. ( 39) can be rewritten using [92,93] h The coefficients s n,α (m) correspond to the function h n,α (x) evaluated in certain simple points: It follows that In Fig. 4 we compare the prediction (80) for the logarithm of the charged moments log Z n (α) with the exact lattice computations done using the correlation matrix in Eq. ( 74).The dependence of Z n (α) in t and α is perfectly reproduced by the exact result.By doing the Fourier transform in Eq. ( 24), we can compute We recognise that Z q (σ z 0 (t)) = 0.This has been checked in the left panel of Fig. 5. Finally, we deduce that the SROE behaves  Left panel: probability of finding q as an outcome of the measurement of Q A .Right panel: SROE as function of t at fixed q.From Eq. ( 70) only the q = 0 charge sector has a non-vanishing entropy which saturates to the constant value log 2 .
as reported in Eq. ( 70): The fact that the OE is different from 0 only in one charge sector is a strong violation of the equipartition that we have found in Section 3.This lack of equipartition is a direct consequence of the fact that the operator entanglement stays always finite.

OE of a Jordan-Wigner string
We now turn to the calculation of the SROE of a Jordan-Wigner string in the tight-binding chain with periodic boundary conditions, i.e.
We place the endpoint of the string at the origin without loss of generality.Therefore, we are interested in the correlation matrix built from by replacing c(t) with c † (t) which anticommutes with all the c's.It is convenient to extract the time dependence as and to perform a Bogoliubov transformation (exactly as done for the total OE in [22]), Thus, the initial state becomes q (JW 0 (t)) with periodic boundary conditions.The solid line is Eq. ( 88), which also takes into account the exact knowledge of the non-universal constants for this system.The SROE of this operator is studied in the right panel: the comparison between numerics and the analytical prediction in Eq. ( 89) is quite good.
in Eq. ( 37), which already were diagonal.However, due to this simplification on the initial state, we can compute the correlation matrix as where C nm (t) = 0 y=−L A i m−n J m−y (t)J n−y (t) [94].Rewriting the initial state as in Eq. ( 86) also allows us to have an analytical insight on the time evolution of the SROE, following the observations done in [22].Indeed, Eq. ( 86 Let us remark that the knowledge of the non-universal constants from [69] allows us to benchmark our analytical prediction against the exact lattice computations without fitting any parameter, as we have done in the left panel of Fig. 6.By doing a Fourier transform and applying the saddle point approximation in the limit t → ∞, we get (89) where the constants h n , δ n , κ n have been defined in Eq. ( 59) and after Eq. ( 54), and This is just the double of the entanglement entropy following a quench from a domain wall state [95].Such a result is not surprising since it was already observed in [22] that the total OE of the JW string is exactly twice the entanglement entropy of a domain wall [96] and we find that the same is true in each symmetry sector.Hence, the equipartition of OE in the symmetry sectors is asymptotically restored with an asymptotic correction of the form q 2 / log 2 (t), with a non-trivial prefactor depending on non-universal quantities.The results of this and of the former subsection suggest that the breaking of the equipartition is a feature directly related to the finiteness of operator entanglement.The latter takes place for operators with a finite support (like σ z , i.e. c † x c x ) while extended operators, like the Jordan-Wigner strings, have diverging (in t) operator entanglement and restore equipartition for large times.

Conclusions
In the light of the examples studied so far, we can now draw our general conclusions on symmetry resolved operator entanglement.We have used the definition of the SROE introduced in [41], which quantifies the OE of a U (1) symmetric operator in a given charge sector and we have analysed in detail three operators.The thermal density has been studied with the twist field formalism and we found that it satisfies the operator area law, i.e. it has a bounded SROE in every charge sector and it displays equipartition, meaning that at leading order in the subsystem size it does not depend on the charge.For free fermion Hamiltonians, to evaluate the SROE of a local operator (in terms of the creation and annihilation operators) evolving in Heisenberg picture, we have exploited the knowledge of the SROE in terms of the correlation functions.We found that the local density operator obeys the area law at any time step and strongly violates the equipartition, while the SROE a Jordan-Wigner string grows logarithmically in time and obeys the equipartition of the entanglement.This remarkable difference between the three operators might be traced back to the fact that the SROE of local operators remains always finite, contrarily to the Jordan-Wigner string (that diverges for large t) or the thermal density matrix (that diverges for large β).The rationale appears to be that operators with a local spatial support have an operator entanglement that does not diverge and cannot satisfy equipartition.
We can now think of possible future directions that one could investigate about the OE and its eventual symmetry resolution.In a previous work about the same subject [41], the authors have studied the SROE of a reduced density matrix after a quantum quench [22,62,63], and it would be interesting to study how the entanglement barrier is affected by a non-unitary time evolution, due, for example, to the effect of local measurements in the dynamics.Another possible direction is studying the connection between the OE and the reflected entropy [97,98], by applying the techniques described in this manuscript to analyse the symmetry resolution also in this context (see also [40] for the first steps towards this direction).This example would also gain insights about the holographic dual of the SROE when O = √ ρ.Similarly one would like to explore the connection of our results with symmetry resolution of the computable cross norm negativity [40,99] which is another measure of entanglement in mixed states.A last interesting point concerns the interplay of OE and symmetries when the latter are broken; this problem can be studied generalising to the operatorial level the recently introduced entanglement asymmetry [100-103].

r i=1 λ 2 i = 1 .
The operators O A,i ∈ End(H A ) (same for O B,i ) obey the orthonormality condition Tr[O †

1 <
B G L k l 9 + S B P H p 3 3 p P 3 7 L 1 + r g 5 5 A 8 8 i + T H e + w d l b K e B < / l a t e x i t > T ↵ l a t e x i t s h a 1 _ b a s e 6 4 = " h P 1 w q

1 <
B G L k l 9 + S B P H p 3 3 p P 3 7 L 1 + r g 5 5 A 8 8 i + T H e + w d l b K e B < / l a t e x i t > T ↵ l a t e x i t s h a 1 _ b a s e 6 4 = " + Q / e s w C m b u r 7 h k v x 4 e R 2 t c v A

1 <
B G L k l 9 + S B P H p 3 3 p P 3 7 L 1 + r g 5 5 A 8 8 i + T H e + w d l b K e B < / l a t e x i t > T ↵ l a t e x i t s h a 1 _ b a s e 6 4 = " + Q / e s w C m b u r 7 h k v x 4 e R 2 t c v A e a X s 7 Z W r p 5 X S w e G w n h m y R t b J F v H I P j k g x + S E N A g j 9 + S R P J F n 5 8 F 5 c V 6 d t 8 H q m D P 0 r J I f c D 6 + A C n M p W Q = < / l a t e x i t > T ↵ < l a t e x i t s h a 1 _ b a s e 6 4 = " + Q / e s w C m b u r 7 h k v x 4 e R 2 t c v A

Figure 2 :
Figure 2: Numerical check of formula (57) for free fermions on the lattice with dispersion ε(k) = − cos k, for two different values of n.The Fermi velocity is v = 1.We take a chain of L sites at inverse temperature β and cut an interval of length L A .The black symbols represent the extrapolated data for fixed L A /β = 0.5.

Figure 4 :
Figure 4: Logarithm of the charged moments for the operator σ z 0 (t) in the tightbinding model (see Eq. (40)) for a bipartition [−L A , 0] ∪ [1, L A ], L A = 30, with periodic boundary conditions.The plots are at fixed α as function of time (left) and at fixed time as function of α (right).The symbols are the numerical data coming from the evaluation of the eigenvalues of the correlation matrix in Eq. (74) while the solid line represents Eq. (80).

Figure 5 :
Figure5: Symmetry resolution of the OE for the operator σ z 0 (t).Left panel: probability of finding q as an outcome of the measurement of Q A .Right panel: SROE as function of t at fixed q.From Eq. (70) only the q = 0 charge sector has a non-vanishing entropy which saturates to the constant value log 2 .
where now |0 d,b ⟩ is the vacuum for the d and b modes.This transformation does not do anything to the Hamiltonian H ⊗ 1 − 1 ⊗ H and to the charge operator Q ⊗ 1 − 1 ⊗ Q T

Figure 6 :
Figure 6: Left panel: logarithm of the charged moments for the Jordan-Wigner string JW 0 (t) in the tight-binding model for a bipartition [−L A + 1, 0] ∪ [1, L A ], L A = 30,with periodic boundary conditions.The solid line is Eq.(88), which also takes into account the exact knowledge of the non-universal constants for this system.The SROE of this operator is studied in the right panel: the comparison between numerics and the analytical prediction in Eq. (89) is quite good.