Entanglement and quantum transport in integrable systems

TL;DR

This paper investigates entanglement dynamics after a quench in integrable systems, revealing how quasiparticles determine entanglement production and steady-state entropy, with results supported by numerical simulations.

Contribution

It introduces a conjectured formula for entanglement production rate and steady-state entropy in integrable systems after joining two reservoirs, based on the quasiparticle picture.

Findings

Entanglement production rate is linked to quasiparticles in the NESS.

Steady-state entanglement entropy matches the thermodynamic entropy of the NESS.

Numerical tDMRG simulations confirm the theoretical predictions.

Abstract

Understanding the entanglement structure of out-of-equilibrium many-body systems is a challenging yet revealing task. Here we investigate the entanglement dynamics after a quench from a piecewise homogeneous initial state in integrable systems. This is the prototypical setup for studying quantum transport, and it consists in the sudden junction of two macroscopically different and homogeneous states. By exploiting the recently developed integrable hydrodynamic approach and the quasiparticle picture for the entanglement dynamics, we conjecture a formula for the entanglement production rate after joining two semi-infinite reservoirs, as well as the steady-state entanglement entropy of a finite subregion. We show that both quantities are determined by the quasiparticles created in the Non Equilibrium steady State (NESS) appearing at large times at the interface between the two reservoirs. Specifically, the steady-state entropy coincides with the thermodynamic entropy of the NESS, whereas the entropy production rate reflects its spreading into the bulk of the two reservoirs. Our results are numerically corroborated using time-dependent Density Matrix Renormalization Group (tDMRG) simulations in the paradigmatic XXZ spin-1/2 chain.

Figures

Equations

$$S(t)\propto 2t\underset{2|v|t<\ell }{\int }d\lambda |v(\lambda )|f(\lambda )+\ell \underset{2|v|t>\ell }{\int }d\lambda f(\lambda ),$$

$$\begin{gathered} S_{YY}=s_{YY}L=L\sum _{n=1}^{\infty }\int d\lambda [{\rho }_n^{(t)}\ln {\rho }_n^{(t)} \\ -{\rho }_n\ln {\rho }_n-{\rho }_n^{(h)}\ln {\rho }_n^{(h)}] \\ \equiv L\sum _n\int d\lambda s_{YY}^{(n)}[{\rho }_n,{\rho }_n^{(h)}]. \end{gathered}$$

$$[\zeta -v_{\zeta ,n}(\lambda )]{\partial }_{\zeta }{\vartheta }_{\zeta ,n}(\lambda )=0,$$

$${\rho }_{\zeta ,n}^{(t)}(\lambda )=a_n(\lambda )+\sum _{n=1}^{\infty }(a_{n,m}\star {\rho }_{\zeta ,m})(\lambda ).$$

$${\vartheta }_{\zeta ,n}(\lambda )={\theta }_H(v_{\zeta ,n}(\lambda )-\zeta )({\vartheta }_n^B(\lambda )-{\vartheta }_n^A(\lambda ))+{\vartheta }_n^A(\lambda ),$$

$$s_{\infty }=\sum _{n=1}^{\infty }\int d\lambda s_{YY}^{(n)}[{\rho }_{\zeta =0,n},{\rho }_{\zeta =0,n}^{(h)}],$$

$$s^{\prime }=\sum _{n=1}^{\infty }\int d\lambda |v_{\zeta =0,n}|s_{YY}^{(n)}[{\rho }_{\zeta =0,n},{\rho }_{\zeta =0,n}^{(h)}].$$

$$H_{XXZ}=J\sum _{i=1}^L[S_i^x{S}_{i+1}^x+S_i^y{S}_{i+1}^y+\Delta S_i^z{S}_{i+1}^z].$$

$${\left[\frac{\sin ({\lambda }_j+i\frac{\eta }{2})}{\sin ({\lambda }_j-i\frac{\eta }{2})}\right]}^L=-\prod _{k=1}^M\frac{\sin ({\lambda }_j-{\lambda }_k+i\eta )}{\sin ({\lambda }_j-{\lambda }_k-i\eta )},$$

$${\lambda }_{n,\gamma }^j={\lambda }_{n,\gamma }+i\frac{\eta }{2}(n+1-2j)+{\delta }_{n,\gamma }^j,$$

$$L{\theta }_n({\lambda }_{n,\alpha })=2\pi I_{n,\alpha }+\sum _{(n,\alpha )\ne (m,\beta )}{\Theta }_{n,m}({\lambda }_{n,\alpha }-{\lambda }_{m,\beta }).$$

$$\begin{gathered} {\Theta }_{n,m}(\lambda )\equiv (1-{\delta }_{n,m}){\theta }_{|n-m|}(\lambda )+2{\theta }_{|n-m|+2}(\lambda ) \\ +\cdots +{\theta }_{n+m-2}(\lambda )+{\theta }_{n+m}(\lambda ). \end{gathered}$$

$${ϵ}_n(\lambda )\equiv -\frac{\sinh (\eta )\sinh (n\eta )}{\cosh (n\eta )-\cos (2\lambda )},z_n({\lambda }_{n,\alpha })=\frac{2\pi I_{n,\alpha }}{L}.$$

$${\rho }_{n,p}(\lambda )\equiv \underset{L\to \infty }{\lim }\frac{1}{L({\lambda }_{n,\alpha +1}-{\lambda }_{n,\alpha })}.$$

$${\rho }_n^{(t)}(\lambda )=a_n(\lambda )+\sum _{n=1}^{\infty }(a_{n,m}\star {\rho }_m)(\lambda ),$$

$$a_{nm}(\lambda )=$$

$$+\dots +2a_{n+m-2}(\lambda )+a_{n+m}(\lambda ),$$

$$a_n(\lambda )=\frac{1}{\pi }\frac{\sinh \left(n\eta \right)}{\cosh (n\eta )-\cos (2\lambda )}.$$

$${\eta }_n(\lambda )=\frac{{\eta }_{n-1}(\lambda +i\frac{\eta }{2}){\eta }_{n-1}(\lambda -i\frac{\eta }{2})}{1+{\eta }_{n-2}(\lambda )}-1,$$

$${\rho }_{n,h}(\lambda )={\rho }_{n,t}(\lambda +i\frac{\eta }{2})+{\rho }_{n,t}(\lambda )-{\rho }_{n-1,h}(\lambda ),$$

$${\eta }_1=\frac{2[2\cosh (\eta )+2\cosh (3\eta )-3\cos (2\lambda ){\sin }^2(\lambda )]}{[\cosh (\eta )-\cos (2\lambda )][\cosh (4\eta )-\cos (4\lambda )]}$$

$$\displaystyle\rho_{1,h}=a_{1}\Big{(}1-\frac{\cosh^{2}(\eta)}{a_{1}\pi^{2}\sin^{2}(2\lambda)+\cosh^{2}(\eta)}\Big{)},$$

$${\eta }_1=\frac{\cos (4\lambda )-\cosh (2\eta )}{{\cos }^2\lambda (\cos (2\lambda )-\cosh (2\eta ))}-1,$$

$${\rho }_{1,h}=a_1(\lambda )+\frac{1}{2\pi }(\omega (\lambda -i\eta /2)+\omega (\lambda +i\eta /2)),$$

$$\begin{gathered} \omega (\lambda )=-\frac{\sinh (\eta )(-2+\cosh (\eta )+2\cosh (2\eta )}{4(\cos (2\lambda )-\cosh (2\eta ){)}^2} \\ +\frac{3\cosh (3\eta )+4\cos (2\lambda )(-\cosh (\eta )+{\sinh }^2(\eta ))}{4(\cos (2\lambda )-\cosh (2\eta ){)}^2}. \end{gathered}$$

$$\delta E_n=e_n({\lambda }_{n,p}^{\ast })-e_n({\lambda }_{n,h}^{\ast }).$$

$$\delta P=z_n({\lambda }_{n,p}^{\ast })-z_n({\lambda }_{n,h}^{\ast }).$$

$$v_n^{\ast }(\lambda )\equiv \frac{\partial e_n}{\partial z_n}=\frac{e_n^{\prime }(\lambda )}{2\pi {\rho }_n(1+{\eta }_n(\lambda ))},$$

$$e_n^{\prime }(\lambda )+\frac{1}{2\pi }\sum _{m=1}^{\infty }\int d\mu e_m^{\prime }(\mu )\frac{{\Theta }_{m,n}^{\prime }(\mu -\lambda )}{1+{\eta }_m^{\ast }(\mu )}={ϵ}_n^{\prime }(\lambda ),$$

$$S^z\equiv S_i^z$$

$$E\equiv S_i^x{S}_{i+1}^x+S_i^y{S}_{i+1}^y+\Delta S_i^z{S}_{i+1}^z-\frac{\Delta }{4},$$

$$J_E\equiv S_{i-1}^x{S}_i^z{S}_{i+1}^y-S_{i-1}^y{S}_i^z{S}_{i+1}^x$$

$$-\Delta S_{i-1}^z{S}_i^x{S}_{i+1}^y+\Delta S_{i-1}^z{S}_i^y{S}_{i+1}^x$$

Full text

I Introduction

The quest for a complete understanding of entanglement spreading in out-of-equilibrium many-body system is a fruitful research theme. Recently, this became experimentally relevant, as it is now possible to measure the entanglement dynamics with cold atoms islam-2015 ; kaufman-2016 . The best-known entanglement diagnostic tool is the von Neumann (entanglement) entropy $S\equiv-\textrm{Tr}\rho_{A}\ln\rho_{A}$, where $\rho_{A}$ is the reduced density matrix of a subsystem $A$ (see Figure 1 (a) for a one-dimensional setup).

A prominent out-of-equilibrium situation is that of the quench from two homogeneous initial states. This is a particular instance of inhomogeneous quantum quenches. The prototypical setup for spin chains is depicted in Figure 1 (a), and it consists in the sudden junction of two chains $A$ and $B$ that are prepared in two homogeneous macroscopically different quantum states $|\Psi_{A}\rangle$ and $|\Psi_{B}\rangle$. Subsequently, the state $|\Psi_{A}\rangle\otimes|\Psi_{B}\rangle$ is evolved in real time using a quantum many-body Hamiltonian $H$. Typically, a non-zero current arises between the two chains. The main focus so far has been on transport of local quantities, such as the local energy and local magnetization. Several techniques have been used, such as Conformal Field Theory spyros-2008 ; bernard-2012 ; bhaseen-2015 ; allegra-2016 ; dubail-2017 ; dubail-2017a (CFT), free-fermion methods eisler-2009 ; de-luca-2013 ; sabetta-2013 ; eisler-2013 ; alba-2014 ; collura-2014 ; de-luca-2016 ; eisler-2016 ; viti-2016 ; kormos-2017 ; perfetto-2017 ; vidmar-2017 , field theory methods de-luca-2014 ; olalla-2014 ; biella-2016 , integrability zotos-1997 ; zotos-1999 ; prosen-2011 ; prosen-2013 ; ilievski-2015 , and numerical techniques fabian-2003 ; gobert-2005 ; langer-2009 ; karrasch-2012 ; karrasch-2013 ; vidmar-2017 . For integrable models a recent breakthrough olalla-2016 ; bertini-2016 allows for an analytic treatment of transport problems using Thermodynamic Bethe Ansatz (TBA) techniques doyon-2016 ; yoshimura-2016 ; de-luca-2016a ; doyon-2017 ; doyon-2017a ; doyon-2017b ; doyon-2017c ; bulchandani-2017 ; bulchandani-2017a ; ilievski-2017 ; doyon-2017d .

Surprisingly, as of now there are no exact results for the entanglement dynamics after inhomogeneous quenches. Notable exceptions are systems that can be mapped to CFTs in curved spacetime allegra-2016 ; dubail-2017 ; dubail-2017a . In this case it is established that $S$ grows logarithmically after the quench dubail-2017 .

In contrast, for homogeneous quenches a well-known quasiparticle picture calabrese-2005 allows one to understand the entanglement dynamics in terms of quasiparticles traveling ballistically, with opposite quasimomenta, through the system. The underlying idea is that only quasiparticles created at the same point in space are entangled. At a generic time $t$ the entanglement entropy $S(t)$ between a subregion $A$ and the rest is proportional to the number of quasiparticles emitted at the same point in space at $t=0$, and being at time $t$ one in subsystem $A$ and the other in $B$. More quantitatively,

[TABLE]

where $f(\lambda)$ depends on the cross section for creating quasiparticles with quasimomentum $\lambda$, and $v(\lambda)$ is their velocity. Eq. (1) holds in the space-time scaling limit, i.e., $\ell,t\to\infty$ with fixed $t/\ell$. In many physical situations a maximum velocity $v_{M}$ exists, for instance, due to the Lieb-Robinson bound lieb-1972 . Then, Eq. (1) predicts a linear growth for $t\leq\ell/(2v_{M})$, followed by a volume-law $S\propto\ell$ at longer times. The validity of (1) has been proven rigorously only for free-fermion models fagotti-2008 ; ep-08 ; nr-14 ; coser-2014 ; cotler-2016 ; buyskikh-2016 , for which $f(\lambda)$ and $v(\lambda)$ can be determined analytically. The quasiparticle picture has been also confirmed in several numerical studies de-chiara-2006 ; lauchli-2008 ; kim-2013 , and using holographic methods hrt-07 ; aal-10 ; aj-10 ; allais-2012 ; callan-2012 ; ls-14 ; bala-2011 ; liu-2014 . Violations of the quasiparticle picture have been observed in CFT with large central charge bala-2011 ; ab-13 ; asplund-2015 ; lm-15 ; kundu-2017 ; sagar-2017 . Remarkably, a framework to render (1) predictive for generic integrable systems has been developed in Ref. alba-2016, (see Ref. mestyan-2017, for an application to the Hubbard chain). Specifically, for integrable systems $f(\lambda)$ coincides with the thermodynamic entropy associated with the post-quench steady state, while the quasiparticle velocities $v(\lambda)$ are those of the low-lying excitations around it.

Here by combining a recent hydrodynamic approach for integrable systems olalla-2016 ; bertini-2016 with the quasiparticle picture (1), we provide the first exact results for the entanglement dynamics after a global quench from a piecewise homogeneous initial state in generic integrable systems. Specifically, we focus on the steady-state entanglement entropy of a finite region $A^{\prime}$ (see Figure 1 (b)), as well as the entanglement production rate between two semi-infinite reservoirs. Our main result is that both are determined by the physics of the Non-Equilibrium-Steady-State (NESS) (ray with $\zeta=0$ in Figure 1 (a)) appearing at the interface between the two reservoirs. The steady-state entanglement entropy coincides with the NESS thermodynamic entropy, whereas the entanglement production rate reflects the spreading of the NESS from the interface into the two reservoirs.

The latter is relevant for quantifying the rate at which quantum information spreads. The study of the quantum information spreading is attracting enormous attention recently. For instane, in the holographic community, the focus has been on the so-called entanglement velocity $v_{E}$, which is defined as hol-vel $dS(t)/dt=v_{E}s_{th}\partial A$, with $s_{th}$ the thermodynamic entropy of the steady state, and $\partial A$ the length of the boundary between the subsystem of interest and the environment. For homogeneous quenches, $v_{E}$ has been investigated using holographic methods hol-vel , exact calculations in free models cotler-2016 ; ff-vel , and recently in the dynamics obtained from random unitary gates ru-vel . For inhomogeneous quenches, the information spreading started to be investigated only recently erdmenger-2017 .

II Entanglement via integrable hydrodynamics

The first key ingredient to derive our results is that the spectrum of integrable models exhibits families of stable quasiparticles. Typically, these are composite objects of elementary excitations. For spin chains, they correspond to bound states of magnons. The possible set of quasimomenta (rapidities) $\lambda$ that can be assigned to the quasiparticles are obtained by solving the so-called Bethe equations taka-book (see also Appendix A). In the thermodynamic limit the rapidities form a continuum. Thermodynamic properties of integrable models are described by the particle densities $\rho_{n}(\lambda)$ and the hole densities $\rho^{\scriptscriptstyle(h)}_{n}(\lambda)$. The latter is the density of unoccupied rapidities. The total density is defined as $\rho^{\scriptscriptstyle(t)}_{n}\equiv\rho_{n}+\rho^{\scriptscriptstyle(h)}_{n}$. Here the subscript $n$ labels different quasiparticle families, and for spin chains is the size of the bound states. For free models in terms of quasimomenta $\rho_{n}^{\scriptscriptstyle(t)}=const$, reflecting that the quasimomenta are equally-spaced. Every set of $\rho_{n},\rho^{\scriptscriptstyle(h)}_{n}$ can be interpreted as a thermodynamic macrostate, which corresponds to an exponentially large number of microscopic eigenstates of the model. Their number is given as $e^{S_{YY}}$, with $S_{YY}$ the so-called Yang-Yang entropy yang-1969

[TABLE]

Similar to free models, $S_{YY}$ counts the number of ways of assigning the different rapidities $\lambda$ to the quasiparticles, compatibly with the densities $\rho_{n},\rho_{n}^{\scriptscriptstyle(h)}$.

The second key ingredient in our approach is the integrable hydrodynamics framework olalla-2016 ; bertini-2016 . Due to integrability, information spreads ballistically from the interface between $A$ and $B$. As a consequence, physical observables depend only on the combination $\zeta\equiv x/t$ (see Figure 1 (a)), with $x$ the distance from the interface between $A$ and $B$. For each fixed $\zeta$, dynamical properties of local and quasilocal observables are described by a thermodynamic macrostate, i.e., a set of densities $\rho_{\zeta,n},\rho_{\zeta,n}^{\scriptscriptstyle(h)}$. Each macrostate identifies a different Generalized Gibbs Ensemble rigol-2008 ; polkovnikov-2011 ; calabrese-2005 ; rigol-2007 ; cazalilla-2006 ; barthel-2008 ; cramer-2008 ; cramer-2010 ; calabrese-2011 ; cazalilla-2012a ; calabrese-2012 ; sotiriadis-2012 ; collura-2013 ; collura-2013a ; fagotti-2013 ; fagotti-2014 ; kcc14 ; delfino-2014 ; sotiriadis-2014 ; ilievski-2015a ; alba-2015 ; essler-2015 ; cardy-2015 ; cardy-2015 ; sotiriadis-2016 ; bastianello-2016 ; vernier-2016 ; vidmar-2016 ; gogolin-2015 ; essler-2016 ; calabrese-2016 ; piroli-2017 ; piroli-2017a . The macrostate with $\zeta=0$ is known as Non Equilibrium Steady State (NESS). The central result of Ref. olalla-2016, ; bertini-2016, is that $\rho_{\zeta,n}$ satisfy the continuity equation

[TABLE]

where $\vartheta_{\zeta,n}\equiv\rho_{\zeta,n}/\rho_{\zeta,n}^{\scriptscriptstyle(t)}$, together with the standard TBA equations taka-book

[TABLE]

Here $a_{n}$ and $a_{n,m}$ are known functions of $\lambda$ (see Appendix A for their expression for the XXZ chain), and the star symbols denotes the convolution $f\star g\equiv\int d\mu f(\lambda-\mu)g(\mu)$. Crucially, in (3) $v_{\zeta,n}$ are the velocities of the low-lying excitations around the macrostate $\rho_{\zeta,n}$, and fully encode the interactions (scatterings) between quasiparticles bonnes-2014 . They can be calculated using standard TBA techniques bonnes-2014 (see Appendix C). The solutions of (3) can be conveniently written as bertini-2016

[TABLE]

with $\theta_{H}(x)$ the Heaviside function, and $\vartheta_{\zeta,n}^{A(B)}$ the densities describing the steady states arising after the homogeneous quenches with initial states $|\Psi_{A}\rangle$ and $|\Psi_{B}\rangle$, respectively. Exact results for $\vartheta_{n}^{\scriptscriptstyle A(B)}$ are available for several quenches caux-2013 ; wouters-2014A ; pozsgay-2014A ; bse-14 ; dwbc-14 ; bucciantini-15 ; ac-16 ; pce-16 ; piroli-2016 ; mestyan-2015 ; brockmann-2014 ; iqdb-15 ; piroli-2016a ; bpc-16 ; caux-2016 (see also Appendix B).

We now present our main results. We start discussing the steady-state entanglement entropy of a finite region $A^{\prime}$ of length $\ell$ embedded in part $A$ and placed next to the interface with $B$ (see Figure 1 (b)). The steady-state entropy is obtained in the limit $\ell/t\to 0$, which identifies the macrostate with $\zeta=0$ (NESS). Since the spatial extension of the region described by this macrostate increases linearly with time, for $t\gg\ell$ the $\zeta=0$ macrostate is expected to describe the entire subsystem $A^{\prime}$. Similar to Ref. alba-2016, , it is natural to conjecture that the entanglement entropy density $s_{\infty}\equiv\lim_{t\to\infty}S(t)/\ell$ becomes that of the Yang-Yang entropy of the NESS. The latter is the thermodynamic entropy of the GGE describing (quasi) local observables in space-time regions with $x/t\to 0$. Using (2), this implies

[TABLE]

where $s^{\scriptscriptstyle(n)}_{YY}$ is calculated using the solutions of (3)(4).

We now turn to the entanglement production rate. For homogeneous quenches, this corresponds to the slope of the linear term in  (1). Physically, Eq. (1) means that $S(t)$ is determined by the total number of quasiparticles that at time $t$ crossed the interface between $A$ and $B$. It is natural to assume that the same applies to the inhomogeneous case. After the quench, the different lightcones associated with different $\zeta$s start spreading from the interface between $A$ and $B$. For short times, all the types of quasiparticles remain confined within $A$, the only ones crossing the boundary between $A$ and the rest are the ones described by the $\zeta=0$ macrostate. At generic time $t$ their number is proportional to the width of the associated lightcone (see Figure 1 (b)). Equivalently, the entanglement growth at short times reflects the spreading of the NESS. Since the lightcone width increases linearly with time, one should expect the linear behavior $S(t)\propto s^{\prime}t$, with $s^{\prime}\equiv dS(t)/dt$ the entanglement production rate given as

[TABLE]

Formally, Eq. (7) is the same as that for the homogeneous quench conjectured in Ref. alba-2016, . However, here the quasiparticle velocities are not odd under parity, i.e., $v_{n}(\lambda)\neq-v_{n}(-\lambda)$, implying that the lightcone is not symmetric (in Figure 1(b) $v_{n}^{\scriptscriptstyle(A)}$ and $v_{n}^{\scriptscriptstyle(B)}$ denote the different quasiparticle velocities in the two reservoirs). Similar to the homogeneous case the velocity of the entangling particles depends only on the local equilibrium state at large times alba-2016 (here the NESS). Finally, we should mention that obtaining the full-time dynamics of the entanglement entropy is an arduous task, unlike for homogeneous quenches. The reason is that it requires reconstructing the quasiparticles trajectories. Precisely, since the quasiparticles velocity $v$ is finite, the typical time $t^{*}$ for the quasiparticles to travel through the subsystem is $t^{*}\sim\ell$, with $\ell$ the subsystem size. This implies that these quasiparticles should be described by $\zeta=\ell/t^{*}={\mathcal{O}}(1)$. A possible direction to determine the quasiparticle trajectories is to use the approach of Ref. doyon-corr, .

In the following we provide numerical evidence for (6) and (7) considering the anisotropic spin-$1/2$ Heisenberg chain (XXZ chain). The Hamiltonian reads

[TABLE]

Here $S_{i}^{x,y,z}\equiv\sigma_{i}^{x,y,z}/2$ are spin-$1/2$ operators, $L$ the chain length and $\Delta$ the anisotropy parameter. We set $J=1$ in (8), restricting ourselves to $\Delta>1$. The XXZ chain is the paradigm of Bethe ansatz solvable models taka-book (see Appendix A).

We consider the inhomogeneous quenches in which parts $A$ or $B$ are prepared in the Néel state $\left|N\right\rangle\equiv(\left|\uparrow\downarrow\uparrow\cdots\right\rangle+\left|\downarrow\uparrow\downarrow\cdots\right\rangle)/\sqrt{2}$, the Majumdar-Ghosh or dimer state $\left|MG\right\rangle\equiv[(\left|\uparrow\right\rangle+\left|\downarrow\right\rangle)/\sqrt{2}]^{L/2}$, and the ferromagnetic state $|F\rangle\equiv\left|\uparrow\uparrow\cdots\right\rangle$. For all these cases the bulk densities $\vartheta_{n}^{\scriptscriptstyle A(B)}$ (cf. (5)) are known analytically (see Appendix B). Also, for $\Delta>1$ the dynamics from states obtained by joining states with opposite magnetization piroli-2017a leads to diffusive or subdiffusive transport ljubotina-2017 ; misguich-2017 ; stephan-2017 . In contrast, all the initial states considered here lead to ballistic transport (see Appendix D).

The strategy to use (6) and (7) is to first solve the coupled systems of integral equations (3) and (4) for $\zeta=0$. The obtained densities $\rho_{\zeta=0,n},\rho_{\zeta=0,n}^{\scriptscriptstyle(h)}$ (some numerical results are shown in Appendix E) are then substituted in (2) (6) and (7). Our results are summarized in Figure 2. The Figure shows the steady-state entropy $s_{\infty}$ (dotted lines) and the entanglement production rate $s^{\prime}$ as a function of $\Delta$. For the quench with initial state $|N\rangle\otimes|F\rangle$, both $s_{\infty}$ and $s^{\prime}$ vanish for $\Delta\to\infty$, reflecting that the Néel state is the ground state of the XXZ chain in that limit. Interestingly, for both $|N\rangle\otimes|F\rangle$ and $|MG\rangle\otimes|F\rangle$, since the ferromagnet is an exact eigenstate of the XXZ chain at any $\Delta$, no quasiparticle production happens in subsystem $A$. As a consequence, $S(t)$ is fully determined by the quasiparticle transport from $B$ to $A$. However, $s_{\infty}$ and $s^{\prime}(t)$ are smaller than the corresponding values for the homogeneous quench from the Néel state and the dimer state alba-2016 . It is also interesting to observe that for these quenches only the quasiparticles with $n=1$ contribute in (6) and (7), i.e., the bound states contribution to the entanglement vanishes (see Appendix E). This is not the case for the quench from $|MG\rangle\otimes|F\rangle$, for which all bound states contribute and quasiparticles are generated in both reservoirs. In this case both $s_{\infty}$ and $s^{\prime}$ exhibit a rather weak dependence on $\Delta$ (see Figure 2).

III Numerical checks in the XXZ chain

We now turn to verify the theoretical predictions presented in Figure 2 using tDMRG simulations white-2004 ; daley-2004 ; uli-2005 ; uli-2011 ; itensor . We first focus on the steady-state entropy density $s_{\infty}$ of a finite block $A^{\prime}$ of length $\ell$ (Figure 1 (b)). We consider only the quenches from the states $|N\rangle\otimes|F\rangle$ and $|MG\rangle\otimes|F\rangle$, as they are easier to simulate. In fact, in our tDMRG simulations we use the comparatively small value of the bond dimension $\chi\approx 75$. This allows us to obtain reliable results, provided that a space-time average of the data is performaed This surprising efficiency of tDMRG for transport-related quantities has been also observed in Ref. leviatan-2017, .

Our tDMRG results in the regime $t\gg\ell$ are reported in Figure 3 (a) and (b), plotting $s_{\infty}$ versus $1/\ell$. The raw tDMRG data at any time after the quench are reported in Appendix F. To minimize the effects of oscillating (with the block size) scaling corrections we averaged the data for $t\gg\ell$. The results are for chains with $L=40$ sites, severals $\Delta$s (different symbols), and $\ell\leq 10$. The star symbols are the Bethe ansatz results (6). The dash-dotted lines are fits to $S/\ell=s_{\infty}+a/\ell+b/\ell^{2}$, with $s_{\infty}$ the thermodynamic limit results, and $a,b$ fitting parameters. Finite-size corrections are clearly visible for small $\ell$. However, for both initial states the numerical data are compatible in the thermodynamic limit with the Bethe ansatz results. Note that for the quench from $|MG\rangle\otimes|F\rangle$ the dependence on $\Delta$ for large $\ell$ is not visible within the numerical precision, as expected from Figure 2.

We now focus on the entanglement production rate in Figure 4, again considering the quench from $|N\rangle\otimes|F\rangle$ (panels (a)-(c)) and $|MG\rangle\otimes|F\rangle$ (panels (d)-(e)), plotting $S/\ell$ versus $t/\ell$. The data are now for chains with $L\leq 128$. We always consider the half-chain entropy, i.e., $A^{\prime}=A$ and $A$ half of the chain (see Figure 1). A crucial consequence is that only one boundary is present between $A^{\prime}$ and the rest. The dash-dotted lines are the Bethe ansatz predictions (Figure 2 and (7)). After an initial transient, in all cases $S(t)$ exhibits linear behavior. The slope of this linear increase is in agreement with the Bethe ansatz results, already for $\ell\approx 16$. A more systematic analysis of finite-size corrections is presented in Appendix G. For the quench from $|MG\rangle\otimes|F\rangle$, again, one should observe the weak dependence on $\Delta$. Finally, in Figure 5 we show the results for the initial state $|N\rangle\otimes|MG\rangle$. Due to the large amount of entanglement, we can only provide reliable tDMRG data for $\ell\leq 24$. tDMRG simulations are performed with $\chi=400$. To highlight finite-size and finite-time corrections we show data for $\ell=8$, which exhibit large deviations from the Bethe ansatz predictions. On the other hand, for $\ell=24$ the data are clearly compatible with (7).

IV Conclusions

We investigated the entanglement dynamics after the quenche from a piecewise initial state in integrable models. We conjectured an analytic formula for the entanglement production rate after joining two semi-infinite reservoirs as well as for the steady-state entropy of a finite subregion ((6) and (7), respectively). Our work opens several promising new directions. For instance, it would be interesting to consider the inhomogeneous quench in which an integrability-breaking term acting at the interface between $A$ and $B$ is added to the Hamiltonian. Already in the case of a defect that preserves integrability this gives non trivial effects bertini-2016 . It would be also enlightening to consider a time-dependent bipartition in which subsystem $A^{\prime}$ moves away from the boundary. By appropriately tuning the speed at which $A^{\prime}$ changes its position, it should be possible to change the thermodynamic macrostate governing the entanglement production. It should be possible to extend the method to treat the mutual information, as done for free fermions eisler-2014 ; kormos-2017a . Finally, it would be interesting to extend the quasiparticle picture to treat multipartite systems as in Ref. savona .

Acknowledgements.

I am very grateful to Bruno Bertini for clarifying discussions on several aspects of the integrable hydrodynamic approach. I also acknowledge very fruitful discussions with Maurizio Fagotti and Pasquale Calabrese. This work was supported by the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowoska-Curie grant agreement No 702612 OEMBS.

Appendix A Bethe ansatz solution of the XXZ chain

Due to the conservation of the total magnetization $S_{z}$, the eigenstates of the XXZ can be classified according to the total magnetization $S_{T}=\sum_{i}S_{i}^{z}$. Equivalently, one can use the total number $M$ of down spins as a good quantum number for the eigenstates.

In the Bethe ansatz taka-book solution of the $XXZ$ chain, the eigenstates in the sector with $M$ down spins (particles) are identified by a set of $M$ rapidities $\lambda_{j}$, which are solutions of the so-called Bethe equations taka-book

[TABLE]

where $\eta\equiv\textrm{arccosh}(\Delta)$. Here we are interested in the thermodynamic limit $L,M\to\infty$ with $M/L$ fixed. Then, the solutions of the Bethe equations (9) form string patterns in the complex plane. Rapidities forming a $n$-string can be written as

[TABLE]

where $j=1,\dots,n$ denotes the different string components, $\lambda_{n,\gamma}$ is the “string center”, and $\delta_{n,\gamma}^{j}$ are the string deviations. For the majority of the eigenstates of the XXZ chain, the string deviations are exponentially small, i.e., $\delta_{n,\gamma}^{j}={\mathcal{O}}(e^{-L})$ (string hypothesis). The $n$-strings describe bound states of $n$ down spins. The string centers $\lambda_{n,\gamma}$ are obtained by solving the Bethe-Gaudin-Takahashi (BGT) equations taka-book

[TABLE]

For $\Delta>1$, one has $\lambda_{n,\gamma}\in[-\pi/2,\pi/2)$. Here we used that $\theta_{n}(\lambda)\equiv 2\arctan[\tan(\lambda)/\tanh(n\eta/2)]$. The scattering phases $\Theta_{n,m}(\lambda)$ are given as

[TABLE]

Each choice of the so-called BGT quantum numbers $I_{n,\alpha}\in\frac{1}{2}\mathbb{Z}$ corresponds to a different set of solutions of (11), i.e., to a different eigenstate of the chain. The corresponding eigenstate energy $E$ and total momentum $P$ are obtained by summing over all the rapidities as $E=\sum_{n,\alpha}\epsilon_{n}(\lambda_{n,\alpha})$, and $P=\sum_{n,\alpha}z_{n}(\lambda_{n,\alpha})$ with

[TABLE]

In the thermodynamic limit, one works with the rapidity densities. The root densities $\rho_{n,p}$ are formally defined as

[TABLE]

It is also convenient to define the associated hole densities $\rho_{n,h}$, i.e., the density of unoccupied rapidities, and the total densities $\rho_{n,t}=\rho_{n}+\rho_{n,h}$. The BGT equations (11) in the thermodynamic limit become a system of integral equations

[TABLE]

where we defined

[TABLE]

with

[TABLE]

Appendix B Macrostates for homogeneous quenches

Here we report the analytical results for the root densities $\vartheta^{\scriptscriptstyle(A)}_{n}(\lambda)$ and $\vartheta^{\scriptscriptstyle(B)}_{n}(\lambda)$ (see the main text). We consider the homogeneous quenches from the Néel state and the Majumdar-Ghosh state. These densities characterize the thermodynamic macrostate in the bulk of the two systems $A$ and $B$, i.e., for $|\zeta|\to\infty$. We first define the densities $\eta_{n}\equiv\rho_{n,h}/\rho_{n}$. In terms of $\eta_{n}$ one has $\theta_{n}=1/(1+\eta_{n})$.

For both the Néel state and the Majumdar-Ghosh state, the $\eta_{n}$ obey the recursive equation ilievski-2015a ; brockmann-2014

[TABLE]

where $\eta_{0}=0$ and $\rho_{0,h}=0$.

For the Néel state one has brockmann-2014

[TABLE]

For the Majumdar-Ghosh one has pozsgay-2014A ; mestyan-2015

[TABLE]

and

[TABLE]

with

[TABLE]

Appendix C Velocities of entangling quasiparticles

A crucial ingredient in the integrable hydrodynamic approach olalla-2016 ; bertini-2016 and in the derivation of our results is the velocity of the low-lying excitations around the TBA macrostates $\rho_{n},\rho_{n}^{\scriptscriptstyle(h)}$. Here we outline its derivation following the approach of Ref. bonnes-2014, . Given a generic thermodynamic macrostate identified by some densities $\rho_{n},\rho_{n}^{\scriptscriptstyle(h)}$, one can imagine of choosing among the eigenstates of the XXZ chain a representative of the macrostate. This would be identified by some BGT quantum numbers $I_{n,\alpha}^{*}$. Low-lying excitations around it can be constructed as particle-hole excitations. In each $n$-string sector, these correspond to the change $I^{*}_{n,h}\to I_{n,p}$, where $I_{n,p}(I^{*}_{n,h})$ is the BGT number of the added particle (hole). Since the model is interacting, this local change in quantum numbers affects all the rapidities obtained by solving the new set of Bethe equations. The excess energy of the particle-hole excitation can be written as bonnes-2014

[TABLE]

Note that (24) is the same as for free models apart from the dressing of the single particle energy. The change in the total momentum is obtained from (13) as

[TABLE]

The group velocity associated with the particle-hole excitations is by definition

[TABLE]

where we used that taka-book $dz_{n}(\lambda)/d\lambda=2\pi\rho_{n,t}$, and we defined $e^{\prime}_{n}(\lambda)=de(\lambda)/d\lambda$. Here $e^{\prime}_{n}(\lambda)$ is obtained by solving an infinite system of Fredholm integral equations of the second kind bonnes-2014

[TABLE]

where $\Theta^{\prime}_{n,m}(\lambda)\equiv d\Theta_{n,m}(\lambda)/d\lambda$ and $\epsilon^{\prime}_{n}(\lambda)\equiv d\epsilon_{n}(\lambda)/d\lambda$ (cf. also (12) and (13)).

Appendix D Integrable hydrodynamics: Conserved charges and currents

In this section we provide some numerical checks of the integrable hydrodynamics approach olalla-2016 ; bertini-2016 , focusing on the XXZ chain in the region with $\Delta>1$. We consider the quenches with initial states $|N\rangle\otimes|F\rangle$ and $|MG\rangle\otimes|F\rangle$. Using tDMRG simulations we investigate the dynamics of the local magnetization $S^{z}$, the local energy density $E$, and the energy current $J_{E}$. These are defined as bertini-2016

[TABLE]

In the framework of the Bethe ansatz, these are written in terms of the root densities $\rho_{\zeta,n}(\lambda)$ as

[TABLE]

Note that all the observables are functions of $\zeta\equiv x/t$. The dependence on $\zeta$ is encoded in the densities $\rho_{\zeta,n}$ bertini-2016 .

The Bethe ansatz results (31)(32)(33) are compared with tDMRG results in Figure 6. The upper and lower panels show tDMRG data for the quench with initial states $|N\rangle\otimes|F\rangle$ and $|MG\rangle\otimes|F\rangle$, respectively. The data are for the XXZ chain with $L=300$ sites, up to times $t\lesssim 100$. The results in panels (a-c) are obtained using bond dimension $\chi\sim 20$. For the quench from the state $|MG\rangle\otimes|F\rangle$ (panels (e-f)) we used $\chi\sim 75$. In all panels the different symbols correspond to different values of $\Delta$. In order to remove spatial and temporal oscillations, we performed a spatio-temporal average. Specifically, for each fixed $\zeta$ the results in the Figure are obtained by averaging the data in a window $\zeta\pm\epsilon$ with $\epsilon\approx 0.1$. The dash-dotted lines are the theoretical results (31)(32)(33). For $|\zeta|\gg 1$ all the observables become $\zeta$ independent. This happens in spatio-temporal regions with $|x|\gg v_{M}t$, with $v_{M}$ the maximum velocity in the system. Remarkably, for both quenches, and for all values of $\Delta$, the tDMRG data are in very good agreement with the Bethe ansatz. Note, however, that the theoretical predictions exhibits some cusp-like behaviors (see for instance the arrow in panel (b)). These are expected and have been discussed in Ref. bb, , although they are not well reproduced by the tDMRG results. It is natural to expect that these deviations should be attributed to the finite $\chi$ and to the finite-time and finite-size effects.

Appendix E Neel-Ferro quench: Bethe ansatz results

In this section we provide some details on the Bethe ansatz results for the quench from the state $|N\rangle\otimes|F\rangle$. Our results are summarized in Figure 9.

Specifically, panels (a) and (b) show the densities $\vartheta_{\zeta,n}$ and $\rho_{\zeta,n}$, respectively. The results are for $\zeta=0$, which identifies the relevant macrostate to describe the entanglement dynamics after the quench (the subscript $\zeta$ in $\vartheta_{\zeta,n}$ and $\rho_{\zeta,n}$ is omitted in the Figure). First, it is interesting to observe that only $\vartheta_{\zeta,1}$ and $\rho_{\zeta,1}$ are nonzero. Moreover, $\vartheta_{1}$ is nonzero only for $\lambda>0$. Similar behavior is observed for $\rho_{1}$ (see panel (b)).

The group velocities of the low-lying excitations around the macrostate with $\zeta=0$ are reported in Figure 9 (c). These are obtained by numerically solving (27) and using (26). As anticipated in the main text, one has $v_{n}(\lambda)\neq-v_{n}(-\lambda)$. This is in stark contrast with homogeneous quenches, where one has alba-2016 $v_{n}(\lambda)=-v_{n}(-\lambda)$. As a consequence the lightcone of the entangling quasiparticles is not symmetric. It is also interesting to observe that for any $\lambda$, $v_{n}<0$ for $n>1$, which implies that there is no transport of bound states from $B$ to $A$ after the quench.

Finally, in Figure 9 (e) we show the contributions to the Yang-Yang entropy density $s_{YY}^{\scriptscriptstyle(n)}$ (see its definition in the main text) of the quasiparticles. Clearly, only quasiparticles with $n=1$ and $\lambda>0$ contribute to the entropy. Together with the results in panel (c), this implies that the entanglement between the two subsystems is generated by the transport of particles from $B$ to $A$.

Appendix F Steady-state entropy: DMRG data

In this section we discuss in more detail the numerical data for the steady-state entanglement entropy presented in Figure $3$ in the manuscript, i.e., for the quench from the state $|N\rangle\otimes|F\rangle$. We consider the entanglement entropy of a subsystem $A^{\prime}$ of length $\ell$ placed next to the boundary between $A$ and $B$. To avoid boundary effects $A^{\prime}$ is embedded in $A$, which is chosen larger than $A^{\prime}$.

Our tDMRG results are presented in Figure 8. The data are obtained using standard tDMRG simulations with bond dimension $\chi\sim 300$ for the XXZ chain with $L=40$ sites and $t\sim 16$. To calculate the von Neumann entropy we employed the techniques described in Ref. ruggiero-2016, . The different panels correspond to different values of $1.5\leq\Delta\leq 10$. The figure shows $S/\ell$ plotted versus the time after the quench. In each panel the different symbols are for subsystems of different length $1\leq\ell\leq 9$. For each $\ell$ a linear increase, followed by a saturation, is observed. The saturation value decreases upon increasing $\Delta$, in agreement with the theoretical predictions. For $t\to\infty$, in the limit $\ell\to\infty$, the entropy density $S/\ell$ should converge to the Bethe ansatz predictions. However, strong finite-size effects are visible in the Figure, due to the small $\ell$ considered. Moreover, it is interesting to observe that upon increasing $\Delta$, the data exhibit strong oscillations with time. The data presented in Figure $3$ correspond to the time average of the tDMRG data in the time window $13\leq t\leq 16$. Finally, we should mention that similar results are obtained for the quench from the state $|MG\rangle\otimes|F\rangle$.

Appendix G Scaling corrections

In this section we discuss in more detail the finite-size corrections of the tDMRG data in Figure $4$ in the manuscript. We focus on the quench from the state $|N\rangle\otimes|F\rangle$, for which the largest system subsystem sizes up to $\ell=64$ are available. We consider the scaling corrections to the entanglement production rate. Figure 9 shows the deviations $\delta S/\ell$ plotted versus the rescaled time $t/\ell$. The theory predictions are obtained using formula (7). The different curves in Figure 9 correspond to different sizes of subsystem $A$ (see Figure $1$ in the manuscript). Clearly, corrections are vanishing in the scaling limit $t,\ell\to\infty$ with fixed ratio $t/\ell$. Interestingly, the dependence of the corrections on $t/\ell$ becomes weaker upon increasing $\ell$. By comparing the data for $\ell=32$ and $\ell=64$, it is clear that corrections decay as $1/\ell$ in the thermodynamic limit. This is supported in the inset of Figure 9 showing the data for fixed $t/\ell\approx 0.1$ and $t/\ell\approx 0.27$ plotted versus $1/\ell$. We should mention that a similar behavior as $1/\ell$ has been observed for homogeneous quenches alba-2016 .