Entanglement asymmetry and quantum Mpemba effect in the XY spin chain

Entanglement asymmetry is a quantity recently introduced to measure how much a symmetry is broken in a part of an extended quantum system. It has been employed to analyze the non-equilibrium dynamics of a broken symmetry after a global quantum quench with a Hamiltonian that preserves it. In this work, we carry out a comprehensive analysis of the entanglement asymmetry at equilibrium taking the ground state of the XY spin chain, which breaks the U(1) particle number symmetry, and provide a physical interpretation of it in terms of superconducting Cooper pairs. We also consider quenches from this ground state to the XX spin chain, which preserves the U(1) symmetry. In this case, the entanglement asymmetry reveals that the more the symmetry is initially broken, the faster it may be restored in a subsystem, a surprising and counter-intuitive phenomenon that is a type of a quantum Mpemba effect. We obtain a quasi-particle picture for the entanglement asymmetry in terms of Cooper pairs, from which we derive the microscopic conditions to observe the quantum Mpemba effect in this system, giving further support to the criteria recently proposed for arbitrary integrable quantum systems. In addition, we find that the power law governing symmetry restoration depends discontinuously on whether the initial state is critical or not, leading to new forms of strong and weak Mpemba effects.


Introduction 3
Charged moments and XY spin chain 5 Entanglement asymmetry in the ground state of the XY spin chain Hot water may freeze faster than cold water: this counter-intuitive statement describes the Mpemba effect.Such phenomenon was already known to Aristotle and was neglected until 1963 when a student called E. B. Mpemba observed it preparing an ice-cream [1].This observation has opened a new research activity devoted to understanding the mechanism and conditions behind the Mpemba effect.Indeed, it has been observed not only in a solution of milk and sugar or in water but in a wide variety of systems, including clathrate hydrates [2], polymers [3], magnetic alloys [4], carbon nanotube resonators [5], granular gases [6], or dilute atomic gases [7] to cite some of them.Today, the Mpemba effect is more generally rephrased as an anomalous relaxation phenomenon where a system initially further out of equilibrium relaxes faster than a system initially closer to equilibrium.Recently, a theoretical framework for the Mpemba effect was developed in Ref. [8,9], followed by a demonstration of the effect in a controlled experimental setting consisting of a colloidal system that is suddenly quenched by placing it in a thermal bath at a lower temperature [10].Further aspects of this framework have been studied, e.g., in [11][12][13][14].We emphasize that one important aspect of these works is the introduction of a distance between the state of the system and the final equilibrium state to characterize the Mpemba effect.Despite considerable effort to understand this phenomenon at a classical level, there are only a few investigations in the quantum realm.Most of them study the relaxation of quantum systems after a quench of the temperature or are subject to non-unitary dynamics [15][16][17][18][19][20].However, a version of the Mpemba effect in a closed many-body quantum system at zero temperature has been recently reported in Ref. [21].In particular, if we prepare a spin-1/2 chain in a state that breaks a U (1) symmetry and we evolve the system unitarily with a Hamiltonian that preserves it, the symmetry may be dynamically restored in a subsystem of the chain and, furthermore, the more the symmetry is initially broken, the faster it may be restored.The situation here is slightly different from the standard classical Mpemba effect: since the system is isolated, the local equilibrium state does depend on the initial state.Then, what defines the quantum Mpemba effect in this context is not the distance from a common asymptotic state but rather the amount of symmetry breaking.
Thus, in order to study the quantum Mpemba effect, we have to use a quantity that does a similar job as the distance considered in Refs.[8,10] to probe the classical counterpart, but at the level of symmetry breaking.To this aim in Ref. [21] the entanglement asymmetry was introduced to measure how much a symmetry is broken in a part of an extended quantum system.So far, the entanglement asymmetry has been studied for the U (1) symmetry associated with transverse magnetization (particle number) in global quantum quenches to the XX spin chain from both the tilted ferromagnetic and Néel states, see Refs.[21] and [22] respectively.While in the first case, the quantum Mpemba effect can be observed, for the tilted Néel state the symmetry is not restored after the quench since the reduced density matrix relaxes to a non-Abelian Generalized Gibbs ensemble.In this case, the asymmetry tends at late times to different non-zero values depending on the initial state, and one cannot define a quantum Mpemba effect.The entanglement asymmetry and the quantum Mpemba effect have also been analyzed in quenches from different initial states to interacting integrable Hamiltonians in the recent Ref. [23], in particular, the Lieb-Liniger model and the rule 54 quantum cellular automaton, using the space-time duality approach developed in Ref. [24].In addition, a general explanation of the microscopic origin of the quantum Mpemba effect in free and interacting integrable systems has also been proposed in Ref. [23].Furthermore, experimental confirmations of this effect have been reported in a trapped-ion setup [25].Entanglement asymmetry has also been employed to analyze the breaking of discrete symmetries in the XY spin-chain [26] and the massive Ising field theory [27], and of compact groups in matrix product states [28].
The goal of the present paper is twofold.On the one hand, we perform a comprehensive analysis of the entanglement asymmetry in the ground state of the XY spin chain, which is the most paradigmatic free integrable system that breaks a U (1) symmetry.On the other hand, taking this ground state, we investigate the time evolution of the entanglement asymmetry after a sudden global quench to the XX spin chain Hamiltonian, which respects the U (1) symmetry.This framework provides the ideal setup to further study the quantum Mpemba effect discovered in Ref. [21] in free fermionic systems and give support to the general mechanism presented in Ref. [23] for integrable models.
Entanglement asymmetry: Before summarizing our main results, let us first define the entanglement asymmetry.We consider an extended quantum system in a pure state |Ψ⟩, which we divide into two spatial regions A and B. The state of A is given by the reduced density matrix ρ A obtained as ρ A = Tr B (|Ψ⟩ ⟨Ψ|), where Tr B denotes the partial trace in the subsystem B. Let us denote by Q the charge operator with integer eigenvalues that generates a U (1) symmetry group.We require that Q is the sum of the charge in each region, If |Ψ⟩ has a defined charge, i.e. it is an eigenstate of Q, then it respects the corresponding symmetry and [ρ A , Q A ] = 0.The latter implies that ρ A is blockdiagonal in the eigenbasis of Q A and each block corresponds to a particular charge sector.This situation has recently been intensively studied in the context of entanglement since entanglement entropy [29][30][31][32] and other entanglement measures [33][34][35][36] admit a decomposition in the charge sectors of the theory, which provides a much better understanding of numerous features of quantum many-body systems [37][38][39][40][41][42][43][44][45][46][47][48].
On the other hand, if |Ψ⟩ is not eigenstate of Q, then it breaks the U (1) symmetry generated by Q and [ρ A , Q A ] ̸ = 0. Therefore, ρ A is not block-diagonal.In this case, a proper measure of how much the symmetry is broken in the subsystem A is the entanglement asymmetry, denoted by ∆S A , and defined as In this definition, the density matrix ρ A,Q is the result of projecting ρ A over all the charge sectors of Q A ; that is, ρ A,Q = q∈Z Π q ρ A Π q , where Π q denotes the projector onto the eigenspace of Q A with charge q ∈ Z.The matrix ρ A,Q is therefore block-diagonal in the eigenbasis of Q A .One can check that, due to the form of ρ A,Q , the entanglement asymmetry ∆S A is equal to the relative entropy between ρ A and ρ . This identity implies that the entanglement asymmetry is non-negative, ∆S A ≥ 0. The other important property as measure of symmetry breaking is that ∆S A vanishes if and only if [ρ A , Q A ] = 0; that is, when the state of A respects the symmetry associated to Q A .The entanglement asymmetry ∆S A can be computed from the moments of the density matrices ρ A and ρ A,Q by applying the well-known replica trick for the entanglement entropy [50,51].If we define the Rényi entanglement asymmetry as one has that lim n→1 ∆S A is easier to calculate for positive integer n values, for which it can be measured in ion trap experiments using protocols based on randomized shadows [25,[52][53][54][55].Moreover, ∆S

(n)
A satisfies the two crucial properties to be a measure of symmetry breaking: it is non-negative [56] and is zero if and only if Main results: As we have already mentioned, the goal of this work is to expand the analysis done in [21] for the tilted ferromagnetic state.We study the entanglement asymmetry in the ground state of the XY spin chain, described by the Hamiltonian where the σ β j are the Pauli matrices at the site j, γ is the anisotropy parameter between the couplings in the x and y directions of the spin and h is the value of the external transverse magnetic field.When the anisotropy parameter γ is not zero, H breaks the U (1) symmetry generated by the transverse magnetization Therefore, in this work, we are interested in the region γ ̸ = 0 for which the ground state entanglement asymmetry associated to Q is non-zero while for γ = 0 it vanishes.For γ ̸ = 0, the XY spin chain is critical along the lines |h| = 1, which belong to the Ising universality class.The tilted ferromagnetic states considered in Ref. [21] are only a subset of the ground states of the Hamiltonian (3) along the curve γ 2 + h 2 = 1 [57,58].In this more general setup, we can compute the entanglement asymmetry for any ground state of (3) and we find that, for a subsystem A of contiguous spins of length ℓ, it reads where g(γ, h) is a function depending on the two parameters γ and h of the Hamiltonian (3).We find that this term is related to the density of Cooper pairs, which are responsible for the breaking of the conservation of the number of particles.We remark that ∆S (n) A increases logarithmically with the subsystem size, both if the system is critical (|h| = 1) or not.
If we choose the ground state of the Hamiltonian in Eq. ( 3) for arbitrary γ and h and we let it evolve with the XX spin chain, that is taking γ = 0 and h = 0 in Eq. ( 3), which commutes with the charge (4), the symmetry is dynamically restored.We derive a quasi-particle picture for the entanglement asymmetry at large times after the quench based on the initial density of Cooper pairs.From it, we find that the Rényi entanglement asymmetry vanishes for large times as t −3 for any initial value of γ when |h| ̸ = 1 and we predict under which conditions for the parameters (γ, h) we observe the Mpemba effect.It turns out that, if the density of Cooper pairs around the slowest modes of the post-quench Hamiltonian is larger for the state that initially breaks less the symmetry, the quantum Mpemba effect occurs, in agreement with the general findings of [23] for integrable systems.On the other hand, when the system is prepared initially in the critical ground state, i.e. |h| = 1, the Rényi asymmetry vanishes as t −1 for any value of γ.Therefore, for critical systems, we can define a strong version of the Mpemba effect for which the relaxation happens algebraically slower regardless of the initial condition for the non-critical state.
Outline: In Section 2, we provide a recipe to evaluate the entanglement asymmetry for Gaussian fermionic operators such as the reduced density matrix of the ground state of the XY spin chain.Section 3 is devoted to the analysis of the entanglement asymmetry in the ground state of the XY model, while Section 4 studies the time-evolution of the asymmetry after a quench to the XX spin chain and the origin of the Mpemba effect.Finally, we draw our conclusions in Section 5 and we include three appendices with additional results and technical details.

Charged moments and XY spin chain
As we have seen in the previous section, by applying the replica trick, the entanglement asymmetry ∆S A can be computed from the Rényi version ∆S (n) A defined in Eq. (2).The advantage of doing this is that, using the Fourier representation of the projector Π q , the projected density matrix ρ A,Q can be rewritten as Therefore, its moments are given by where α = {α 1 , . . ., α n } and Z n (α) are the (generalized) charged moments with α ij ≡ α i − α j and α n+1 = α 1 .Since in general [ρ A , Q A ] ̸ = 0, the order in which these operators enter in the expression of In this manuscript, we are particularly interested in calculating the charged moments Z n (α) and, from them using Eq. ( 7), the Rényi entanglement asymmetry ∆S (n) A in the ground state of the XY spin chain (3).As well-known, this Hamiltonian is easily diagonalizable as follows [59].We can first map it to the fermionic operators c j = (c † j , c j ) via a Jordan-Wigner transformation, namely By performing now a Fourier transformation to momentum space d k = j∈Z e −ikj c j and then the Bogoliubov transformation the XY spin chain is diagonal in terms of the Bogoliubov modes η k , where ϵ k is the single-particle dispersion relation Thus the ground state is the Bogoliubov vacuum |0⟩ that is annihilated by all the operators η k , i.e. η k |0⟩ = 0 for all k.For γ ̸ = 0, this state breaks the U (1) symmetry associated to the conservation of the total transverse magnetization (4), i.e. [ρ, Q] ̸ = 0 with ρ = |0⟩ ⟨0|, and the asymmetry ∆S (n) A is non-zero.On the other hand, for γ = 0, |0⟩ is an eigenstate of Q and ∆S (n) A vanishes.Therefore, |0⟩ is an ideal state to explore ∆S (n) A .The ground state of the XY spin chain is a Slater determinant and, consequently, the reduced density matrix ρ A is a Gaussian operator in terms of c j [60].This simplifies the calculation of ∆S (n) A since, due to the Wick's theorem, ρ A is univocally determined by the two-point correlation matrix with j, j ′ ∈ A. If A is an interval of contiguous sites of length ℓ, then Γ is a 2ℓ × 2ℓ block Toeplitz matrix; that is, their entries are the Fourier coefficients [61] Γ of the 2 × 2 symbol Under the Jordan-Wigner transformation, the transverse magnetization Q in Eq. ( 4) is mapped to the fermion number operator Q = j (c † j c j − 1/2) and e iαQ A turns out to be Gaussian, too.Therefore, Eq. ( 8) is the trace of the product of Gaussian fermionic operators, ρ A and e iα j,j+1 Q A .As explicitly shown in Appendix B of Ref. [22], using the special properties of Gaussian operators [62,63], the trace of Eq. ( 8) can be re-expressed as a determinant involving the two-point correlation matrix Γ, with Eq. ( 17) allows to exactly compute numerically ∆S and is the starting point to derive analytic expressions for Z n (α) and ∆S

(n)
A for large subsystem sizes.

Entanglement asymmetry in the ground state of the XY spin chain
In this section, we study the entanglement asymmetry in the ground state of the XY spin chain.As we have previously shown, this state is the vacuum |0⟩ of the Bogoliubov modes that diagonalize the Hamiltonian (3) of the chain.Since the reduced density matrix ρ A = Tr B (|0⟩ ⟨0|) is Gaussian, we can apply Eq. ( 17) to study both numerically and analytically the charged moments Z n (α), from which the Rényi entanglement asymmetry can be derived using Eqs.( 7) and (2).

Charged moments
For simplicity, let us first consider the case n = 2 and afterwards we will generalize the results to any n.Observe that, for n = 2, the expression (17) of the charged moments Z n (α) in terms of the two-point correlation function Γ simplifies, after a change of variable α 12 = α, as The matrix Γ α ≡ Γe iαn A is block Toeplitz with symbol Therefore, in Eq. ( 18), we have the product Γ α Γ −α of two block Toeplitz matrices, which in general is not block Toeplitz, and the well-known results on the determinant of this kind of matrices cannot in principle be applied.However, in Ref. [22], we found the following result for the asymptotic behavior of determinants that contain a product of block Toeplitz matrices like the one in Eq. ( 18).If we denote as T ℓ [g] the (ℓd)×(ℓd) dimensional block Toeplitz matrix with symbol the d × d matrix g, then for large ℓ det where the coefficient A is given by If we apply Eq. ( 21) in Eq. ( 18), then we obtain that the n = 2 charged moments behave for large subsystem size ℓ as and In Fig. 1, we numerically test this result.We plot the logarithm of the ground state charged moment Z 2 (α)/Z 2 (0) as a function of the angle α for a fixed subsystem of length ℓ = 40 and two different sets of values for h and γ; in the left panel, we consider h = γ = 0.5 while in the right one we take h = 2 and γ = 0.5.The dots are the exact value of Z 2 (α) calculated using Eq. ( 18) and the solid lines correspond to the asymptotic analytic prediction of Eq. ( 23).As evident in the plot, for |h| ≤ 1, log(Z 2 (α)/Z 2 (0)) presents a cusp at α = ±π/2 while, for |h| > 1, this non-analiticity disappears.In the inset of the right panel, we check that the discrepancy between the analytic prediction and the exact points around α = π/2 is due to subleading corrections in ℓ, see the caption for details.The result of Eq. ( 23) for n = 2 can be rewritten in a more appealing form that straightforwardly suggests its generalization to any integer n ≥ 2. In fact, observe that the coefficient A 2 (α) of Eq. ( 24) can be recast in the following factorized expression where As we show in Appendix A, this result can be extended to any integer n ≥ 2. The charged moments behave for large ℓ similarly to the case n = 2, cf.Eq. ( 23),  27) while the points are the exact values obtained using directly Eq. ( 18).The bottom panel represents the extrapolated data with extrapolation form a+b/ℓ to check that the discrepancy observed in the right panel is only a finite-size effect.To extrapolate the data, we have taken into account the numerical data for ℓ = 40, 50, 60.Interestingly, the extrapolated points are exactly equal to our analytical prediction of the linear growth of log Z 2 (α) in Eq. ( 27) (solid curve).
where the coefficient A n (α) admits the following factorization in the replica space, In Fig. 2, we check numerically Eq. ( 27) for the case n = 3.
It is important to remark that, along the critical lines |h| = 1, we have numerically observed that the expression (8) for the charged moments Z n (α) includes an additional subleading term Z n (α) ∼ Z n (0)e An(α)ℓ ℓ mn(α) .Unfortunately, the explicit form of m n (α) cannot be obtained with the methods employed in this manuscript.However, since the factor ℓ mn(α) produces a Figure 2: Logarithm of the n = 3 charged moment Z 3 (α 1 , α 2 , α 3 )/Z 3 (0, 0, 0) in the ground state of the XY spin chain as a function of α 2 for α 1 and α 3 constant, different values of the couplings h and γ and subsystem size ℓ = 40.In the upper panels, we take its real part while the corresponding imaginary part is represented in the plots of the lower row.The points correspond to the exact numerical values calculated using Eq.(17).The solid lines correspond to the asymptotic expression of Eq. ( 27) employing as coefficient A 3 (α) the prediction of Eq. (28).
subleading term in the entanglement asymmetry, we can safely neglect it in the rest of the paper.

Asymptotic behavior of the entanglement asymmetry
As we explain in Sec. 2, once we have the charged moments (8), the Rényi entanglement asymmetry ∆S

(n)
A can be determined by plugging them into the the n-dimensional integral of Eq. ( 7) and then using Eq.(2).In general, this integral can only be calculated by numerical means but, employing a saddle point approximation, we can derive analytically the asymptotic behavior of ∆S (n) A for large subsystems.
To do so, we can follow the same strategy applied in Ref. [22].By taking into account that the phases α jj+1 satisfy n j=1 α jj+1 = 0, we can reduce the n-fold integral (7) to an (n − 1)-fold one after the change of variables αj = α jj+1 , If we insert in this expression the prediction of Eq. ( 27) for the charged moments at large ℓ, the integral takes the form where we have explicitly used the factorization in the replica space found in Eq. ( 28) for the coefficient A n (α).One can check that there are 2 n−1 points in the region of integration [−π, π] ×(n−1) that satisfy the saddle point condition Around all the saddle points, the integrand of Eq. ( 30) has the same behavior at quadratic order in αj and, therefore, their leading contribution to the integral is the same.Hence, if we expand the exponent in Eq. ( 30) around αj = 0 and we properly count the number of saddle points, then Eq. ( 30) can be approximated by the Gaussian integral with The integral of Eq. ( 32) is solvable using the standard formulae, Finally, plugging this result in Eq. ( 2), we obtain that for the ground state of the XY spin chain, the Rényi entanglement asymmetry behaves as The integral of Eq. ( 33) that gives the term g(γ, h) can be computed explicitly.In fact, if we perform the change of variables z = e ik , it can be rewritten as a contour integral in the complex z-plane.Using then the residue theorem, we find In Figs. 3 and 4, we investigate the validity of Eq. ( 35) for n = 2 and n = 3 respectively.In these plots, we represent the ground state entanglement asymmetry as a function of the subsystem size taking different couplings h and γ.The points are the exact numerical values of ∆S (n) A calculated with Eq. (17).The dashed lines correspond to assume the prediction of Eq. ( 27) for the charged moments and then calculate numerically its exact Fourier transform (7) to get ∆S (n) A .In this case, we obtain a good agreement with the numerical points for all the values of h and γ considered.The solid lines represent the asymptotic behavior obtained in Eq. ( 34) using the saddle point approximation.Observe that, for the range of subsystem sizes considered, Eq. ( 35) describes well the exact numerical results for |h| ≤ 1 and any γ, both at n = 2 and n = 3.The same occurs for |h| > 1 and γ > 1.However, for |h| > 1 and γ < 1, the saddle point approximation requires to consider larger subsystems.
In Fig. 5, we plot the saddle point approximation of Eq. ( 35) for ∆S A as a function of γ and several fixed values of h (left panel) and viceversa (right panel) taking as susbsystem size ℓ = 1000 in both cases.Observe in the left panel that ∆S A grows monotonically with the anisotropy parameter γ.Therefore, by varying γ, we can tune how much the U (1) symmetry generated by Q is broken.In particular, as we already pointed out, at γ = 0, the Hamiltonian (3) corresponds to the XX spin chain which commutes with Q. Hence the ground state respects the corresponding U (1) symmetry and the entanglement asymmetry is expected to vanish.However, according to the asymptotic expression (35), ∆S A → −∞ when γ → 0. The reason of this apparent discrepancy is that the limits ℓ → ∞ and γ → 0 do not commute.The other remarkable property of the ground state entanglement asymmetry can be seen in the right panel.As evident also from Eq. ( 36), for large ℓ, the entanglement asymmetry is independent of the transverse magnetic field h in the ferromagnetic phase (|h| < 1) while, in the paramagnetic phase (|h| > 1), it monotonically decreases with h.In fact, at h → ±∞, the ground state of the XY spin chain is |↑↑ • • • ↑⟩ and |↓↓ • • • ↓⟩ respectively, which are eigenstates of Q, and ∆S (n) A = 0.When we take this limit in the asymptotic expression (35), the entanglement asymmetry diverges ∆S (n) A → −∞ since the limits ℓ → ∞ and h → ±∞ do not commute, similarly to the case γ → 0.
Finally, it is interesting to note that the asymptotic result (35) for ∆S A admits an interpretation in terms of the density of Cooper pairs in the ground state |0⟩.Observe that the factor g(γ, h) that enters in Eq. ( 35) only depends, as an integral in momentum space, on the quantity sin 2 ∆ k , see Eq. (33).Using the two-point correlation matrix Γ of Eq. ( 15), it is easy to see that sin ∆ k is related to the correlator ⟨d † A is proportional to the logarithm of ℓg(γ, h) according to Eq. (35), it monotonically increases with the density of Cooper pairs present in the state |0⟩ and the U (1) symmetry associated to particle conservation is more broken.In fact, this symmetry is respected if and only if the correlations ⟨d † 2π−k d † k ⟩ vanish, i.e. in the absence of Cooper pairs.This interpretation of Cooper pairs as the excitations responsible of how much the particle number symmetry is broken will be further supported in the next section, where we elaborate a quasi-particle picture for the entanglement asymmetry after a quench in terms of them.

Entanglement asymmetry out-of-equilibrium
In this section, we study the global quantum quench from the ground state of the XY spin chain (3) with γ ̸ = 0, |Ψ(0)⟩ = |0⟩, which breaks the particle number symmetry generated by Q, to the XX spin chain Hamiltonian H XX , which corresponds to take γ = 0 and h = 0 in Eq. (3) and, therefore, it commutes with Q and the U (1) symmetry is expected to be dynamically restored in the subsystem A, i.e. lim t→∞ [ρ A (t), Q A ] = 0. Thus the time-evolved state is In order to evaluate the time evolution of the entanglement asymmetry in this quench protocol, we first derive a quasi-particle description for the dynamics of the charged moments defined in Eq. ( 8).

∆S
(2) A as a function of the subsystem length ℓ for different values of h and γ.The dots are the exact numerical values of the asymmetry computed using Eq.(17).The solid lines correspond to the asymptotic result of Eq. ( 35) while the dashed ones correspond to the evaluation of ∆S (2) A without the saddle point approximation in the Fourier transformation (7) of the charged moments (27).

Time evolution of the charged moments
In Sec. 3, we have exploited the fact that the reduced density matrix ρ A of the ground state of the XY spin chain is Gaussian and, in virtue of Wick theorem, the charged moments Z n (α, t = 0) are univocally determined by the two-point correlation matrix Γ of Eq. ( 17).Since the XX Hamiltonian is quadratic in terms of the fermionic operators c j , Eq. ( 17) also applies for the reduced density matrix ρ A (t) of subsystem A after the quench.Furthermore, given that the post-quench Hamiltonian preserves the translational invariance of the system, the time-evolved two-point correlation matrix Γ(t) is still block Toeplitz and reads [61] where the symbol G(k, t) is now

∆S
(3) A for the ground state of the XY spin chain as a function of ℓ for different values of h and γ.The points represent the exact numerical value obtained with Eq. ( 17).The solid lines are the result of Eq. ( 35) for large subsystem sizes while the dashed ones have been obtained calculating exactly the Fourier transformation (7) of the charged moments Z n (α) using their analytic expression (27).
with cos ∆ k and sin ∆ k defined in Eqs.(11) and ϵ XX (k) = − cos(k) is the one-particle dispersion relation of the post-quench Hamiltonian H XX .
In order to find the analytic expression that describes the charged moments Z n (α, t) in the ballistic regime t, ℓ → ∞ with ζ = t/ℓ fixed, we first determine their stationary value at large times.It can be obtained by averaging the time dependent terms in the symbol G(k, t) of Eq. ( 39).As t → ∞, the terms e ±2itϵ XX (k) average to zero and the symbol reduces to Observe that the correlators ⟨Ψ(t)|c j c j ′ |Ψ(t)⟩ and ⟨Ψ(t)|c † j c † j ′ |Ψ(t)⟩ vanish in the stationary regime.This is the first signature of the dynamical restoration of the particle number symmetry in the subsystem A.
For n = 2, the stationary behavior of Z 2 (α, t) can be determined by applying the conjecture of Eq. ( 21), as we did in Eq. ( 23) for the charged moments of the ground state.In this case, and, using the time-averaged symbol of Eq. ( 40), we find where we have introduced and n(k is the density of occupied modes with momentum k.This result implies that Z 2 (α, t → ∞) ∼ Z 2 (0, t → ∞); in fact, we recover the result predicted in Ref. [61] for the stationary value of the entanglement entropy in this quench protocol.
For n > 2, we cannot employ the conjecture of Eq. ( 21) to derive the stationary value of Z n (α, t) at large times.In general, the expression (17) for the charged moments does not simplify as for the case n = 2, cf.Eq. ( 18), and it contains the inverse matrix (I − Γ(t)) −1 .Nevertheless, in Eq. ( 66) of Appendix A, we report a formula that predicts the asymptotic behavior of a determinant like the one in Eq. ( 17), with a product of block Toeplitz matrices that also includes the inverse of block Toeplitz matrices.Since the time-averaged symbol I − G(k, t → ∞) of the matrix I − Γ(t) is invertible, we can directly apply Eq. ( 66) to (17) in the large time limit, where W j (k) = (I +G(k, t → ∞))(I −G(k, t → ∞)) −1 e iα j,j+1 σz .Using Eq. ( 40) and calculating directly the determinant, we find that is, At this point, we know both the charged moments Z n (α, t) at the initial time from Eq. ( 27) and its asymptotic behavior at t → ∞ in Eq. (45).These two ingredients are enough to reconstruct the dynamics of Z n (α, t) for any finite time t by exploiting the quasi-particle picture of entanglement.The underlying idea is that the pre-quench initial state has very high energy with respect to the ground state of the Hamiltonian governing the post-quench dynamics; hence, it can be seen as a source of quasi-particle excitations at t = 0. We assume that quasi-particles are uniformly created in pairs with momenta ±k and velocity v(k) = dϵ XX (k)/dk.At a generic time t, the entanglement between a subsystem A and B is proportional to the total number of quasi-particles that were created at the same spatial point and are shared between A and B at that moment, which is given by the function min(2t|v(k)|, ℓ).This idea has been firstly proposed to compute the entanglement dynamics after a global quantum quench in [64][65][66].However, we can also apply it here to determine the time evolution of the charged moments Z n (α, t), in the same way as it was done in Refs.[21] and [22] for the tilted ferromagnetic and Néel states respectively.If we subtract from the stationary value (45) of Z n (α, t) its initial asymptotic behavior, obtained in Eq. ( 27), we get the contribution to Z n (α, t) at t → ∞ of the pairs of entangled quasi-particle generated in  (37) for n = 2 (left panel) and n = 3 (right panel).We plot it as a function of t/ℓ taking several initial ground states with different couplings γ, h and various values of the subsystem size ℓ and the phases α j,j+1 .The symbols were obtained numerically using Eq. ( 17) and the continuous lines correspond to the analytic prediction of Eq. ( 47).
the quench and shared between A and B, This expression can be extended to finite times by properly counting the number of entangled excitations that A and B share at each moment.This can be done by simply inserting the function min(2ζ|v(k)|, 1) in the momentum integrals of the right hand side of Eq. ( 46).We then obtain the exact time evolution after the quench of the charged moments (8) in the scaling limit t, ℓ → ∞ with ζ = t/ℓ fixed, where Z n (0, t) and B n (α, ζ) read respectively and The coefficient A n (α) is given in Eq. ( 28).The expression (47) is the main result of this section, and we benchmark it against exact numerical calculations in Fig. 6 taking as initial configuration the ground state of the XY spin chain for different values of the couplings γ and h: the symbols have been obtained using Eq. ( 17), while the solid lines are Eq.( 47).This expression is valid in the limit ℓ → ∞, and we observe that the agreement improves as ℓ increases.= 80  A (t) after the quench (37).The symbols are the exact numerical results for a subsystem of length ℓ = 100 (n = 2) and ℓ = 80 (n = 3), and different initial conditions for γ, h.The continuous lines are our prediction obtained by plugging the charged moments reported in Eq. ( 47) into Eqs.( 7) and ( 2).Right panels: Square of the density of the Cooper pairs at time t = 0 in Eq. (11).The crossing of two densities is a necessary condition for the presence of quantum Mpemba effect, according to the criterion explained in the main text.
We can observe the quantum Mpemba effect if an only if conditions (51) and ( 52) are satisfied.
Using the asymptotic expression (35) for the ground state of the XY Hamiltonian, the condition (51) at t = 0 can be rewritten in terms of the density of Cooper pairs of the two initial configurations as On the other hand, according to Eq. ( 50), the inequality ( 52) is satisfied if and only if It is clear that it is sufficient to enforce this second condition only for large times.Let us then study more carefully the behavior of the function b(ζ, γ 1 , h 1 ) in the limit t → ∞ or, equivalently, the limit ζ → ∞.In this case, it is useful to apply the identity, where Θ is the Heaviside Theta function, such that Θ(x) = 1 when x > 0. Plugging this result in Eq. ( 50), we firstly observe that is non-vanishing for the modes Outside the critical lines |h| ̸ = 1, sin 2 ∆ k vanishes around k = 0 and π and, therefore, we can take the approximation log If we perform the change of variables k ′ = k − π in the second integral of the expression above, we then find where Therefore, the condition (52 for large ζ, to observe the quantum Mpemba effect can be re-expressed in terms of the densities of Cooper pairs in the initial states as Given the form of sin 2 ∆ k , Υ k (γ, h) is a definite positive, even function of k that vanishes at k = 0 for any value of γ > 0 and |h| ̸ = 1.Therefore, there always exists a large enough time t I for which the integral condition of Eq. ( 59) can be replaced by Eqs. ( 53) and ( 60) are the necessary and sufficient microscopic conditions to observe the quantum Mpemba effect between a pair of ground states of the XY spin chain after a quench to the XX spin chain.According to them, the quantum Mpemba effect occurs when the state that initially breaks less the symmetry, and therefore contains a smaller net number of Cooper pairs (condition ( 53)), has instead a larger density of Cooper pairs around the modes with the slowest velocity v(k) (condition ( 60)), which correspond to the momenta k = 0 and k = π.This is a very natural condition since the entanglement asymmetry satisfies the quasi-particle picture of Eq. ( 50) and, therefore, its leading behavior at large times is determined by the slowest excitations.In the right panels of Fig. 7, we plot the function Υ k (γ, h) that enters in the condition (60) for some of the initial states studied in the left panels of that figure: observe that, whenever the inequality ( 60) is met for a pair of couplings (γ, h) that also satisfy (53), the curves that describe their asymmetries intersect at certain time and Eq. ( 52) is fulfilled.
Notice that the simultaneous validity of Eqs. ( 51) and ( 52) then requires that the density of Cooper pairs corresponding to two different quenches should cross, as made explicit in Fig. 7.In addition, we observe that the conditions ( 53) and ( 59) are valid for any value of the Rényi index n.For the condition at t = 0, the reason is that all the dependence on γ and h in Eq. ( 35) is in the term g(γ, h), which is independent of n.For the large time condition, the starting point ( 50) from which it is derived does not depend on n.Many of the former considerations are valid generically in integrable systems [23].Specializing on our quench, we can obtain a set of conditions for the quantum Mpemba effect equivalent to the microscopic ones but only involving the couplings γ, h of the initial states.For the inequality (51) at t = 0, this can be straightforwardly done using the asymptotic expression (35), together with Eq. ( 36) for the term g(γ, h).In the case of the condition (52) at large times, we need to determine explicitly the leading behavior of ∆S (n) For |h| ̸ = 1, this can be done from Eq. ( 57) by expanding the functions v(k) and sin 2 ∆ k around k = 0 or k = π in each integral and k * (ζ) around ζ = ∞.We find that, at leading order in large ζ, i.e. it vanishes for large times as t −3 for any value of γ and h.The fact that the prefactor in Eq. ( 61) monotonically increases as a function of γ and it depends non-trivially on h reflects that it is not enough starting from a state with larger γ to reach before ∆S (n) A (t) → 0, but the dependence on h is crucial to observe the Mpemba effect.Fixing γ, we notice that, for |h| < 1, Eq. ( 61) is a monotonically increasing function of h; since the initial asymmetry grows with γ and does not depend on h in this region, then it is necessary that γ 2 > γ 1 and h 2 < h 1 to satisfy the Mpemba conditions (51) and (52).In particular, they are always met by any pair of ground states with couplings belonging to the curve h 2 + xγ 2 = 1 for a fixed parameter x > 0, which describes an ellipse in the (h, γ)-plane.In fact, for any initial state on this curve, In this case, the prefactor of the t −3 decay is a monotonously decreasing function of γ, and, therefore, we always observe that the more the symmetry is broken, the faster it is restored.Interestingly, for large subsystems, the spectrum of the correlation matrix Γ is the same for all the ground states along a curve h 2 +xγ 2 = 1 and, consequently, they have equal entanglement entropy [69][70][71].On the other hand, the discussion in the region |h| > 1 is more involved because both the initial entanglement asymmetry (35) and its large time behavior (61) are monotonic decreasing functions of h.The replica limit n → 1 in Eq. ( 61) is not well defined.In Appendix B, we carefully A (t) after the quench (37) from a critical state, h = 1 and different γ's.The symbols are the exact numerical results for a subsystem of length ℓ = 80.The continuous lines are our prediction obtained by plugging the charged moments reported in Eq. ( 47) into Eqs.( 7) and (2).As explained in the main text, whatever is the initial value of γ, for large time t the lines collapse and, eventually, the symmetry is restored almost simultaneously, a weak version of Mpemba effect.perform it, starting from the Fourier transform of the charged moments in Eq. (7).The final result reads Observe that, while the Rényi entanglement asymmetry in Eq. ( 61) decays to zero at large times as ℓ 4 /t 3 , in the limit n → 1 it behaves as ℓ 4 log(t)/t 3 , being the logarithmic correction log(t) a particular feature of this case.This also happens for the von Neumann entanglement entropy, as it has been found in [61].
When |h| = 1, we can find an expression similar to Eq. (61).In this case, sin 2 ∆ k ̸ = 0 at k = 0 and the approximation of Eq. ( 57) is not valid.If we take Eq. ( 56) instead and we expand at leading order the integrands around the modes k = 0 and k = π respectively and the function k * (ζ) around ζ = ∞, then we obtain We observe that the behavior of the entanglement asymmetry as a function of ζ is different if the initial configuration is the ground state of a critical Hamiltonian or not: in the former case, it decreases as 1/ζ, while if we start outside the critical line the decay to zero is algebraically faster, as 1/ζ 3 .Therefore, if we consider a critical state and a non-critical one that breaks more the symmetry, the symmetry is always restored faster in the latter.This can be seen as a strong quantum Mpemba effect.In fact, in the classical Mpemba effect, the system relaxes exponentially to the equilibrium state, but in certain particular situations, the decay is exponentially faster, a phenomenon dubbed as strong Mpemba effect [9].By analogy, in our quantum setup, the asymmetry reaches the equilibrium always following a power law but with a smaller exponent in the case of critical states, so in a much slower fashion.In addition, note that the prefactor of Eq. ( 64) does not depend on γ, while the initial entanglement asymmetry along the lines |h| = 1 grows monotonically with γ according to Eq. ( 35).This means that, independently of how much the symmetry is initially broken, for critical states, it is restored (almost) at the same time, as we show in Fig. 8.We can call this phenomenon weak quantum Mpemba effect.
As occurs in the non-critical region, the limit n → 1 in Eq. ( 64) is also not well defined.Repeating the same steps as in Appendix B for the gapped phase, we find which differs from the result of Eq. ( 64) for the Rényi entanglement asymmetry in the logarithmic correction log(t).

Conclusions
In this manuscript, we have investigated the U (1) symmetry breaking in the XY spin chain using the entanglement asymmetry, completing the analysis initiated in Ref. [21] for the tilted ferromagnetic state and specializing the general discussion on the quantum Mpemba effect for integrable systems done in Ref. [23].We have first studied the behavior of the entanglement asymmetry in the ground state of this model, finding that, at leading order, it grows logarithmically with the subsystem size whether (|h| ̸ = 1) or not the system is gapped.We remark that this is quite different with respect to what happens for the total entanglement entropy, quantity from which the entanglement asymmetry is defined: when |h| ̸ = 1, the entanglement entropy saturates to a constant value for large subsystems [72,73], while a violation of the area law occurs only along the critical lines, where the entropy scales logarithmically with the subsystem size [51].Another important result of this work is that we find that the entanglement asymmetry depends on the density of the Cooper pairs, |⟨d 2π−k d k ⟩|, of the ground state.This is a natural result if we take into account that the breaking of the U (1) particle number symmetry in the XY spin chain can be traced back to the presence of superconducting pairing terms in the corresponding fermionic Hamiltonian.
In addition, we have investigated the evolution of the entanglement asymmetry in a global quantum quench, starting from the ground state of the XY spin chain and letting the system evolve with the XX Hamiltonian that preserves the particle number such that the symmetry is dynamically restored in the subsystem.With the help of the quasi-particle picture of entanglement, we have derived a closed-form analytic expression for the asymmetry at large times, from which we have deduced the necessary and sufficient conditions to observe the quantum Mpemba effect in terms of the density of Cooper pairs of the initial states.Essentially, if the density of the slowest Cooper pairs is larger for the state that breaks less the symmetry, then the Mpemba physics shows up, meaning that the more the symmetry is broken, the faster it is restored.The set of microscopic conditions that we obtain here are in agreement with the criteria derived in Ref. [23] for an arbitrary integrable quantum system.
It would be interesting to investigate several questions in the future.The first one is an explanation of the mechanism of the quantum Mpemba effect when the symmetry is restored by a non-integrable Hamiltonian, or the evolution is non-unitary.This analysis has been initiated in Ref. [21] for systems of a few sites, showing the robustness of this phenomenon also if the evolution Hamiltonian is non-integrable.So far, only the breaking of Abelian symmetries has been investigated, but we would like to use the entanglement asymmetry to explore the symmetry breaking of non-Abelian groups.Finally, the analysis done in this manuscript has revealed that the critical lines of the XY model, |h| = 1, are peculiar since an extra term appears in the charged moments.It would be interesting to find its exact expression and determine its (subleading) contribution to the entanglement asymmetry, not only to have a more accurate prediction of it but also to understand if it contains information about the underlying conformal field theory that describes these critical lines.with where W β,j (k) stands for the 2 and take the limit β → ∞, we find Eq.( 27) with By calculating explicitly the determinant in the integrand of Eq. ( 70) for different integer values of n, one can check that this limit actually yields the factorized form of Eq. ( 28) for the coefficient A n (α).

B Large time behavior of the von Neumann entanglement asymmetry
In Eq. ( 61), we notice that the replica limit n → 1 is not well-defined.Therefore, in this Appendix, we carefully derive the asymptotic expression in the limit t → ∞ of the entanglement asymmetry (1).As already observed in Ref. [23], in this regime we can explicitly compute the Fourier transform in Eq. ( 7).Indeed, by expanding the charged moments for small values of (1 − min(2ζ|v(k)|, 1)) (i.e.large ζ), we find The integral above has been done in Eq. [sm-54] of [23], and, by identifying max , we can report here the final result in our case, We can deduce the replica limit n → 1 after doing an analytic continuation of the result above to any complex value of n and, in the large time regime, we find and, finally, we get Eq.( 63) of the main text.
Along the critical lines |h| = 1, close to k = 0, we find Therefore, the main difference with respect to the non-critical case is that the leading term at large ζ in the series of Eq. ( 75) is not j = −1 but we have now to consider all of them.By taking into account that −1 j=−∞ (−1) j /j = log 2, we obtain at leading order in ζ from which Eq. ( 65) is derived.
C Comparison between the charged moments Z n (α) and the FCS The expression for the charged moments in Eq. ( 8) when n = 1 is also known as full counting statistics (FCS), χ(α) = Tr(ρ A e iαQ A ), see Refs.[74][75][76][77][78] for different studies of it in the XY spin chain.Given the result for generic n in Eq. ( 27) for the ground state, one might be tempted to deduce that, if the U (1) symmetry is broken, the charged moments Z n (α) factorize into the product of the FCS with different phases α j,j+1 .However, using the results for the FCS obtained in [74,78], we will show in the following that this is not always true.The FCS can be cast as the determinant of a Toeplitz matrix with symbol f (e i∆ k , α/2) [74,78], where the function f is given in Eq. (26).Thus one can use the theorems on the asymptotic behavior of Toeplitz determinants to analyze χ(α) for ℓ ≫ 1.For |h| > 1 and any value of α or when h < 1 and α ∈ (−π/2, π/2), the symbol f (e i∆ k , α/2) is a non-zero continuous function in k and the Szegő theorem holds, We observe that the integral satisfies the following equality which implies that, in this regime of the parameters, the result in Eq. ( 27) is a factorization of the charged moments into the FCS.However, when |h| < 1 and α ∈ [−π, −π/2] ∪ [π/2, π], the symbol f (e i∆ k , α) acquires winding number +1.In this case, the prediction in Eq. ( 79) is not valid and it must be modified as log χ(α) ∼ ℓ If we consider the analytic continuation of f (e i∆ k , α) from the unit circle z = e ik to the complex plane, then z 0 denotes the zero of such analytic continuation with |z 0 | < 1 and closest to the unit circle z = e ik .This point can be either The presence of this winding number is the responsible that the charged moments Z n (α) do not exactly factorize when ℓ → ∞ into the FCS Tr(ρ A e iα j,j+1 Q A ).In other words, if we compare Eq. ( 27) with ( 79) and (81), the factorization only works in principle when α j,j+1 ∈ (−π/2, π/2) for all j.But taking into account the periodicity properties in α j,j+1 of the charged moments, it can be extended to α j,j+1 ∈ [−π, π] by introducing the parameter σ j , which vanishes if |α j,j+1 | ≤ π/2 and σ j = π otherwise, i.e. we can write e iσ j /2 Tr(ρ A e i(α j,j+1 −σ j )Q A ).
The term σ j = π ensures that we are always in the regime where Eq. ( 79) is valid.The Fourier transform of the FCS yields the probability distribution p(q) for the transverse magnetization Q A (or particle number) to take the value q.We can make a comparison between our final result in Eq. ( 35) and the Rényi-Shannon entropy for the distribution p(q), or Rényi number entropy, where p(q) is the probability for the observable Q A to take the value q.The result for H n reads where the O(1) term does depend on (γ, h) and, in general, it is different with respect to what we find in Eq. (35).In fact, it is clear from that expression that the entanglement asymmetry only takes into account the number of Cooper pairs as the O(1) term only depends on sin ∆ k , and not on the total number of fermions which contribute to H n .

Figure 1 :
Figure 1: Logarithm of the n = 2 charged moment Z 2 (α)/Z 2 (0) in the ground state of the XY spin chain as a function of α for different values of h and γ and subsystem size ℓ = 40.The solid lines correspond to the analytic prediction of Eq. (27) while the points are the exact values obtained using directly Eq.(18).The bottom panel represents the extrapolated data with extrapolation form a+b/ℓ to check that the discrepancy observed in the right panel is only a finite-size effect.To extrapolate the data, we have taken into account the numerical data for ℓ = 40, 50, 60.Interestingly, the extrapolated points are exactly equal to our analytical prediction of the linear growth of log Z 2 (α) in Eq. (27) (solid curve).
2π−k d † k ⟩ by the equality ⟨d † 2π−k d † k ⟩ = i sin ∆ k /2.The modulus |⟨d † 2π−k d † k ⟩| can be thought as the density of Cooper pairs of momentum k that the state |0⟩ contains.Therefore, since ∆S (n)

Figure 5 :
Figure5: Plot of the asymptotic expression(35) in the limit n → 1 of the entanglement asymmetry for the ground state of the XY spin chain as a function of the anisotropy parameter γ and several values of the external magnetic field h (left panel) and viceversa (right panel).In both plots we take as subsystem length ℓ = 10 3 .

Figure 6 :
Figure6: Time evolution of Z n (α, 0) after the quench(37) for n = 2 (left panel) and n = 3 (right panel).We plot it as a function of t/ℓ taking several initial ground states with different couplings γ, h and various values of the subsystem size ℓ and the phases α j,j+1 .The symbols were obtained numerically using Eq.(17) and the continuous lines correspond to the analytic prediction of Eq. (47).

Figure 7 :
Figure 7: Left panels: Time evolution of the Rényi entanglement asymmetry ∆S

1 Figure 8 :
Figure 8: Time evolution of the Rényi entanglement asymmetry ∆S