Symmetry-resolved entanglement in fermionic systems with dissipation

Let us introduce the symmetry-resolved Rényi and von Neumann entropies. The latter is the entropy calculated from the different charge sectors of the reduced density matrix as

$$\begin{equation} S_{1}\left(q\right): = -\mathrm{Tr} \left[\rho_A\left(q\right)\log\left(\rho_A\left(q\right)\right)\right], \end{equation} \tag{ 21 }$$

where, again, $\rho_A(q)$ is the normalised block of ρ_A corresponding to the charge q. To proceed, we plug equation (20) into the definition of the von Neumann entropy, to obtain the decomposition as

$$\begin{equation} S_{1} = \sum_q p\left(q\right)S_{1}\left(q\right)-\sum_q p\left(q\right)\log \left(p\left(q\right)\right): = S^{\left(c\right)}+S^{\left(\mathrm{num}\right)}, \end{equation} \tag{ 22 }$$

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^{(\mathrm{num})}$ measures the charge fluctuations between subsystem A and its complement, and it is called number entropy [14, 59–62].

Notice that from $\rho_A(q)$ one can define the symmetry-resolved Rényi entropies [46] as

$$\begin{equation} S_r\left(q\right): = \frac{1}{1-r}\log \mathrm{Tr} \left[\rho^r_A\left(q\right)\right]. \end{equation} \tag{ 23 }$$

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 $\Pi_q$ is a nonlocal operator. Here we follow a different strategy, which was put forward in reference [18, 19] and which we now outline. The central objects of this approach are the charged moments, $Z_r(\alpha)$. They are defined as

$$\begin{equation} Z_r\left(\alpha\right) = \mathrm{Tr}\left[\rho_A^r e^{i\alpha Q_A}\right]. \end{equation} \tag{ 24 }$$

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.

$$\begin{equation} \mathcal{Z}_r\left(q\right): = \mathrm{Tr}\left(\Pi_q\rho^r\right) = \int_{-\pi}^{\pi}\frac{\mathrm d\alpha}{2\pi}e^{-iq\alpha}Z_r\left(\alpha\right),\quad p\left(q\right) = \mathcal{Z}_1\left(q\right), \end{equation} \tag{ 25 }$$

where $\Pi_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

$$\begin{equation} S_r\left(q\right) = \frac{1}{1-r}\log\left( \frac{\mathcal{Z}_r\left(q\right)}{\mathcal{Z}_1^r\left(q\right)}\right),\quad S_1\left(q\right) = \lim_{r \to 1}S_r\left(q\right), \end{equation} \tag{ 26 }$$

where the symmetry-resolved von Neumann entropy S₁ 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]

$$\begin{equation} \log\left(Z_r\left(\alpha\right)\right) = \mathrm{Tr}\log\left[C_A^re^{i\alpha}+\left(1-C_A\right)^r\right], \end{equation} \tag{ 27 }$$

where C_A (cf equation (8)) is the fermionic correlation matrix restricted to part A of the chain.

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 \cup A_2$ (see figure 1), where A₁ (A₂) is an interval of length $\ell_1$ ($\ell_2$). We can write

$$\begin{equation} \rho_A = \sum_{ijkl}\langle e^1_i,e^2_j|\rho_A|e^1_k, e^2_l\rangle\vert e_i^1,e_j^2\rangle\langle e^1_k,e^2_l\vert, \end{equation} \tag{ 28 }$$

where $\vert e_j^1\rangle$ and $\vert e_k^2\rangle$ are orthonormal bases for the Hilbert spaces of A₁ and A₂, respectively. The partial transpose $\rho^{T_1}_A$ of ρ_A with respect to A₁ is defined as

$$\begin{equation} \left(\vert e_i^1,e_j^2\rangle\langle e^1_k,e^2_l\vert\right)^{T_1} : = \vert e_k^1,e_j^2\rangle\langle e^1_i,e^2_l\vert. \end{equation} \tag{ 29 }$$

The negativity is written in terms of the negative eigenvalues of $\rho_A^{T_1}$, and it is essentially the trace norm of $\rho_A^{T_1}$. In terms of its eigenvalues λ_i , the trace norm of $\rho_A^{T_1}$ can be rewritten as

$$\begin{equation} \mathrm{Tr}|\rho_A^{T_1}| = \sum_i |\lambda_i| = \sum_{\lambda_i >0}|\lambda_i|+\sum_{\lambda_i <0}|\lambda_i| = 1+2\sum_{\lambda_i<0}|\lambda_i|, \end{equation} \tag{ 30 }$$

where in the last equality we use the normalisation $\sum_i \lambda_i = 1$. The negativity is given by

$$\begin{equation} \mathcal{N}: = \frac{\mathrm{Tr}|\rho_A^{T_1}|-1}{2}, \end{equation} \tag{ 31 }$$

and, as it is clear by its definition, in the absence of negative eigenvalues, $\mathrm{Tr}|\rho_A^{T_1}| = 1$ and it 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].

Zoom In Zoom Out Reset image size

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 TR partial transpose, can be computed. To introduce the fermionic TR partial transpose it is convenient to introduce Majorana fermion operators w_j as

$$\begin{equation} w_{2j-1} = c_j+c^\dagger_j, \quad w_{2j} = i\left(c_j-c^\dagger_j\right), \end{equation} \tag{ 32 }$$

with c_j the original fermions. The density matrix in the Majorana representation takes the form

$$\begin{equation} \rho_A = \sum_{\substack{\kappa,\tau \\ |\kappa|+|\tau| = \mathrm{even}}} M_{\kappa,\tau} w_{m_1}^{\kappa_{m_1}} \dots w_{2m_{\ell_1}}^{\kappa_{2m_{\ell_1}}} w_{m^{^{\prime}}_1}^{\tau_{m^{^{\prime}}_1}} \dots w_{2m^{^{\prime}}_{\ell_2}}^{\tau_{2m^{^{\prime}}_{\ell_2}}}. \end{equation} \tag{ 33 }$$

Here, $w^0_x = \unicode{x1D7D9}$, $w^1_x = w_x$, $\kappa_i, \tau_j \in \{0,1\}$, and κ (τ) is a $2m_{\ell_1}$-component vector ($2m^{^{\prime}}_{\ell_2}$) with entries $0,1$ and with norm $|\kappa| = \sum_j \kappa_j$ ($|\tau| = \sum_j \tau_j$). The constraint on the parity of $|\kappa|+|\tau|$ is due to the fact that we require the density matrix to commute with the total fermion-number parity operator. In (33) $M_{\kappa,\tau}$ are complex numbers. Using equation (33), the partial TR transpose $\rho^{R_1}_A$ with respect to the subsystem A₁ is defined by

$$\begin{equation} \rho_A^{R_1} = \sum_{\substack{\kappa,\tau \\ |\kappa|+|\tau| = \mathrm{even}}} i^{|\kappa|}M_{\kappa,\tau} w_{m_1}^{\kappa_{m_1}} \dots w_{2m_{\ell_1}}^{\kappa_{2m_{\ell_1}}} w_{m^{^{\prime}}_1}^{\tau_{m^{^{\prime}}_1}} \dots w_{2m^{^{\prime}}_{\ell_2}}^{\tau_{2m^{^{\prime}}_{\ell_2}}}. \end{equation} \tag{ 34 }$$

The matrix $\rho_A^{R_1}$ is not necessarily Hermitian, and may have complex eigenvalues, although $\mathrm{Tr}\rho_A^{R_1} = 1$. Still the eigenvalues of the operator $[\rho_A^ {R_1} (\rho_A^{R_1})^{\dagger} ]$ are all real. Following references [102–104] we define the fermionic logarithmic negativity as

$$\begin{equation} {\cal E}: = \log\mathrm{Tr}\sqrt{\rho_A^{R_1}\left(\rho_A^{R_1}\right)^\dagger}. \end{equation} \tag{ 35 }$$

One can also define the (fermionic) Rényi logarithmic negativities ${\cal E}_r = \log(R_r)$ with R_r given by

$$\begin{equation} R_r = \begin{cases} \mathrm{Tr}\left(\rho_A^{R_1}\left(\rho_A^{R_1}\right)^{\dagger}\dots \rho_A^{R_1}\left(\rho_A^{R_1}\right)^{\dagger}\right), &\quad r \quad \mathrm{even},\\[12pt] \mathrm{Tr}\left(\rho_A^{R_1}\left(\rho_A^{R_1}\right)^{\dagger}\dots \rho_A^{R_1}\right),& \quad r \quad \mathrm{odd}. \end{cases} \end{equation} \tag{ 36 }$$

The negativity can be obtained as $\displaystyle \mathcal{E} = \lim_{r_e \to 1} \log(R_{r_e})$, 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 TR

$$\begin{equation} N_r\left(\alpha\right) = \begin{cases} \mathrm{Tr}\left(\rho_A^{R_1}\left(\rho_A^{R_1}\right)^{\dagger}\dots \rho_A^{R_1}\left(\rho_A^{R_1}\right)^{\dagger}e^{i{Q}_A\alpha}\right), &\quad r \quad \mathrm{even},\\[12pt] \mathrm{Tr}\left(\rho_A^{R_1}\left(\rho_A^{R_1}\right)^{\dagger}\dots \rho_A^{R_1}e^{i{Q}_A\alpha}\right), & \quad r \quad \mathrm{odd}. \end{cases} \end{equation} \tag{ 37 }$$

A Fourier transform allows us to derive the moments of the TR partial transpose ${\mathcal{Z}}_{R_1,r}$ as

$$\begin{equation} \mathcal{Z}_{R_1,r}\left(q\right) = \displaystyle \int_{-\pi}^{\pi}\frac{\mathrm d\alpha}{2\pi }e^{-i\alpha q}N_{r}\left(\alpha\right), \quad p\left(q\right) = \displaystyle \int_{-\pi}^{\pi}\frac{\mathrm d\alpha}{2\pi }e^{-i\alpha q}N_{1}\left(\alpha\right), \end{equation} \tag{ 38 }$$

from which we define the ratios $R_r(q)$ and the charge-imbalance-resolved negativity ${\mathcal E}(q)$ as

$$\begin{equation} R_{r}\left(q\right) = \frac{\mathcal{Z}_{R_1,r}\left(q\right)}{p^r\left(q\right)}, \qquad \mathcal{E}\left(q\right) = \lim_{r_e\to 1}\log\left(\frac{{\mathcal Z}_{R_1,r_e}\left(q\right)}{p^{r_e}\left(q\right)}\right). \end{equation} \tag{ 39 }$$

Let us stress that the replica limit for $R_{r_e}(q)$ requires an analytic continuation $r_e \to 1$ from the even sequences $R_{r_e}$, whereas in the denominator we can set r → 1.

Let us now discuss how to compute the charged moments $N_r(\alpha)$ in free-fermion systems. The covariance matrix Γ is the central object to compute the charged moments and the symmetry-resolved negativity. Γ is defined as

$$\begin{equation} \Gamma_{jl} = \unicode{x1D7D9}-2C_{jl}, \quad j,l\in A, \end{equation} \tag{ 40 }$$

where $\unicode{x1D7D9}$ is the identity matrix restricted to subsystem $A = A_1\cup A_2$, and C_jl is the fermionic correlation matrix (17). If the subsystem A is made of two intervals as $A = A_1\cup A_2$ (see figure 1), we can decompose Γ as

$$\begin{equation} \Gamma = \begin{pmatrix} \Gamma_{11}& \Gamma_{12} \\ \Gamma_{21} & \Gamma_{22} \end{pmatrix}, \end{equation} \tag{ 41 }$$

where Γ₁₁ and Γ₂₂ are the reduced covariance matrices of the two subsystems A₁ and A₂, respectively, while Γ₁₂ and Γ₂₁ contain the cross correlations between them. By simple Gaussian states manipulations, the covariance matrices $\Gamma_\pm$ associated with $\rho^{R_1}_A$ and $(\rho^{R_1}_A)^{\dagger}$ (cf (34) for their definitions) are obtained as

$$\begin{equation} \Gamma_{\pm} = \begin{pmatrix} -\Gamma_{11}& \pm i \Gamma_{12} \\ \pm i\Gamma_{21} & \Gamma_{22} \end{pmatrix}. \end{equation} \tag{ 42 }$$

Here we are interested in $N_r(\alpha)$ defined as

$$\begin{equation} N_{r}\left(\alpha\right) = \mathrm{Tr}\left[\left(\rho_A^{R_1} \left(\rho_A^{R_1}\right)^{\dagger}\right)^{r/2}e^{iQ_A\alpha}\right], \end{equation} \tag{ 43 }$$

where r is an even integer. We also consider $N_1(\alpha)$ defined as

$$\begin{equation} N_{1}\left(\alpha\right) = \mathrm{Tr}\left[\rho_A^{R_1} e^{i Q_A\alpha}\right]. \end{equation} \tag{ 44 }$$

To proceed we need the covariance matrix $\Gamma_\mathrm{x}$ associated with the normalised composite density operator $\rho_{\mathrm{x}} = \rho_A^{R_1} (\rho_A^{R_1})^{\dagger}/Z_{\mathrm{x}}$. This is given as [79]

$$\begin{equation} \Gamma_{\mathrm{x}} = \left(\unicode{x1D7D9}+\Gamma_+\Gamma_-\right)^{-1}\left(\Gamma_++\Gamma_-\right), \end{equation} \tag{ 45 }$$

where the normalisation factor is $Z_{\mathrm{x}} = \mathrm{Tr}(\Gamma_{\mathrm{x}}) = \mathrm{Tr}(\rho_A^2)$.

To proceed, let us define the charged Rényi negativities ${\cal E}_r(\alpha)$ (for even r) as

$$\begin{align} {\cal E}_r\left(\alpha\right): = \log\left(N_r\left(\alpha\right)\right)& = -2i\alpha\ell +\mathrm{Tr}\log\left[ \left(\frac{\unicode{x1D7D9}-\Gamma_\mathrm{x}}{2}\right)^{r/2}+ e^{i\alpha}\left(\frac{\unicode{x1D7D9}+\Gamma_\mathrm{x}}{2}\right)^{r/2}\right] \nonumber\\ & \quad +\frac{r}{2}\mathrm{Tr}\log\left[C_A^2+\left(\unicode{x1D7D9}-C_A\right)^2\right], \end{align} \tag{ 46 }$$

where $N_r(\alpha)$ are defined in (43) and we have set $\ell_1 = \ell_2 = \ell$. Here C_A is the restricted fermionic correlation matrix, i.e. C_jl with $j,l\in A$. Now, equation (46) can be rewritten as

$$\begin{align} {\cal E}_r\left(\alpha\right)& = -2i\alpha \ell+\sum_{j = 1}^{2\ell}\log \left[ \left(\frac{1-\nu_j^{\mathrm{x}}}{2} \right)^{r/2}+ e^{i\alpha} \left(\frac{1+\nu_j^{\mathrm{x}}}{2} \right)^{r/2} \right]\nonumber\\ & \quad +\frac{r}{2}\sum_{j = 1}^{2\ell}\log \left[ \zeta_j^{2}+ \left(1-\zeta_j\right)^{2} \right], \end{align} \tag{ 47 }$$

where $\nu_j^{\mathrm{x}}$ are eigenvalues of the matrices $\Gamma_{\mathrm{x}}$ (45), and ζ_j are the eigenvalues of the fermion correlation matrix C_A introduced in the section 2. The charged negativity ${\cal E}(\alpha)$ is defined as

$$\begin{equation} {\cal E}\left(\alpha\right): = \lim_{r\to1}{\cal E}_r\left(\alpha\right),\quad\mathrm{with}\,\, r\,\,\mathrm{even}. \end{equation} \tag{ 48 }$$

Finally, let us observe that $N_1(\alpha)$ is written in terms of the eigenvalues $\nu^+_j$ of $\Gamma_+$ (42)

$$\begin{equation} \log\left(N_{1}\left(\alpha\right)\right) = -2i \alpha \ell+\sum_{j = 1}^{2\ell}\log \left[ \left(\frac{1-\nu^+_j}{2} \right)+e^{i\alpha} \left(\frac{1+\nu^+_j}{2} \right) \right]. \end{equation} \tag{ 49 }$$

The strategy to compute the charge-imbalance-resolved negativity is to perform the Fourier transform of the quasiparticle prediction for $\log(N_{r}(\alpha))$ and $\log(N_{1}(\alpha))$.

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 $\ell,t\to\infty$, $\gamma^\pm\to0$, at fixed and arbitrary $t/\ell$ and $\gamma^\pm \ell$, with $\ell$ 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^{\dagger}c^{\dagger}$ 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 $\widetilde 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

$$\begin{equation} \log\left(Z_r\left(\alpha\right)\right) = \sum_{j = 0}^\infty z_j \mathrm{Tr}\,C_A^j, \end{equation} \tag{ 50 }$$

where C_A is C_jl (cf (17)) with $j,l\in A$, and the coefficients z_j are the coefficients of the Taylor expansion of $\log[x^r e^{i\alpha}+(1-x)^r]$ around x = 0.

By following the same steps as [86], we can find the quasiparticle picture for the charged moments. Specifically, we first consider the expansion (50). Thus, the problem is reduced to that of determining the hydrodynamic behaviour of $\mathrm{Tr}\,C_A^j$ 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

$$\begin{align} \log\left(Z_r\left(\alpha\right)\right) & = \ell\int \frac{\mathrm dk}{2\pi}\left[\mathrm{max}\left(1-2|v\left(k\right)|t/\ell,0\right) h^\mathrm{mix}_r\left(\alpha,k\right)\right.\nonumber\\ & \left.\quad +\mathrm{min}\left(2|v\left(k\right)|t/\ell,1\right) h^{\mathrm{q}}_r\left(\alpha,k\right)\right], \end{align} \tag{ 51 }$$

where $v(k) = \sin(k)$, i.e. the same group velocity of the fermions in the non-dissipative case. In (51) we introduced the two contributions $h_r^\mathrm{mix}(\alpha,k)$ and $h_r^\mathrm{q}(\alpha,k)$ defined as

$$\begin{align} h_r^\mathrm{mix} & = \mathrm{Tr}\, h_r\left(\alpha,\hat t^{^{\prime}}\left(k\right)\right), \end{align} \tag{ 52 }$$

$$\begin{align} h_r^\mathrm{q} & = h_r\left(\alpha,\rho_t\left(k\right)\right). \end{align} \tag{ 53 }$$

Here the function $h_r(\alpha,x)$ is defined as

$$\begin{equation} h_r\left(\alpha,x\right) = \log\left(e^{i\alpha}x^r+\left(1-x\right)^r\right). \end{equation} \tag{ 54 }$$

The term (52) is of purely dissipative origin and it vanishes in the limit $\gamma^\pm\to0$. Precisely, $h^\mathrm{mix}_r(\alpha,k)$ corresponds to the charged moments obtained from (27) by replacing C_A with the full-system correlator. Notice that in (52) the trace is over the two-by-two matrix $\hat t^{^{\prime}}(k)$ (cf (18)). The term $h_r^\mathrm{q}$ in (51) is a 'quantum' one, meaning that it is determined by both the dissipative and the unitary part of the evolution (15). In the limit $\gamma^\pm\to0$ this term gives the result of the quasiparticle picture in the non dissipative case. Importantly, $h^\mathrm{q}_r(\alpha,k)$ depends only on the fermion density $\rho_t(k)$. In the absence of dissipation $\rho_t(k)$ does not depend on time, which is not the case at nonzero $\gamma^\pm$. For quenches in the tight-binding chain in the presence of gain/loss dissipation the time-evolved density $\rho_t(k)$ is obtained as [83]

$$\begin{equation} \rho_t\left(k\right) = e^{-\gamma\left(k\right) t}\rho_0\left(k\right)+\frac{\alpha\left(k\right)}{\gamma\left(k\right)}\left(1-e^{-\gamma\left(k\right) t}\right). \end{equation} \tag{ 55 }$$

Here we have

$$\begin{equation} \gamma\left(k\right) = \gamma^++\gamma^-,\quad \alpha\left(k\right) = \gamma^+. \end{equation} \tag{ 56 }$$

As is clear from (56), for the quench from the Néel state and the dimer state $\gamma(k)$ and $\alpha(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 [83, 84]. In (55), $\rho_0(k)$ is the density of fermionic quasiparticles that fully characterises the generalised Gibbs ensemble (GGE) in the absence of dissipation. Specifically, $\rho_0(k)$ is computed as $\rho_0(k) = \langle\psi_0|c^\dagger_k c_k|\psi_0\rangle$. For the quench from the Néel and the dimer state, a straightforward calculation gives the densities $\rho^\mathrm{N}_0(k)$ and $\rho^\mathrm{D}_0(k)$ as

$$\begin{equation} \rho^\mathrm{N}_0\left(k\right) = \frac{1}{2},\quad \rho_0^\mathrm{D}\left(k\right) = |\cos\left(k\right)|. \end{equation} \tag{ 57 }$$

Finally, by using the explicit form of $\hat t^{^{\prime}}(k)$ (cf (18)) for the Néel and the dimer quench, we obtain that $h^{\mathrm{mix}}_r$ is the same for both of them, and it reads as

$$\begin{equation} h^\mathrm{mix}_r\left(\alpha,k\right) = \frac{h_r\left(\alpha,n_0\right)+h_r\left(\alpha,n_1\right)}{2}. \end{equation} \tag{ 58 }$$

Here n₀ and n₁ are given as

$$\begin{equation} n_0 = \frac{\alpha\left(k\right)}{\gamma\left(k\right)}\left(1-e^{-\gamma\left(k\right)t}\right),\quad n_1 = e^{-\gamma\left(k\right)t}+\frac{\alpha\left(k\right)}{\gamma\left(k\right)}\left(1-e^{-\gamma\left(k\right)t}\right), \end{equation} \tag{ 59 }$$

where $\gamma(k)$ and $\alpha(k)$ are given in (56). Notice that n₀ is obtained from (55) by setting $\rho_0(k) = 0$ and n₁ by fixing $\rho_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 (55). It is also interesting to observe that we can rewrite (51) as

$$\begin{equation} \log\left(Z_r\left(\alpha\right)\right) = \int\frac{\mathrm dk}{2\pi} \left\{\ell h_r^\mathrm{mix}\left(\alpha,k\right)+ \min\left(2|v\left(k\right)|t,\ell\right)\left[h_r^\mathrm{q}\left(\alpha,k\right)-h_r^\mathrm{mix}\left(\alpha,k\right)\right]\right\}. \end{equation} \tag{ 60 }$$

From (60) it is clear that $h_r^\mathrm{mix}$ plays a twofold role. Specifically, $h_r^\mathrm{mix}$ gives a volume-law contribution (first term in (60)) which does not contain information about the dynamics of the quasiparticles. On the other hand, $h_r^\mathrm{mix}$ also appears with a minus sign in the second term in (60), which is the 'quantum' one because it is the only one surviving in the limit $\gamma^\pm\to0$.

Let us now discuss some limits of $\log(Z_r(\alpha))$. First, at t = 0, since $h_r(\alpha,0) = 0$ and $h_r(\alpha,1) = i\alpha$, one has that $\log(Z_r(\alpha)) = i\alpha\ell/2$. For $t\to\infty$, one has that $n_0,n_1\to \alpha(k)/\gamma(k) = n_\infty$ (cf (17)) and $\rho_t\to n_\infty$. This implies that the second term in (60) vanishes and finally one obtains $\log(Z_r(\alpha))\to \ell h_r^\mathrm{mix}(\alpha,n_\infty)$.

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(\alpha))$ in (25) and perform the Fourier transform. This is, in general, a challenging task, although closed-form expressions have been obtained in the nondissipative case [68]. Clearly, the Fourier transform in (25) can be performed numerically and we report the analytical predictions for the symmetry-resolved entropies in figure 2. However, this computation suffers from a numerical instability in solving the integral, and for this reason we omit the result for short time.

Zoom In Zoom Out Reset image size

Another strategy is to use a saddle point approximation. Indeed, in (25) one has to assume $q\propto \ell$ in order to have nontrivial results. We have

$$\begin{equation} {\mathcal Z}_r\left(q\right) = \int_{-\pi}^\pi\frac{\mathrm d\alpha}{2\pi} e^{\ell h_r\left(\alpha\right)}, \end{equation} \tag{ 61 }$$

where we defined $h_r(\alpha)$ as

$$\begin{equation} h_r\left(\alpha\right) = -i\frac{q}{\ell}\alpha+\int_{-\pi}^\pi \frac{\mathrm d k}{2\pi}\left\{h^\mathrm{mix}_r\left(\alpha,k\right) +\mathrm{min}\left(2|v|t/\ell,1\right)\left[h^\mathrm{q}_r\left(\alpha,k\right)- h^\mathrm{mix}\left(\alpha,k\right)\right]\right\}. \end{equation} \tag{ 62 }$$

The saddle point condition for α reads as

$$\begin{equation} \partial_\alpha h_r\left(\alpha\right)\Big|_{\alpha = \alpha^*} = 0. \end{equation} \tag{ 63 }$$

Since the integration over k in (62) can be performed analytically, the resulting saddle point condition for $\alpha^*$ can be effectively solved numerically. Moreover, analytic solutions of (62) are also possible, although the obtained expressions are cumbersome, and we do not report them. Instead, we now discuss some features of the solution $\alpha^*$.

In the absence of dissipation, $\alpha^*$ 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 start to be nonzero for $t\gt t_d$, with t_d the delay time [68]. Specifically, for each value of the charge q in the quench from the Néel state we have that $\alpha^*$ in the nondissipative situation is given as

$$\begin{equation} \alpha^* = i\log\left(\frac{\mathcal{J}\left(t\right)-2\Delta q}{\mathcal{J}\left(t\right)+2\Delta q}\right),\quad\mathrm{with}\,\, \Delta q: = q-\frac{\ell}{2}. \end{equation} \tag{ 64 }$$

Here we define

$$\begin{equation} \mathcal{J}\left(t\right): = \int_{-\pi}^\pi\frac{\mathrm dk}{2\pi}\min\left(2|v\left(k\right)|t,\ell\right). \end{equation} \tag{ 65 }$$

In the presence of gain and loss dissipation the delay effect is suppressed, as was already shown in [21].

For simplicity, we focus on a quench from the Néel state $|\mathrm{N}\rangle$. First, the saddle point condition (63) 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 $\gamma^+\to 0$ ($\gamma^-\to 0$), i.e. pure loss or pure gain, we still find a delay time for $q\gt\ell/2$ ($q\lt\ell/2$). Indeed, in this case the saddle point condition has two solutions and, at short times $t \unicode{x2A7D} 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 nondissipative case, one finds that the delay time is obtained by solving the equation

$$\begin{equation} {\mathcal J}\left(t\right)-2|\Delta q| = 0,\quad\Delta q = q-\frac{\ell}{2}. \end{equation} \tag{ 66 }$$

In figure 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\lt\ell/2$ it displays a delay, which is also confirmed by the exact lattice computations.

Zoom In Zoom Out Reset image size

Having solved, for instance numerically, the saddle-point constraint (63), ${\mathcal Z}_r(q)$ for $t\gt t_d$ is given as

$$\begin{equation} {\mathcal Z}_r\left(q\right) = e^{\ell h_r\left(\alpha^*\right)}\sqrt{\frac{1}{2\pi\ell|h_r^{^{\prime\prime}}\left(\alpha^*\right)|}}, \end{equation} \tag{ 67 }$$

where $\alpha^*$ is the solution of (63) with $\mathrm{Re}(\alpha^*) = 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

$$\begin{equation} \log\left(Z_r\left(\alpha\right)\right)\simeq\log\left(Z_r\left(0\right)\right)+i\alpha\mathcal{B}^{\left(1\right)}_r-\frac{\alpha^2}{2}\mathcal{B}^{\left(2\right)}_r, \end{equation} \tag{ 68 }$$

where

$$\begin{equation} \mathcal{B}^{\left(1\right)}_r = \partial_{\alpha}\log\left(Z_r\left(\alpha\right)\right) \Big |_{\alpha = 0},\quad \mathcal{B}^{\left(2\right)}_r = \partial^2_{\alpha}\log\left(Z_r\left(\alpha\right)\right) \Big |_{\alpha = 0}. \end{equation} \tag{ 69 }$$

With this quadratic expansion, we can perform the Fourier transform in (25). We obtain

$$\begin{equation} \mathcal{Z}_r\left(q\right) = Z_r\left(0\right)e^{-\frac{\left(q-\mathcal{B}^{\left(1\right)}_r\right)^2}{2\mathcal{B}^{\left(2\right)}_r}}\frac{1}{\sqrt{2\pi\mathcal{B}^{\left(2\right)}_r }}. \end{equation} \tag{ 70 }$$

This approximation is valid only until $\Delta q = q-\mathcal{B}^{(1)}_r$ is ${\mathcal O}(\sqrt{\mathcal{B}^{(2)}_r})$, or, in other words, for small fluctuations $\Delta q$. Within this approximation, the corresponding symmetry-resolved entropies read

$$\begin{equation} S_r\left(q\right) = S_r-\frac{1}{2}\log\left(2\pi\right)-\frac{1}{2\left(1-r\right)}\log\frac{\mathcal{B}^{\left(2\right)}_r}{\left(\mathcal{B}^{\left(2\right)}_1\right)^r}-\frac{\left(q-\mathcal{B}^{\left(1\right)}_r\right)^2}{2\left(1-r\right)}\left[\frac{1}{\mathcal{B}^{\left(2\right)}_r}- \frac{r}{\mathcal{B}^{\left(2\right)}_1}\right]. \end{equation} \tag{ 71 }$$

From equation (71), we can also predict the asymptotic value of the symmetry-resolved entropies for $t \to \infty$ at leading order in $\ell$. For a Néel state, we obtain

$$\begin{align} S_1\left(q\right)& = -\ell n_{\infty}\log \left(n_{\infty}\right)-\ell \left(1-n_{\infty}\right)\log \left(1-n_{\infty}\right)-1/2 \log\left(2\pi \ell\right)-1/2 \log\left(\left(1-n_{\infty}\right)n_{\infty} \right) \nonumber\\ & \quad +\frac{\tanh ^{-1}\left(1-2 n_{\infty}\right) \left(-2 \ell n_{\infty}^2 q+\ell n_{\infty} \left(n_{\infty} \left(\ell+2 n_{\infty}-3\right)+1\right)+\left(2 n_{\infty}-1\right) q^2\right)}{\ell \left(n_{\infty}-1\right) n_{\infty}}\nonumber\\ & \quad -\frac{\left(\ell n_{\infty}-q\right)^2}{2\ell n_{\infty}\left(1-n_{\infty}\right)}. \end{align} \tag{ 72 }$$

Now, equation (72) predicts a volume-law scaling $\propto\ell$, and a subleading logarithmic one $\propto\log(\ell)$. Importantly, equation (72) relies on the validity of the quadratic approximation (71).

It is also interesting to derive the number entropy $S^\mathrm{(\mathrm{num})}$ (cf (22)) within the quadratic approximation. It is obtained from $p(q) = \mathcal{Z}_1(q)$, and it is defined as

$$\begin{equation} S^{\left(\mathrm{num}\right)} = -\sum_q p\left(q\right) \log\left(p\left(q\right)\right). \end{equation} \tag{ 73 }$$

By using the quadratic approximation for $\mathcal{Z}_1(q)$ (cf (70)), we obtain that

$$\begin{align} p\left(q\right) = \frac{1}{\sqrt{\pi}}\frac{1}{\sqrt{2\ell \rho_t \left(1-\rho_t\right)-2\mathcal{J}^{^{\prime}}\left(t\right)\left(\rho_t-n_0\right)^2}} \exp\left({-\frac{\left(q-\ell \rho_t\right)^2}{2\ell \rho_t\left(1-\rho_t\right)-2\mathcal{J}^{^{\prime}}\left(t\right) \left(\rho_t-n_0\right)^2}}\right), \nonumber\\ \end{align} \tag{ 74 }$$

where we define $\mathcal{J}^{^{\prime}}(t)$ as

$$\begin{equation} \mathcal{J}^{^{\prime}}\left(t\right) = \ell \int \frac{\mathrm d\lambda}{2\pi}\mathrm{max}\left(1-2|v|t/\ell,0\right). \end{equation} \tag{ 75 }$$

In (74) n₀ and n₁ are the same as in (59), and ρ_t is given in (55) with $\rho_0(k) = 1/2$ for a Néel state. The change of variables $q \to q -\ell \rho_t$ gives

$$\begin{equation} S^{\left(\mathrm{num}\right)}\simeq \int_{-\infty}^{\infty}\mathrm dq p\left(q\right) \log\left(p\left(q\right)\right) = \frac{1}{2}\left[1+\log\left(\pi\left(2\ell \rho_t\left(1-\rho_t\right)- 2\mathcal{J}^{^{\prime}}\left(t\right)\left(\rho_t-n_0\right)^2\right) \right)\right]. \end{equation} \tag{ 76 }$$

We anticipate that (76) correctly describes the exact value of the number entropy. In the large time limit, from (76) we find

$$\begin{equation} \lim_{t \to \infty}S^{\left(\mathrm{num}\right)} = \frac{1}{2}+\frac{1}{2}\log \left[2\pi \ell n_{\infty}\left(1-n_{\infty}\right) \right], \end{equation} \tag{ 77 }$$

We will provide numerical benchmarks of (77) in section 7.2.

Let us start discussing the hydrodynamic limit of the charged moments $Z_r(\alpha)$ (cf (27) for their definitions in the lattice model). We report numerical results for the charged moments in figures 6 and 7 for the quench from the Néel state and the Majumdar–Ghosh state, respectively. In the figures we show $\log(Z_r(\alpha))/\ell$ plotted versus $t/\ell$. The results are for a finite interval of length $\ell$ embedded in the infinite chain. We show results for several values of the gain/loss rates $\gamma^\pm$ (symbols in the figure). Importantly, although we are interested in the weak-dissipative limit in which $\gamma^\pm\to0$ with $\gamma^\pm \ell$ fixed, here we show results for finite $\gamma^\pm$. The continuous lines in the figures is the theoretical prediction (51). Notice that at t = 0 one has that

$$\begin{equation} \left.\log\left(Z_r\left(\alpha\right)\right)\right|_{t = 0} = \ell\int_{-\pi}^\pi\frac{\mathrm dk}{2\pi}h^{\mathrm{mix}}_r \left(\alpha\right). \end{equation} \tag{ 90 }$$

On the other hand, in the limit $t\to\infty$, one has that (51) reduces to

$$\begin{equation} \left.\log\left(Z_r\left(\alpha\right)\right)\right|_{t\to\infty} = \ell\int_{-\pi}^\pi\frac{\mathrm dk}{2\pi}h^{\mathrm{q}}_r\left(\alpha\right). \end{equation} \tag{ 91 }$$

We should stress that since both $h^\mathrm{mix}_r(\alpha)$ and $h^\mathrm{q}_r(\alpha)$ depend on time, one has to take the limits t = 0 and $t\to\infty$ also within the integrand in (90) and (91). Now, as it is clear from figures 6 and 7, the agreement between the numerical data (symbols in the figures) and the quasiparticle picture is perfect, even away from the weak dissipation limit.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Having confirmed our theoretical results for the charged moments (51), we now consider the symmetry-resolved von Neumann entropy. As discussed in section 5.2, the strategy is to plug the quasiparticle prediction (51) for the charged moments $Z_r(\alpha)$ in (25), then performing the Fourier transform numerically. In figure 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 sectors are populated as soon as the state evolves, in contrast with what happens without dissipation [68] or in the presence of pure loss/gain.

Zoom In Zoom Out Reset image size

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(\alpha)$ (cf (27)) around α = 0. Then, one can perform the Fourier transform in (25) analytically. The result is reported in figure 8 (left) as dashed lines. Interestingly, the quadratic approximation works well for the case of balanced gain and loss dissipation, i.e. for $\gamma^+ = \gamma^-$. 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) = \mathcal{Z}_1(q)$. We show results for p(q) in figure 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 $\gamma^+ \ll \gamma^-$ and $\gamma^- = \gamma^+$. 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 $\gamma^- = \gamma^+$ (i.e. $n_{\infty} = 1/2$), p(q) remains peaked around $q = \ell/2$, for $\gamma^+\ll\gamma^-$, since $n_{\infty} \ll 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 figure 8 are for fixed value of q, for $\gamma^+\ll\gamma^-$ we have that $\Delta q = q-\bar q$ drifts as time passes to larger and larger values because of the drift of $\bar q$. Now the quadratic approximation is valid only in the neighbourhood of $\Delta q = 0$ (i.e. for $\Delta q \ll \mathcal{B}^{(2)}_r$). Working at fixed q, such condition breaks down unless $\gamma^- = \gamma^+$, for which $\bar q$ stays constant. This observation physically explains the failure of the quadratic approximation observed in figure 8. Moreover, from figure 8 it is clear that, for $\gamma_+ = 0.01, \gamma_- = 0.1$, there is only a tiny time window (corresponding to $\Delta q\simeq 0$) where the quadratic approximation is valid.

Zoom In Zoom Out Reset image size

It is important to investigate explicitly the scaling behaviour in the weak-dissipative hydrodynamic limit. This is discussed in figure 10. We plot exact numerical data for the symmetry-resolved von Neumann entropy for $\ell,t\to\infty$ at fixed $t/\ell$ and $\gamma^\pm\to0$ with fixed $\gamma^\pm\ell$. Moreover, we also fix $q/\ell$. Clearly, figure 10 shows that there are strong finite $\ell$ corrections. Indeed, in the weak-dissipative hydrodynamic limit all the data in figure 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_1(q)$. This is due to the fact that the result for the symmetry-resolved entropies should behave like (for $t \to \infty$)

$$\begin{equation} S_1\left(q\right) = a\ell-\frac{1}{2}\log \ell+O\left(1/\ell\right) .\end{equation} \tag{ 92 }$$

Zoom In Zoom Out Reset image size

The logarithmic correction in (92) is captured already by the quadratic approximation (cf (72)). 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 figure 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^{(\mathrm{num})}$, whereas the continuous line is the analytic result obtained by using the quasiparticle picture. The result is obtained by using (60) (see also (51)), 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 $\gamma^\pm$. In the figure we also show the long time limit (77) (horizontal line), as obtained from the quadratic approximation, which describes quite well the lattice results.

Zoom In Zoom Out Reset image size

We now numerically investigate the validity of our results for the fermionic negativity (see section 6). First, we consider the charged negativity ${\cal E}(\alpha)$ defined in (87). In figure 12 we show exact numerical data for ${\cal E}_1(\alpha)$ and ${\cal E}(\alpha)$ (left and right panels, respectively) in the tight-binding 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 ${\cal E}(\alpha)/\ell$ versus $t/\ell$. To reach the weak-dissipation limit we rescale the dissipation rates $\gamma^\pm$ as $1/\ell$. Specifically, we choose $\gamma^+ = 1/(2\ell)$ and $\gamma^- = \gamma^+/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 = \ell/2$. Our theoretical predictions are (85) and (88). Now, the symbols in figure 12 are exact numerical data for ${\cal E}(\alpha)$ and ${\cal E}_1(\alpha)$. Our theoretical predictions are shown as continuous lines in the figure, and are in perfect agreement with the numerical data.

Zoom In Zoom Out Reset image size

Finally, we discuss the charge-imbalance-resolved negativity ${\mathcal{E}}(q)$. This is obtained by using (85) and (88) in (38), performing the Fourier transform numerically, and using (39). We show numerical results in figure 13. The top panels in the figure are for fixed $\gamma^- = 0.1$ and $\gamma^+ = 0.01$, $\ell = 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 = \ell/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 figure 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 $\mathcal{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 figure 13 in the bottom panel. We now show results in the limit $\ell,t,q\to\infty$ with $t/\ell,q/\ell$ fixed. Moreover, we consider $\gamma^\pm\to0$ at fixed $\gamma^\pm\ell$. Now, although for finite $\ell$ some deviations are present, upon increasing $\ell$ the numerical results approach the analytic predictions (lines in the figure).

Zoom In Zoom Out Reset image size