Symmetry-resolved entanglement in fermionic systems with dissipation

We investigate symmetry-resolved entanglement in out-of-equilibrium fermionic systems subject to gain and loss dissipation, which preserves the block-diagonal structure of the reduced density matrix. We derive a hydrodynamic description of the dynamics of several entanglement-related quantities, such as the symmetry-resolved von Neumann entropy and the charge-imbalance-resolved fermionic negativity. We show that all these quantities admit a hydrodynamic description in terms of entangled quasiparticles. While the entropy is dominated by dissipative processes, the resolved negativity is sensitive to the presence of entangled quasiparticles, and it shows the typical ‘rise and fall’ dynamics. Our results hold in the weak-dissipative hydrodynamic limit of large intervals, long times and weak dissipation rates.


Introduction
The study of entanglement dynamics is of fundamental importance to understand the process of thermalisation or the lack thereof in out-of-equilibrium quantum many-body systems [1][2][3][4][5].In one-dimensional integrable systems the entanglement dynamics after a quantum quench is amenable of a hydrodynamic description within the framework of the famous quasiparticle picture [6][7][8][9][10][11][12].Specifically, the entanglement dynamics is determined by the ballistic propagation of entangled pairs of quasiparticles.Importantly, recent progress with cold-atom and trapped-ion systems allows us to experimentally measure entanglement-related quantities [13][14][15][16] and their dynamics.
The starting point for all entanglement-related quantities is the reduced density matrix ρ A of a subsystem A. For example, the entanglement Rényi entropies are defined as The limit n → 1 gives the von Neumann entropy as The Rényi entropies and the von Neumann entropy are proper entanglement measures for bipartite pure states only.If the system is in a bipartite mixed state, neither the von Neumann entropy nor the mutual information are bona fide entanglement measures.In these situations, for fermionic systems, one can use the logarithmic negativity or the fermionic logarithmic negativity (see section 4 for precise definitions).
Still, most of the effort so far focussed on closed quantum many-body systems.Since in generic systems the interaction with an environment is unavoidable, it is crucial to understand the behaviour of entanglement and symmetry-resolved entanglement in out-of-equilibrium open quantum many-body systems.This is in general a challenging task, and results are available only in few situations, for instance, within the simplified framework of the Lindblad master equations [82].For quadratic Lindblad master equations, i.e., for free-fermion and free-boson models with linear Lindblad operators, the out-of-equilibrium dynamics of the Rényi entropies, and of the mutual information can be understood within the quasiparticle picture [83][84][85].This result holds in the weak-dissipative hydrodynamic limit of long times and large subsystems, and weak dissipation.Moreover, some results are available for the fermionic logarithmic negativity as well.For instance, the dynamics of the negativity in free-fermion systems subject to gain and loss dissipation was investigated in Ref. [86].Interestingly, due to the mixedness of the global state, the negativity is not related to the Rényi mutual information, in contrast with the unitary case [87].The dynamics of both the von Neumann entropy and the negativity in the presence of localised losses has been studied recently [88,89].Relatedly, it has been shown that in free-fermion models in contact with localised thermal baths, which are treated within the Lindblad formalism, the mutual information exhibits a logarithmic scaling [90].Finally, symmetry-resolved entanglement in open quantum many-body systems has received comparatively little attention, although some numerical and perturbative results are available [21].
Here we build a quasiparticle picture for several symmetry-resolved entanglement measures in free fermions subject to gain and loss dissipation.Such dissipation shall violate the particle number conservation, i.e. the U (1) symmetry.However, as we will discuss below, the local Lindblad jump operators, which define the dissipative mechanism, preserve the blockdiagonal structure of the reduced density matrix.We can refer to this scenario as a weak U (1) symmetry [21,[91][92][93][94][95].Specifically, we provide analytic results for both the symmetry-resolved von Neumann entropy, and the charge-imbalance-resolved fermionic negativity.Similar to the non-resolved von Neumann entropy, the symmetry-resolved one is dominated by volume-law terms, which originate from dissipation.These are accompanied by terms that retain information about entangled quasiparticles.Precisely, these terms have the same structure as in the non-dissipative case.However, the entanglement content of the correlated quasiparticles is renormalised by the dissipative processes, and it vanishes exponentially at long times, reflecting how dissipation destroys genuine quantum correlation.While the entropies are dominated by dissipative contributions, we show that the fermionic symmetry-resolved negativity is determined by the entangled quasiparticles.Indeed, the resolved negativity exhibits the typical Figure 1: Schematic view of our setup.A fermionic chain is subject to gain and loss dissipation: Fermions are added or removed incoherently at rates γ ± .The incoherent processes act on each site of the chain independently.The chain is initially prepared in a low-entangled state |Ψ 0 .Here we choose as initial states the fermionic Néel state |N and the Majumdar-Ghosh state |D .The chain evolves under the action of the tight-binding fermionic Hamiltonian and the Lindblad operators.We are interested in the symmetry-resolved entanglement entropy of a subregion A with the rest.We also consider the charge-imbalance-resolved negativity between two equal intervals A 1 and A 2 of length and placed at distance d.
"rise and fall" dynamics, i.e., it grows at short times t 1/γ ± , with γ ± the dissipation rates, and it vanishes at long times.Interestingly, for generic gain/loss processes the symmetryresolved entropies do not exhibit a time delay, meaning that they are different from zero for any t > 0, as already numerically observed in [21].This reflects that different charge sectors are immediately "populated" due to the presence of the dissipation.This represents a crucial difference with respect to the case without dissipation [69].
The manuscript is organised as follows.We start in section 2 reviewing quantum quenches in free-fermion models without dissipation.In section 3 we introduce the dissipation and its treatment within the framework of the Lindblad master equation.In section 4 we introduce the symmetry-resolved von Neumann entropy and the resolved negativity, discussing also their calculation in free-fermion models.Section 5 is devoted to the presentation of our main results.Specifically, in section 5.1 we discuss the charged moments of the reduced density matrix, which are essential ingredients to compute the symmetry-resolved entropy.In section 5.2 we discuss how to obtain the symmetry-resolved entropy and several useful approximations.In section 6 we present our results for the charge-imbalance-resolved negativity.In particular, we discuss the charged moments of the partial time-reversed reduced density matrix, and the resolved negativity.In section 7 we provide numerical benchmarks.Specifically, in section 7.1 and section 7.2 we discuss the charged moments and the symmetry-resolved entropies, respectively.In section 7.3 we address the validity of the quadratic approximation for the charged moments.Finally, in section 7.4 we discuss the scaling in the weak-dissipative hydrodynamic limit.In particular, we highlight the presence of logarithmic corrections, which arise from the charge fluctuations between the subsystem and the rest.These are revealed in the scaling of the number entropy.In section 7.5 we benchmark our results for the resolved negativity.We present our conclusions in section 8.

Quantum quenches in free-fermion models
We consider the fermionic chain defined by the Hamiltonian Here L is the chain size, c † j and c j are canonical fermionic creation and annihilation operators with anticommutation relations {c † j , c l } = δ jl and {c j , c l } = 0.The Hamiltonian (3) is diagonalised by going to Fourier space defining the fermionic operators b k := 1 √ L j e ikj c j , with the quasimomentum k = 2πp/L and p = 0, 1, . . ., L − 1.The Hamiltonian (3) becomes diagonal as Here we defined the single-particle energy levels ε(k).In the following we are considering the thermodynamic limit L → ∞.For later convenience we define the group velocity of the fermionic excitations v(k) as We focus on the nonequilibrium dynamics after the quench from the fermionic Néel state |N and Majumdar-Ghosh (dimer) state |D , defined as with |0 the fermionic vacuum state.
Let us now discuss the quench protocol.At time t = 0 we prepare the system in |N or |D .At t > 0 the chain undergoes unitary dynamics under the Hamiltonian (3).The fermionic correlation function C jl is the central object to address entanglement related quantities in free-fermion systems [96].The matrix C jl is defined as Here • denotes the expectation value calculated with the time-dependent wavefunction.The dynamics of the correlation function C jl (t) (cf.( 8)) after the quench from both the Néel state and the Majumdar-Ghosh state can be derived analytically (see for instance, Ref. [80]).For the Néel state it is straightforward to obtain where J x (t) is the Bessel function of the first kind.In ( 9) is to stress that the dynamics is unitary.
It is useful to exploit the invariance of the Néel state and the Majumdar-Ghosh state under translation by two sites.Thus, we rewrite (8) as where now j, l label the position of a two-site unit cell and t(k) is a two-by-two matrix.The factor 2 in the exponent in the integral in (10) reflects translation invariance by two sites.From (9), one obtains the matrix t(k) for the Néel state as For the Majumdar-Ghosh state we use tD (k) given by [69,97] with where ε(k) is given in (4).
3 Lindblad evolution in the presence of gain and loss dissipation Let us now discuss the out-of-equilibrium dynamics in the tight-binding chain (cf.( 3)) with fermionic gain and loss processes.To describe these dissipative processes, we employ the formalism of quantum master equations [82].Precisely, the Lindblad equation describes the time-dependent density matrix ρ t of the full system as Here, {x, y} is the anticommutator.The dissipation is modelled by L j,α , which are known as Lindblad jump operators.For gain and loss processes they are defined as L j,− := γ − c j and L j,+ := γ + c † j , with γ ± the gain and loss rates.As it is clear from Eq. ( 15), incoherent absorption and emission of fermions happen independently at each site of the chain.This choice of the jump operators, L j,± , falls into the category of weak symmetries, i.e. the generator of the master equation as a whole commutes with the symmetry operation, while each of the jump operator does not commute individually with the symmetry [91].Their main feature is that they preserve the block-diagonal form of the reduced density matrix.Indeed, if the system is prepared in an initial symmetric state (like |N or |D ), the system state remains supported in the same symmetric eigenspace (with a given charge) when no jumps take place.On the other hand, when the first jump L j 1 ,± occurs at time t 1 , the system is transformed to another eigenspace with different charge.If we interpret Eq. ( 15) as a set of quantum trajectories, each one corresponding to the solution of a stochastic Schrödinger equation, the quantum trajectories are symmetric at all times: a single quantum jump only changes the number of particles by an integer number, but the total number of particles is conserved along each trajectory [94].
In (16), Γ ± are L × L matrices defined as Γ ± jl = γ ± δ jl .Since Γ ± jl are diagonal, C jl can be rewritten in terms of the correlation matrix C jl (t) describing the dynamics in the absence of dissipation, i.e., with γ ± = 0 (cf.(8)).Precisely, one has where 1 jl is the identity matrix element.Moreover, Eq. ( 17) suggests that it is convenient to modify the matrix t(k) (cf.( 10)) defining t (k) as where 1 2 is the 2 × 2 identity matrix, and t(k) is the original matrix defined in (10).The correlator C jl is then obtained as Here t (k) is given in (18), with t(k) defined in (11) and ( 12).

Symmetry-resolved entanglement
In this section we overview the definitions of our quantities of interest, namely the symmetryresolved entanglement entropies and the charge-imbalance-resolved negativity.We focus on quantum systems with an internal U (1) symmetry generated by a local operator Q.For the tight-binding chain in (3) this corresponds to the conservation of the number of fermions.The U (1) symmetry implies that [ρ, Q] = 0.If the system is bipartite in two complementary subsystems A and B, since the operator Q is sum of local terms, Q is decomposed as the sum of operators acting separately on the degrees of freedom of A and B, i. e., , implying that the reduced density matrix ρ A is block diagonal.The different blocks correspond to a different charge sector with eigenvalue q of Q A .We have where Π q projects over the eigenspace with eigenvalue q and we defined p(q) = Tr (Π q ρ A ) as the probability that Q A has value q.Here we normalise ρ A (q) such that Trρ A (q) = 1.In the following section we introduce several entanglement-related quantities that can be defined from ρ A (q).This allows to understand how the different charge sectors of the reduced density matrix contribute to the entanglement.We stress again that this block-diagonal structure of ρ A is preserved by the dissipative evolution in Eq. ( 15) within our choice of L j,± [21,91,93].

Symmetry-resolved entropies
Let us introduce the symmetry-resolved Rényi and von Neumann entropy.The latter is the entropy calculated from the different charge sectors of the reduced density matrix as where, again, ρ A (q) is the normalised block of ρ A corresponding to the charge q.To proceed, we plug Eq. ( 20) into the definition of the von Neumann entropy, to obtain the decomposition as where p(q) is the probability of measuring charge q in the subsystem.The first term in (22) quantifies the average of the von Neumann entropy of each charge sector, and it is called configurational entropy [14,48,50].On the other hand, S (num) measures the charge fluctuations between subsystem A and its complement, and it is called number entropy [14,[59][60][61][62][63].
Notice that from ρ A (q) one can define the symmetry-resolved Rényi entropies [46] as Here, however, we only consider the von Neumann entropy.We now discuss how to compute the symmetry-resolved entropy.In general, its computation, for instance in a field-theory setup, requires the spectrum of the symmetry-resolved reduced density matrix.This is challenging to extract because the projector Π q is a nonlocal operator.Here we follow a different strategy, which was put forward in Ref. [18,19], and which we now outline.The central objects of this approach are the charged moments, Z r (α).They are defined as By performing the Fourier transform of the charged moments, we get the moments of the reduced density matrix in a given charge sector, i.e., where Π q , again, is the projector, which selects the block of the reduced density matrix that corresponds to value of the charge q.The symmetry-resolved Rényi entropies S r (q) with Rényi index r are defined as where the symmetry-resolved von Neumann entropy S 1 is obtained as usual by taking the limit r → 1 of S r .Finally, for the free-fermion Hamiltonian (3) one can write the charged moments in terms of the two-point fermionic correlation function as [18] log where C A (cf. Eq. ( 8)) is the fermionic correlation matrix restricted to part A of the chain.

Charge-imbalance-resolved fermionic negativity
Let us now consider the logarithmic negativity.We start by reviewing the definition of the partial transpose and its relation to the time-reversal transformation.Let us consider a density matrix ρ A in which A is further partitioned as A = A 1 ∪ A 2 (see Fig. 1), where A 1 (A 2 ) is an interval of length 1 ( 2 ).We can write where |e 1 j and |e 2 k are orthonormal bases for the Hilbert spaces of A 1 and A 2 , respectively.The partial transpose ρ T 1 A of ρ A with respect to A 1 is defined as The negativity is written in terms of the negative eigenvalues of ρ T 1 A , and it is essentially the trace norm of ρ T 1 A .In terms of its eigenvalues λ i , the trace norm of ρ T 1 A can be rewritten as where in the last equality we used the normalisation i λ i = 1.Moreover, in the absence of negative eigenvalues, Tr|ρ T 1 A | = 1 and the negativity vanishes.For a bosonic system, it is known [98] that the partial transpose is equivalent to the time-reversal partial transpose.This allows one to calculate the negativity in terms of the covariance matrix in bosonic systems, such as the harmonic chain [99].
For fermionic systems, the standard negativity is not easily computed in terms of the two-point correlation function, because the partial transpose is not a Gaussian operator [100,101].On the other hand, the fermionic negativity [102], which is defined from the timereversed partial transpose, can be computed.To introduce the fermionic time-reversed partial transpose it is convenient to introduce Majorana fermion operators w j as with c j the original fermions.The density matrix in the Majorana representation takes the form Here, w 0 x = 1, w 1 x = w x , κ i , τ j ∈ {0, 1}, and κ (τ ) is a 2m 1 -component vector ( 2m2 ) with entries 0, 1 and with norm |κ| = j κ j (|τ | = j τ j ).The constraint on the parity of |κ| + |τ | is due to the fact that we require the density matrix to commute with the total fermionnumber parity operator.In (32) M κ,τ are complex numbers.Using Eq. ( 32), the partial time-reversed (TR) transpose ρ R 1 A with respect to the subsystem A 1 is defined by The matrix ρ R 1 A is not necessarily Hermitian, and may have complex eigenvalues, although ] are all real.Following Refs.[102][103][104] we define the fermionic negativity as One can also define the (fermionic) Rényi logarithmic negativities E r = log(R r ) with R r given as The negativity can be obtained as E = lim re→1 log(R re ), where r e denotes an even integer [102].
We can adapt the approach discussed in section 4.1 to the computation of the charged moments of the time-reversed (TR) partial transpose obtaining A Fourier transform allows us to derive the moments of the TR partial transpose Z R 1 ,r as from which we define the ratios R r (q) and the charge-imbalance-resolved negativity E(q) as Let us stress that the replica limit for R re (q) requires an analytic continuation r e → 1 from the even sequences R re , whereas in the denominator we can set r → 1.
Let us now discuss how to compute the charged moments N r (α) in free-fermion systems.The covariance matrix Γ is the central object to compute the charged moments and the symmetry-resolved negativity.Γ is defined as where 1 is the identity matrix restricted to subsystem A = A 1 ∪ A 2 , and C jl is the fermionic correlation matrix (17).If the subsystem A is made of two intervals as A = A 1 ∪ A 2 (see Fig. 1), we can decompose Γ as where Γ 11 and Γ 22 are the reduced covariance matrices of the two subsystems A 1 and A 2 , respectively, while Γ 12 and Γ † 21 contain the cross correlations between them.By simple Gaussian states manipulations, the covariance matrices Γ ± associated with ρ R 1 A and (ρ R 1 A ) † (cf.(33) for their definitions) are obtained as Here we are interested in N r (α) defined as where r is an even integer.We also consider N 1 (α) defined as To proceed we need the covariance matrix Γ x associated with the normalised composite density operator x .This is given as [80] where the normalisation factor is Z x = Tr(Γ x ) = Tr(ρ 2 A ).To proceed, let us define the charged Rényi negativities E r (α) (for even r) as where N r (α) are defined in (42) and we have set 1 = 2 = .Here C A is the restricted fermionic correlation matrix, i.e., C jl with j, l ∈ A. Now, Eq. ( 45) can be rewritten as where ν x j are eigenvalues of the matrices Γ x (44), and ζ j are the eigenvalues of the fermion correlation matrix C A introduced in the section 2. The charged negativity E(α) is defined as .We show results for the quench from the Néel state in the tight-binding chain.We plot the density of entropy S(q)/ versus rescaled time t/ for different charge sectors q.Here we fix q/ = 1/8 and q/ = 1/4.To ensure the weak-dissipative hydrodynamic limit the dissipation rates γ ± are fixed as γ − = 1/ , γ + = 5/(100 ).
Finally, let us observe that N 1 (α) is written in terms of the eigenvalues ν + j of Γ + (41) The strategy to compute the charge-imbalance-resolved negativity is to perform the Fourier transform of the quasiparticle prediction for log(N r (α)) and log(N 1 (α)).
5 Quasiparticle picture for symmetry-resolved entropies In this section we derive the quasiparticle picture for the symmetry-resolved entropies.In section 5.1 we focus on the charged moments of the reduced density matrix, whereas in section 5.2 we present our results for the symmetry-resolved von Neumann entropy.In section 5.2 we also discuss some useful approximations that are needed to derive analytically the behaviour of the resolved entropies.

Charged moments of the reduced density matrix
Here we discuss the quasiparticle picture for the charged moments (27) for the quench from the Néel state (6) and dimer state (7).We consider the weakly dissipative hydrodynamic limit, which corresponds to , t → ∞, γ ± → 0, at fixed and arbitrary t/ and γ ± , with the length of the subsystem A.
Since we are interested in the dynamics of Gaussian states, the system is fully characterised at any time after the quench by its two-point correlation functions.In the absence of pairing terms like c † c † or cc in the Hamiltonian at time t = 0, the relevant covariance matrix is C jl (cf.(17)).Again, the correlation matrix in the presence of gain/loss dissipation is simply related to the dynamical correlation function C jl (cf.(10)) for the quench without dissipation.To determine the quasiparticle picture for the charged moments, it is convenient to employ the representation of the fermionic correlator in (19) and (10) (11).To proceed, we formally expand (27) as where C A is C jl (cf.( 17)) with j, l ∈ A, and the coefficients c j are the coefficients of the Taylor expansion of log[x r e iα + (1 − x) r ] around x = 0.By following the same steps of Ref. [86], we can find the quasiparticle picture for the charged moments.Specifically, we first consider the expansion (49).Thus, the problem is reduced to that of determining the hydrodynamic behaviour of Tr C j A for arbitrary j.This last step can be performed by using the multidimensional version of the stationary phase approximation [105].By closely following the same steps as in [86], we obtain for both the quench from the Néel state and the dimer state that where v(k) = sin(k), i.e., the same group velocity of the fermions in the non-dissipative case.
In (50) we introduced the two contributions h mix r (α, k) and h q r (α, k) defined as Here the function h r (α, x) is defined as The term ( 51) is of purely dissipative origin and it vanishes in the limit γ ± → 0. Precisely, h mix r (α, k) corresponds to the charged moments obtained from ( 27) by replacing C A with the full-system correlator.Notice that in (51) the trace is over the two-by-two matrix t (k) (cf.( 18)).The term h q r in ( 50) is a "quantum" one, meaning that it is determined by both the dissipative and the unitary part of the evolution (15).In the limit γ ± → 0 this term gives the result of the quasiparticle picture in the non dissipative case.Importantly, h q r (α, k) depends only on the fermion density ρ t (k).In the absence of dissipation ρ t (k) does not depend on time, which is not the case at nonzero γ ± .For quenches in the tight-binding chain in the presence of gain/loss dissipation the time-evolved density ρ t (k) is obtained as [84] Here we have As it is clear from (55) for the quench from the Néel state and the dimer state γ(k) and α(k) do not depend on the quasimomentum k.This simplification holds for the gain/loss dissipation but it is not expected for the generic linear dissipation [84,85].In (54) ρ 0 (k) is the density of fermionic quasiparticles that fully characterises the Generalised Gibbs Ensemble (GGE) in the absence of dissipation.Specifically, ρ 0 (k) is computed as For the quench from the Néel and the dimer state, a straightforward calculation gives the densities ρ N 0 (k) and ρ D 0 (k) as Finally, by using the explicit form of t (k) (cf.( 18)) for the Néel and the dimer quench, we obtain that h mix r is the same for both of them, and it reads as Here n 0 and n 1 are given as where γ(k) and α(k) are given in (55).Notice that n 0 is obtained from (54) by setting ρ 0 (k) = 0 and n 1 by fixing ρ 0 (k) = 1.This reflects that at t = 0 for both the Néel state and the dimer state half of the eigenvalues of the correlation matrix C jl (cf.( 8)) are zero and half are one.They do not depend on time in the absence of dissipation.For gain/loss dissipation, the eigenvalues evolve independently according to (54).It is also interesting to observe that we can rewrite (50) as From ( 59) it is clear that h mix r plays a twofold role.Specifically, h mix r gives a volume-law contribution (first term in (59)) which does not contain information about the dynamics of the quasiparticles.On the other hand, h mix r also appears with a minus sign in the second term in (59), which is the "quantum" one because it is the only one surviving in the limit γ ± → 0.

Computing the symmetry-resolved entropies
To determine the quasiparticle picture for the symmetry-resolved entropies S r (q) (cf.( 26)), one has to employ the quasiparticle picture for the charged moments log(Z r (α)) in ( 25) and perform the Fourier transform.This is in general a challenging task, although closedform expressions have been obtained in the nondissipative case [69].Clearly, the Fourier transform in (25) can be performed numerically and we report the analytical predictions for the symmetry-resolved entropies in Fig. 2.However, this computation suffers from a numerical instability in solving the integral, and for this reason we omit the result for short time.
Another strategy is to use a saddle point approximation.Indeed, in (25) one has to assume q ∝ in order to have non trivial results.We have where we defined h r (α) as The saddle point condition for α reads as Since the integration over k in (61) can be performed analytically, the resulting saddle point condition for α * can be effectively solved numerically.Moreover, analytic solutions of (61) are also possible, although the obtained expressions are cumbersome, and we do not report them.Instead, we now discuss some features of the solution α * .In the absence of dissipation, α * is a complex number and the presence of a nonzero real part leads to a delay in the evolution of the symmetry-resolved entropies, i.e., the entropies ).We plot S 1 (q)/ versus the rescaled time t/ , and fixed = 80, for the quench from the fermionic Néel state in the tight binding chain.The different lines are the results obtained from the quasiparticle picture for different values of the charge q and dissipation γ + , using the saddle point approximation.For some values of t/ we report the numerical value of the symmetry-resolved entropy obtained through exact simulations of the dynamics (symbols).We clearly observe a delay time in agreement with Eq. ( 65), which does not depend on the dissipative term γ + .
start to be nonzero for t > t d , with t d the delay time [69].Specifically, for each value of the charge q we have that α * in the non-dissipative situation is given as Here we defined In the presence of gain and loss dissipation the delay effect is suppressed, as it was already shown in [21].
For simplicity, we focus on a quench from the Néel state |N .First, the saddle point condition (62) has three solutions.By a standard saddle point analysis [105], we observe that the only physical solution describing the exact lattice computations are the ones with zero real part for any time t.In the limit γ + → 0 (γ − → 0), i.e. pure loss or pure gain, we still find a delay time for q > /2 (q < /2).Indeed, in this case the saddle point condition has two solutions and, at short times t ≤ t d , the dominant solution of the saddle point condition gives a complex entropy which physically means that it averages to zero.As in the non-dissipative case, one finds that the delay time is obtained by solving the equation In Fig. 3 we report the behaviour of the symmetry-resolved entropy using the saddle point approximation in the presence of pure gain and we show that, for q < /2 it displays a delay, which is also confirmed by the exact lattice computations.
Having solved, for instance numerically, the saddle-point constraint (62), Z r (q) for t > t d is given as where α * is the solution of ( 62) with Re(α * ) = 0.
Another useful strategy to gain insights into the symmetry-resolved entropies is to use a quadratic approximation, by expanding the charged moments around α = 0.One obtains log(Z r (α)) log(Z r (0)) + iαB (1)  r − where With this quadratic expansion, we can perform the Fourier transform in (25).We obtain This approximation is valid only until ∆q = q −B (1) r ), or, in other words, for small fluctuations ∆q.Within this approximation, the corresponding symmetry-resolved entropies read From Eq. ( 70), we can also predict the asymptotic value of the symmetry-resolved entropies for t → ∞ at leading order in .For a Néel state, we obtain Now, Eq. ( 71) predicts a volume-law scaling ∝ , and a subleading logarithmic one ∝ log( ).Importantly, Eq. ( 71) relies on the validity of the quadratic approximation (70).
It is also interesting to derive the number entropy S (num) (cf.( 22)) within the quadratic approximation.It is obtained from p(q) = Z 1 (q), and it is defined as By using the quadratic approximation for Z 1 (q) (cf.( 69)), we obtain that where we defined J (t) as In ( 73) n 0 and n 1 are the same as in (58), and ρ t is given in (54) with ρ 0 (k) = 1/2 for a Néel state.The change of variables q → q − ρ t gives We anticipate that (75) correctly describes the exact value of the number entropy.In the large time limit, from (75) we find We will provide numerical benchmarks of ( 76) in section 7.2.
6 Quasiparticle picture for the charge-imbalance resolved fermionic negativity We now discuss the quasiparticle prediction for the charge-imbalance-resolved fermionic negativity.Here we focus on two intervals A 1 and A 2 embedded in an infinite chain.The two intervals are of the same length , and are at distance d (see Fig. 1).We first determine the scaling behaviour of the charged moments of the partial time-reversed transpose N r (α) (cf.( 42)) after the quench from the Néel state.Here we only consider the case of the Néel quench because this allows us to exploit the results of Ref. [86].To derive the quasiparticle prediction for N r (α) and N 1 (α) we proceed as in section 5. To illustrate the idea, let us first consider the case of log(N 1 (α)) (cf.( 48)).The main simplification in the treatment of N 1 is that it depends only on the eigenvalues of the matrix Γ + (cf.( 41)).We can write N 1 (α) as where the coefficients c 1,j (α) are obtained from with r = 1.Here Γ ± are the matrices defined in (41).Now, to determine the quasiparticle prediction for log(N 1 (α)) one has to compute the hydrodynamic limit of Tr(Γ j + ) for generic j.This has been computed in [86], and it reads We also have In ( 79) and ( 80), b = e −(γ + +γ − )t (cf.( 17)) and we defined a as Here the functions Θ 1 and Θ 2 are defined as One has Θ 1 = 1 at t = 0, and it decreases linearly up to t = /(2v(k)).At longer times Θ 1 is identically zero.On the other hand, Θ 2 is zero for t ≤ d/(2v(k)).At larger times Θ 2 grows linearly with time, up to t = (d + )/(2v(k)), where it reaches a maximum.Then, Figure 4: Effect of gain and loss dissipation in the charge-imbalance-resolved negativity E(q): E(q) is plotted versus t/ .Here we consider the negativity between two equal adjacent intervals of size .The data are for the quench from the Néel state in the tight-binding chain.We consider several values of the charge imbalance q.Our results are in the weak-dissipative hydrodynamic limit, i.e., at fixed t/ , q/ and γ + = 1/(2 ) and γ − = 1/(4 ), with = 80 (continuous line) and = 40 (solid line).
Θ 2 exhibits a linear decrease up to t = (d + 2 )/(2v(k)).At longer times Θ 2 is zero.The same function Θ 2 describes the dynamics of the mutual information after quantum quenches without dissipation [87,106].By using ( 79), ( 80) in ( 77), we obtain Here we introduced h r,α as To compute the charge-imbalance-resolved negativity, let us first consider the charged negativity E(α) as where N r (α) is defined in (42).The strategy, again, is to employ the results of [86].Unlike the case in (84), to evaluate log(N r (α)) one has to obtain the behaviour of Γ r x (cf.( 44)) for arbitrary r.These terms, however, can be obtained by using the results of Ref. [86].The calculation is cumbersome and, since it is similar to Ref. [86], we do not report it.For even r, the result reads as where one has to sum over the ± signs.The functions a and b are the same as in (81).In (87) the function h r,α is the same as in (78), whereas Θ 1 and Θ 2 are defined in (82) and (83).
In (87) we defined the new function Θ 3 as Now, the first integral in (87) originates from the contribution of Γ x (see the first term in ( 46)), whereas the second one is due to the second term in (46).Notice that Θ 3 is the same function describing the dynamics of the entropies of A 1 ∪ A 2 (see Fig. 1).The charged negativity is obtained by setting r = 1 in (87).
From the expression for the charged negativity E r (α) (cf.( 87)) we obtain the chargeimbalance-resolved negativity by first deriving the moments of the TR partial transpose Z R 1 ,re (q) via the Fourier transform in (37).This step has to be performed numerically.Finally, the charge-imbalance-resolved negativity E(q) is obtained from (38).
Our theoretical results for E(q) are shown in Fig. 4 and Fig. 5 for two equal-length adjacent and disjoint intervals, respectively.For two adjacent intervals we plot the density of negativity E(q)/ versus t/ .Our results hold in the weak-dissipative hydrodynamic limit.This corresponds to the usual hydrodynamic limit t, → ∞ with their ratio t/ fixed.We also consider the limit q → ∞ with fixed q/ .Finally, we consider the weak dissipation limit γ ± → 0 with fixed γ ± .The charge-imbalance-resolved negativity exhibits a typical "rise and fall" dynamics.This behaviour is observed also in the absence of dissipation [86].However, the short-time dynamics is not linear, in contrast with the non dissipative case.For each value of q/ in Fig. 4 we show results for = 40 and = 80 (continuous and dashed line, respectively).The dependence on , which could arise because of the Fourier tranform, is quite weak.Finally, we observe that at very short times E(q) suffers from numerical instability in the evaluation of the Fourier transform for this geometry, therefore we do not show the plot for short times.While it is likely that the presence of generic dissipation suppresses the time delay, as it happens for the entropies, clarifying this issue would require considerable numerical effort, and we leave it as an open problem.In Fig. 5 we consider two disjoint intervals at distance d = /4.The behaviour of E(q) is similar to the case of adjacent intervals, the only difference being that now the negativity is zero up to t d/(2v max ), with v max the maximum velocity in the system.

Numerical benchmarks
In this section we provide numerical benchmarks for the results derived in sections 5 and 6.Precisely, in section 7.1 we discuss the charged moments, and in section 7.2 and section 7.3 we focus on the symmetry-resolved entropy.In section 7.4 we investigate the weak dissipative limit and logarithmic scaling corrections.In section 7.5 we consider the charge-imbalanceresolved negativity.

Charged moments
Let us start discussing the hydrodynamic limit of the charged moments Z r (α) (cf.(27) for their definitions in the lattice model).We report numerical results for the charged moments in Fig. 6 and Fig. 7 for the quench from the Néel state and the Majumdar-Ghosh state, respectively.In the figures we show log(Z r (α))/ plotted versus t/ .The results are for a finite interval of length embedded in the infinite chain.We show results for several values of the gain/loss rates γ ± (symbols in the figure).Importantly, although we are interested in the weak-dissipative limit in which γ ± → 0 with γ ± fixed, here we show results for finite γ ± .The continuous lines in the figures is the theoretical prediction (50).Notice that at t = 0 one has that log(Z r (α On the other hand, in the limit t → ∞, one has that (50) reduces to We should stress that since both h mix r (α) and h q r (α) depend on time, one has to take the limits t = 0 and t → ∞ also within the integrand in ( 89) and (90).Now, as it is clear from Fig. 6 and Fig. 7, the agreement between the numerical data (symbols in the figures) and the quasiparticle picture is perfect, even away from the weak dissipation limit.

Symmetry-resolved entropies
Having confirmed our theoretical results for the charged moments (50), we now consider the symmetry-resolved von Neumann entropy.As discussed in section 5.2, the strategy is to plug the quasiparticle prediction (50) for the charged moments Z r (α) in (25), then performing the Fourier transform numerically.In Fig. 8 we compare our analytic prediction against exact numerical data for the symmetry-resolved entropy.We show results for both the quench from the Néel state and the Majumdar-Ghosh state (left and right panel, respectively).We stress that the fact that the symmetry-resolved entropies are non-vanishing for t > 0 for any charge sector is due to the presence of dissipative terms: They suppress the delay and all the charge   25) as Fourier transform of the quasiparticle prediction (50), while the dashed line for the Néel quench is obtained from the quadratic approximation (cf.Eq. ( 73)).
sectors are populated as soon as the state evolves, in contrast with what happens without dissipation [69] or in the presence of pure loss/gain.

Quadratic approximation
Let us now discuss the regime of validity of the quadratic approximation discussed in section 5.2.This is obtained by expanding the charged moments Z r (α) (cf.( 27)) around α = 0.Then, one can perform the Fourier transform in (25) analytically.The result is reported in Fig. 8 (left) as dashed lines.Interestingly, the quadratic approximation works well for the case of balanced gain and loss dissipation, i.e., for γ + = γ − .However, away from the case of balanced dissipation the quadratic approximation fails to describe the dynamics of the entropy.
The reason of this failure is understood by monitoring the behaviour of the charge probability distribution p(q) = Z 1 (q).We show results for p(q) in Fig. 9.The left and right panels in the figure are for the quench for the Néel and the Majumdar-Ghosh state, respectively, while top and bottom panels show results for γ + γ − and γ − = γ + .As we can notice from the figure, the evolution of p(q) is affected by the gain/loss terms.The main difference is that while for γ − = γ + (i.e.n ∞ = 1/2), p(q) remains peaked around q = /2, for γ + γ − , since n ∞ 1 (cf.( 17)), as t increases, the peak of p(q) moves toward smaller and smaller values of q.Hence, since the symmetry resolved entropies in Fig. 8 are for fixed value of q, for γ + γ − we have that ∆q = q − q drifts as time passes to larger and larger values because of the drift of q.Now the quadratic approximation is valid only in the neighbourhood of ∆q = 0 (i.e. for ∆q B (2) r ).Working at fixed q, such condition breaks down unless γ − = γ + , for which q Figure 10: Symmetry-resolved entanglement entropy for the tight-binding chain after the Néel quench in the presence of gain and loss dissipation.Weak-dissipative hydrodynamic scaling limit.We plot S(q)γ versus the rescaled time tγ, with γ = γ + + γ − .To reach the scaling limit we rescale also the dissipation rates as γ + = 1/(2 ), γ − = 2/ .Moreover, we fix q/ = 2/5.Sizeable corrections are visible and prevent from observing a data collapse.In the inset we show data at fixed t/ = 1.8 as a function of 1/ .The continuous line is a fit to A − 1/2γ log( ) + B/ , with A, B fitting parameters.
stays constant.This observation physically explains the failure of the quadratic approximation observed in Fig. 8.Moreover, from Fig. 8 it is clear that, for γ + = 0.01, γ − = 0.1, there is only a tiny time window (corresponding to ∆q 0) where the quadratic approximation is valid.

Logarithmic scaling corrections and the number entropy
It is important to investigate explicitly the scaling behaviour in the weak-dissipative hydrodynamic limit.This is discussed in Fig. 10.We plot exact numerical data for the symmetryresolved von Neumann entropy for , t → ∞ at fixed t/ and γ ± → 0 with fixed γ ± .Moreover, we also fix q/ .Clearly, Fig. 10 shows that there are strong finite corrections.Indeed, in the weak-dissipative hydrodynamic limit all the data in Fig. 10 are expected to collapse on the same curve.Deviations from the expected behaviour can be attributed to the presence of logarithmic terms in S(q).This is due to the fact that the result for the symmetry-resolved entropies should behave like (for t → ∞) The logarithmic correction in ( 91) is captured already by the quadratic approximation (cf.( 71)).In the sum rule (22), this logarithmic term cancels with the number entropy and so it can be better understood by studying the latter which is shown in Fig. 11.The figure shows that at t = 0 the number entropy is zero, reflecting that the initial state is a product state.At later times the entropy grows, and it exhibits an exponential saturation at long times.The symbols in the figures are exact numerical data for S (num) , whereas the continuous line is the analytic result obtained by using the quasiparticle picture.The result is obtained by using (59) (see also (50)), performing the Fourier transform (25) numerically.The predictions of the quasiparticle picture are in perfect agreement with the lattice results, even away from the weak dissipative limit, i.e., for finite γ ± .In the figure we also show the long time limit (76) (horizontal line), as obtained from the quadratic approximation, which describes quite well the lattice results.75)).The horizontal line is the result in the long time limit (76).

Charge-imbalance-resolved negativity
We now numerically investigate the validity of our results for the fermionic negativity (see section 6).First, we consider the charged negativity E(α) defined in (86).In Fig. 12 we show exact numerical data for E 1 (α) and E(α) (left and right panels, respectively) in the tightbinding chain with gain and loss of fermions.We consider the quench from the fermionic Néel state.Since we are interested in the weak-dissipative hydrodynamic limit, we plot E(α)/ versus t/ .To reach the weak-dissipation limit we rescale the dissipation rates γ ± as 1/ .Specifically, we choose γ + = 1/(2 ) and γ − = γ + /2.We show results for several values of α.The top panels are for two adjacent intervals, whereas the bottom ones are for two disjoint ones at d = /2.Our theoretical predictions are (84) and (87).Now, the symbols in Fig. 12 are exact numerical data for E(α) and E 1 (α).Our theoretical predictions are shown as continuous lines in the figure, and are in perfect agreement with the numerical data.Finally, we discuss the charge-imbalance-resolved negativity E(q).This is obtained by using ( 84) and ( 87) in (37), performing the Fourier transform numerically, and using (38).We show numerical results in Fig. 13.The top panels in the figure are for fixed γ − = 0.1 and γ + = 0.01, = 40 and q = 0, 5, 10.The top-right panel shows results for two adjacent intervals, i.e., for d = 0, whereas the top-left one is for two disjoint intervals at d = /2.The continuous lines in the figures are the theoretical results obtained from the quasiparticle picture.The agreement between lattice results and hydrodynamic limit is perfect for two disjoint intervals.Deviations from the scaling limit are visible for two adjacent intervals (top-right panel in Fig. 13).These are due to the fact that our results are valid only in the weak-dissipative hydrodynamic limit.Similar corrections are observed in the absence of dissipation [87].The right plot also clearly shows that there is no delay in the presence of dissipation, since E(q) is different from zero as soon as t > 0, at least for q = 0, 5.The scaling behaviour for two adjacent intervals is better investigated in Fig. 13 in the bottom panel.We now show results in the limit , t, q → ∞ with t/ , q/ fixed.Moreover, we consider γ ± → 0 at fixed γ ± .Now, although for finite some deviations are present, upon increasing the numerical results approach the analytic predictions (lines in the figure).We plot the densities E(α)/ and E 1 (α)/ as a function of t/ .We fix γ + = 1/(2 ), γ − = γ + /2.The data for different at fixed α collapse on the same curve, which is in perfect agreement with (87) (left panels), and (84) (right panels).

Conclusions
We investigated the effects of gain and loss dissipation in the dynamics of the symmetryresolved entropies and charge-imbalance-resolved negativity after a quantum quench in the tight-binding chain.This choice of dissipative terms preserves the block-diagonal structure of the reduced density matrix, and it is therefore suitable to study the symmetry resolution of entanglement in each charge sector.We derived quasiparticle picture formulas for the dynamics of charged moments of the reduced density matrix, the symmetry-resolved von Neumann entropy, and the charge-imbalance-resolved negativity.Our results hold in the hydrodynamic limit of large intervals and long times, with their ratio fixed.Moreover, to ensure a nontrivial scaling behaviour we also considered weak dissipation.We showed that while the symmetryresolved entropies are dominated by dissipative processes, the resolved negativity exhibits the typical rise and fall dynamics, reflecting that it is sensitive to entangled quasiparticles.Furthermore, we observed that the entropy does not exhibit a time delay (unless we consider a pure gain/loss evolution), contrarily to what happens in the non-dissipative case.Let us now mention some possible directions for future research.First, it would be interesting to extend our results to other types of dissipation, for instance, for generic quadratic Lindblad master equations [84,85].While this is feasible for the symmetry-resolved entropies, it is a challenging task for the negativity.Indeed, even for the non-resolved negativity only the case of free fermions with gain and loss dissipation was investigated so far [86].Clearly, it would be important to go beyond the case of quadratic Lindblad equations, considering more complicated dissipators, such as dephasing [107].For instance, it would be interesting  to employ recent results obtained for this dissipator in Ref. [108].Finally, one could study the effects of dissipative terms which do not preserve the block-diagonal structure of the reduced density matrix and quantify how much the symmetry is broken in this setup by using the entanglement asymmetry studied in [73,75].Similarly, one could consider initial states that explicitly breaks the internal symmetry and study the evolution of the asymmetry subject to gain and loss.

Figure 2 :
Figure2: Effect of gain and loss dissipation in the symmetry-resolved von Neumann entropy S(q).We show results for the quench from the Néel state in the tight-binding chain.We plot the density of entropy S(q)/ versus rescaled time t/ for different charge sectors q.Here we fix q/ = 1/8 and q/ = 1/4.To ensure the weak-dissipative hydrodynamic limit the dissipation rates γ ± are fixed as γ − = 1/ , γ + = 5/(100 ).

Figure 3 :
Figure 3: Time-delay in the dynamics of the symmetry-resolved von Neumann entropy with pure gain (γ − = 0).We plot S 1 (q)/ versus the rescaled time t/ , and fixed = 80, for the quench from the fermionic Néel state in the tight binding chain.The different lines are the results obtained from the quasiparticle picture for different values of the charge q and dissipation γ + , using the saddle point approximation.For some values of t/ we report the numerical value of the symmetry-resolved entropy obtained through exact simulations of the dynamics (symbols).We clearly observe a delay time in agreement with Eq. (65), which does not depend on the dissipative term γ + .

Figure 5 :
Figure 5: Same as in Fig. 4 for two disjoint intervals at d = /4.

Figure 6 :
Figure 6: Effect of gain/loss dissipation in the tight-binding chain after the quench from the fermionic Néel state: Real part of the logarithm of the charged moments Re[log(Z r (α))]/ plotted versus t/ .We show results for r = 1 and r = 2 in the left and right panel, respectively.The symbols are exact lattice data for several values of α and gain/loss rates γ ± .The data are for fixed = 80.The continuous lines are the quasiparticle predictions (cf.(50)).

Figure 7 :Figure 8 :
Figure 7: Same as in Fig. 6 for the quench from the Majumdar-Ghosh state.

Figure 9 :
Figure 9: Charge probability distribution of the tight-binding chain after the Néel (left) or dimer quench (right) as a function of q for fixed values of , tγ ± .The solid line is obtained by performing numerically the integral in Eq. (25) as Fourier transform of the quasiparticle prediction (50), while the dashed line for the Néel quench is obtained from the quadratic approximation (cf.Eq. (73)).

Figure 11 :
Figure 11: Dynamics of the number entropy (cf.(22)) S (num) after the quench from the Néel state in the tight-binding chain with gain/loss dissipation.The symbols denote exact lattice results, whereas the solid line is the prediction of the quasiparticle picture (cf.(75)).The horizontal line is the result in the long time limit(76).

Figure 12 :
Figure 12: Time evolution of E(α) (left panels) and E 1 (α) (right panels) after a quench from the Néel state in the tight-biding model.The top and bottom panels are for two adjacent and two disjoint intervals at distance d = /2, respectively (see Fig.1).We plot the densities E(α)/ and E 1 (α)/ as a function of t/ .We fix γ + = 1/(2 ), γ − = γ + /2.The data for different at fixed α collapse on the same curve, which is in perfect agreement with (87) (left panels), and (84) (right panels).