Entangled multiplets, asymmetry, and quantum Mpemba effect in dissipative systems

Recently, the entanglement asymmetry emerged as an informative tool to understand dynamical symmetry restoration in out-of-equilibrium quantum many-body systems after a quantum quench. For integrable systems the asymmetry can be understood in the space-time scaling limit via the quasiparticle picture, as it was pointed out in Ares et al (2023 Nat. Commun. 14 2036) . However, a quasiparticle picture for quantum quenches from generic initial states was still lacking. Here we conjecture a full-fledged quasiparticle picture for the charged moments of the reduced density matrix, which are the main ingredients to construct the asymmetry. Our formula works for quenches producing entangled multiplets of an arbitrary number of excitations. We benchmark our results in the XX spin chain. First, by using an elementary approach based on the multidimensional stationary phase approximation we provide an ab initio rigorous derivation of the dynamics of the charged moments for the quench treated in Ares et al (2023 SciPost Phys. 15 089). Then, we show that the same results can be straightforwardly obtained within our quasiparticle picture. As a byproduct of our analysis, we obtain a general criterion ensuring a vanishing entanglement asymmetry at long times. Next, by using the Lindblad master equation, we study the effect of gain and loss dissipation on the entanglement asymmetry. Specifically, we investigate the fate of the so-called quantum Mpemba effect (QME) in the presence of dissipation. We show that dissipation can induce QME even if unitary dynamics does not show it, and we provide a quasiparticle-based interpretation of the condition for the QME.


Introduction
Symmetries are often used by physicists as Lego blocks upon which they build our understanding of Nature and its fundamental laws.Indeed, the concepts of symmetry and symmetrybreaking are ubiquitous in Physics.In recent years there has been renewed interest in studying symmetry-breaking both in and out-of-equilibrium systems [1,2,3,4,5,6,7,8,9,10,11]. Let us consider the prototypical setup of the quantum quench [12] in which a system, here a one-dimensional one, is prepared in an initial state and let to evolve under a many-body Hamiltonian H. Let us also assume that the Hamiltonian commutes with a charge operator Q, which generates an Abelian symmetry U (1).One can consider the situation in which the initial state breaks explicitly the symmetry.An intriguing consequence of thermalization, or equilibration to a Generalized Gibbs Ensemble [12] (GGE) is that the symmetry is, at least for typical initial states (see Ref. [13]), dynamically restored at long times [14].Very recently, the question how fast a symmetry is restored as a function of the degree of symmetry breaking in the initial state, and how this is reflected in entanglement-related quantities attracted a lot of attention.The so-called entanglement asymmetry [1] emerged as an informative tool.The asymmetry (see section 3) is defined from the reduced density matrix ρ A of a subsystem A (see Fig. 1).Let us consider the case of a U (1) charge operator Q = Q A + Q Ā, with A a finite region, and Ā its complement.Since the initial state breaks the U (1) symmetry of the model, ρ A does not commute with the charge Q A restricted to the subsystem, implying that ρ A does not exhibit a block-diagonal structure.This means that if one insisted on using a basis of eigenstates of Q A , ρ A would have matrix elements between states with different eigenvalues of Q A .The asymmetry is sensitive to these matrix elements.If the symmetry is dynamically restored, matrix element connecting different charge sectors, and hence the entanglement asymmetry, vanish at long times.Surprisingly, it was observed that in several quenches the more the symmetry is broken in the initial state, the faster is restored [1].This provides a counterpart of the classical Mpemba effect in out-of-equilibrium quantum many-body systems.The classical version of the Mpemba effect can be summarised in the sentence "hot water freezes faster than cold water", since it was primarily observed in water [15].Classical Mpemba effect was observed also in several systems, such as clathrate hydrates [16], magnetic alloys [17], carbon nanotube resonators [18], granular gases [19], colloidal systems [20] or dilute atomic gases [21].A microscopic explanation of the quantum Mpemba effect in both interacting and free integrable spin chains has been provided in [5] (see also [9] for a further study of the asymmetry in free integrable quantum systems).
Very recently, the asymmetry was measured experimentally by using a trapped-ion quantum simulator, and the quantum Mpemba effect was confirmed [11] (see also Ref. [22,23,24,25] for the observation of another version of the quantum Mpemba effect in a single qubit).The dynamics of the entanglement asymmetry has been studied for a quench from the tilted ferromagnetic state [1], and, more in general, from the ground state of the XY spin chain [9], which does not preserve the transverse magnetization.In these situations, the dynamics determined by a free U (1)-invariant Hamiltonian leads to the symmetry restoration of the U (1) charge in the subsystem.Crucially, the dynamics of the entanglement asymmetry, and the onset of the quantum Mpemba effect, can be understood in the space-time scaling limit (hydrodynamic limit) ℓ, t → ∞ with the ratio t/ℓ fixed.Indeed, the main formula obtained in Refs.[1,9] is reminiscent of the so-called quasiparticle picture [26,27,28,29] for the entanglement entropy.Still, a full-fledged quasiparticle picture was not available so far.For instance, a general prescription to determine how different configurations of entangled quasiparticles contribute to the asymmetry was not completely understood.
Here we first review the results of Refs.[1,2,9].We focus on a dynamics governed by the XX chain, considering quantum quenches from the ground state of the XY chain and the tilted Néel state.In this setup, the U (1) charge corresponds to the total magnetization, We are interested in the dynamics of the entanglement asymmetry of a finite region A embedded in an infinite chain.(b) Quasiparticle picture for the asymmetry: the case of entangled pairs.If entangled pairs of quasiparticles are produced after the quench, only the pairs that are fully contained in the subsystem contribute.Here the pair contributes for t ≤ t 1 .At longer times the pair is not fully contained in A, and the asymmetry vanishes.(c) For the tilted Néel state quadruplets of entangled quasiparticles are produced with velocities v 1 = v 2 and v 3 = v 4 .In contrast with (b), in the absence of dissipation the asymmetry does not vanish at long times because there are always paired fermions in A.
which commutes with the Hamiltonian of the XX chain.On the other hand, both the initial states that we consider do not have a well-defined total magnetization, and hence break the U (1) symmetry.The dynamics after the quench from the ground state of the XY chain is understood in terms of the ballistic propagation of pairs of entangled quasiparticles, whereas the quench from the tilted Néel state gives rise to entangled quadruplets of quasiparticles.Interestingly, while in the former case the symmetry is restored at long times, and the entanglement asymmetry vanishes, in the latter this does not happen.This is due to the fact that the XX chain possesses some non-Abelian conservation laws [2,13].Specifically, we show that the framework of entangled multiplets developed in Refs.[30,31,32] is the natural one to build the quasiparticle picture for the entanglement asymmetry.First, by using multidimensional stationary phase methods we provide an ab initio derivation of the results of Ref. [2].In contrast with Ref. [1] (see also Ref. [2]) our derivation does not rely heavily on the Toeplitz structure of the two-point fermionic correlation function.Based on these results, we conjecture the quasiparticle picture for the charged moments of the reduced density matrix, which are the main ingredients to obtain the entanglement asymmetry.Our formula is by construction readily generalized to quantum quenches generating entangled multiplets of arbitrary size.As in the standard quasiparticle picture, the dynamics of the charged moments and the asymmetry depend on the way the members of the entangled multiplet are shared between the subsystem A and the rest.Crucially, in our generalized quasiparticle picture the entanglement asymmetry receives nonzero contributions only from configurations in which at least two entangled quasiparticles are in the subsystem (see Fig. 1) (b).This suggests that in generic quenches governed by the XX Hamiltonian producing multiplets of arbitrary size, the symmetry is always restored at long times.This happens because at long times any two quasiparticles in the multiplet cannot be in A at the same time, provided that the quasiparticles forming the multiplet have different velocities.Interestingly, the last condition is violated in the quench from the tilted Néel state, which gives rise to quadruplets of quasiparticles with velocities v 1 (k) = v 2 (k) and v 3 (k) = v 4 (k), with k the quasimomentum.This means that at long time there is always a pair of entangled quasiparticles in A (see Fig. 1 (c)).This prevents the symmetry restoration, which is reflected in a nonzero entanglement asymmetry.
Next we study the effect of dissipation on the quantum Mpemba effect.Indeed, the question as to whether the Mpemba effect persists in the presence of an environment has attracted recent attention.For instance, in Ref. [11] it was shown that Mpemba effect does not occur under pure dephasing dynamics.On the other hand, it has been shown that Mpemba effect can occur in Markovian dynamics [33,34,23].We focus on the dynamics starting from the tilted Néel state in the XX chain in the presence of gain and loss dissipation, which describe incoherent emission or absorption of fermions from the chain.We employ the framework of the Lindblad master equation [35].Crucially, for the gain and loss dissipation the Lindblad operators are linear, so that the Gaussianity of the initial state is preserved and we can study the quantum Mpemba effect analytically.In contrast with the unitary dynamics, which does not support symmetry restoration and Mpemba effect, we show that in the presence of gain/loss dissipation both symmetry restoration and quantum Mpemba effect can occur.We also derive the conditions that give rise to the Mpemba effect, comparing them with the case without dissipation.
The manuscript is organised as follows.In section 2 we introduce the XY and XX chains.We also introduce the framework of the Lindblad equation and gain/loss dissipation.In section 3, we introduce the entanglement asymmetry, and the charged moments of the reduced density matrix.In section 4, we provide an ab initio derivation of the charged moments in the hydrodynamic limit after the quench from the tilted Néel state in the XX chain.In Section 4.1 we consider the unitary dynamics, whereas in Section 4.2 we discuss the effect of gain and loss dissipation.In Section 5 we conjecture a generalized quasiparticle picture for the charged moments.In particular, in Section 5.1 we introduce the framework of entangled multiplets and the quasiparticle picture for the von Neumann entropy, discussing the generalization to the charged moments in Section 5.2.In Section 5.2.1, by using the generalized quasiparticle picture, we rederive the formula describing the dynamics of the charged moments after the quench from the tilted Néel state in the XX chain.In Section 5.2.2 we obtain the charged moments after the quench from the ground state of the XY chain that was derived in Ref. [9].In Section 6 we discuss the effect of gain and loss dissipation on the quantum Mpemba effect, focusing on the quench from the tilted Néel state in the XX chain.In particular, we discuss the condition ensuring the onset of Mpemba effect.In Section 6.1 we compare that condition with the case without dissipation.In Section 7 we discuss numerical benchmarks for the results of Section 6.Specifically, in Section 7.1 we discuss the dynamics of the charged moments, whereas in Section 7.2 we consider the entanglement asymmetry.We conclude in Section 8. Finally, in Appendix A we detail the derivation of some of the results of Section 4. Some generalizations to the dissipative case are discussed in Appendix A.1.In Appendix B we provide the proof of some of the formulas of Section 5.

Quantum quenches in the XY chain with dissipation
Here we deal with the out-of-equilibrium dynamics after a quantum quench in one-dimensional integrable systems [12].Given an initial state |Ψ(0)⟩ and a Hamiltonian H, the time-evolved state Ψ(t) is given as We focus on the XY spin chain given in terms of the Hamiltonian H as where σ x,y,z j are Pauli matrices acting on site j of the chain, and h is a magnetic field.The chain is infinite.After a Jordan-Wigner transformation, Eq. ( 2) is mapped to a quadratic fermionic Hamiltonian as where c j , c † j are standard fermionic operators.Eq. ( 3) is diagonalized by first rewriting it in terms of the Fourier transformed fermionic operators c k := ∞ j=−∞ e −ikj c j , and then by applying a Bogoliubov transformation as Here, η k , η † k are new fermionic operators, and the Bogoliubov angle ∆ k is defined as The Hamiltonian (3) in terms of η k is diagonal as with ε XY (k) := (h − cos k) 2 + γ 2 sin 2 k the single-particle dispersion.The ground state (not normalized) of the XY spin chain is where |0⟩ is the state annihilated by the Dirac operators c j , and satisfies η k |Ω⟩ = 0 ∀k.For γ = 0, the Hamiltonian (2) becomes that of the XX spin chain given as The fermionic Hamiltonian (3) for γ = h = 0 becomes the tight-binding free-fermion chain Eq. ( 9) can be diagonalised in terms of the Fourier transformed operators c k (cf.( 4)) as where the one-particle dispersion relation is ε XX (k) = − cos(k).For the following, it is important to stress that H XX , unlike the XY chain (2), preserves the transverse global magnetization, i.e., the (8) satisfies In the fermionic language Eq. ( 11) can be written as Q = j (c † j c j − 1/2).For the following, it is also crucial to observe that the Hamiltonian (6) of the XY chain can be written in terms of m different species of quasiparticles, each having a dispersion relation ε j (k) (j = 1, . . ., m), as Here d j (k) are the fermionic operators associated to the different species of quasiparticles and satisfy the usual anticommutation relations {d i (k), d † j (p)} = δ ij δ(k − p), while B is the reduced Brillouin zone B ⊆ [0, 2π].The operators d j (k) are obtained from the operators η(k) diagonalizing (2).The number of species is determined by the prequench initial state.The reason for the rewriting ( 12) is because the prequench initial state |0⟩ can induce non-diagonal correlations in quasimomentum space, i.e., ⟨0|η This means that different species of quasiparticles share nontrivial correlations.Indeed, typically the quasiparticles of different species are entangled, forming an entangled multiplet.These entangled multiplets are essential to accurately describe the dynamics of entanglement-related quantities [31,32,36].
Here we consider the dynamics starting from the initial Gaussian state Ψ(0), defined as where |θ⟩ is the tilted Néel state Crucially, for θ ̸ = mπ, m ∈ Z, it breaks the U (1) symmetry associated with the conservation of the transverse magnetisation (11).
A crucial quantity to determine entanglement properties [37] of the XY chain is the twopoint correlation matrix Γ jj ′ (t) defined as where c j := (c j , c † j ), and c j are the fermionic operators appearing in (9).In this work we are interested in understanding the effect of dissipation on the quantum Mpemba effect.To account for dissipation, we employ the framework of the Lindblad master equation [35].The Lindblad equation describes the dynamics of the full-system density matrix ρ(t) as Here L j,α are the Lindblad jump operators encoding dissipation.We restrict ourselves to the paradigmatic case of gain and loss dissipation.They encode incoherent processes in which fermions are injected or removed in the chain at rates γ + and γ − , respectively.The Lindblad jump operators, L j,α , for gain and loss processes are defined as L j,− := γ − c j and L j,+ := γ + c † j .The generic Lindblad equation is not exactly solvable, not even for integrable Hamiltonians.Remarkably, if H is quadratic in the fermion operators, and the Lindblad operators are linear in the fermions, Eq. ( 16) can be solved analytically via the "third quantization" approach [38].Finally, it is crucial to observe that for gain and loss processes the Gaussianity of ρ is preserved during the time evolution under Eq.( 16).

Entanglement asymmetry & charged moments
Let us review the formal definition of the entanglement asymmetry [1], which is the main topic of this manuscript.We consider a pure state |ψ⟩ living in the bipartition A ∪ Ā (see Fig. 1 (a)).From here, one can build the reduced density matrix ρ A = Tr Ā(|ψ⟩ ⟨ψ|) describing the subsystem A. Let us consider a local charge operator Q, which generates a U (1) symmetry for the full system, and it can be decomposed as operator acting in subsystem A and Ā, respectively.If |ψ⟩ is an eigenstate of Q, [ρ A , Q A ] = 0 and ρ A has a block-diagonal structure in the eigenbasis of Q A , with each block labelled by an integer q ∈ Z.If [ρ A , Q A ] ̸ = 0, ρ A exhibits non-zero entries outside the blocks.Thus, we can define a new density matrix ρ A,Q := q Π q ρ A Π q by using the projectors on the eigenspace of q, Π q , such that is obtained from ρ A by removing the matrix elements that connect blocks of ρ A with different values of Q A .In analogy with the definitions of the Rényi entropies S (n) A (ρ) = 1/(1 − n) ln Trρ n , we can introduce the Rényi entanglement asymmetry as ∆S (n) i.e., the difference between the Rényi entropy obtained from ρ A,Q and the one obtained from ρ A .In the replica limit n → 1, the Rényi entropies provide the von Neuman entropy, S(ρ) = −Tr(ρ ln ρ), and ∆S (n) A reduces to The entanglement asymmetry satisfies two crucial properties.First, ∆S A ≥ 0. Second, it vanishes if and only if the symmetry is not broken, i.e. [ρ A , Q A ] = 0.
Before digging into the physics of the asymmetry, we introduce the charged moments, which play a fundamental role to compute ∆S (n) A .By exploiting the Fourier representation of the projector operator, Π q , the reduced density matrix ρ A,Q reads From ( 19) one obtains where α = {α 1 , . . ., α n } and with α ij := α i − α j and α n+1 ≡ α 1 .The quantities Z n (α) are known as charged moments and, by using (20), allow to obtain the entanglement asymmetry as Crucially, the asymmetry depends only on the ratio Zn(α) Zn(0) .Since at any time, both ρ A (t) and Q A = j∈A c † j c j are Gaussian operators in terms of the fermionic operators c j = (c j , c † j ), the charged moments can be obtained from the two-point correlation function Γ(t) (cf.( 15)) restricted to the subsystem A [1].By using the composition rules for the trace of a product of Gaussian operators [13], the charged moments (21) can be written as [1] Z where W j (t) = (1 + Γ(t))/(1 − Γ(t))e i(α j −α j+1 )N A , with N A a diagonal matrix taking into account the presence of the charge operator Q A in (21) and such that (N A ) 2j−1,2j−1 = −1 and (N A ) 2j,2j = 1.
4 Hydrodynamic limit of the charged moments after generic quenches in the XX chain: A rigorous derivation As we have explained in the previous section, the charged moments in Eq. ( 21) are the starting point to compute the asymmetry.In the following we focus on the dynamics of the charged moments after a quantum quench.Specifically, we provide a rigorous derivation of their leading-order behaviour in the hydrodynamic limit t → ∞, ℓ → ∞ with fixed ratio t/ℓ.We first revisit the non-dissipative case [2] in section 4.1 considering the quench from the tilted Néel state (cf.( 14)) in the XX chain.Specifically, we provide a thorough analytical proof of the results reported in [2].Our proofs do not rely on Toeplitz structure of the fermionic correlators, unlike the results of Ref. [2], and hence is straightforwardly generalizable to other situations.Then we extend our results accounting for the effect of gain/loss dissipation in section 4.2.

Non-dissipative case
We start by rewriting Eq. ( 23) as The correlation matrix Γ(t) (cf.( 15)) can be cast as a block Toeplitz matrix, Here the so-called symbol G t (k) is a 4 × 4 matrix For the quench from the tilted Néel in the XX chain one has [2] C t (k) = cos(2tε(k/2))g 11 (k, θ) e ik/2 (g 12 (k, θ) The entries of ( 27) and ( 28) are given in terms of [2] The main result of this section is that at the leading order in the hydrodynamic limit Eq. ( 24) gives where n,∞ (k) are 4 × 4 matrices (cf.( 45), ( 46), ( 49)) built, respectively, from the symbol G 0 (k) at t = 0 (cf.(25)) and from the non-oscillating part G ∞ (k) (cf.(38)).In (30) we defined v(k (10)).In the following paragraph, we report the ab initio derivation of (30).Notice that Eq. ( 30) was conjecture in Ref. [2] based on robust analytical and numerical evidence.Finally, we remark that the dispersion relation appearing in Eqs. ( 27) and ( 28) is ε(k/2), rather than ε(k).The reason is that the Brillouin zone of the post-quench Hamiltonian has been halved to rearrange the entries of Γ(t) as a block Toeplitz matrix, taking into account the two-site periodicity of the initial state (13) [2].
Proof of Eq. (30).Let us start the proof by first evaluating the moments M a,b defined as where N A originates from W j (t) (cf.( 24)).The moments M a,b appear in the Taylor expansion of (24).Analytical knowledge of M a,b as a function of a, b, α, in the hydrodynamic limit allows one to obtain ln(Z n (α, t)).We exploit the structure for Γ(t) in Eq. ( 25) to write M a,b as where M := a 1 +...+a n is the total number of quasimomentum variables, n A := diag(−1, −1, 1, 1), and we identify k M +1 ≡ k 1 .To proceed, we transform the sums in (32) into integrals, by using the identity We then decompose G t (k) in Eq. ( 26) according to the different temporal dependences, i.e., where In (34) we have factored out the matrix U (k) and its inverse, which cancel out when we take the trace in (32).Now, after substituting the decomposition ( 34) in (32), one obtains 3 M terms.Each term is the trace of a string of powers of G ± (k), G ∞ (k), n A multiplied by a time-oscillating phase.We observe that most of the traces of these strings are zero.Indeed, it is not hard to explicitly check that, because of the symmetries of G ± (k) and G ∞ (k), we have: , where P is any product of matrices G ∞ (k) and n A ; • Tr [G ± (k)P ] = 0, with P as before.
Thus, the only strings of matrices that give a nonzero contribution to the moments (32) are those containing substrings of alternated G + (k) and G − (k), possibly interspersed with factors G ∞ (k) and n A .
Let us now perform the integration in the 2M −2 variables k 2 , ..., k M , ξ 2 , ..., ξ M in Eq. ( 32), applying the stationary phase approximation.In the limit ℓ → ∞, the stationary phase states that where x j ∈ D are the stationary points of f (x), i.e., such that ∇f (x j ) = 0, D is the integration domain, H is the Hessian matrix of f and σ its signature (precisely, the difference between the number of positive and negative eigenvalues).Let us study the stationarity conditions for the generic term in (32).The derivatives of the phase with respect to the ξ j variables give the conditions k j = k j+1 , j = 2, ..., M , which imply k j = k 1 ∀j.By taking the derivative of the phase with respect to the k j variables, we obtain with j = 2, . . ., M .Precisely, the stationarity condition with respect to k j gives (41) if a term G ± (k j ) is present in the string of operators appearing in (32).The ± in ( 41) is the same as in G ± (k j ).If k j is not present in the string, the stationarity condition gives (42).Keeping in mind that only strings with alternate occurrences of G ± (k) can appear in (32), the stationarity conditions (41) imply that either except for the only string containing only G ∞ (k) and n A , which gives ξ j = ξ 1 ∀j.
Since any stationary point has to be inside the integration domain, one must have This condition puts a constraint on ξ 1 , except for the string without any G ± (k).Thus, by integrating over ξ 1 , the terms originating from strings with at least a G ± (k) yield a factor max 0, 2 )), whereas the strings with only G ∞ and n A give a factor 2. We can also compute the determinant of the Hessian matrix evaluated at the stationary points (cf.( 40)), which is 4 2−2M , and its signature, which is simply 0. Substituting everything back into Eq.( 32), and evaluating the functions in the stationary point, we obtain where ).
Notice that now the terms inside the trace in (43) do not depend on time, unlike in (32).The only time dependence is in the kinematic factor min(2v k t/ℓ, 1).Notice also that at t = 0 and t → ∞ only the terms with G 0 and G ∞ remain, respectively.We find it useful to rewrite Eq. ( 43) as Now, let us define ln(Z (24).Analogously, we define ln(Z (24).Thus, analiticity of (24) as a function of Γ and N A implies that the leading order of the charged moments in the hydrodynamic limit is Eq.(30).Specifically, to obtain (30) one can Taylor expand (24), apply (44) to each term of the series, and resum the obtained expression.
To proceed, one has to determine ln(det Z α n,0 ) and ln(det Z α n,∞ ).Crucially, this involves computing functions of 4 × 4 matrices, unlike (24), which requires dealing with 2ℓ × 2ℓ matrices.The derivation is quite technical and we report it in Appendix A. The final result for ln(det where and f ij are defined in (29).Moreover, one obtains ln(det with and Again, Eq. ( 45) and Eq. ( 46) were conjectured in Ref. [2].

Dissipative case
Let us now focus on the time evolution of the charged moments (21) in the presence of fermionic gain and loss processes (see section 2).
In the presence of gain and loss dissipation, the symbol of Γ(t) in Eq. ( 26) must be modified.In particular, the correlations C t and F t (cf.( 26)) are affected differently by the dissipation.A straightforward calculation using the formalism of Ref. [38] shows that in the presence of gain and loss one has that where n ∞ := γ + /(γ + + γ − ), with γ ± the dissipation rates.In (50) the symbol ⊗ denotes the Kronecker product.The term σ z ⊗1 implies that while the correlator C t is shifted and rescaled due to the dissipation, F t is only rescaled.In order to compute the entanglement asymmetry, we start by studying the time evolution of the charged moments in the weakly-dissipative hydrodynamic limit t, ℓ → ∞, γ + , γ − → 0 with fixed t/ℓ, γ + t, γ − t.Similar to the nondissipative case, we anticipate that the leading order contribution in the weakly-dissipative hydrodynamic limit of Eq. ( 24) is given as which is one of the main results of this section.Now ), although they are built from the symbol G ′ t (k) (cf.( 50)).Notice that unlike the non-dissipative case, there is a residual time dependence in (51) due to the terms exp(− (γ + + γ − ) t), which, however, are constant in the weaklydissipative hydrodynamic limit.In the following paragraph, we derive Eq. ( 51).We anticipate that, in the presence of gain and loss, we are able to analytically compute Z n (α, t) only in some particular cases.
Proof of Eq. (51): The derivation of Eq. ( 30) can be adapted to the dissipative case if we split the symbol G ′ t (k) (see Eq. ( 50)) as with Here ε(k) is the energy dispersion (10), and G ± (k) are defined in (38).Notice that in defining G ′ 0 we performed the t → 0 limit in (52) keeping the terms exp(−(γ + + γ − )t) because they are constant in the weakly-dissipative hydrodynamic limit.Clearly, Eq. ( 51) is more involved than in the non-dissipative case.Still, the derivation of (51) proceeds as for the non-dissipative (see section 4.1).Precisely, we can still define the moments M a,b as in (31).In (32), one has to replace G t with G ′ t (cf.( 52)).The conditions on the sequences of G ′ yielding zero trace are the same as for the G, thus Eq. ( 44) remains the same after replacing 53)-( 56)).Finally, this allows us to define Z ′ (α) n,0 (k, t) and Z ′ (α) n,∞ (k, t), obtaining (51).Now, by using (51), one can obtain ln(Z n (α, t)) for any n in the weakly-dissipative hydrodynamic limit.The computation requires diagonalizing 4 × 4 matrices.Crucially, to understand the behavior of the asymmetry ∆S (n) A (cf. ( 17)), one has to know the analytic dependence on n of the charged moments.Unfortunately, we are not able to find a closed-form expression for ln(det Z ′(α) n,0 ) and ln(det Z ′(α) n,∞ ) for generic γ ± .Analytic manipulations are somewhat easier for balanced gain and losses, i.e., for γ + = γ − , which we now discuss.We have n ∞ = 1/2 and Eq. ( 50) becomes , where λ(t) := e −(γ + +γ − )t .Again, even for balanced gains and losses, although one can numerically obtain G n,0 for any integer n, a closed-form analytic expression cannot be obtained.We refer to Appendix A.1 for more details.Nevertheless, we obtain a closed formula for ln(det Z ′ (α) n,∞ ) for generic n.To do that, one has to multiply f 11 (k, θ) and g 12 (k, θ) in Eq. ( 29) by λ(t) in Eq. (46).The derivation is discussed in Appendix A.1.The final result reads as where (49).The difficult term to evaluate for generic n is the short-time one ln det Z ′ (α) 2,0 , but we are, at least, able to find a closed formula for n = 2, reading If we relax the condition γ + = γ − , the formulas become even more cumbersome, and we report them only for small n.For n = 2 we find (59) Similarly, we can derive the longtime term ln det Z ′ (α) n,∞ only for small values of n.We give some detail in Appendix A.1.For instance, for n = 2 we obtain where h′ n (α, k, θ, t) is obtained by substituting The results above allow us to obtain an analytic expression for the leading-order term of ln Z 2 (α, t) for generic γ ± , by plugging Eqs. ( 59) and (60) in Eq. (51).We should stress that the case with n = 2 is sufficient to probe the Mpemba effect, as it has been recently shown [11].
5 Generalized quasiparticle picture for the charged moments Both Eqs. ( 30) and ( 51) have a form that is reminiscent of the quasiparticle picture of the entanglement dynamics [26,27,28,29].The main ingredient of the quasiparticle picture is the presence of entangled multiplets of excitations, which are formed after the quench.In the most basic incarnation, the multiplets are formed by two quasiparticles, i.e., by entangled pairs of quasiparticles.The entangled quasiparticles travel ballistically, as free excitations.The entanglement entropy between two complementary subsystems is proportional to the number of entangled pairs shared between the intervals.The entanglement content shared by each pair is related to the thermodynamic entropy of the GGE describing the post-quench steady state [27,28].For free-fermion and free-boson models, the quasiparticle picture describes the dynamics of both von Neumann and Rényi entropies, and the associated mutual information [39].The quasiparticle picture can be also used to describe the dynamics of the logarithmic negativity [40].The dynamics of quantum information after quenches from inhomogeneous initial states has been addressed [30,41,42,43].Very recently, the quasiparticle picture has been adapted to describe entanglement dynamics in systems described by quadratic fermionic and bosonic Lindblad master equations [44,45,46,47].The generalization of the quasiparticle picture for interacting integrable models was conjectured in Ref. [28] for quenches producing entangled pairs only.Interestingly, while a quasiparticle picture works for the von Neumann entropy [28,48,49], it has been shown in Ref. [50] that in the presence of interactions it fails for the Rényi entropies, although a hydrodynamic description of the growth of the Rényi entropies is possible.The steady-state Rényi entropies have been obtained in [51,52,53] also for interacting systems.Very recently, the validity of the quasiparticle picture for two-dimensional free systems has been explored in Refs.[54,55].Remarkably, the quasiparticle picture, at least for free models, can be extended to quantum quenches producing entangled multiplets, i.e., going beyond the entangled pair paradigm [31,32,36].
Here we show that the quasiparticle picture for entangled multiplets [31,32] provides the natural framework to describe the dynamics of the charged moments in free-fermion systems.Notice that although some elements of the quasiparticle picture were presented already in Refs.[1,5,9], a full-fledged quasiparticle picture was still lacking.Specifically, even in the case of quenches producing entangled pairs, a framework to determine the contribution of the quasiparticles to the charged moments and the asymmetry for generic initial states was not completely specified.Therefore, we conjecture a formula for ln(Z n (α, t)), which allows to describe their evolution in the hydrodynamic limit.We describe our framework for the non-dissipative case, showing that it is consistent with the results of Sec.4.1 (and also with the results of [2]).For the dissipative case, we do not have a well-established framework, even though the results of Ref. [43,44,45] suggest that it should be possible to generalize our approach to that case.

Entangled multiplets and von Neumann entropy
Let us briefly review the quasiparticle picture for the von Neumann entropy S A (ρ A ) in the presence of entangled multiparticle excitations.As anticipated, this generalized quasiparticle picture explains the spreading of quantum correlations as a consequence of the propagation of entangled multiplets of quasiparticles.Let us assume that the Hamiltonian of the system can be recast in the form (12), with multiplets formed by m quasiparticles.Let us assume that in the basis of the post-quench fermionic modes d j (k) (cf ( 12)), the two-point correlation function on the initial state is block-diagonal in the quasimomentum space, i.e., if we define The correlator G(k) := 2 G(k) − 1 contains the correlations among the various species of particles constituting the multiplet with quasimomentum k, and it can be usefully decomposed into m × m blocks as where . .m labelling the different species, and Let us now focus on the dynamics of the entangled multiplets.One has to be careful in identifying the propagation velocity, that in general is not v j (k) = dε j (k)/dk, with k the quasimomentum in the folded Brillouin zone B = [0, 2π/m].This is due to the fact that the quasimomentum k ∈ [0, 2π/m] of the quasiparticles of species j is a nontrivial function k(k ′ ) of the quasimomentum k ′ ∈ [0, 2π] before the folding.In general, k is related to k ′ by a shift and a sign change.On the other hand, the velocity of the quasiparticles has to be obtained from the dispersion ε(k ′ ) before the folding, because ε(k ′ ) is the dispersion of the physical excitations of the system.Thus, we have v j (k) = dε j (k ′ )/dk ′ = dε j (k)/dk • dk/dk ′ , where dk/dk ′ accounts for the minus sign originating from the folding.
Finally, the contribution s A (x, t, k) to the von Neumann entropy of a subsystem A at time t originating from a multiplet produced initially at position x is conjectured to be [30,31] where G A(x,t) is the restriction of the matrix G(k) (cf.(62)) to the indices describing the subset A(x, t) of quasiparticles that are contained in the subsystem A at time t.The Rényi entropies can be found in a similar way.To understand Eq. ( 63), it should be noticed that the correlator in quasimomentum-space (cf.(61)) and the real-space correlator (Γ + 1)/2, with Γ defined in (15), have the same spectrum, because the change of basis to go to momentum space is unitary (the canonical anticommutation rules are preserved).Clearly, this does not hold for the spectrum of the correlator restricted to a subsystem A, because the bipartition breaks the translation symmetry.This renders the computation of entanglement-related quantities challenging.However, one can imagine of dividing the chain into "mesoscopic" cells of lengtt dx, such that a ≪ dx ≪ |A|, where a is the lattice spacing.Each cell is characterized by a density matrix ρ x , which, since the system is in a Gaussian state, is fully described by the two-point function G x , whose spectrum, since dx ≫ a, will be approximately the same as that of full system, i.e., k∈B G(k).Thus, we can write ρ x ≈ k∈B ρ x,k , with the state of the whole chain ρ = x ρ x , and ρ x,k the density matrix identified by G x (k).To connect with the standard quasiparticle picture, we view each cell as a source of entangled multiplets of quasiparticles, the quasiparticles being the single-particle eigenmodes d † 1 (k), ..., d † m (k) of the Hamiltonian.The correlations between them are encoded in G x (k).The multiplets spread because of the ballistic motion of the quasiparticles.Suggestively, we could associate to each multiplet k produced in the cell at x a non-local correlation matrix G x,k,t .This matrix describes the global state of the m cells at positions x + v j (k)t⟩} j=1,...,m .The only nonzero entries of G x,k,t are the same as the entries of G x,k with the caveat that the mode d † j is that of the cell at x + v j (k)t, meaning that the different modes are in different cells.So G x,k,t describes the way in which different cells are correlated with each other by the modes originating from the cell at x. Since the state is Gaussian, one could associated to G x,k,t a density matrix ρ x,k,t describing the global state of the cells.The state of the whole chain is, then, ρ t = x ρ x,t ≈ x k∈B ρ x,k,t .Now the state of a subsystem A is obtained by tracing over the cells that are not in A in each ρ x,k,t , which corresponds to tracing over the modes outside of A in the matrix G x,k,t .Finally, from the state of A one can compute the entanglement entropy, and the result is (63).
Notice also that in the case m = 2 this prescription yields the well-known formula for the case of entangled pairs (see e.g.[27]), as shown in [32].It has been checked in [31] that our procedure gives the correct hydrodynamic limit for the von Neumann and Rényi entropies, and also for the tripartite mutual information in [32].Moreover, this procedure allows to obtain the dynamics of the von Neumann entropy starting from inhomogeneous initial states [30].

Entangled multiplets and charged moments
Having explained the quasiparticle picture for the von Neumann entropy in the presence of entangled multiplets, we now show how it can be adapted to describe the charged moments (cf.Eq. ( 23)).Inspired by Eq. ( 24), since ln Z n (α, t) is a function of the real-space correlator Γ restricted to the considered subsystem, we can generalize the argument leading to (63).We conjecture that the leading order behaviour of ln Z n (α, t), in the hydrodynamic limit and for any subsystem A, is obtained as a sum of the contributions z n (x, t, k) from each multiplet, reading where G A(x,t) is the same as in (63).Here, we defined w A j (x, t, k) as and, importantly, the restriction n A of the charge operator N A (cf. Eq. ( 23)) to the block A has to be written in the same basis as G A , i.e., in the basis of the fermionic modes d j (cf.( 12)).For the models and quenches discussed here, n A remains diagonal as n A = diag(−1, . . ., −1, 1, . . ., 1), although this is not true in general.We now show that (64) yields Eqs. ( 30) with ( 45) and ( 46) for the XX chain, and it correctly reproduces known results for the dynamics of the asymmetry after quenches in the XX chain [2,9].

5.2.1
First example: quench from the tilted Néel state in the XX chain.
As a first example of the correctness of (64) here we discuss the quench from the tilted Néel state in Eq. ( 13) with the XX Hamiltonian (10).To apply the multiplet framework, we first define the fermionic operators associated with the different species d with The number of quasiparticles forming the entangled multiplet is dictated by the structure of the initial state.For the tilted Néel state the quadruplet structure arises from the fact that the expectation values of both the correlators c † x c y and c x c y are nonzero, and because of the two-site periodicity of the initial state.According to our discussion in section 5.2, the velocities of the four species are v k)/dk and v 4 (k) = −dε 4 (k)/dk because for the species with j = 2, one has k ′ = π − k and for j = 4, one has k ′ = 2π − k.The two-point fermionic correlation functions calculated over the tilted Néel state are block-diagonal in terms of d j (k), which is the reason for introducing them.Indeed, the only nonzero correlators in the k-space are and The functions f ij and g ij appearing in the correlators in (67), ( 68) and ( 69), (70) are the same as in Eq. ( 29).Notice the rescaling k → 2k as compared with (29), which reflects that the Brillouin zone is halved for the tilted Néel state.Again, only c k , c π±k , c 2π−k appear in (67) (68), which makes apparent the reason for our definitions of the different species.We can now write down G(k) (cf.( 62)) in block diagonal form.The blocks C(k) and and We stress that C(k) and F (k) are different from C(k) and F(k) in Eqs. ( 27) and ( 28), as can be seen from the definition of G(k) based on Eq. (61) (see [2]).Clearly, C(k) and F (k) are not diagonal, signaling nontrivial correlation and entanglement between the different species of quasiparticles.To obtain the quasiparticle picture for the charged moments one has to determine the kinematics of the quasiparticles.Since we observed that v 1 = v 2 and v 3 = v 4 , the contribution of the quadruplet to the charged moment can originate from the configurations in which two quasiparticles forming the multiplet are in A (either 1, 2 or 3, 4).The other possibility is that all four quasiparticles are in A (indeed, the configuration with no quasiparticles in A immediately yields zero).We anticipate that this configuration gives a non-trivial contribution, differently from the quasiparticle picture for the entanglement entropy, for which only quadruplets that are shared between the subsystem and the rest contribute to the entanglement between them.
The fact that the quadruplets d 1 , d 2 , d 3 , d 4 are formed by two pairs of quasiparticles with the same velocity, prevents the symmetry restoration by contributing non-trivially to the charged moments and to the integrals in Eq. (107) at any time t.Physically, the fact that v 1 = v 2 and v 3 = v 4 implies that even at infinite time there are always pairs of "correlated" fermions in A, which prevent the Cooper pairs correlator from vanishing (see also [2] for an explanation of the lack of restoration of the U (1) symmetry in terms of a non-Abelian generalised Gibbs ensemble).Let us now observe that the total number of multiplets with the pair 1 − 2 or 3 − 4 of quasiparticles in A are min(ℓ, 2|v 1 (k)|t) and min(ℓ, 2|v 3 (k)|), respectively.On the other hand, the number of multiplets with all the four quasiparticles inside A is max (ℓ − 2|v 1 (k)|t, 0).
Putting everything together and using (64), we find the predicted leading order behaviour of the charged moments as where z (4)  n (k) is the multiplet contribution (64) for all the four quasiparticles inside the subsystem, while z {1,2} n (k), z {3,4} n (k) are the contributions for the case with quasiparticles 1, 2 and 3, 4 in A.
We now explicitly show that Eq. ( 73) coincides with the quasiparticle picture (30) for the charged moments reported in Sec.4.1.Let us give a closer look at Eq. ( 73).First, z (4) is given as (cf.(64)) where G(k) is defined in (62), and C(k) and F (k) are given in ( 71) and (72).We can reorder rows and columns of the 8 × 8 matrix G(k) into a block-diagonal form with 4 × 4 blocks, such that Eq. (74) becomes with and are defined in (65), where now n A is a 4 × 4 matrix.One can check that G (1)  (38)).This allows us to obtain (see the proof of Eq. ( 45) in Appendix A) The contribution z {1,2} n corresponding to two of the quasiparticles forming the quadruplet being in A is, instead, with the restricted two-point function G {1,2} (k) reading The matrix G {1,2} is obtained from (62) by selecting rows and columns corresponding to the quasiparticles in A. The contribution of the configuration with quasiparticles 3, 4 in A reads As for z (4)  n , it is easy to show that G {1,2} (k) is isospectral to G ∞ (2k) (cf.(38)) and that, in the eigenbasis of G {1,2} (k), we can rewrite n A = σ x ⊗ 1.Thus, by using the same steps in the proof of Eq. ( 46) in Appendix A, we obtain If we finally substitute Eqs. ( 77) and (80) in Eq. ( 73), and we change the integration variable as k ′ = 2k, we recover Eq. ( 30), proving the equivalence of ( 73) and (30).

Second example: Quench from XY spin chain
Here we consider the quench in the XX chain starting from the ground state of the XY chain.
We employ (30) to recover the quasiparticle picture for the dynamics of the charged moments, which was derived in [9].To employ our generalized quasiparticle picture, we first observe that the quench produces only pairs of entangled quasiparticles.Let us define 6)).This yields opposite velocities The two-point fermionic correlation function on the initial state, i.e., the ground state |Ω⟩ of the XY chain, is block diagonal with respect to these two species of quasiparticles.Specifically, the only two nonzero correlators in the quasimomentum space are Thus, the correlation matrix G(k) (cf.( 62)) reads with where the row and column indices of the matrices G(k) and F (k) identify the two species.
The leading order of the charged moments in the hydrodynamic limit for an interval A of size ℓ was obtained in Ref. [9] as where We now rederive (84) by using the general formula (64).Since the quench produces only entangled pairs, one has configurations with no quasiparticles in A, or one or two quasiparticles in A. If no quasiparticles are present in A, the contribution to the charged moments is trivially zero.If one or two quasiparticles are in A there are nonzero contributions that we denote as z {1} n (k, α) and z {2} n (k, α), respectively.The number of configurations corresponding to one or both quasiparticles in A is min(ℓ, 2|v 1 (k)|t) and max(ℓ − 2|v 1 (k)|t), respectively.Notice that they are the same as in (73).In conclusion, we obtain Let us now consider the matrix G {1} = G {2} , which is obtained by restricting the row and column indices to 1 (or 2) in (83) and using (82).It reads as Since G {1} and G {2} commute with the matrix n A = −σ z , the one-particle contributions z {2} n and z {2} n are independent from α.It is straightforward to obtain that where h n (x) is defined in (48) and n(k) in (86).On the other hand, the contribution z {1,2} n due to both members of the pair being in A can be derived following the same steps to obtain Eq. ( 45) (see Appendix A).The result reads as where f k is defined in (85).We report the derivation of (90) in Appendix B. Finally, after substituting Eqs. ( 89) and (90) in Eq. (87), and using that the integrand in (84) is invariant under k → 2π − k, we obtain formula (84).The results derived above allow us to draw some general conclusions on the dynamics of the charged moments and the asymmetry after generic quenches producing entangled pairs of quasiparticles.As anticipated, in contrast with the von Neumann and Rényi entropies, when an entangled pair is shared between the subsystem A and its complement A, typically it does not contribute to the entanglement asymmetry ∆S (n) A .For the quench from the ground state of the XY chain this is clear because the asymmetry depends on the ratio (Z n (α)/Z n (0)) (cf.( 22)).Now, since (89) does not depend on α, Eq. ( 89) gives no contribution to the asymmetry.Interestingly, this happens for any free-fermion quench producing entangled pairs travelling with opposite velocities and such that the restricted 2 × 2 correlation matrix G {j} commutes with the charge operator n A .Since both G {j} and n A are 2 × 2 matrices, it can be easily shown that the commutation condition is equivalent to n A being diagonal; it ensures that z {j} n (k, α) = ln h n (n j (k)), where n j (k) is the density of quasiparticles of type j produced after the quench.This implies that the only nontrivial contribution to the entanglement asymmetry originates from pairs fully contained in the subsystem.At infinite time there are no pairs with both particles contained in any finite subsystem A, implying that for quenches producing entangled pairs, the symmetry is always asymptotically restored, provided that [G {j} , n A ] = 0. Notice that this reasoning can be straightforwardly generalized to quenches producing larger multiplets, if all the species of the multiplet propagate with different velocities: also in these cases the symmetry is asymptotically restored.

Dissipation and the quantum Mpemba effect
Here we discuss the effect of dissipation on the symmetry restoration and the quantum Mpemba effect.We focus on the quench from the tilted Néel state (cf.( 14)) in the XX chain in the presence of gain and loss dissipation.Crucially, in the absence of dissipation, the U (1) symmetry of the XX chain is not restored.This can be interpreted as the consequence of non-Abelian conservation laws [13] of the Hamiltonian [2].Here we show that in the presence of dissipation the U (1) symmetry of the XX chain is restored.Concomitantly, we show that quantum Mpemba effect is robust against dissipation.
Let us first remark that by using Eq. ( 51) for the charged moments in the weaklydissipative hydrodynamic limit and the definitions ( 20) and (17), it is straightforward to compute, at least numerically, the entanglement asymmetry for generic dissipation.
In the following, we first analytically derive the leading order of the entanglement asymmetry ∆S (n) A .Specifically, for balanced gain and losses we provide the formula for ∆S (n) A in the limit t → ∞ for any n ̸ = 1.For γ + ̸ = γ − we provide the result for ∆S (2) A in the limit t → ∞.Moreover, by comparing the short-and long-time behaviour of the asymmetry, we obtain a condition for the Mpemba effect in the presence of gain/loss dissipation.In Section 6.1 we compare this condition with the case without dissipation.
When t = 0, i.e., before the quench, the entanglement asymmetry is the same as in the case without dissipation, which was already derived in [2] for large subsystem size ℓ.The result reads as [2] with f ij defined in (29).To determine the behaviour of the asymmetry in the limit t → ∞, let us start from the ratio of the charged moments in the hydrodynamic limit As it is clear from the definitions (18) (19) (20) (21), Eq. (92) allows one to obtain the entanglement asymmetry.In the large-time limit t/ℓ → ∞, only the term in (92) multiplying min(1, 2|v(k)|t/ℓ) survives.Thus, by plugging the result in Eq. ( 20), we find that the entanglement asymmetry is To proceed, let us focus on the case with n = 2.If we plug Eq. (60) in Eq. ( 93) and we expand at the leading order in the small parameter λ(t), we find where α = α 1 − α 2 , and s = (γ + − γ − )/(γ + + γ − ).If we define the remaining integration in α yields

∆S
(2) where I 0 is the modified Bessel function.Finally, we obtain Let us now consider the case γ + = γ − .We can evaluate the leading order behaviour of the Rényi entanglement asymmetry for any n ̸ = 1, by replacing Eq. ( 57) into Eq.( 93) and expanding at the leading order in λ(t), such that we obtain We stress that the result above holds for n ̸ = 1 because we are not able to analytically continue in n Eqs. ( 57) and ( 49).Thus, we cannot prove that the limits n → 1 and t → ∞ commute (see for example [9]).Using that λ(t) = e −2γ + t → 0 in the large time limit, we can expand the exponential in the integral above.The resulting n-fold integral in α can be computed by observing that dxdy sin 2 (x − y) = 2π 2 and, since the sum over p j 's involves n 2 terms, we get This allows us rewrite Eq. (98) as As we anticipated, our goal is to determine the condition under which the Mpemba effect occurs.Given two initial tilted Néel states with tilting angles θ 1 > θ 2 , in formulas, this condition reads For all the cases in which we were able to explicitly compute the entanglement asymmetry, by using Eqs.(91), ( 97) and (100), the condition (101) becomes We notice that, for θ ∈ [0, π/2], the first condition is always satisfied as far as θ 1 > θ 2 .One expects the criterion (102) to hold also for generic gain/loss processes, although we were not able to prove it.In the following, we compare the condition (102) with the one obtained in other (non-dissipative) free-fermion quench protocols [9].

Multiplet picture and onset of Mpemba: Comparison with the nondissipative case
Here we show that the generalized quasiparticle picture for the charged moments introduced in section 5 allows to understand the condition (102) giving rise to the Mpemba effect both in the non-dissipative case, as well as in the presence of dissipation.
We can now revisit the condition (102) under which the Mpemba effect occurs.We start discussing the condition for Mpemba in the absence of dissipation.We first consider the quench from the ground state of the XY spin chain in the XX chain.Then we focus on the quench in the XX chain starting from the tilted Néel state.For the latter (cf.( 13)) in the absence of dissipation the symmetry is not asymptotically restored (unless θ = π/2 [2]).Therefore, quantum Mpemba effect does not occur.Nevertheless, even though the symmetry is not restored, a crossing between the curves for the asymmetry corresponding to different tilting angles can occur.We can still write the condition (101) which ensures that there is a crossing in time between ∆S (n) The main goal here is to show that the condition for Mpemba to happen in the presence of dissipation is similar to the crossing condition in the absence of symmetry restoration.Specifically, in these conditions all the quasimomenta contribute, in contrast with the condition for Mpemba found in [5,9].
Let us first discuss the condition for the standard Mpemba effect, focussing on the quench from the ground state of the XY chain (see section 2 and 5).The post-quench dynamics is driven by the XX Hamiltonian.The condition to observe the Mpemba effect discussed in Ref. [9] can be interpreted in terms of multiplets.For the quench from the ground state of the XY chain only pairs (k, 2π − k) are produced.First, at time t = 0, the leading order contribution to the asymmetry is [9] where ∆ k (γ, h) is the Bogoliubov angle in Eq. ( 5), while at long time it becomes Thus, the quantum Mpemba effect occurs if there exists a time, t I , such that where The first inequality does not depend on whether the dynamics is unitary or dissipative.On the other hand, the second condition in (105) involves only the modes in the interval The reason is that unlike the quench from the tilted Néel only pairs of quasiparticles traveling with opposite velocity are present.This means that at long times the asymmetry vanishes because the number of pairs in A vanishes.We now discuss the quench from the tilted Néel state.Now, neither restoration of symmetry nor Mpemba effect occur in the absence of dissipation, although curves for different values of the tilting angle can cross at intermediate time.To write the crossing condition, let us notice that the leading order of the Rényi entanglement asymmetry at long time reads [2] where p n is, in general, a cumbersome function of n ± (k, θ) (cf.(47)), that also depends on the Rényi index n.For example, p 2 (n + , n − ) = h 2 (n + )h 2 (n − ) with h n defined in (48).For the following we do not need the explicit expression for p n (see Ref. [2]).By using ( 106) and (91), the crossing condition in Eq. ( 101) becomes . (107) The conditions in (107) are quite similar to the ones in the presence of dissipation (cf.( 102)), in contrast with the standard Mpemba conditions (105).The first condition is the same, since at t = 0 the dissipation has no effect, whereas the second one differs for the denominator p n in the integrand, which is absent in the dissipative case.The generalized quasiparticle picture discussed in Section 5 allows to interpret the crossing condition in Eq. (107).To establish the condition (107) one has to compare the asymmetry in the initial state with that at long times.
In the initial state, the asymmetry is proportional to the correlators ⟨c k c 2π−k ⟩ and ⟨c k c π−k ⟩.Again, the former correlator is due to the members of the quadruplet traveling with opposite velocity.The latter one corresponds to the quasiparticles traveling at the same velocity, which are responsible for the lack of symmetry restoration at long times.The integral in the first inequality in Eq. ( 107) is the sum of the contributions coming from ⟨c k c π−k ⟩ and ⟨c k c 2π−k ⟩ as it is clear from Eq. (69).At the initial time all the quasiparticles in the quadruplet contribute to the asymmetry.On the other hand, the integral in the second inequality of Eq. ( 107) involves only ⟨c k c π−k ⟩ because at long times, only the correlator ⟨c k c π−k ⟩ appears.The correlator ⟨c k c 2π−k ⟩ vanishes at long times because the associated quasiparticles have opposite velocities, and cannot be both in A at the same time.The integrand in the second condition in (107) can be somehow identified as the "asymmetry content" of the pair (k, π −k).
Notice that the normalization p n (n + , n − ) depends on both the correlator ⟨c k c 2π−k ⟩ and the correlator ⟨c † k c k ⟩, via the functions f 11 (2k, θ) and g 12 (2k, θ) in Eq. (79).Summarizing, the crossing in the entanglement asymmetry occurs if initial states with larger asymmetry lead to pairs c k c π−k with smaller asymmetry.
The interpretation of condition (102) for the quench in the presence of dissipation remains qualitatively the same.First, dissipation suppresses the correlations carried by the quadruplets at a rate independent of k.This means that the correlator ⟨c k c π−k ⟩ vanishes at long times, implying symmetry restoration and the possibility of having Mpemba effect.Specifically, the condition to observe the Mpemba effect is, again, that the initial state with larger total asymmetry leads to pairs (k, π − k) with smaller asymmetry content, that at long times are destroyed by the dissipation at the same rate.Indeed, Eq. (107) and Eq. ( 102) are the same, except for the factor, p n (n + (k, θ), n − (k, θ)), which is dropped in Eq. ( 102).The reason is that it is subleading in λ(t).Notice that in the presence of dissipation, the entanglement asymmetry ∆S (n) A (θ, t) decays exponentially with time with a rate 2(γ + + γ − ), as can be seen from Eq. (100).This is different in the unitary integrable cases studied so far, where the asymmetry displays an algebraic decay in time (see Ref. [5]).On the other hand, the classical Mpemba effect, which occurs, for instance, in colloidal systems [56], exhibits a similar exponential behaviour.Finally, another similarity with the classical Mpemba effect is that in the presence of dissipation, the steady state arising at long times does not depend on the initial state, in contrast with the quantum Mpemba effect.

Numerical results
In this section we provide numerical evidence supporting our results.We focus on the quench from the tilted Néel state in the XX chain because, as it was discussed, in the absence of dissipation the U (1) symmetry of the XX chain, which is broken by the initial state, is not restored at long times.Our results confirm that in the presence of dissipation the symmetry can be restored, and Mpemba effect can occur.The section is structured as follows.In subsection 7.1 we benchmark our results for the charged moments in the presence of gain and loss dissipation.In section 7.2 we focus on the Mpemba effect discussing the Rényi entanglement asymmetry.

Charged moments
In Fig. 2 we show numerical results for the charged moments ln(Z n (α, t)/Z n (0, t))/ℓ for n = 2 and n = 3 (left and right panel, respectively).The charged moments are plotted against the rescaled time t/ℓ, with ℓ the size of subsystem A. We show results in the weakly-dissipative hydrodynamic limit t, ℓ → ∞ with the ratio t/ℓ fixed and γ ± → 0 with fixed γ + ℓ, γ − ℓ.At fixed γ ± , in the limit t → ∞ the charged moments exhibit an exponential relaxation to the steady state.In Fig. 2 we consider several tilting angles θ, and values of α ij := α i − α j .The symbols in Fig. 2 are exact lattice results.The data are obtained by applying directly Eq. ( 24) to the real-space correlator Γ(t) computed using the symbol in Eq. ( 50).The charged moments are always negative.At θ = 0 they vanish, whereas upon increasing θ they become larger in absolute value.The continuous lines in the Figure are obtained by using the quasiparticle picture (cf.( 51)).The agreement with the numerics is perfect.

Entanglement asymmetry
We now discuss the Rényi entanglement asymmetry.The entanglement asymmetry is obtained by plugging Eq. ( 51) in (20) and performing the Fourier transform numerically.The analytic prediction obtained by using the quasiparticle picture is reported with the solid lines in Fig. 3.We plot it for different tilting angles θ, dissipation parameters γ + , γ − and for n = 2, 3. Specifically, the left panel shows results for n = 2 and γ + = 4γ − .Fig. 3 (right panel) plots results for n = 3 and balanced gain and losses, i.e., for γ + = γ − .For n = 2 we can also take advantage of the general analytical expression for the charged moments given in Eqs. ( 59) and (60).The symbols in the plots are exact numerical data.The agreement with the quasi-particle prediction is remarkable.
Fig. 3 shows that there exist states that have larger asymmetry at time t = 0 but faster restoration during the time evolution (for instance, the blue and red lines in the left panel).This is a clear signature of the quantum Mpemba effect.However, this behaviour is not generic, meaning that having a larger tilting angle θ ∈ [0, π/2] does not automatically imply Mpemba.This is in contrast with what happens for the quench from the tilted ferromagnetic state [1].The scenario depends dramatically on the dissipation.For instance, for balanced gain and losses (right panel in Fig. 3) the Mpemba effect seems stronger.In the insets in Fig. 3 we show the condition (107) for Mpemba effect.We plot f 2 11 (k, θ) as a function of k.Clearly, when Mpemba occurs we observe that larger values of θ, i.e., stronger symmetry breaking in the initial state give rise to multiplets with smaller asymmetry content, and smaller integral dkf 2 11 (k, θ).This confirms the validity of (107).We finally check the asymptotic expansion of the entanglement asymmetry in Eq. (100).In Fig. 4 we show the exact quasiparticle picture obtained by performing the Fourier transform of Eq. ( 51) plugged in Eqs. ( 20) and ( 17) (solid lines), and we prove that they exponentially decay as ∼ e −2t(γ + +γ − ) for n = 2, 3 (dashed black lines).

Conclusions
We investigated the quantum Mpemba effect in one-dimensional fermionic chains subject to dissipation.We employed the framework of Lindblad master equation [35].First, by using the multidimensional stationary phase approximation we provided a fully rigorous derivation of quantum Mpemba effect for several quenches in free-fermion models in the absence of dissipation (see Ref. [1] and [2]).Based on these results we conjectured a formula for the charged moments in the hydrodynamic limit of long times and large intervals with their ratio fixed.Our formula exploits the fact that the generic quench produces entangled multiplets of excitations, rather than entangled pairs.Our results suggest that the framework of entangled multiplets is the natural one to formulate the quasiparticle picture interpretation for the dynamics of the entanglement asymmetry, and provides valuable insights also for quenches producing entangled pairs.Then, we focused on the effects of gain and loss dissipation on the dynamics of the entanglement asymmetry.Specifically, we considered the quench from the tilted Néel state in the XX chain.In the absence of dissipation the U (1) symmetry of the model is not restored by the dynamics, and the quantum Mpemba effect does not occur.This is due to the presence of peculiar non-Abelian conservation laws of the XX chain [2,13].
Here we show that gain/loss dissipation can restore the symmetry and the quantum Mpemba effect can occur.Interestingly, the condition that explains the effect is different with respect to the one studied in [5,9].The reason can be traced back to the presence of entangled quadruplets of quasiparticles and not just entangled pairs.The correlator ⟨c k c π−k ⟩ vanishes in the presence of dissipation, which explains why the symmetry is now restored.However, at long times the correlator gets contributions only from quasiparticles traveling at the same velocity, since they can be both in A at the same time.This is in contrast with quenches producing entangled pairs [9].As these paired quasiparticles with the same velocity exist for any k, all the modes contribute to the asymmetry at large time (cf. the second line of (102)), not just the "slow" ones, in contrast with the standard Mpemba effect [9].Let us now discuss some possible future directions.First, the quasiparticle picture that we proposed in Sec. 5 works for the charged moments.However, our results suggest that it may be possible to extend the approach to obtain directly the entanglement asymmetry.To do that, one has to first determine the relationship between the integrands in Eqs. ( 91)-(100) and the fermionic correlation matrix in momentum space (cf.(64)).It would be also interesting to generalize our results for the charged moments and for the entanglement asymmetry to other types of Lindblad master equations, by exploiting the results of Ref. [45].This would allow to better clarify the conditions under which "dissipative" quantum Mpemba effect occurs.Furthermore, our results showed that the Mpemba effect, or the absence thereof, is intertwined with the production of nontrivial entangled multiplets after the quench.It would be interesting to explore the effect of the size of the multiplets and of dissipation on the asymmetry, and on the quantum Mpemba effect.To this purpose, one could exploit the results of Ref. [32].Clearly, it would be interesting to go beyond the case of quadratic Lindblad equations, to study the effect of interactions.While it is challenging to determine the entanglement asymmetry, by exploiting integrability it should be possible to study the effect of interactions and dissipation on the quantum Mpemba effect at the level of correlation functions [57].

A.1 Considerations on the dissipative case
In the presence of gain and loss dissipation, it is challenging to evaluate the two determinants in (51), in contrast with the non dissipative case (see Appendix A).The reason is that the Fourier transform G ′ (k) of the fermionic correlator, defined in (50), loses some crucial properties that we used in the non-dissipative case.Indeed, it is no longer true that (G ′ 0 (k, θ, t)) 2 = 1 (cf.( 53)), which we heavily exploited to calculate the short-time term ln(det Z (α) n,0 (k, t)).Moreover, the generalization of (119) exhibits a block-diagonal structure in the eigenbasis of G ′ ∞ , although it is more complicated (cf.( 136)).In the following we show that for generic dissipation rates γ ± we are able to determine analytically ln(det Z ′ (α) n,0 ) and ln(det Z ′ (α) n,∞ ) only for small values of n.For balanced gain and losses, i.e., γ + = γ − we provide analytic results for ln(det Z ′ (α) n,∞ ) for arbitrary n, whereas we provide a closed-form expression for ln(det Z ′ (α) n,0 ) only for n = 2.
By using (140) and ( 135) in ( 51) we obtain the leading order of ln Z 2 (α, t) in the weaklydissipative hydrodynamic limit.
B Proof of Eq. (90) Here we derive (90).We consider the XX chain (see section 2).Th pre-quench initial state is the ground state of the XY chain.As explained in section 2, only entangled pairs of quasiparticles are produced after the quench.Eq. ( 90) is the contribution to the charged moments ln(Z (α) n ) due to the configuration in which both the quasiparticles forming the pair are in subsystem A. To prove (90) one has to compute where W {1,2} (k) := n j=1 (1+G {1,2} (k)))(1−G {1,2} (k)) −1 e i(α j −α j+1 )n A , with G {1,2} (k) defined in (82) and n A = diag(−1, −1, 1, 1).The subscript {1, 2} means that one has to restrict the correlator to the species of quasiparticles that are in A. Notice that since there are only two species of quasiparticles one has that G {1,2} (k) = G(k), i.e., the restricted correlator is the full correlator.Now, we can adapt the strategy used in Appendix A, obtaining the equivalent of Eq. (111) as where W {1,2} (k) := n j=1 (1 + G(k))e i(α j −α j+1 )n A .Again, we proceed by Taylor-expanding the exponentials e i(α j −α j+1 )n A .Thus, one defines the moments M Eq. ( 143) can be rewritten as To perform the trace in (144) we diagonalize n A (1 + G(k)).The only doubly-degenerate eigenvalue is 2 cos ∆ k .Thus, we obtain We then sum over h j , reconstructing the exponentials of n A .Each term yields a factor cos α j,j+1 + i cos ∆ k sin α j,j+1 = f k (α j,j+1 ).We obtain Finally, by substituting (146) and ( 147) into (142), and taking the (halved) logarithm, we obtain the desired Eq. (90).

Figure 1 :
Figure 1: Setup considered in this work.(a) Quench from the tilted Néel state, with θ the tilting angle, in the XX chain in the presence of global gain and loss dissipation.We are interested in the dynamics of the entanglement asymmetry of a finite region A embedded in an infinite chain.(b) Quasiparticle picture for the asymmetry: the case of entangled pairs.If entangled pairs of quasiparticles are produced after the quench, only the pairs that are fully contained in the subsystem contribute.Here the pair contributes for t ≤ t 1 .At longer times the pair is not fully contained in A, and the asymmetry vanishes.(c) For the tilted Néel state quadruplets of entangled quasiparticles are produced with velocities v 1 = v 2 and v 3 = v 4 .In contrast with (b), in the absence of dissipation the asymmetry does not vanish at long times because there are always paired fermions in A.

Figure 3 :
Figure 3: Rényi entanglement asymmetry after the quench from the tilted Néel state in the XX chain as a function of time t for a subsystem size ℓ = 80, dissipation γ + = 4γ − = 1/200 (left panel) γ + = γ − = 1/200 (right panel) and for different initial tilting angle θ.The curves have been obtained by performing the numerical Fourier transform of Eq. (51), while the symbols are the exact numerical values.The inset shows f 2 11 (k, θ) as a function of quasimomentum k.The vanishing behaviour at long times and the crossing at short one signal the quantum Mpemba effect.
Thus, the only term to evaluate in (121) is Tr[W (↑) ∞ (k)].Again, it is useful to first calculate the moments M