Non-equilibrium entanglement asymmetry for discrete groups: the example of the XY spin chain

Entanglement asymmetry is a novel quantity that, using entanglement methods, measures how much a symmetry is broken in a part of an extended quantum system. So far, it has only been used to characterise the breaking of continuous Abelian symmetries. In this paper, we extend the concept to cyclic ZN groups. As an application, we consider the XY spin chain, in which the ground state spontaneously breaks the Z2 spin parity symmetry in the ferromagnetic phase. We thoroughly investigate the non-equilibrium dynamics of this symmetry after a global quantum quench, generalising known results for the standard order parameter.


Introduction
In the last few years, the interplay between entanglement and symmetries has become the centre of an intense research activity, yielding up a more refined vision of the behaviour of many-body quantum systems [1][2][3][4][5][6] and creating a new framework to study entanglement [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21].For example, novel quantities, such as the symmetry-resolved entanglement entropy [7][8][9], have been conceived to analyse how entanglement distributes in the symmetry sectors of a theory.On the other hand, it has recently arisen the idea of using entanglement tools to study symmetry breaking.To this end, entanglement asymmetry has been proposed in Ref. [22] as a measure of how much a symmetry is broken in a part of an extended quantum system.So far, the entanglement asymmetry has been applied to examine the time evolution of an initially broken U (1) symmetry after a quench with a Hamiltonian that preserves it, both in free [22,23] and interacting integrable systems [24].Generally, it is expected that the symmetry is dynamically restored in a subsystem after the quench.However, unexpectedly, the entanglement asymmetry reveals that the symmetry can be restored earlier when it is initially more broken [22].This surprising phenomenon can be seen as a quantum version of the very counter-intuitive Mpemba effect: the more a system is initially out of equilibrium, the faster it relaxes.Entanglement asymmetry has been also employed in Ref. [23] to analyse the special case in which symmetry is not restored and the subsystem relaxes to a non-Abelian Generalised Gibbs ensemble.
The previous works restrict to the continuous U (1) group.The goal of the present manuscript is to extend the study of the entanglement asymmetry to the finite, discrete cyclic Z N groups.In particular, we consider the Z 2 spin flip parity symmetry of the XY spin chain which, in the thermodynamic limit, is spontaneously broken by the ground state in the ferromagnetic phase [25].Usually this symmetry breaking is detected using the one-site longitudinal magnetisation as order parameter.This order parameter presents the inconvenience that, despite a non-zero value indicating that the symmetry is broken, the converse is not always true [26].On the contrary, the entanglement asymmetry unambiguously determines if the symmetry is or not respected.Moreover, although for a single spin the order parameter can be used to measure symmetry breaking, for larger subsystems there are longer range correlations that can break the symmetry and are only taken into account by the entanglement asymmetry.Using it, we study the non-equilibrium time evolution of the Z 2 spin parity symmetry after a global quantum quench in the XY spin chain Hamiltonian, extending the study performed in Refs.[27][28][29] employing the order parameter.Global quenches in this system have been studied from many perspectives, see e.g.[30][31][32][33][34][35][36][37][38][39][40][41][42][43][44].We distinguish between chains with open and periodic boundary conditions, since the approach to calculate the entanglement asymmetry and the results obtained differ: while in the periodic case this Z 2 symmetry is always restored, the same does not always occur in the open chain due the presence of boundary modes.Moreover, in the periodic case, by comparing with the analytic results obtained in Refs.[27,28] for the order parameter in the thermodynamic limit, we conjecture an effective expression that describes the evolution of the entanglement asymmetry after the quench.
The paper is organised as follows: in Sec. 2, we define the entanglement asymmetry, we review its main known features for the U (1) group and we extend it to cyclic Z N groups, describing some of its principal properties in this case.In Sec. 3, we introduce the XY spin chain and we present the behaviour of the entanglement asymmetry of the spin parity in the ground state.Sec. 4 is devoted to calculate and discuss the entanglement asymmetry after a quantum quench of the XY spin chain with OBC, while in Sec. 5 we carry out a similar analysis but taking PBC.In Sec. 6, we draw our conclusions.The manuscript also includes several appendices that describe in detail those points of the main text that are more technical.

Entanglement asymmetry
Let us consider an extended quantum system that can be divided into two spatial parts A and B.
If we assume that the total system is in a pure state |Ψ⟩, then the state of one of the regions, for example A, is described by the reduced density matrix ρ A = Tr B (|Ψ⟩ ⟨Ψ|), which is generally a mixed state.We further take a local charge operator Q given by the sum of the charge in A and We assume that Q generates an Abelian group with elements of the form e iαQ where α is a real number.This includes the unitary group U (1) or the cyclic groups Z N .
It is well-known that, when |Ψ⟩ respects the symmetry associated to Q, i.e. it is an eigenstate of Q, the reduced density matrix of any subsystem A displays a block-diagonal structure in the eigenbasis of Q A , and the blocks correspond to the sectors of a given charge.In this case, the entanglement entropy of subsystem A can be decomposed into the contributions of each charge sector, known as symmetry-resolved entanglement entropies [7][8][9].
In the case of symmetry-breaking states, entanglement entropy can be used to quantify how much the symmetry generated by Q is broken in a subsystem.In fact, if the state |Ψ⟩ locally breaks the symmetry in the region A, then ρ A does not commute with Q A .This means that, in the eigenbasis of Q A , ρ A presents non-zero entries outside the sectors of fixed charged ρ q j , , here represented by * .From this state, it is always possible to define another density matrix ρ A,Q that respects the symmetry, i.e. [ρ A,Q , Q A ] = 0.If we take the projection of ρ A on the different charge sectors of Q A and we sum all them, we obtain a block diagonal matrix where Π q j is the projector on the subspace of eigenvalue q j .The density matrix ρ A,Q can be physically interpreted as the result of performing a non-selective measurement of Q A in ρ A .From ρ A and its projection ρ A,Q , we can define the entanglement asymmetry [22] in terms of their von Neumann entropies The entanglement asymmetry measures how much the symmetry generated by Q is broken in the susbsystem A. In fact, ∆S A satisfies two crucial properties: it is always non-negative, ∆S A ≥ 0, and it vanishes if and only if [ρ A , Q A ] = 0 [45].These properties are satisfied whether the total system is in a pure or mixed state.Therefore, the entanglement asymmetry can also be applied to the study of symmetry breaking in mixed states and, in particular, to a system at a certain temperature.
Since the direct calculation of the von Neumann entropy is normally a difficult task, the usual strategy [46,47] is to replace it in Eq. ( 1) by the Rényi entropies Then we define the Rényi entanglement asymmetry as Note that, if we take in this expression the limit n → 1, we recover Eq. ( 1).Moreover, Rényi entanglement asymmetries are also non-negative quantities, ∆S A ≥ 0, and they cancel if and only if the symmetry is respected in subsystem A, i.e. [ρ A , Q A ] = 0 [48].It is also important to remark that the Rényi entanglement asymmetries for integer n ≥ 2 are experimentally measurable in ion traps through randomised measurements [49][50][51][52][53].

Review of entanglement asymmetry for the U (1) group
In the literature, the entanglement asymmetry has only been analysed for the Abelian unitary group U (1).In this case, the paradigmatic example to understand how it works is the tilted ferromagnetic state, For θ ̸ = 0, this state breaks the U (1) symmetry corresponding to rotations around the z axis, which are generated by the transverse magnetisation Q = 1 2 j σ z j .It has been shown that, if the subsystem A is an interval of contiguous spins of length ℓ, the entanglement asymmetry behaves as [22] ∆S for large ℓ.To better understand this result, let us focus on the n = 2 Rényi asymmetry ∆S A .Using the Fourier representation of the projectors Π q j , the projected matrix ρ A,Q can be written as Since |θ⟩ is a separable state, ρ A = |θ, ℓ⟩ ⟨θ, ℓ|, and |θ, ℓ⟩ is the restriction of |θ⟩ to ℓ contiguous spins.We can see Eq. ( 7) as the state of A after applying a random element of the U (1) group to ρ A with a uniform probability density.Consequently, ∆S A can be interpreted as the quantum information loss about the subsystem A after such an operation.Despite the infinite number of elements in the U (1) group, this information loss is finite because two states |θ, ℓ⟩ and e iαQ A |θ, ℓ⟩ are not orthogonal for α ̸ = 0 and finite ℓ but have a non-zero overlap.In fact, the n = 2 Rényi entanglement asymmetry can be written as

Entanglement asymmetry for cyclic groups Z N
Let us now consider the finite cyclic group Z N of rotations around the z axis by angles multiple of 2π/N , which is a subgroup of the U (1) symmetry generated by the transverse magnetisation Q.The density matrix ρ A,Q N of this group is obtained by projecting on the different charge sectors of the transverse magnetisation Q modulo N .Using the Fourier representation of the projectors, it can be written as This matrix can be interpreted either as the state of subsystem A after a non-selective measurement of the transverse magnetisation Q modulo N or after randomly applying to ρ A one of the elements of the Z N group with equal probabilities.The second interpretation gives the intuition that the maximal information loss after such an operation is log N since N is the dimension of the group and one can always recover all the information on the state of A by learning which element of the group has been applied.In the left panel of Fig. 1, we study the entanglement asymmetry of the tilted ferromagnetic state |θ⟩ as a function of the subsystem length ℓ for different subgroups Z N of the U (1) group of rotations around the z axis seen before.For small subsystem size, the entanglement asymmetry presents the same behaviour as in the U (1) case.However, for large ℓ, it precisely saturates to log N .This saturation occurs at larger values of ℓ as the dimension of the subgroup increases, recovering the result of Eq. ( 6) for the U (1) group in the limit N → ∞.In particular, for N = 2, we have ∆S where are the functions for the Rényi-n entropies of one bit of information.Their limit n → 1 reads   5) with θ = π/4.We represent them as a function of the subsystem length.For large subsystems they tend to log N .The solid line is the asymptotic entanglement asymmetry for the U (1) group, cf.Eq. ( 6), which diverges logarithmically with ℓ.Right panel: the exact analytic expression (11) for the subgroup Z 2 as a function of the tilting angle θ of the ferromagnetic state and different subsystem sizes.
In the right panel of Fig. 1, we plot Eq. (11) (i.e. for N = 2) as a function of the tilting angle θ for n = 1 and different subsystem sizes.We can see that ∆S A grows monotonically with θ until θ = π/2, value at which the U (1) symmetry is maximally broken since all the spins point in the x direction.Observe that, as the subsystem size increases, the asymmetry saturates faster to its maximum value, log 2, when the tilted angle θ turns on.In Fig. 2, we study the Rényi entanglement asymmetry for the subgroups Z 2 (left panel) and Z 3 (right panel) as a function of the tilting angle θ for different Rényi indices n and a fixed subsystem length.As we can see, for the subgroup Z 2 , ∆S A is monotonic in θ for n → ∞.On the other hand, for Z 3 , it is non monotonic in θ when n → ∞, presenting the same peaked behaviour as in the U (1) case [22].This is expected since, as we show in the left panel of Fig. 1

, ∆S (n)
A for the subgroups Z N tends to the U (1) asymmetry in the limit N → ∞.
We can generally prove that, for Z N groups, the entanglement asymmetry ∆S A is bounded by the dimension of the group.This is actually a consequence of the following inequality for the von Neumann entropy.If we consider a set of density matrices {ρ j } and of positive coefficients {λ j } with j λ j = 1, then [54] Therefore, if in the definition (1) of the entanglement asymmetry we take into account that for a Z N group the projected density matrix ρ A,Q N can be written in the form of Eq. ( 10), and we apply the previous inequality identifying λ j = 1/N and ρ j = e iα j Q A ρ A e −iα j Q A , then we can easily conclude that ∆S A ≤ log N.
The equality in Eq. ( 14) is attained if and only if the states ρ j have support in orthogonal subspaces.This implies that the entanglement asymmetry for a Z N group saturates to its maximal  value, ∆S A = log N , if and only if all the states e iα j Q A ρ A e −iα j Q A have orthogonal support.Therefore, the Z N symmetry is maximally broken if the action of any element of the group, except the identity, maps ρ A to different orthogonal subspaces.In the tilted ferromagnetic state, this is equivalent to requiring that |θ, ℓ⟩ and e iα j Q A |θ, ℓ⟩ are orthogonal, similarly to the discussion done in Sec.2.1 for the U (1) case, but now the average ( 8) is over a finite number of states and, therefore, is bounded.In Appendix A, we prove that these properties are also satisfied by the Rényi entanglement asymmetry (4) for any index n.Interestingly, the equality in Eq. ( 14) corresponds to a particular instance of Holevo theorem [54], a result of quantum information theory that states that the maximum information that one can extract from a mixed state such as ρ A,Q N is given by the difference between S( j λ j ρ j ) and j λ j S(ρ j ), which in our case is precisely log N , as we already announced.

The XY spin chain
In the rest of the paper, we study the non-equilibrium dynamics of a Z 2 broken symmetry in the XY spin chain.This section is devoted to introducing this model and discuss some of its general features that will be useful afterwards.The Hamiltonian of the XY spin chain reads where is a boundary term closing the chain with Periodic Boundary Conditions (PBC).If this term is not present, then the chain has Open Boundary Conditions (OBC).Observe that the Hamiltonian contains nearest neighbour interactions between the x and y components of the spins.The parameter γ > 0 tunes the anisotropy between the couplings in these two directions.There is also an external magnetic field of intensity h > 0 in the z direction.

∆S
(n) )-parameter plane of the XY spin chain (16).For γ > 0, the system is critical (zero mass gap) along the line h = 1, which separates the ferromagnetic (h < 1) and the paramagnetic (h > 1) phases.The line γ = 1 corresponds to the quantum Ising chain.In the thermodynamic limit, the ground state spontaneously breaks Z 2 spin parity symmetry in the ferromagnetic phase and the corresponding entanglement asymmetry quickly saturates to log 2 in all the region when the subsystem is large enough.A special line is h 2 + γ 2 = 1, in which the ground state is the tilted ferromagnetic state |θ⟩ analysed in Secs.2.1 and 2.2.
The reason why we are interested in this system is that it exhibits a Z 2 symmetry corresponding to rotations of π rad around the z axis, a subgroup of the U (1) symmetry discussed in Sec.2.1 and in Refs.[22,23].It is associated to the parity operator which maps the spins σ x,y j to −σ x,y j leaving σ z j invariant.This symmetry is spontaneously broken by the ground state of Eq. ( 16) in the thermodynamic limit L → ∞ when h < 1 (ferromagnetic phase) while, for h > 1 (paramagnetic phase), the symmetry is respected.Both phases are separated by the critical line h = 1.
These two phases can be distinguished using the one-site longitudinal magnetisation ⟨σ x j ⟩ as order parameter.In fact, for h < 1, we have that ⟨σ x j ⟩ ̸ = 0 while, for h > 1, ⟨σ x j ⟩ = 0. Since a non-zero magnetisation in the x direction obviously breaks the Z 2 symmetry of π rad rotations about the z axis, then the condition ⟨σ x j ⟩ ̸ = 0 implies that such symmetry is broken in the ground state.
We can also analyse this spontaneous breaking of the Z 2 parity symmetry with the Rényi entanglement asymmetry (4).In Fig. 3, we represent the (h, γ)-parameter plane of the XY spin chain, indicating the large subsystem value of the asymmetry in each phase.We find that, in the ferromagnetic phase h < 1, ∆S A quickly saturates for any Rényi index n to its maximal value, log 2, as the length of the subsystem ℓ increases.On the contrary, in the paramagnetic region h > 1, the asymmetry vanishes, signaling that the ground state respects the symmetry.Therefore, for infinitely large subsystems ∆S (n) A acts somehow as a topological indicator of symmetry-breaking.It is also worth mentioning that the tilted ferromagnetic state |θ⟩, for which we obtained its exact Z 2 entanglement asymmetry in Sec.2.2, is the ground state of the XY spin chain along the circle γ 2 + h 2 = 1 [55,56].
In this work, we are interested in studying the time evolution of this Z 2 symmetry in a global quantum quench from the ground state of the Hamiltonian with couplings (h 0 , γ 0 ) in the ferromagnetic phase, i.e. h 0 < 1, to another XY Hamiltonian with different couplings (h f , γ f ).One has [H, P ] = 0 for any post-quench XY Hamiltonian H. Hence, after this quench, the state that describes a subsystem A is expected to relax to a Generalised Gibbs Ensemble (GGE) that respects the Z 2 symmetry [28,29,38].This symmetry restoration has actually been studied using the order parameter ⟨σ x j ⟩ in Refs.[27,28].It was found to exhibit an exponential decay in time, both for quenches to the ferromagnetic and paramagnetic phases.In our analysis, we distinguish two cases.In Sec. 4, we consider the entanglement asymmetry of an interval attached to one of the extremes of an open chain, situation that in the thermodynamic limit L → ∞ corresponds to a semi-infinite chain.In Sec. 5, we examine the entanglement asymmetry in the thermodynamic limit of a chain with periodic boundary conditions.In each case, we apply a different approach and we find that, in some situations, the evolution of the asymmetry after the quench differs.

Jordan-Wigner transformation
One of the most important features of the XY spin chain ( 16) is that it is free, i.e. it can be mapped to a Hamiltonian describing non-interacting spinless fermions [57].This can be done by means of the Jordan-Wigner transformation where c j and c † j are creation and annihilation fermionic operators, satisfying the anticommutation relations Under this mapping, the parity operator (18) becomes which corresponds to the parity of the number of fermions N .The Hamiltonian ( 16) is mapped to the Kitaev fermionic chain which is local in terms of the fermionic operators in spite of the non-locality of the Jordan-Wigner transformation.Observe that the transverse magnetic field h in Eq. ( 16) translates into a chemical potential and the anisotropy parameter γ controls the superconducting pairing terms c † j c † j+1 + h.c.. Importantly, these pairing terms break the fermion number symmetry, but not their parity since they create/destroy particles by pairs.
The boundary term H b reads after the Jordan-Wigner transformation as Observe that it includes the parity operator P and, therefore, depending on the parity sector, yields periodic or antiperiodic boundary conditions in the Kitaev chain.

Entanglement asymmetry in an open chain
For OBC, there is no boundary term H b in the XY spin chain Hamiltonian (16), which is a quadratic form in terms of the fermionic operators c j and c † j , as we have seen in Eq. (22).In order to simplify the calculations and to better understand some physical results later, it is very useful to introduce the Majorana fermionic operators [58] ǎ2j−1 = c † j + c j and ǎ2j = i(c † j − c j ), (24) by doubling the sites of the lattice.These operators are real valued, ǎ † m = ǎm , and satisfy the anticommutation relation Since the Hamiltonian ( 22) is quadratic, following Ref.[57], it can be diagonalised by performing a Bogoliubov rotation to a new set of fermionic operators η k and η † k .If we organise the Majorana ǎm and Bogoliubov operators η k , η † k in single vectors, then there exists a matrix R h,γ of dimension 2L × 2L that relates η and ǎ such that the Hamiltonian ( 22) is diagonal in terms of the Bogoliubov modes, where is the single-particle dispersion relation and the momenta k j satisfy a particular quantisation condition.In actual applications, we will restrict for simplicity to the line γ = 1 (quantum Ising chain).
In Appendix B, we report in detail how we obtain the matrix R h,γ in that case.Therefore, the eigenstates of H are the different configurations K of occupied Bogoliubov modes where |∅⟩ is the Bogoliubov vacuum defined as η k |∅⟩ = 0 for all k satisfying the quantisation condition.In particular, for h < and the coefficients β m are determined by the Bogoliubov rotation R h,γ of Eq. ( 27) that relates η and ǎ.One can check that the energy of this boundary mode b tends to zero as L → ∞.Therefore, in the thermodynamic limit, the states |∅⟩ and b |∅⟩ are degenerate and the ground state spontaneously breaks the Z 2 parity symmetry since ⟨σ x j ⟩ ̸ = 0.The state (32) is a linear combination of Slater determinants that in general does not satisfy Wick theorem.This means that its reduced density matrix ρ A is not Gaussian and the usual methods to univocally determine it in terms of the two-point fermionic correlations [59], see also Appendix C, cannot be applied in this case.For example, let us take h = 0 and γ = 1.One of the two degenerate symmetry-broken ground states (32) is In this configuration, it is easy to see that a non-zero order parameter implies a non-zero odd-point fermionic correlation function, and, therefore, Wick theorem is not satisfied.

Generalised Wick theorem
For symmetry-breaking states of the form (32), it is actually possible to obtain a generalised version of Wick theorem that allows to express any correlator in terms of the one-point and two-point fermionic correlation functions.The rules that we must apply depend on the behaviour of the operators under the parity transformation (18).For even operators [O e , P ] = 0, which correspond to even-point fermionic correlators, we have where we took into account that the odd-point functions ⟨∅| bO e |∅⟩ = ⟨∅| O e b |∅⟩ = 0 since |∅⟩ is a Slater determinant and we can apply Wick theorem.Using the definition (31) of the boundary mode b, we find Observe that, in the correlators that appear in the right hand side of this equation, we can apply the usual Wick theorem and they can be expressed in terms of the two-point fermionic correlation matrix The only extra ingredient are the coefficients β m , which determine the boundary mode b of Eq. ( 31) and actually correspond to the one-point correlation functions For odd operators, {O o , P } = 0, which correspond to odd-point fermionic correlations, one can check that the presence of the boundary mode makes them non-zero, where {O o , ǎm } is an even operator and its expectation value in the vacuum |∅⟩ can be computed from the usual Wick theorem.

Reduced density matrix
In the light of the previous discussion, for states of the form (32), it is in principle possible to use the generalised Wick theorem of Eqs.(36) and (39) for even and odd correlators to reconstruct the reduced density matrix of any subsystem A from the observables with support on it through the expression The problem of this formula is that, since the number of operators in a subsystem grows exponentially with its size, one cannot use it to study large subsystems.However, in Ref. [60], it is proposed an effective expression of ρ A for certain states that break parity symmetry.That ansatz can be generalised to the ground state of Eq. ( 32) for which we find provided A is an interval of contiguous spins starting from the edge as in Fig. 5.In this ansatz, ρ A,e is the even part of ρ A under parity P , which is given by the sum of two Gaussian density matrices where ρ 0 and ρ0 are the reduced density matrices to A of the states |∅⟩ and b |∅⟩ respectively.We have numerically checked that, when we take the thermodynamic limit L → ∞ (semi-infinite chain), ρ 0 ≈ ρ0 ; therefore, in that case, we assume that both ρ 0 and ρ0 are given by the restriction Figure 5: To study the entanglement asymmetry in a quench of the XY spin chain with OBC from the ground state (32), we consider an interval starting from the edge where the boundary mode b is localised.This mode has a characteristic length ℓb ∼ −(log h) −1 .The asymmetry of (32) depends on the fraction of the mode that is contained in the interval, tending to log 2 when ℓ ≫ ℓb.
to A, Γ A 0 , of the correlation matrix Γ 0 in Eq. (37).The operator bA in Eq. ( 41) is the (renormalised) restriction of the boundary mode b of Eq. ( 31) to the subsystem A, One can check that Eq. ( 41) for ρ A correctly reproduces all the expectation values of even and odd operators predicted by the generalised Wick theorem (36) and (39) and, importantly, it allows to efficiently compute the entanglement asymmetry of the parity, as we will see in Sec.4.4.

One and two-point correlation functions
In the ground state: In order to compute the one and two-point correlation functions of Eqs.(38) and (37) in the ground state (32), the knowledge of the matrix R h,γ of Eq. ( 27) that diagonalises the fermionic Hamiltonian ( 22) is sufficient.In fact, the vacuum |∅⟩ is characterised by the absence of excitations, i.e. ⟨η † k η k ⟩ = 0, and therefore Consequently, by expressing the Majorana fermions ǎm in terms of the Bogoliubov operators η k , η † k using the matrix R h,γ , one can compute the total system two-point correlation matrix (37) as On the other hand, since b = η † k 0 +η k 0 , the one-point fermionic functions β m of Eq. ( 38) are given by the sum of the two first lines of the matrix R h,γ which correspond to the boundary mode according to Eq. ( 27) or, alternatively, After a quench: Once we let the ground state (32) for (h 0 , γ 0 ), with h 0 < 1, unitarily evolve with the XY Hamiltonian (16) of another set of couplings (h f , γ f ), both the one and two-point correlation functions ⟨ǎ m (t)⟩ and ⟨ǎ m (t)ǎ m ′ (t)⟩ can be computed numerically using the Heisenberg picture.The Majorana operators ǎm can be evolved in time using the change of basis that relates the Bogoliubov modes that diagonalise the pre and post-quench Hamiltonians.Since the time evolution in the Heisenberg picture of the post-quech Bogoliubov operators ηk and η † k is determined by we can then apply successive changes of basis to compute the dynamics of the correlation matrix (37), where R h 0 ,γ0;h f ,γ f is the matrix that expresses the modes of the pre-quench Hamiltonian in terms of the new ones, R h 0 ,γ 0 ;h f ,γ f = R h 0 ,γ 0 R −1 h f ,γ f .Similarly, we can get the one-point functions β m (t) after the quench with the equation In this way, we have access to all the information about the dynamics of the symmetry-broken state (32) after the quench.

Entanglement asymmetry
With the results of the previous subsections, we can compute the entanglement asymmetry of the spin parity after a quench from the state (32) for a subsystem attached to the boundary of the chain, as in Fig. 5.
In general, according to its definition (4), the calculation of the Rényi entanglement asymmetry requires to determine the moments of ρ A and its projection ρ A,P over the charged sectors of the parity operator P in Eq. ( 18).In the eigenbasis of the parity in the subsystem A, P A , the reduced density matrix ρ A takes the form where the charge sectors ρ ++ and ρ −− are even, i.e. commute with the parity operator P A , while the off-diagonal blocks ρ +− and ρ −+ are odd, i.e. anticommute with it.Thus we can always write Therefore, for the parity operator, the moments of ρ A,P are equal to the ones of the even part of ρ A , Tr ρ n A,P = Tr ρ n A,e , which is a Gaussian operator determined by the restriction to A of the correlation matrix Γ 0 .On the other hand, since the trace of odd operators is zero, the moments of ρ A can be obtained from those of the even part of ρ n A , which we denote as Therefore, the definition of the Rényi entanglement asymmetry in Eq. ( 4) can be re-expressed in the form In Appendix D.1, we employ this result together with the ansatz of Eq. ( 41) for ρ A to obtain an expression for the entanglement asymmetry of an interval attached to one of the boundaries in terms of the one and two-point correlation functions β m and Γ 0 seen before.The final formula reads  A in terms of one and two-point correlation functions after a quench in OBC from the state (32).We consider two quenches to the ferromagnetic and paramagnetic phases and different values of the Rényi index n.The continuous lines correspond to Eq. ( 53) whereas the points have been computed using exact diagonalisation.
where Γ A 0 denotes the restriction to the subsystem A of the matrix Γ 0 , Γ A 1 is a 2ℓ × 2ℓ matrix with entries and The brackets {Γ τ 1 , . . ., Γ τn } stand for the composition of the correlations matrices of Gaussian operators, see Appendix C for its precise definition.While Eq. ( 53) becomes very cumbersome as the Rényi index n increases, it simplifies to a rather simple expression for the small ones; for example, for n = 2, we find ∆S We can check that Eqs. ( 53) and ( 56) are correct by comparing them with the results obtained using exact diagonalisation, as we do in Fig. 6 for Rényi indexes n = 2 and n = 5 in two different quenches.The points are the exact numerical values of ∆S (n) A calculated using the reduced density matrix ρ A determined from the exact diagonalisation of the XY spin chain Hamiltonian (16) with OBC.The continuous lines correspond to Eq. ( 53) for n = 5 and Eq. ( 56) in the case n = 2.
Of course, the exact diagonalisation approach is only efficient for small systems, as the ones considered in Fig. 6.For larger sizes, one must resort to the expression obtained in Eq. (53) in terms of the one and two-point correlations.In what follows, we discuss the results that we find for the entanglement asymmetry using it and taking L → ∞ (semi-infinite line) in quenches from the symmetry-breaking ground states (32) in the ferromagnetic region to different XY Hamiltonians both in the paramagnetic and in the ferromagnetic phases.Before proceeding, it is interesting to analyse with Eq. ( 53) the behaviour of the entanglement asymmetry in the initial ground state (32) as a function of the subsystem size.As Fig. 5 illustrates, A is an interval of contiguous spins of length ℓ attached to the edge of the chain where the boundary mode b of Eq. ( 31) is localised.This mode has characteristic length ℓ b = −1/(log h).Therefore, we expect that the asymmetry of the interval depends on which fraction of the boundary mode is included in the interval: if ℓ < ℓ b, then 0 < ∆S This is confirmed by Fig. 7, where we plot the entanglement asymmetry (53) as a function of ℓ/ℓ b for the ground state (32) with different values of h, obtaining an excellent overlap between them, signalling that the entanglement asymmetry is a scaling function of ℓ/ℓ b.
Quenches to the paramagnetic phase: When we quench to the critical point h = 1 (left panel of Fig. 8) or to the paramagnetic phase h > 1 (right panel of Fig. 8), we observe that the entanglement asymmetry ∆S (n) A tends to zero at large times.This implies that the Z 2 spin parity symmetry, initially broken, is dynamically restored in the subsystem A due to the melting of the boundary mode.As evident in Fig. 8, the main feature is that the entanglement asymmetry does not instantaneously decrease after the quench, but it describes a plateau until times proportional to the subsystem length ℓ.
This result has a very simple interpretation: when quenching to the paramagnetic phase, the boundary mode of the pre-quench Hamiltonian can be expressed as a sum of the quasi-particle excitations of the post-quench Hamiltonian, which all have a non-zero speed of propagation.Consequently, all these modes end up escaping the subsystem.Therefore, the time at which ∆S (n) A starts to drop is roughly the time necessary for the fastest excitations to reach the end of the subsystem, that is t ∼ ℓ/v max where v max = max k |ε ′ k | is the maximal velocity of propagation of the post-quench Bogoliubov modes, which for γ = 1 and h ≥ 1 is one.
In Fig. 9, we study the dependence of the results on the Rényi index n.From that plot, it is clear that it has not great relevance as the entanglement asymmetry essentially shows the same behaviour when varying its value.
Quenches to the ferromagnetic phase: If we quench to the ferromagnetic phase, h < 1, the   32) to the critical line h = 1 (left panel) and to the paramagnetic phase h > 1 (right panel).We take as subsystem an interval of length ℓ attached to the boundary.The points have been calculated using Eq. ( 53).As we discuss in the main text, the asymmetry remains constant after the quench until a time t ∼ ℓ/v max , where v max is the maximum velocity of the post-quench excitations, and then it drops to zero, signaling the restoration of the spin parity symmetry at large times.A for different values of n in a quench from the state (32) to the paramagnetic region.We take as subsystem an interval of fixed length ℓ = 20 starting from the boundary of the chain.The points have been obtained using Eq. ( 53) spin parity symmetry is only partially restored, as shown in Fig. 10.In the left panel, we observe that the farther the initial magnetic field h 0 is from the final one h f , the smaller the stationary value of ∆S (n) A is, i.e. the more the symmetry is restored.This phenomenon is due to the existence of a boundary mode in the post-quench Hamiltonian.In fact, the pre-quench boundary mode b, defined in Eq. ( 31), can be expressed in terms of the post-quench Bogoliubov excitations as Figure 10: Dynamics of the n = 2 entanglement asymmetry of an interval of length ℓ attached to the boundary of a semi-infinite chain after quenches from states of the form (32) with different values of h 0 to the same point in the ferromagnetic region (left panel) and from the same initial state (32) to different points in that phase (right panel).Contrary to quenches to the paramagnetic phase, cf.Fig. 8, the entanglement asymmetry does not vanishes at large times, and the spin parity symmetry is not restored, due to the existence of a boundary mode in the post-quench Hamiltonian, see also the main text.
Therefore, it has a non-zero overlap α 0 with the boundary mode bf of the post-quench Hamiltonian, which has zero velocity.Consequently, part of the initial boundary mode, which is responsible for the parity symmetry breaking at t = 0, stays in the subsystem and the symmetry is only partially restored.In the right panel of Fig. 10, we consider different quenches from the same initial state.We see that the symmetry is more restored in A as we quench closer to the critical point, where it is fully restored.This can be explained by the fact that the length of the post-quench boundary mode diverges as we quench closer to h f = 1 and, moreover, the overlap α 0 diminishes too.

Entanglement asymmetry in a periodic chain
In this section, we study the time evolution after a quench of the entanglement asymmetry of spin parity in a XY spin chain (16) with Periodic Boundary Conditions (PBC).In this case, the presence of the boundary term H b , given by Eq. ( 23), in the fermionic Hamiltonian (22) makes it non-quadratic.However, using the fact that the Hamiltonian ( 16) is invariant under parity P transformations (18), we can solve this issue by separating the Hilbert space into two sectors, called Neveu-Schwarz (NS) for positive parity and Ramond (R) for the negative one.Then the XY Hamiltonian ( 16) can be written as where H N S and H R are the fermionic Hamiltonian (22) with antiperiodic and periodic boundary conditions respectively according to Eq. ( 23).Note that both commute with the parity P .Therefore, if we restrict to one of these sectors, the Hamiltonian is quadratic and, by performing a Bogoliubov rotation (27) as in the case of open boundary conditions, we can diagonalise them.
For the NS sector, we find where ϵ k is the one-particle dispersion relation introduced in Eq. ( 29).In this case, the momenta k j satisfy the quantization condition For the R sector, we have with the momenta p j satisfying These two sectors define two different vacuums, |∅⟩ N S and |∅⟩ R , such that If the chain has a finite size, then its ground state is the vacuum |∅⟩ N S .In the thermodynamic limit L → ∞, as we already mentioned, the states |∅⟩ N S and η † 0 |∅⟩ R become degenerate in the ferromagnetic phase h < 1 and the ground state is the linear superposition which spontaneously breaks the Z 2 parity symmetry (21).As happens for the ground state (32) with OBC, the state (64) does not satisfy the Wick theorem.However, the main difference is that the odd part in Eq. ( 64) consists of an excitation over a vacuum (R vacuum) different from the the vacuum of the even part (NS vacuum).Thus one cannot use the ansatz (41) for ρ A in OBC.
We have to take an alternative route based on cluster decomposition to calculate the entanglement asymmetry in a quench from the state (64).
In the thermodynamic limit, the correlation functions of even operators O e are equal for |∅⟩ N S and η † 0 |∅⟩ R [28].Therefore, the even-point fermionic correlations functions of the state (64) can be calculated using Wick theorem and, consequently, are fully characterised by the two-point correlation matrix Γ 0 defined in Eq. ( 37), which can be determined using the same procedure as in OBC.In the thermodynamic limit, it takes the following form after the quench [28] where Here ∆ k is the difference ∆ k = θ f k − θ 0 k between the Bogoliubov angles of the post and pre-quench Hamiltonians with couplings (h f , γ f ) and (h 0 , γ 0 ) respectively and εk is the dispersion relation of the post-quench Hamiltonian.The Bogoliubov angle for arbitrary couplings h, γ is given by On the other hand, the calculation of odd-point fermionic correlations in the state (64) is a rather difficult task.However, following the ideas in Ref. [61], we can use the cluster decomposition principle to determine them from the even-point ones, exploiting the fact that the latter are equal in both parity sectors in the thermodynamic limit and are given by the matrix Γ 0 .In fact, let us take an arbitrary odd operator O o with support in the subsystem A. We further consider the auxiliary spin operator σ x r on the site r that does not belong to the interval A. If we take the even correlation function ⟨O o σ x r ⟩, then we expect than in the infinite distance separation, r → ∞, the spin at the site r is uncorrelated from the spins in the subsystem A and Therefore, Using this idea, in Appendix E, we show that it is possible to obtain the odd part of ρ A for the state (64).Applying some of the methods to compute the entanglement entropy of disjoint intervals in the XY spin chain [61], we find that ρ A can be expressed in terms of two Gaussian operators ρ 0 and ρ which are respectively determined by the restriction to A, Γ A 0 , of the correlation matrix (65) and by Γ A 1 , defined as As we show in Appendix E, the latter can be obtained from Γ 0 using cluster decomposition as a Schur complement.

Entanglement asymmetry
From the expression (71) of the reduced density matrix ρ A for the ground state (64), we can derive using Eq. ( 52) an efficient formula to calculate the entanglement asymmetry of the spin parity symmetry for a large subsystem.Its obtention involves some tedious algebraic manipulations, which are described in detail in Sec.D.2.The final result reads where Γ A 0 and Γ A 1 are the restriction to A of the correlation matrices (65) and (72) respectively, s n = n j=1 χ j and P 1 is the diagonal matrix . The Rényi-2 asymmetry has a simpler expression ∆S (2) As we did for OBC, we check that the formula of Eq. (73) produces the correct entanglement asymmetry by comparing with calculations done using exact diagonalisation for small systems.In Fig. 11, we analyse the time evolution of ∆S (n) A for n = 7 in a particular quench from the state (64).The continuous line corresponds to Eq. (73) while the points have been obtained with exact diagonalisation for two different total sizes L of the chain.Observe that, while the curve and the points match just after the quench, they deviate at large times.The reason is that Eq. ( 73  A after a quench in PBC from the state (64) in terms of the one and two-point fermionic correlation functions.The points are the value of ∆S (n) A obtained using exact diagonalisation for two chains of different length L. The continuous line corresponds to Eq. ( 73), which is valid in the thermodynamic limit.This explains the discrepancies with the points at times in which the finite-size effects become relevant.

gives ∆S (n)
A in the thermodynamic limit L → ∞.In fact, recall that this expression relies on the cluster decomposition between the spins of subsystem A and an auxiliary spin.Therefore, it is expected that our method only agrees with the exact diagonalisation results at early times and small correlations lengths, before finite-size effects, as the revivals of Fig. 11, appear.
The rest of this section is devoted to discussing the results that we obtain for the entanglement asymmetry after quenches from the symmetry-breaking ground state (64) in the ferromagnetic phase to different points of the (h, γ)-plane of the XY spin chain.We first consider the case of a subsystem made of a single site in order to compare the entanglement asymmetry with the usual order parameter ⟨σ x 1 ⟩.We will then take larger subsystem sizes and obtain an analytic expression that describes the time evolution of the asymmetry at leading order.
One-site subsystems -comparison with the order parameter: As we already pointed out in Sec. 3, in the thermodynamic limit, the breaking of the Z 2 symmetry associated to spin parity can be detected with the standard order parameter ⟨σ x 1 ⟩.It is then interesting to compare the results for ⟨σ x 1 ⟩ with those of the asymmetry for a subsystem of one site.By cluster decomposition, we can indeed get ⟨σ x 1 ⟩ 2 by computing the expectation value of a string of Majorana fermions as described in Ref. [28].
According to Eq. ( 40), the reduced density matrix of a one site subsystem is simply and its projection on the Z 2 parity sectors is Using now the anticommutation relations of Pauli matrices, one can derive Figure 13: Left panels: Time evolution in an infinite chain of the Rényi entanglement asymmetry of spin parity for different quenches in which the chain is initially prepared in the symmetry breaking ground state (64) for a specific set of couplings (h 0 , γ 0 ) in the ferromagnetic region and then is let evolved with the XY Hamiltonian (16) with different couplings (h f , γ f ) in the paramagnetic (upper panel) and in the ferromagnetic (lower panel) phases.Each panel corresponds to a particular quench, in the upper one we take two subsystem lengths ℓ and fixed Rényi index n and vice versa in the lower one.Right panels: same plots taking logarithmic scale in the y axis.The points have been obtained using Eq.(73) while the solid curves correspond to the analytic conjecture (79) in which the free parameter C h0,γ0;h f ,γ f has been adjusted to obtain the best agreement with the numerical points.
of the plateau in terms of the boundary mode present in the subsystem A. In PBC, even if we do not have the boundary mode picture, by analogy, we can think that the quasi-particle excitations generated after the quench emerge from the two end-points of the subsystem.This would explain that the length of the plateau is twice shorter, ∼ ℓ/(2v max ), than in OBC (recall that, in that case, the interval is attached to the boundary and the quasi-particles can only emerge from one point).We find that the exponential decay is rather universal for all the quenches.The change of the initial and final couplings (h 0 , γ 0 ) and (h f , γ f ) only affects the decay rate, which is independent of the subsystem size ℓ and the Rényi index n.This observation is crucial since, as we have seen above, both the asymmetry of a single site and the order parameter ⟨σ x 1 ⟩ 2 show the same exponential decay at late times, which is described by Eq. (78).We can use this fact to conjecture an effective analytic expression for the asymmetry of a subsystem of length ℓ after the quench.We find that it  73) and the solid curves correspond to the analytic conjecture of Eq. ( 79).As can be seen, the farther the initial state is from the paramagnetic region the faster the symmetry is restored.Center: Initial state asymmetry for different h 0 and γ = 1.The points have been obtained with Eq. (73).Right: same plot as in center but at time t = 10 after a quench to the Hamiltonian (h f , γ f ) = (1.2, 1).Comparing center and right panels, we observe that the points are inverted: at any length scale ℓ, the symmetry is more restored at t = 10 by the states that break it more at t = 0, a reminiscence of the quantum Mpemba effect observed in U (1) symmetries.
is actually well reproduced by the formula where C n = 2 + (n − 2)C h 0 ,γ 0 ;h f ,γ f and C h 0 ,γ 0 ;h f ,γ f is an unknown coefficient that depends on the couplings of the pre and post-quench Hamiltonian.Observe that Eq. (79) satisfies the required properties, for t < l/2v max it simplifies to ∆S (n) A (0) = log 2 while at larges times it exponentially decays at leading order as the order parameter (78), which is independent of the subsystem lenght ℓ and Rényi index n, and in particular is valid in the limit n → 1.The solid curves in all the plots of Fig. 13 correspond to the analytic expression (79).
For each quench, the value of the unknown coefficient C h 0 ,γ 0 ;h f ,γ f has been adjusted by hand to obtain the best agreement with the numerical points calculated using Eq. ( 73).Quantum Mpemba effect?In the case of the U (1) symmetry associated to the transverse magnetisation, a quantum version of Mpemba effect has been observed, in which the symmetry is restored faster after a quench for the initial states that break it more [22].For the Z 2 spin parity symmetry, we cannot in principle observe the same effect since the asymmetry at time t = 0 converges to log 2 very fast as the length ℓ of the subsystem increases.Therefore, for large enough intervals, the spin parity symmetry is maximally broken for any initial state of the form (64) that we can consider.However, when studying a series of quenches to the same point in the paramagnetic region from different initial states such as the ones plotted in the left panel of Fig. 14, it is very tempting to conclude that the symmetry is faster restored for the ground state corresponding to h 0 = 0 than for h 0 = 0.2 or h 0 = 0.4, which are parameters closer to the paramagnetic region.
On the other hand, it is important to remember that a value log 2 for the asymmetry of a ℓ site subsystem means that the symmetry is completely broken at that length scale.If we analyse the entanglement asymmetry of the initial states (64) as a function of the subsystem size, as we do in the middle panel of Fig. 14, we can see that, for any subsystem length, the entanglement asymmetry for h 0 = 0.8 is smaller than for h 0 = 0.5 and h 0 = 0.2.In the right panel of Fig. 14, we plot the value of entanglement asymmetry at time t = 10 after performing a quench to the same point (h f , γ f ) = (1.2, 1) from the initial states considered in the middle panel.As we see, at t = 10, the points are inverted with respect to the ones at initial time, now the symmetry is more broken for the state corresponding h 0 = 0.8 than for h 0 = 0.5 and h 0 = 0.2.This is a phenomenon somehow reminiscent of the quantum Mpemba effect found in the U (1) case in the sense that the more the spin parity symmetry is broken at any length scale, the faster it restores at any length scale.

Conclusions
In this manuscript, we have extended to finite discrete cyclic groups Z N the study of the recently introduced entanglement asymmetry.Contrary to the U (1) symmetry considered in Ref. [22], the entanglement asymmetry is bound by the dimension of the group.As an application, we have examined the entanglement asymmetry of the Z 2 group generated by the spin parity in the XY spin chain.In particular, we have analysed its time evolution after a quench from the ferromagnetic phase, where this symmetry is spontaneously broken by the ground state in the thermodynamic limit.This extends the analysis carried out in Ref. [27,28] using the standard order parameter.While the order parameter could be employed as a measure of symmetry breaking for one-site subsystems, the longer range correlations that break the symmetry in extended subsystems are only taken into account in the asymmetry, which is then a proper quantifier of symmetry breaking for larger subsystems.Furthermore we showed that the asymmetry is more suited to detect symmetry breaking since we found examples in which the order parameter vanishes when the symmetry is still broken.
We were faced with the technical issue that the states that break this Z 2 symmetry do not satisfy the Wick theorem, and the powerful techniques for Gaussian reduced density matrices employed in the case of the U (1) symmetry in Refs.[22] and [23] cannot be directly applied here.We have solved this problem obtaining the reduced density matrix by using a generalised version of Wick theorem that takes into account the one-point functions in open boundary conditions and by applying the cluster decomposition principle when the chain has periodic boundary conditions.
According to our results, the entanglement asymmetry rapidly saturates for large subsystems to its maximal value, log 2, in all the ferromagnetic phase, being zero in the paramagnetic region; therefore, it acts as a topological indicator of symmetry breaking for large subsystems.The entanglement asymmetry also allows to determine when and how the Z 2 parity symmetry is restored in a subsystem after a quench to another set of couplings of the XY spin chain.In particular, we have found that, just after the quench, the entanglement asymmetry describes a plateau of length proportional to the subsystem size and then it suddenly decreases its value.For open boundary conditions, if the quench is to the paramagnetic region, the asymmetry of an interval starting from the boundary drops to zero at large times and the symmetry is fully restored, while in a quench to the ferromagnetic phase, it tends to a non-zero value and the symmetry is only partially restored due to the boundary mode of the post-quench Hamiltonian.For periodic boundary conditions, the symmetry is fully restored in all the cases.Remarkably, we have found that, at large times and any Rényi index, the entanglement asymmetry shows the same exponential decay to zero as the order parameter, for which an analytic expression was derived in Refs.[27,28].From it, we have conjectured an analytic expression for the entanglement asymmetry that effectively reproduces its time evolution after the quench.

B Diagonalisation of the XY spin chain Hamiltonian with OBC
As we discuss in Sec.4.4, in order to diagonalise the Hamiltonian (16) for Open Boundary Conditions, it is sufficient to know the rotation matrix R h,γ between the Majorana and Bogoliubov fermions introduced in Eq. ( 27).In the quantum Ising line, γ = 1, we obtain it following the detailed discussion presented in Ref. [62].First, it is necessary to solve the quantisation condition with k ∈ (0, π) and j ∈ {1, 2, . . ., L}, which has L − 1 real solutions for h < 1 and L solutions for h ≥ 1.For h < 1 the last mode is the boundary mode, given by the complex momentum k 0 = π + iν where ν is obtained by solving tanh(ν(L + 1)) = sinh ν h + cosh ν with ν > 0. (89) Once we have numerically obtained the momenta with Eq. ( 88), we can construct the modes η k that diagonalise the Hamiltonian (22) as where A k is a normalisation constant, Similarly, the boundary modes are given by A k 0 (−1) j (sinh(ν(L + 1 − j))ǎ 2j−1 + i sinh(νj)ǎ 2j ) , (93) where

C Gaussian density matrices and their composition rules
In general, a Gaussian density matrix is generated from a quadratic form such that where the factor Z(W ) ensures that Tr(ρ W ) = 1.A Gaussian density matrix ρ W satisfies the Wick theorem and, therefore, it is univocally determined by the two-point correlation matrix [59] δ mm ′ + iΓ mm ′ = Tr(ρ W ǎm ǎm ′ ), (97) which can be written as We can thus use the notation ρ W ≡ ρ[Γ].For large subsystems, it is in practice impossible to directly use full reduced density matrices because their dimension grows exponentially as 2 ℓ with the subsystem size ℓ, whereas correlation matrices have a size 2ℓ × 2ℓ.The computation of the Rényi entanglement asymmetry (4) involves traces of products of different Gaussian matrices Tr (ρ[Γ 1 ] • • • ρ[Γ n ]).To that end, let us first consider two Gaussian density matrices ρ[Γ] and ρ[Γ ′ ] associated to the correlation matrices Γ and Γ ′ .The special properties of Gaussian operators allow to write their product in the form [61,63] where the composition rule × for correlation matrices is defined as and the trace of two Gaussian density matrices can be obtained from their respective correlation matrices as where µ j are the eigenvalues of ΓΓ ′ with halved degeneracy.Observe that, in principle, there is an ambiguity in the global sign in the last equality of the expression above, which is fixed in our case since the traces we compute always have a positive sign.Using these rules, any trace of a product of Gaussian operators can be expressed as where can be computed by recurrence.

D Derivation of the expressions to calculate ∆S (n) A
In this Appendix, we describe in detail how we obtain the formulas of Eqs. ( 53) and (73) that we use in the main text to compute efficiently the time evolution of entanglement asymmetry after a quench in OBC and PBC respectively.

D.1 In OBC
According to Eq. ( 52), in order to calculate the Rényi entanglement asymmetry we must determine the traces Tr ρ n A,e and Tr(ρ n A | e ).In Sec. 4, we have seen that in OBC the reduced density matrix ρ A is given by the ansatz of Eq. (41).Therefore, in this case, the even and odd parts of ρ A are on reduced density matrices of disjoint intervals found in Ref. [61].Indeed, correlation matrices of the form (136) can be computed as Schur complements.Since the string of Eq. (135) overlaps with subsystem A at m = 1, we need to consider the subsystem A ′ as defined in Fig. 15.Then we introduce the correlation matrix Γ A ′ 1,r with entries δ mm ′ + i(Γ A ′ 1,r ) mm ′ = ⟨ǎ m ǎm ′ P R ⟩, with m, m ′ ∈ A ′ .
Let us now consider the restriction of the two-point correlation matrix Γ 0 to the subsystem R∪A ′ , .
The correlation matrix Γ A ′ 1,r can then be obtained as a Schur complement [61] and it has the following structure: .
By definition (Γ A ′ 1,r ) mn = (Γ A 1 ) mn for the sites m, n ∈ A ∩ A ′ .Using now the translational invariance of the chain and the fact that when r → −∞ we can take a separation of one less site with no consequence, we have Note that all this construction is based on the assumption that ⟨σ x 1 ⟩ ̸ = 0.It is a priori possible to derive ρ A,o from other odd correlators, but we found that it is always sufficient numerically to use this expression.When ⟨σ x 1 ⟩ vanishes because of oscillations after the quench, we simply used the same reasoning applied to σ y 1 and σ y r to cure the divergence.

Figure 1 :
Figure 1: Left panel: the symbols are the entanglement asymmetry for different subgroups Z N of the U (1) group of rotations around the z in the tilted ferromagnetic state (5) with θ = π/4.We represent them as a function of the subsystem length.For large subsystems they tend to log N .The solid line is the asymptotic entanglement asymmetry for the U (1) group, cf.Eq. (6), which diverges logarithmically with ℓ.Right panel: the exact analytic expression(11) for the subgroup Z 2 as a function of the tilting angle θ of the ferromagnetic state and different subsystem sizes.

Figure 2 :
Figure 2: Rényi entanglement asymmetry in the tilted ferromagnetic state for the subgroups Z 2 (left panel) and Z 3 (right panel) of the U (1) group of rotations around the z axis.We take different Rényi indices n and subsystem length ℓ = 5.

Figure 4 :
Figure 4: The initial ground state (32) considered to study the entanglement asymmetry of spin parity in OBC contains a boundary Majorana excitation b localised at the extreme of the chain indicated above.

Figure 6 :
Figure 6: Numerical check of the expression (53) to compute ∆S (n)

Figure 7 :
Figure7: n = 2 entanglement asymmetry of spin parity for the state(32) and for different values of h < 1 and γ = 1 as a function of ℓ/ℓb.Here ℓb is the characteristic length of the edge mode b localised at the boundary of the chain where the subsystem of length ℓ starts.The perfect overlap of the points with different h shows that the entanglement asymmetry of the state (32) is a scaling function of the ratio ℓ/ℓb.

Figure 8 :
Figure8: Evolution of the n = 2 entanglement asymmetry in a semi-infinite chain for a quench from the state (32) to the critical line h = 1 (left panel) and to the paramagnetic phase h > 1 (right panel).We take as subsystem an interval of length ℓ attached to the boundary.The points have been calculated using Eq.(53).As we discuss in the main text, the asymmetry remains constant after the quench until a time t ∼ ℓ/v max , where v max is the maximum velocity of the post-quench excitations, and then it drops to zero, signaling the restoration of the spin parity symmetry at large times.

Figure 9 :
Figure 9: Time evolution in a semi-infinite chain of the Rényi entanglement asymmetry ∆S (n)

A / log 2 PBC, h f = 1 . 2 , γ = 1 , ℓ = 40 hA / log 2 PBC 8 Figure 14 :
Figure14: Time evolution of the n = 2 entanglement asymmetry of spin parity in a quench from ground states (64).Left: different h 0 and γ 0 = 1 taking the same post-quench Hamiltonian with h f = 1.2 and γ f = 1, i.e., in the paramagnetic region.The points have been obtained using Eq.(73) and the solid curves correspond to the analytic conjecture of Eq. (79).As can be seen, the farther the initial state is from the paramagnetic region the faster the symmetry is restored.Center: Initial state asymmetry for different h 0 and γ = 1.The points have been obtained with Eq. (73).Right: same plot as in center but at time t = 10 after a quench to the Hamiltonian (h f , γ f ) = (1.2, 1).Comparing center and right panels, we observe that the points are inverted: at any length scale ℓ, the symmetry is more restored at t = 10 by the states that break it more at t = 0, a reminiscence of the quantum Mpemba effect observed in U (1) symmetries.

Figure 15 :
Figure 15: Setup for the cluster decomposition described in the text to determine the operator ρ 1 .