Dynamics of Charge-Resolved Entanglement after a Local Quench

TL;DR

This paper investigates how charge-resolved entanglement measures evolve over time after a local quench in 1+1D conformal systems, revealing that the total entanglement growth stems from a subset of charge sectors, with results validated by numerical and exact methods.

Contribution

It introduces a conformal field theory approach to analyze charge-resolved entanglement dynamics post-quench, linking sector contributions to overall entanglement growth.

Findings

Total entanglement grows as log(t) after a quench.

Charge sector contributions scale as sqrt(log(t)).

Numerical and exact results agree with theoretical predictions.

Abstract

Quantum entanglement and its main quantitative measures, the entanglement entropy and entanglement negativity, play a central role in many body physics. An interesting twist arises when the system considered has symmetries leading to conserved quantities: Recent studies introduced a way to define, represent in field theory, calculate for 1+1D conformal systems, and measure, the contribution of individual charge sectors to the entanglement measures between different parts of a system in its ground state. In this paper, we apply these ideas to the time evolution of the charge-resolved contributions to the entanglement entropy and negativity after a local quantum quench. We employ conformal field theory techniques and find that the known dependence of the total entanglement on time after a quench, $S_A \sim \log(t)$, results from $\sim\sqrt{\log(t)}$ significant charge sectors, each of which contributes $\sim\sqrt{\log(t)}$ to the entropy. We compare our calculation to numerical results obtained by the time-dependent density matrix renormalization group algorithm and exact solution in the noninteracting limit, finding good agreement between all these methods.

Equations

$$S_A=-\mathrm{Tr}{\rho }_A\log ({\rho }_A).$$

$$S_A^{(n)}=\mathrm{Tr}{\rho }_A^n.$$

$${\mathcal{N}}_{A_1,A_2}=\frac{\left|\left|{\rho }_{A_1\cup A_2}^{T_2}\right|\right|-1}{2},$$

$${\rho }_A=\sum _{i,j,i^{\prime },j^{\prime }}{c}_{i,j,i^{\prime },j^{\prime }}{|i⟩}_{A_1}{|j⟩}_{A_2}{⟨i^{\prime }|}_{A_1}{⟨j^{\prime }|}_{A_2}\to$$

$${\rho }_A^{T_2}=\sum _{i,j,i^{\prime },j^{\prime }}{c}_{i,j,i^{\prime },j^{\prime }}{|i⟩}_{A_1}{|j^{\prime }⟩}_{A_2}{⟨i^{\prime }|}_{A_1}{⟨j|}_{A_2},$$

$${\mathcal{N}}_{A_1,A_2}^{(n)}=\mathrm{Tr}{\left({\rho }_{A_1\cup A_2}^{T_2}\right)}^n,$$

$$\begin{array}{rl}& S_A(N_A)=-\mathrm{Tr}{\rho }_A^{(N_A)}\log {\rho }_A^{(N_A)},\\ & S_A^{(n)}(N_A)=\mathrm{Tr}{\left({\rho }_A^{(N_A)}\right)}^n.\end{array}$$

$$\begin{array}{rl}& {\mathcal{N}}_{A_1,A_2}(Q)=\frac{\left|\left|{\left({\rho }_A^{T_2}\right)}^{(Q)}\right|\right|-1}{2},\\ & {\mathcal{N}}_{A_1,A_2}^{(n)}(Q)=\mathrm{Tr}{\left({\left({\rho }_A^{T_2}\right)}^{(Q)}\right)}^n.\end{array}$$

$$S_A^{(n)}\propto ⟨{\mathcal{T}}_n(x=x_L,\tau =0){\tilde{\mathcal{T}}}_n(x=x_R,\tau =0)⟩,$$

$$d_n=\overline{d_n}=\frac{c}{24}\left(n-n^{-1}\right),$$

$$S_A^{\left(n\right)}=\mathrm{Tr}{\rho }_A^n=c_n{\left(\frac{L_A}{a}\right)}^{-c\left(n-n^{-1}\right)/6},$$

$$S_A=\frac{c}{3}\log \frac{L_A}{a}+c_1^{\prime },$$

$$\omega (z)=\frac{z}{ϵ}+\sqrt{{\left(\frac{z}{ϵ}\right)}^2+1},$$

$$S_A^{\left(n\right)}\left(t\right)\propto {\left(\frac{{ϵ}^2+(vt{)}^2}{aϵ/2}\right)}^{-c\left(n-n^{-1}\right)/12},$$

$$S_A(t)=\frac{c}{3}\log \frac{vt}{a}+k_0,$$

$$\omega (z)=\coth \left(\frac{\pi ϵ}{2L}\right)\frac{1+\zeta (z)}{1-\zeta (z)},$$

$$\text{with}\zeta (z)=\sqrt{\frac{\sinh \frac{\pi }{L}(z+ϵ)}{\sinh \frac{\pi }{L}(z-ϵ)}}.$$

$$S_A^{(n)}(\alpha )=\mathrm{Tr}{\rho }_A^n{e}^{i{\hat{N}}_A\alpha }.$$

$$S_A^{(n)}(N_A)={\int }_{-\pi }^{\pi }\frac{d\alpha }{2\pi }{e}^{-i{N}_A\alpha }{S}_A^{(n)}(\alpha ).$$

$$\begin{array}{rl}& H=\frac{1}{2}v\left(K(\Pi {)}^2+\frac{1}{K}({\partial }_x\varphi {)}^2\right),\end{array}$$

$$\mathcal{V}=e^{i\frac{\alpha }{2\pi }{\varphi }_j},$$

$$d^{(\mathcal{V})}=\frac{K}{2}\left(\frac{{\alpha }^2}{2\pi }\right).$$

$$d_n(\alpha )=\frac{c(n-n^{-1})}{24}+\frac{d^{(\mathcal{V})}}{n}.$$

$$S_A^{(n)}(N_A)\cong S_A^{(n)}(\alpha =0)\sqrt{\frac{\pi n}{2K\ln L_A/a}}{e}^{-\frac{{\pi }^{2}n\Delta N_A^2}{2K\ln L_A/a}},$$

$$\mathcal{V}=\sum _{m=-\infty }^{\infty }{w}_{|m|}{e}^{i{\left(\frac{\alpha +2\pi m}{2\pi }\right)}^2\varphi },$$

$$S_A^{(n)}(\alpha ,t)=S_A^{(n)}(\alpha =0,t)\times$$

$$\sum _{m=-\infty }^{\infty }{w}_m{\left(\frac{L_{\mathrm{eff}}}{a}\right)}^{-\frac{c}{12}(n-n^{-1})-\frac{K}{n}{\left(\frac{\alpha +2\pi m}{2\pi }\right)}^2}.$$

$${\mathcal{N}}_{A_1,A_2}^{(n)}\propto ⟨{\mathcal{T}}_n(x_L,0){\tilde{\mathcal{T}}}_n^2(x_B,0){\mathcal{T}}_n(x_R,0)⟩.$$

$${\mathcal{N}}_{A_1,A_2}^{(n)}=\prod _i{\left|\frac{d\omega }{dz}\right|}_{z_i}^{d_i}⟨{\mathcal{T}}_n({\omega }_L){\tilde{\mathcal{T}}}_n^2({\omega }_B){\mathcal{T}}_n({\omega }_R)⟩,$$

$$⟨{\mathcal{T}}_n({\omega }_1){\tilde{\mathcal{T}}}_n^2({\omega }_2){\mathcal{T}}_n({\omega }_3)⟩=$$

$$\prod _i\frac{1}{|{\omega }_i-\widetilde{{\omega }_i}{|}^{d_{n,i}}}{\left(\frac{{\eta }_{1,3}^{d_n^{\left(2\right)}-2d_n}}{{\eta }_{1,2}^{d_n^{\left(2\right)}}{\eta }_{2,3}^{d_n^{\left(2\right)}}}\right)}^{1/2}\mathcal{F}\left(\{{\eta }_{j,k}\}\right),$$

$${\mathcal{N}}_{A_1,A_2}^{(n)}(\alpha )=\mathrm{Tr}(({\rho }_A^{T_2}{)}^n{e}^{i\hat{Q}\alpha }),$$

$${\mathcal{N}}_{A_1,A_2}^{(n)}(Q)={\int }_{-\pi }^{\pi }\frac{d\alpha }{2\pi }{\mathcal{N}}_{A_1,A_2}^{(n)}(\alpha )e^{-i\alpha Q}.$$

Full text

I Introduction

The discussion of entanglement started in the early days of quantum mechanics by Einstein, Podolsky, and Rosen Einstein et al. (1935) as well as Schrödinger Schrödinger (1935), yet quantum entanglement remains an active topic of research in several fields of quantum theory  Srednicki (1993); Ryu and Takayanagi (2006); Calabrese and Cardy (2009); Nishioka et al. (2009), and specifically in quantum many body systems Amico et al. (2008); Horodecki et al. (2009); Laflorencie (2016). Entanglement in many body systems is used for elucidating their physics Osterloh et al. (2002); Osborne and Nielsen (2002); Vidal et al. (2003); Calabrese and Cardy (2009), for understanding the limits of simulating quantum systems on a classical computer Cirac and Verstraete (2009); Verstraete et al. (2008); Schollwöck (2011); Orús (2014), and to characterize their utility as a resource for various quantum information applications Ekert (1991); Bennett and Wiesner (1992); Bennett et al. (1993); Shor (1997); Bouwmeester et al. (1997); Wootters (1998); Nielsen and Chuang (2000); Gisin et al. (2002); Eisert (2006); Harrow et al. (2009).

In order to define the entanglement measures we study in this paper, we introduce the density matrix (DM) for a state $\left|\psi\right\rangle$, $\rho=\left|\psi\right\rangle\left\langle\psi\right|$. We also define the reduced density matrix (RDM): For two subsystems $A$ and $B$, and a pure state $\left|\psi\right\rangle$ of the combined subsystems, the RDM is defined to be $\rho_{A}=\mathrm{Tr}_{B}\rho$. In this case (pure state of the total system), the basic measure for the entanglement of the subsystem $A$ with its environment $B$ is the von Neumann entanglement entropy (vNEE) Von Neumann (1932):

[TABLE]

We also introduce the moments of the density matrix, which we will refer to as Rényi entropies: The $n$th Rényi entropy (RE) is defined to be:

[TABLE]

We stress that the definition of the REs here is different from the standard definition, $\ln(\mathrm{Tr}\rho_{A}^{n})/(1-n)$. The REs obey $S_{A}=-\partial_{n}S^{(n)}_{A}|_{n\to 1}$. The REs are entanglement monotones, but they do not possess all the useful properties that the vNEE doesWilde (2016). However, they are easier to calculate, and can be measured experimentally more easilyHorodecki and Ekert (2002); Moura Alves and Jaksch (2004); Daley et al. (2012); Abanin and Demler (2012); Pichler et al. (2013); Islam et al. (2015); Banchi et al. (2016); Elben et al. (2018); Vermersch et al. (2018); Linke et al. (2018); Brydges et al. (2019); Cornfeld et al. (2018) (although a protocol for the measurement of the spectrum of the RDM of a bosonic system was proposed in Ref. Pichler et al., 2016).

When the total state of the two considered subsystems is not pure, different entanglement measures are needed. For two subsystems $A_{1}$ and $A_{2}$, coupled to an environment $B$, a popular entanglement measure is the entanglement negativity Peres (1996):

[TABLE]

where $||\cdot||$ denotes the trace norm, and the superscript $T_{2}$ stands for the partial transpose:

[TABLE]

where $(\left|i\right\rangle_{A_{1}},\left|i^{\prime}\right\rangle_{A_{1}})$ and $(\left|j\right\rangle_{A_{2}},\left|j^{\prime}\right\rangle_{A_{2}})$ are orthonormal bases for subsystems $A_{1}$ and $A_{2}$ respectively. Here too we define Rényi negativities (RNs):

[TABLE]

which could be analytically continued from an even integer $n$ to yield the negativity, using $\left|\left|\rho_{A_{1}\cup A_{2}}^{T_{2}}\right|\right|=\lim_{n\rightarrow 1/2}\mathcal{N}_{A_{1},A_{2}}^{(2n)}$. While the RNs are not even entanglement monotones (as opposed to the REs), they are still useful indicators of entanglement since (like the REs), they are experimentally measurable for both bosons Cornfeld et al. (2018); Gray et al. (2018) and fermions Cornfeld et al. (2018), and are easier to calculate.

We study a system that has some conserved charge $\hat{N}$, that obeys $\hat{N}_{A}+\hat{N}_{B}=\hat{N}$, for $\hat{N}_{i}$ the charge on subsystem $i$, such as a spin component or particle number. We assume that the state of the total system has some fixed value of $\hat{N}$, hence $\left[\rho,\hat{N}\right]=0$. Performing a partial trace over the equation above, we get $\left[\rho_{A},\hat{N}_{A}\right]=0$. This implies that the RDM is block diagonal, each block corresponding to some eigenvalue of $\hat{N}_{A}$ and denoted by $\rho_{A}^{(N_{A})}$, $\rho_{A}=\oplus_{N_{A}}\rho_{A}^{(N_{A})}$. This allows to define the charge resolved vNEE and REs Laflorencie and Rachel (2014); Goldstein and Sela (2018); Xavier et al. (2018); Barghathi et al. (2018, 2019):

[TABLE]

The blocks $\rho_{A}^{(N_{A})}$ are not normalized and thus $S_{A}=\sum_{N_{A}}S_{A}(N_{A}),\quad S_{A}^{(n)}=\sum_{N_{A}}S_{A}^{(n)}(N_{A})$. A normalized version has been used in some recent works Xavier et al. (2018a); Bonsignori et al. (2019); Barghathi et al. (2019). The corresponding entropies then characterize the entanglement following a projective measurement of the particle number in subsystem $A$. We prefer to follow our earlier convention Goldstein and Sela (2018); Cornfeld et al. (2018) and not to normalize the individual blocks, since the contributions of the unnormalized blocks to the total entanglement are not only more straightforwardly accessible in the calculation, but are also directly accessible experimentally, using either the protocol introduced in our previous works Goldstein and Sela (2018); Cornfeld et al. (2018) or the method based on random evolution introduced in Refs. Vermersch et al., 2018; Linke et al., 2018, if number-conserving evolution is used, as in those papers.

For the case of two subsystems coupled to an environment, it is useful to define the charge imbalance: $\hat{Q}=\hat{N}_{A_{1}}-\hat{N}_{A_{2}}^{T_{2}}$. By performing a partial trace over the equation shown above we get Cornfeld et al. (2018) $0=\left([\rho_{A},\hat{N}_{A}]\right)^{T_{2}}=[\rho_{A}^{T_{2}},\hat{N}_{A_{1}}-\hat{N}_{A_{2}}^{T_{2}}]=[\rho_{A}^{T_{2}},\hat{Q}]$. Note that in a standard Fock space basis, $\hat{N}_{A_{2}}^{T_{2}}=\hat{N}_{A_{2}}$. We may then define the charge imbalance resolved negativity and RNs Cornfeld et al. (2018):

[TABLE]

By definition, $S_{A}^{(1)}(N_{A})$ is the charge distribution in subsystem A, and $\mathcal{N}_{A_{1},A_{2}}^{(1)}(Q)$ is the charge imbalance distribution, demonstrating an inherent relation between entanglement and charge distribution. The charge resolved entanglement can be used as an instrument to study entanglement properties and gain a better understanding of the charge block structure of the entanglement spectrum. A particular interest in the relation between charge distribution and entanglement has risen for systems undergoing a local quench: We prepare two subsystems $A$ and $B$ in the ground state, couple the two subsystems at $t=0$ and study the evolution of entanglement between them. In Ref. Klich and Levitov, 2009 a relation between charge distribution (quantum noise) and the entanglement has been derived for noninteracting fermions following this type of quench, and motivated further study of the relation between entanglement and charge distribution Hsu et al. (2009); Song et al. (2010, 2011a, 2011b, 2012); Calabrese et al. (2012); Vicari (2012); Eisler and Rácz (2013); Eisler (2013); Chien et al. (2014); Petrescu et al. (2014); Thomas and Flindt (2015); Dasenbrook and Flindt (2015); Sinitsyn and Pershin (2016). With the charge-resolved entanglement measures just discussed it becomes apparent that one should not separately address the dynamics of charge and entanglement, but rather their combined measures. The goal of this work is to addresses this question.

In this paper we combine the methods from Refs. Goldstein and Sela, 2018; Cornfeld et al., 2018 for calculating the time dependent charge-resolved entanglement in 1+1D conformal field theory (CFT) systems. We compare them to exact results for the XX model, as well as to time dependent density matrix renormalization group (tDMRG) results for the XXZ model [Eq. (30)]. For the dynamics of the charge resolved vNEE, we get the contribution of the different charge sectors to the total entropy. We find that the behavior $S_{A}(t)\sim\log(t)$, as extracted in Ref. Calabrese and Cardy, 2009 and decribed below, mostly originates from $\sim\sqrt{\log(t)}$ significant charge sectors, each contributing $\sim\sqrt{\log(t)}$. We get satisfying results comparing the prediction above to numerical results. As for the charge imbalance resolved negativity, the CFT results involve hard-to-calculate conformal blocks, but we can still derive approximate expressions which qualitatively match the numerical results.

The rest of the paper is organized as follows: In Sec. II we present the theoretical background for the calculation of entanglement entropies and negativities for 1+1D CFT, and extend this method for calculating the charge resolved entanglement entropies and charge imbalance resolved negativities. We then use this method to calculate the charge resolved vNEEs and imbalance resolved RNs after a local quench. In Sec. III we compare the CFT predictions to numerical results for the XX and XXZ model. We summarize our findings and outline future directions in Sec. IV. Appendices A and B include the full formulas for the symmetry resolved entanglement entropies and negativities, as these are too cumbersome to appear in full in the main text.

II Conformal Field Theory Analysis

II.1 The Entanglement Entropy

We will first recap the calculation of the total REs of a 1+1D CFT in the ground state and after a local quench Calabrese and Cardy (2009), and then show how to generalize these techniques to the charge-resolved REs. The calculation is based on the replica trick. Space-imaginary time will be represented by the complex plane, with $z=x+iv\tau$, with $v$ the velocity of excitations. Representing the DM as a path integral, the RDM (when the total system is in its ground state) is represented as a path integral over the complex plane with a cut at $\tau=0$ at subsystem $A$’s coordinates. Different boundary conditions on this cut give different matrix elements of the RDM, see Fig 1a for the geometric representation. Sewing $n$ copies of the RDM together is equivalent to the multiplication of the n copies of the RDM, and by doing so one obtains the $n$th RE as a path integral on an $n$-sheet Riemann surface, see Fig. 1b. The transition between the copies is effected by twist field operators $\mathcal{T}_{n}$ and $\tilde{\mathcal{T}}_{n}$. $\mathcal{T}_{n}$ and $\tilde{\mathcal{T}}_{n}$ transfer particles from one copy to the next clockwise and counterclockwise, respectively. Hence, the RE is proportional to the correlation function of these twist fields:

[TABLE]

for $x=x_{L}$ and $x=x_{A}$ the endpoints of subsystem $A$ (taken as a single interval). The scaling dimension of the twist fields was calculated in Ref. Calabrese and Cardy, 2009:

[TABLE]

and the resulting RE and vNEE, again from Ref. Calabrese and Cardy, 2009, are:

[TABLE]

where $L_{A}=|x_{L}-x_{R}|$ is the length of subsystem $A$, $a$ is a cutoff corresponding to, e.g., a lattice spacing, $c$ is the conformal central charge, and $c_{n}$ and $c_{1}^{\prime}$ are constants that cannot be predicted using CFT.

Let us now go on to the time-dependent local quench scenario. We prepare two identical systems in the ground state, and at time $t=0$ we couple them at one contact point. For this time dependent case, the geometry representing the RDM is slightly more complicated: The two halves of the system are not connected before $t=0$. Upon analytical continuation to imaginary time, this is expressed as slits in each sheet, parallel to the imaginary time axis. The slits are separated by a small gap $\epsilon$, which serves as a convergence factor. This is shown in Fig. 2a. We then use the conformal transformation Calabrese and Cardy (2009):

[TABLE]

which takes the Riemann sheet with the slits into the right half plane, as demonstrated in Fig. 2b.

A system living on the right half plane has a boundary on the imaginary time axis, and requires the use of boundary CFT (BCFT) Cardy (2004): We separate our field into a holomorphic part and an anti-holomorphic part, and “unfold” the anti-holomorphic part into the left half plane. This also duplicates the operators in our system, turning an $n$-point function into a $2n$-point function, as demonstrated in Fig. 2c. For an infinite system with one boundary point between the systems $A$ and $B$ we are then left with a two-point function for the composite twist fields. The final result for the entanglement when $A$ and $B$ are the two parts of the system brought together by the local quench, is:

[TABLE]

where the last equation applies for $vt\gg\epsilon$, and $k_{0}$ is a non-universal constant. We can define the effective length $L_{\mathrm{eff}}(t)=\sqrt{\epsilon^{2}+(vt)^{2}}$ and notice that the REs take the form $S^{(n)}_{A}=\tilde{c}_{n}\left(\frac{L_{\mathrm{eff}}(t)}{a}\right)^{d_{n}}$, which is equivalent to Eq. (9) with $L_{A}\rightarrow L_{\mathrm{eff}}(t)$. One can think of the entanglement as carried by quasiparticles with velocity $v$, propagating from the quench point to the rest of the system. $L_{\mathrm{eff}}(t)$ represents the parts of the system reached by these quasiparticles. Even for more complicated geometries and conformal transformation, the result will always be of the form above, with $L_{\mathrm{eff}}(x,t)$ modified according to the geometry.

For finite subsystems $L_{A}=L_{B}=L/2$, one may use the following conformal transformation Stéphan and Dubail (2011):

[TABLE]

The effective lengths in this case, with the quench point at the boundary point between the subsystems or shifted from it, appear in Appendix A.

Let us now proceed to the calculation of the charge-resolved entropy in the ground state. In Ref. Goldstein and Sela, 2018, the twist fields are multiplied by a vertex operator $\mathcal{V}(\alpha)$. For a charge $\hat{N}$ related to an abelian $U(1)$ symmetry, the vertex operator is coupled to the charge such that a particle going counter-clockwise around $\mathcal{V}(\alpha)$ will acquire a phase $e^{i\alpha}$. If the vertex operators are placed at the edges of subsystem $A$, as the twist fields are, the resulting correlation function will give:

[TABLE]

We can think of this vertex operator as introducing a flux $\alpha$ for the particles, and define the measure in (13) to be the flux resolved RE, see Fig. 1c. A Fourier transform will lead us to the charge resolved RE,

[TABLE]

For concreteness, let us consider the case of gapless 1D spin-$\frac{1}{2}$ chain, or, equivalently (via the Jordan-Wigner transformation Jordan and Wigner (1928)), gapless 1D interacting fermion chain. According to the Luttinger liquid paradigm von Delft and Schoeller (1998); Sénéchal (1999), the system can be described both in terms of fermions and in terms of noninteracting bosons. In the spinless case, the relation between the fermion field $\psi$ and the boson field $\phi$ is $\psi\propto e^{i\phi}$. The bosonic field is governed by the gapless non-interacting Hamiltonian density:

[TABLE]

where $\Pi$ is the canonically conjugate momentum to $\phi$, and $v$ and $K$ are model-dependent parameters. $v$ is the field velocity. $K$ is called the Luttinger parameter, $K=1,K<1$ and $K>1$ corresponding to noninteracting, repulsive and attractive fermions, respectively.

The boson-fermion relation above allows us to define the vertex operator. We annotate the bosonic field living on the the $j$th plain (i.e., $j$th copy of the system) by $\phi_{j}$. The vertex operator can be chosen as

[TABLE]

with a scaling dimension derived in Ref. Goldstein and Sela, 2018:

[TABLE]

One may now derive an expression for the scaling function of the composite twist field  Goldstein and Sela (2018):

[TABLE]

This leads to Gaussian dependence of the charge-resolved RE on the charge, provided $\ln(L_{A}/a)\gg 1$:

[TABLE]

where $\Delta N_{A}=N_{A}-\left\langle\hat{N}_{A}\right\rangle$, $\left\langle\hat{N}_{A}\right\rangle$ being the expectation value of $\hat{N}_{A}$ in the ground state. For different geometries, one can replace $L_{A}$ in Eq. (19) by the suitable effective length $L_{\mathrm{eff}}$. The Gaussian dependence on the charge is a result of the unnormalized form of the RDM blocks. In Refs Xavier et al., 2018a; Bonsignori et al., 2019, the blocks are normalized and a the different normalized blocks rise to equal entropies, a property these works refer to as “equipartition”.

These results could be extended to smaller subsystems. We notice that the most general form for $\mathcal{V}$ is in fact:

[TABLE]

where $m$ is integer and $w_{|m|}$ is the corresponding weight. For a large enough $L_{\mathrm{eff}}$, the zeroth order term ($w_{|m|}=0$ for $m\neq 0$) is sufficient, since its correlation function features the slowest decay with subsystem size. In the time dependent case (to be discussed shortly), we often find it necessary to include the next order, $w_{1}\neq 0$, since its contribution is becoming important for short times. For the XX model studied below, the parameter $w_{1}$ can be extracted for the ground state using similar methods to those employed in Ref. Calabrese and Essler, 2010. However, in the time dependent case no such results are available, and we resort to extracting $w_{1}$ from a fit to our numerical results.

Having laid out all the necessary groundwork, we may now derive our new results for the charge resolved entropies following a local quench. Combining Eqs. (11), (17) and (20), we find the following expression for the dynamics of the flux resolved RE after a local quench:

[TABLE]

Plugging this into Eq. (14) will give us the charge resolved RE. For a finite system as studied below, we use the transformation (12) and get the relevant expression $L_{\mathrm{eff}}$ to be substituted in (21). For the case where the boundary between subsystems $A$ and $B$ is shifted from the quench point, we use the same conformal transformation (12), but place the twist fields away from the slits Stéphan and Dubail (2011), and obtain a modified expression for $L_{\mathrm{eff}}(t)$. The final formulas we used appear in Appendix A.

II.2 The Entanglement Negativity

The replica trick for the negativity is derived in Ref. Calabrese et al., 2013, and a geometric representation for it appears in Fig. 3. For two adjacent single-interval subsystems $A_{1}=[x_{L},x_{B}],A_{2}=[x_{B},x_{R}]$, twist field operators are added at the three boundary points $x_{L},x_{B},x_{R}$ such that:

[TABLE]

Calculating the negativity for the local quench case requires introducing these twist fields into $n$-copies of the geometry depicted in Fig. 2a. We are again forced to use BCFT and double the number of operators. We are left with calculating a 6-point function. For a $j$-point function with $j\geq 4$, the result can be predicted by CFT only up to a nonuniversal function $\mathcal{F}$, Francesco et al. (1997) which depends on the full operators content of the studied theory. We thus restrict ourselves to limits in which $\mathcal{F}$ is approximately constant. It is for this reason that we do not study the case where $A_{1}$ and $A_{2}$ are disjoint — following the arguments above, we will need to calculate an 8-point function, for which the effect of nonuniversal function $\mathcal{F}$ is expected to be even more significant.

In Ref. Wen et al., 2015 the RN after a local quench was found to be:

[TABLE]

where $\omega_{i}=\omega(x_{i})$ for $i=L,B,R$, as defined in Eq. (10) or (12). The unfolded 6-point function was found to be:

[TABLE]

where $d_{n}^{(2)}$ is the scaling dimension of $\mathcal{T}_{n}^{2}$ calculated in Ref. Calabrese et al., 2013, and is different for $n_{e}$ and $n_{o}$, even and odd values of $n$, respectively: $d_{n}^{(2)}=\frac{c}{12}\left(\frac{n_{e}}{2}-\frac{2}{n_{e}}\right)\quad d_{n}^{(2)}=\frac{c}{24}\left(n_{o}-\frac{1}{n_{o}}\right)$. In addition, $\tilde{\omega}_{i}=-\omega_{i}^{*}$ and $\eta_{i,j}=\frac{\left(\omega_{i}-\omega_{j}\right)\left(\tilde{\omega}_{i}-\tilde{\omega}_{j}\right)}{\left(\tilde{\omega}_{i}-\omega_{j}\right)\left(\omega_{i}-\tilde{\omega}_{j}\right)}$. $\mathcal{F}$ approaches a constant value for $vt\ll l,vt=l+0^{+}$, and $vt>l$ for two identical subsystems $L_{A_{1}}=L_{A_{2}}=l$, which is the case on which we focus here. We will only consider $l\ll L$, and so neglect corrections due to boundary conditions.

The charge-imbalance resolved RNs are again obtained from the flux resolved RNs,

[TABLE]

[TABLE]

We obtain the flux resolved RNs by adding vertex operators at the boundaries between the subsystemsCornfeld et al. (2018). The additivity of the scaling dimensions, Eq. (18), results in

[TABLE]

We can now present our new results. Combining Eqs. (10) and (27), and using the vertex operators correlation function from Ref. Francesco et al., 1997, we obtain the flux resolved RN for two adjacent systems with an infinite environment (we present the expression for the simple case, $w_{|m|>0}=0$)

[TABLE]

where $i$ and $j$ count the 6 positions of the vertex operators, $\omega^{\prime}_{i}=\omega_{L},\omega_{B},\omega_{R},\tilde{\omega}_{L},\tilde{\omega}_{B},\tilde{\omega}_{R}$ and $\alpha_{i}=\alpha(x_{i})=\alpha,-2\alpha,\alpha,-\alpha,2\alpha,-\alpha$, respectively. To the first order (and taking $w_{m}=0$ for $|m|>1$ in Eq. (20)), we get

[TABLE]

where $\mu=1\cdots 6$ , and $\alpha_{i}^{(\mu)}=\alpha_{i}+2\pi\delta_{i,\mu}-2\pi\delta_{i+3,\mu}$. The explicit expressions for $\omega^{\prime}_{i}-\omega^{\prime}_{j}$ were derived in Ref. Wen et al., 2015, and are reproduced for completeness in Appendix B.

III Numerical Results

III.1 Entanglement Entropy

The XXZ model. We compare our CFT predictions to numerically obtained results for the XXZ spin chain,

[TABLE]

when $\sigma^{+},\sigma^{-},\sigma^{z}$ are the Pauli matrices. Using the Jordan Wigner tranformation Jordan and Wigner (1928)

[TABLE]

the XXZ model can be interpreted as a spinless fermionic chain.

[TABLE]

where the $c_{n}$ annihilation operators obey the fermionic anti-commutation relations. The system is a gapless Luttinger liquid for $-1<\Delta\leq 1$. Its Luttinger parameter and velocity can be extracted from the Bethe ansatz Giamarchi (2003),

[TABLE]

We note that for lower values of $\Delta$, corresponding to higher values of $K$, the Gaussian distribution of the charge resolved entanglement is expected to be wider, and the effect of higher orders of $\alpha$ smaller, as can be seen from Eq. (21).

For $\Delta=0$ Eq. (32) describes spinless noninteracting fermions, allowing an exact calculation of the entanglement. Then, the entanglement Hamiltonian $\hat{H}_{A}$, defined as $\rho_{A}=e^{-\hat{H}_{A}}$, is quadratic. We follow Ref. Goldstein and Sela, 2018 (based on the method introduced in Ref. Peschel, 2003), and obtain exact results for the entanglement entropies,

[TABLE]

where $f_{l}=1/\left(e^{\epsilon_{l}}+1\right)$, and $\epsilon_{l}$ are the eigenvalues of $\hat{H}_{A}$. $f_{l}$ are the eigenvalues of the subsystem correlation matrix $C_{ij}=\left\langle c_{i}^{\dagger}c_{j}\right\rangle,i,j=1\dots L_{A}$, which can be obtained exactly for the noninteracting case by explicit calculation of the matrix $C$ of the full system in the momentum space $\left|k\right\rangle=\sum_{i=1}^{N}\sin\left(\frac{\pi ik}{N+1}\right)\left|i\right\rangle$, $k=1...N$, where the basis states are eigenstates of the energy, $C_{k,k^{\prime}}(t)=C_{k,k^{\prime}}(0)e^{-i(E_{k^{\prime}}-E_{k})t}$, for $E_{k}=2J\cos\left(\frac{\pi k}{N+1}\right)$. $C(0)$ is the correlation matrix in the ground state of the decoupled system, , which is diagonal in the corresponding eigenbasis $\left|q;L\right\rangle=\sum_{i=1}^{N/2}\sin\left(\frac{\pi iq}{N/2+1}\right)\left|i\right\rangle$ and $\left|q;R\right\rangle=\sum_{i=N/2+1}^{N}\sin\left(\frac{\pi(i-N/2)q}{N/2+1}\right)\left|i\right\rangle$, for $q=1...N/2$. In this basis $Q_{q\ell,q^{\prime}\ell^{\prime}}=\delta_{qq^{\prime}}\delta_{\ell\ell^{\prime}}\Theta(-E_{q})$, for $\ell,\ell^{\prime}=L,R$. $C=C^{L}\otimes C^{R}$ and $C^{\ell}_{q,q^{\prime}}=\delta_{q,q^{\prime}}\Theta(-E_{q})$.

Results for the case in which the boundary point between $A$ and $B$ is the quench point are presented in Fig. 4a, and for the case when the boundary point is moved away by $x=l$ from the quench point in Fig. 4b. We use open boundary conditions throughout. In all cases in this study, $a$ and $w_{1}$ from Eq. (21) were used as fitting parameters. The entanglement is periodic in time, as predicted by Eq. (12). This is the result of the entanglement being carried by the quasiparticles moving in velocity $v$ as mentioned above, bumping at the ends of the system and coming back to the other side Stéphan and Dubail (2011). In Fig. 4b, CFT predicts the entanglement to be constant for $\mod(vt,L)<l$ or $\mod(vt,L)>L-l$, when these quasiparticles are allegedly outside of subsystem A. The exact results present some tails in these time regimes, which are caused by excitations moving with velocities smaller than $v$, which are not accounted for by CFT Stéphan and Dubail (2011). We present results both for the zeroth order, calculated from Eq. (21) with $w_{m}=0$ for $m\neq 0$, and for the first order, in which $w_{1}\neq 0$. The zeroth order seems to be a satisfying approximation from Figs. 4a and 4b, but as can be seen in Fig. 4c, it is insufficient for large values of $|\alpha|$.

For nonzero values of $\Delta$, we compare the CFT prediction to numerical results obtained by the tDMRG algorithm Schollwöck (2011); Vidal (2003, 2004); Daley et al. (2004); White and Feiguin (2004), employing the QSpace tensor library Weichselbaum (2012). The extraction of the charged resolved entanglement spectrum is natural in this method, thanks to the block diagonal form of the Matrix Product State (MPS) matrices that are used to represent the state in this algorithm: in the MPS representation, one separates the state of an $N$-sites system into $N$ rank-3 tensors Schollwöck (2011). One entry of the tensor is called the physical index, and ranges over the physical Hilbert space of a single site, whose dimension is 2 in our case. The other two entries connect the site to the rest of the system from left and right, and represent the Hilbert space of all the sites on the left of the site or all the sites on the right of it, truncated to the so called “bond dimension” see Fig. 5a. For a subsystem $A$ extending the left edge of the system to site $i$, the eigenvalues of $\rho_{A}$ are extracted by placing the orthogonality center Schollwöck (2011) at the $i$th tensor and combining its physical index and its left index, thus getting a matrix with one entry representing subsystem $A$ and the second entry representing the rest of the system to the right of $A$. One can exploit the symmetry in the system by separating these tensors into blocks corresponding to different charge sectors. From these charge sector blocks we extract the eigenvalues of each block of the RDM separately.

We used a second order Trotter approximation with a timestep of $\Delta t=10^{-2}J^{-1}$. We set the MPS bond dimension Schollwöck (2011) to 1,024 and the MPS truncation error Schollwöck (2011) to $10^{-8}$ in all tDMRG runs for the entanglement entropy (in practice the truncation error was $~{}10^{-10}$ or less).

In order to stay well within the CFT region, we chose $-0.5\leq\Delta\leq 0.5$. Results for several values of $\Delta$ for open boundary conditions are plotted in Fig. 6. Here too results for both $w_{1}=0$ and $w_{1}\neq 0$ are presented. $w_{1}=0$ appears to generally be a satisfactory fit, as in the noninteracting case. We see oscillations in the numerical data not predicted by CFT: These are oscillations due to finite lattice spacing, which are absent in the CFT approximation and decrease for large $t$ (which is why they were not seen in Fig. 4, where the system length and achievable timescales are much longer). We note that the contribution of the corresponding spatial oscillations was calculated for the ground state case of the XX model in Ref. Calabrese and Essler, 2010, but these results are not straightforward to extend to the current time-dependent case.

III.2 Entanglement Negativity

The numerical method for the exactly solvable XX model is developed in Ref. Cornfeld et al., 2018, and is an extension of the entanglement Hamiltonian method used for the entanglement entropy case. We rely on the result from Refs. Eisler and Zimborás, 2015, 2016 that in this noninteracting case the partially transposed RDM is a sum of two Gaussian matrices: $\rho_{A}^{T_{2}}=\sum_{\sigma=\pm 1}u_{\sigma}\frac{\hat{O}_{\sigma}}{\mathrm{Tr}{\hat{O}_{\sigma}}}$. $u_{\sigma}$ are the coefficients $u_{\sigma}=\frac{1}{\sqrt{2}}e^{-\frac{\pi}{4}\sigma}$, and $\hat{O}_{\sigma}=e^{\sum_{ij}c_{i}^{\dagger}W^{(\sigma)}_{ij}c_{j}}$, where the matrices $W^{(\sigma)}_{ij}$ can be extracted from the correlation matrix $C$ defined in the previous subsection,

[TABLE]

where $C^{(\alpha\beta)}$ are the blocks of $C$ is the regions corresponding to subsystems $A_{\alpha},A_{\beta}$. In Ref. Cornfeld et al., 2018 it was shown that for a quadratic operator $\hat{X}=\sum_{ij}c_{i}^{\dagger}X_{ij}c_{j}$,

[TABLE]

where $u_{\{\sigma\}}=\prod_{i}u_{\sigma_{i}}$. The charge imbalance $\hat{Q}$ is a particular case of such a quadratic operator: $\hat{Q}=\sum_{i,j}c_{i}^{\dagger}(-1)^{\gamma(i)}\mathbb{I}_{ij}c_{j}$, where $\gamma(i)=0$ if $i\in A_{1}$, and $\gamma(i)=1$ if $i\in A_{2}$. Combining the equation above with Eq. (25), which is of the form of the RHS of Eq. (36), we can obtain the exact RNs for the XX model.

In this method, the only RNs that are numerically accessible are $\mathcal{N}_{A_{1},A_{2}}^{(1)},\mathcal{N}_{A_{1},A_{2}}^{(2)}$, and $\mathcal{N}_{A_{1},A_{2}}^{(3)}$, since the matrices in the denominator in Eq. (36) are almost singular, and only for $n\leq 3$ one can explicitly cancel the small denominator against a corresponding factor in the numerator.

For nonzero values of $\Delta$ we used the tDMRG algorithm and extracted the spectrum of $\rho_{A}^{T_{2}}$ following the method described in Ref. Ruggiero et al., 2016. A naive method for extracting the negativity of an MPS goes as follows: We contract all of the tensors representing one subsystem $A_{i}$ into one big tensor, and combine all the physical indices to one physical index of degree $2^{|A_{i}|}$. We then explicitly create $\rho_{A}^{T_{2}}$ by taking two copies of each rank-3 tensor, tracing out the environment, and applying a partial transpose to the physical indices of subsystem $A_{2}$, see Fig. 5. This partial transpose guarantees that the block form of the matrix, formerly corresponding to the charge on both subsystems, will now correspond to the charge imbalance, as illustrated in Fig. 5d. The extraction of the tensor in this method for an MPS bond dimension $D$ is $D^{8}$. In Ref. Ruggiero et al., 2016 a more sophisticated way of constructing an equivalent tensor is presented in complexity $D^{6}$, which has the same charge imbalance block structure. We use the latter method in our calculations.

In this method all negativities are accessible, and we can use the block structure just explained for extracting the charge imbalance resolved negativities. However, since the dependence of the runtime on the bond dimension is stronger for the extraction of the negativity as compared to the entropy, we reduced the bond dimension to 256, leading to a truncation error of $\sim 10^{-5}$.

By definition, $\mathcal{N}_{A_{1},A_{2}}^{(1)}(Q)$ is simply the probability distribution of $Q$. However, it is worthwhile to fit it to the CFT prediction. The reason is that for $n=1$ we do not apply any twist fields. In this case, our vertex operator 6-point function is just an ordinary free-boson 6-point function, for which the nonuniversal function $\mathcal{F}$ should be equal to unity Francesco et al. (1997). Fig. 7 presents the fit to Eq. (29) for the XX model for a very large system using the transformation (10), and for the interacting case (using tDMRG) for a finite system using the transformation (12). The fits are seen to work well in the current $n=1$ case.

The second RN is by definition simply the second RE, or purity, of $A=A_{1}\cup A_{2}$, as mentioned in Ref. Calabrese et al., 2013. We thus skip to the first non-trivial RN, $\mathcal{N}_{A_{1},A_{2}}^{(3)}$. Fig. 8 presents comparisons of the CFT prediction to the exact ($\Delta=0$) and tDMRG ($\Delta\neq 0$) results. Here things get complicated: It looks like the function $\mathcal{F}$ modifies some of the qualitative behavior of the prediction as well. We see that while the overall size and trends of the charge-resolved negativities is generally correct for $vt\ll l$, the qualitative behavior is not always fully reproduced. We are left with the conclusion that the effect of the non-universal contribution for the vertex operators 6-point function (imbalance-resolved RNs) has a more pronounced effect than it does for the twist fields 6-point function (total RNs), similarly to the static case Wen et al. (2015).

IV Conclusions and Future Outlook

We expanded the understanding of the charge-resolved entanglement entropies, defined in Ref. Goldstein and Sela, 2018, and charge-imbalance-resolved entanglement negativities, defined in Ref. Cornfeld et al., 2018, by studying for the first time the time-dependence of these quantities following a local quench. We started by analytically calculating the time dependent charge resolved entanglement for 1+1D CFT: We added flux-like vertex operators to the time dependent replica trick following Ref. Goldstein and Sela, 2018. We then compared the prediction to exact results for the noninteracting XX model, as well as tDMRG results for the critical range of the XXZ model. While the total entanglement scales as the log of the protocol-dependent effective length $L_{\mathrm{eff}}(t)$ (Eq. (11)), we find that the contribution comes from $\propto\sqrt{\log(L_{\mathrm{eff}}(t))}$ significant charge sectors, each contributing an amount proportional to $\sqrt{\log(L_{\mathrm{eff}}(t))}$ to the entanglement (Eqs. (19, 21)). For example, this could be used to understand the number of states one needs to keep in each conserved number block in tDMRG Pollmann et al. (2009).

The dependence of the flux-resolved entanglement on the flux appears to be more complicated for smaller systems (whether the entire system is small in size, or just time is short, making the effective length small), and we needed to amend our expression by expanding the vertex operator as in Eq. (20). However, the CFT predictions for the vNEE agrees nicely with the numerical results for Luttinger liquids, both with and without taking the full expansion for the vertex operator.

For the charge imbalance resolved negativity, the time-dependent calculation forced the use of boundary CFT. In this case, high-order correlation functions appear, which are non-universal. General CFT arguments are sufficient for obtaining the qualitative behavior, but the calculation of the full nonuniversal confromal block remains an interesting question for the future.

Our results pave the way towards studying additional models, such as 1+1D conformal systems with non-abelian symmetries Goldstein and Sela (2018) or systems undergoing global quenches (based on the the appropriate conformal transformation from Ref. Calabrese and Cardy, 2009). It will also be interesting to study the behavior of charge resolved entanglement in non-CFT systems (e.g., topological ones  Kitaev and Preskill (2006); Levin and Wen (2006); Cornfeld et al. (2019)) or conformal systems of higher dimension. One possible motivation is that the charge-resolved entanglement measures are much more numerically accessible for high values of $\Delta N_{A}$ or $\Delta Q$ — these charge or charge imbalance sectors carry less entanglement and hence their corresponding blocks in $\rho_{A},\rho_{A}^{T_{2}}$ respectively are smaller and easier to diagonalize Laflorencie and Rachel (2014). Thus, if the expected $\Delta N_{A}$ or $\Delta Q$ is known, one can use calculations in these smaller sectors to characterize all sectors.

Acknowledgements.

We would like to thank G. Cohen, E. Cornfeld, E. Grosfeld, and E. Sela for useful discussions. Support by the Israel Science Foundation (Grant No. 227/15), the German Israeli Foundation (Grant No. I-1259-303.10), the US-Israel Binational Science Foundation (Grant No. 2016224), and the Israel Ministry of Science and Technology (Contract No. 3-12419) is gratefully acknowledged.

Appendix A Full Formulas for the Symmetry-Resolved Entanglement Entropies

The final formula used for fitting the numerical results is a combination of Eqs. (12), (14) and (21) from the main text. These include the effective length $L_{\mathrm{eff}}(t)$. For the entanglement entropy, we treated two cases. In the first case, where $L_{A}=L_{B}=L/2$ and the quench point is the boundary point between the two subsystems, we get from Eq. (12)

[TABLE]

When the quench point is in the middle of the full system ($x=L/2$), but the boundary point between the two subsystems is shifted from the quench point by $l$, so that $L_{A}=L/2-l$, $L_{B}=L/2+l$, we obtain, using the conformal transformation (12) again but now adding the twist field at $z=l+i\tau$,

[TABLE]

Both terms can be calculated straightforwardly by performing the conformal transformation (12) and combining it with Eq. (29) in Ref. Calabrese and Cardy, 2009.

We now substitute $L_{\mathrm{eff}}$ into Eqs. (18) and (20), and get:

[TABLE]

We can now perform a derivation by $n$ in order to get the flux-resolved vNEE:

[TABLE]

following a Fourier transform (Eq. (14)) which will gives us the charge-resolved vNEE:

[TABLE]

where $\text{erf}=\frac{2}{\sqrt{\pi}}\int_{0}^{z}e^{-t^{2}}dt$ is the error functionAbramowitz and Stegun (1964).

Appendix B Final Formulas for the Entanglement Negativity

The final formula for the flux resolved entanglement negativities, Eq. (29), involves the expressions $|\omega^{\prime}_{i}-\omega^{\prime}_{j}|$ defined in the main text. These expressions for the case of two adjoint subsystems of size $l$ (for an infinite total system) at time $t$ after a local quench were derived in Ref. Wen et al., 2015, and are reproduced here for completeness:

[TABLE]