Molecular dynamics simulation of entanglement spreading in generalized hydrodynamics

TL;DR

This paper uses a molecular dynamics simulation method called flea gas within generalized hydrodynamics to numerically study entanglement spreading in an integrable quantum system, specifically the XXZ spin chain, after inhomogeneous quenches.

Contribution

It applies flea gas simulations to the XXZ model for entanglement growth, providing accurate numerical results and full-time mutual information dynamics where analytical solutions are lacking.

Findings

Flea gas method approximates entanglement growth well in XXZ model.

Numerical results match analytical predictions in certain cases.

Full-time mutual information dynamics are obtained for inhomogeneous quenches.

Abstract

We consider a molecular dynamics method, the so-called flea gas for computing the evolution of entanglement after inhomogeneous quantum quenches in an integrable quantum system. In such systems the evolution of local observables is described at large space-time scales by the Generalized Hydrodynamics approach, which is based on the presence of stable, ballistically propagating quasiparticles. Recently it was shown that the GHD approach can be joined with the quasiparticle picture of entanglement evolution, providing results for entanglement growth after inhomogeneous quenches. Here we apply the flea gas simulation of GHD to obtain numerical results for entanglement growth. We implement the flea gas dynamics for the gapped anisotropic Heisenberg XXZ spin chain, considering quenches from globally homogeneous and piecewise homogeneous initial states. While the flea gas method applied to the XXZ chain is not exact even in the scaling limit (in contrast to the Lieb--Liniger model), it yields a very good approximation of analytical results for entanglement growth in the cases considered. Furthermore, we obtain the {\it full-time} dynamics of the mutual information after quenches from inhomogeneous settings, for which no analytical results are available.

Figures

Equations

$$S=-\text{Tr}{\rho }_A\ln {\rho }_A,$$

$$S(t)=\sum_{\alpha}\Big{[}2t\!\!\!\!\!\!\int\limits_{\!2|v_{\alpha,\lambda}|t<\ell}\!\!\!\!\!\!d\lambda|v_{\alpha,\lambda}|s_{\alpha,\lambda}+\ell\!\!\!\!\!\!\int\limits_{2|v_{\alpha,\lambda}|t>\ell}\!\!\!\!\!\!d\lambda s_{\alpha,\lambda}\Big{]}.$$

$$H=\sum _{i=1}^L\frac{1}{2}(S_i^{+}{S}_{i+1}^{-}+S_i^{-}{S}_{i+1}^{+})+\Delta \sum _{i=1}^L{S}_i^z{S}_{i+1}^z.$$

$$\begin{split}\quad|\mathrm{N},\theta\rangle=\left(\frac{1+\mathcal{T}}{\sqrt{2}}\right)\ &\,\bigg{\{}[\cos(\theta/2)\left|\uparrow\right\rangle+i\sin(\theta/2)\left|\downarrow\right\rangle]\,\otimes\\ &\,\otimes[\sin(\theta/2)\left|\uparrow\right\rangle-i\cos(\theta/2)\left|\downarrow\right\rangle]\bigg{\}}^{\otimes L/2}.\end{split}$$

$$|\mathrm{D}⟩=\left(\frac{1+\mathcal{T}}{\sqrt{2}}\right){\left(\frac{|\uparrow \downarrow ⟩-|\downarrow \uparrow ⟩}{\sqrt{2}}\right)}^{\otimes L/2}.$$

$${\vartheta }_{\alpha ,\lambda }\equiv \frac{{\rho }_{\alpha ,\lambda }}{{\rho }_{\alpha ,\lambda }^{(\mathrm{t})}},{\eta }_{\alpha ,\lambda }\equiv \frac{{\rho }_{\alpha ,\lambda }^{(\mathrm{h})}}{{\rho }_{\alpha ,\lambda }}.$$

$${\rho }_{\alpha ,\lambda }+{\rho }_{\alpha ,\lambda }^{(\mathrm{h})}=a_{\alpha ,\lambda }-\sum _{\beta =1}^{\infty }{\int }_{-\pi /2}^{\pi /2}d\mu T_{\alpha \beta }(\lambda -\mu ){\rho }_{\beta ,\mu },$$

$$a_{\alpha ,\lambda }=\frac{1}{\pi }\frac{\sinh (\alpha \eta )}{\cosh (\alpha \eta )-\cos (2\lambda )},$$

$$T_{\alpha ,\beta }(\lambda )=(1-{\delta }_{\alpha \beta })a_{|\alpha -\beta |,\lambda }+2a_{|\alpha -\beta |+2,\lambda }+\cdots +2a_{\alpha +\beta -2,\lambda }+a_{\alpha +\beta ,\lambda },$$

$$\frac{E}{L}=\sum _{\alpha }{\int }_{-\pi /2}^{\pi /2}d\lambda {ϵ}_{\alpha ,\lambda }{\rho }_{\alpha ,\lambda },\mathrm{with}{ϵ}_{\alpha ,\lambda }=-\frac{\sinh \eta \sinh (\alpha \eta )}{\cosh (\alpha \eta )-\cos (2\lambda )}.$$

$$\frac{⟨\hat{Q}⟩}{L}=\sum _{\alpha }{\int }_{-\pi /2}^{\pi /2}d\lambda q_{\alpha ,\lambda }{\rho }_{\alpha ,\lambda },$$

$$v_{\alpha ,\lambda }^{\mathrm{bare}}\equiv \frac{{ϵ}_{\alpha ,\lambda }^{\prime }}{p_{\alpha ,\lambda }^{\prime }}\mathrm{with}{ϵ}_{\alpha ,\lambda }^{\prime }\equiv \frac{d{ϵ}_{\alpha ,\lambda }}{d\lambda }\text{ and }{p}_{\alpha ,\lambda }^{\prime }=\frac{d{p}_{\alpha ,\lambda }}{d\lambda },$$

$$\frac{1}{L}\ln (\text{\# of eigenstates})=s_{\mathrm{Y}Y}=\sum _{\alpha =1}^{\infty }{\int }_{-\pi /2}^{\pi /2}d\lambda s_{\alpha ,\lambda },$$

$$s_{\alpha ,\lambda }={\rho }_{\alpha ,\lambda }[\log (1+{\eta }_{\alpha ,\lambda })+{\eta }_{\alpha ,\lambda }\log (1+{\eta }_{\alpha ,\lambda }^{-1})].$$

$${\eta }_{\alpha ,\lambda }=\frac{{\eta }_{\alpha -1,\lambda +i\eta /2}{\eta }_{\alpha -1,\lambda -i\eta /2}}{1+{\eta }_{\alpha -2,\lambda }}-1\mathrm{with}\alpha \ge 2,$$

$$1+{\eta }_{1,\lambda }=\frac{T_1(\lambda +i\frac{\eta }{2})}{\varphi (\lambda +i\frac{\eta }{2})}\frac{T_1(\lambda -i\frac{\eta }{2})}{\bar{\varphi }(\lambda -i\frac{\eta }{2})},$$

$$T_1(\lambda )$$

$$+[\cos (2\theta )+3][2\cos (2\lambda )-1]+2{\sin }^2(\theta )\cos (4\lambda )\},$$

$$\varphi (\lambda )$$

$$\bar{\varphi }(\lambda )$$

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

$${\partial }_t{\rho }_{\alpha ,\lambda }(\zeta )+{\partial }_x(v_{\alpha ,\lambda }(\zeta ){\rho }_{\alpha ,\lambda }(\zeta ))=0.$$

$$v_{\alpha ,\lambda }(\zeta )=v_{\alpha ,\lambda }^{\text{bare}}(\zeta )+\sum _{\beta }{\int }_{-\pi /2}^{\pi /2}d\mu \frac{T_{\alpha \beta }(\lambda -\mu )}{a_{\alpha ,\lambda }}{\rho }_{\beta ,\mu }(\zeta )(v_{\alpha ,\lambda }(\zeta )-v_{\beta ,\mu }(\zeta )),$$

$${\partial }_t\rho (v_{\mathrm{b}};x,t)+{\partial }_x(v(v_{\mathrm{b}};x,t)\rho (v_{\mathrm{b}};x,t))=0.$$

$$v(v_{\mathrm{b}};x,t)=v_{\mathrm{b}}+d\int dw\rho (w;x,t)(v(v_{\mathrm{b}};x,t)-v(w;x,t)).$$

$$d_{\alpha ,\beta }(\lambda ,\mu )=\frac{T_{\alpha ,\beta }(\lambda -\mu )}{a_{\alpha ,\lambda }}.$$

$$N_p=L\sum _{\alpha }\int d\lambda {\rho }_{\alpha ,\lambda }$$

$$r_{\alpha }=\frac{L}{N_{\mathrm{p}}}\int d\lambda {\rho }_{\alpha ,\lambda }.$$

$$v_{\alpha ,\lambda }^{\text{bare}}=\frac{\sinh (\eta )a_{\alpha ,\lambda }^{\prime }}{2a_{\alpha ,\lambda }}.$$

$$S(t)=\Big{\langle}\sum_{\textrm{\tiny shared pairs $(\lambda,-\lambda)$}}\frac{s_{\alpha,\lambda}}{\rho_{\alpha,\lambda}}\Big{\rangle}_{t}.$$

$${\rho }_0=\frac{1}{Z}{e}^{-\beta H+(\beta h)S^z},$$

Full text

Molecular dynamics simulation of entanglement spreading in generalized hydrodynamics

Márton Mestyán1*, Vincenzo Alba2

1 International School for Advanced Studies (SISSA) and INFN,

Via Bonomea 265, 34136, Trieste, Italy

2 Delta Institute for Theoretical Physics, University of Amsterdam,

Science Park 904, 1098 XH Amsterdam, the Netherlands

• mestyan@fmf.uni-lj.si

Abstract

**We consider a molecular dynamics method, the so-called flea gas for computing the evolution of entanglement after inhomogeneous quantum quenches in an integrable quantum system. In such systems the evolution of local observables is described at large space-time scales by the Generalized Hydrodynamics approach, which is based on the presence of stable, ballistically propagating quasiparticles. Recently it was shown that the GHD approach can be joined with the quasiparticle picture of entanglement evolution, providing results for entanglement growth after inhomogeneous quenches. Here we apply the flea gas simulation of GHD to obtain numerical results for entanglement growth. We implement the flea gas dynamics for the gapped anisotropic Heisenberg XXZ spin chain, considering quenches from globally homogeneous and piecewise homogeneous initial states. While the flea gas method applied to the XXZ chain is not exact even in the scaling limit (in contrast to the Lieb–Liniger model), it yields a very good approximation of analytical results for entanglement growth in the cases considered. Furthermore, we obtain the full-time dynamics of the mutual information after quenches from inhomogeneous settings, for which no analytical results are available. **

1 Introduction

In the last decade, the study of isolated quantum many-body systems out of equilibrium provided new insights on the deep interplay between entanglement and thermodynamics, shedding new light on the fundamental question how statistical ensembles emerge from the out-of-equilibrium dynamics after a quantum quench [1, 2, 3, 4, 5, 6].

Integrable models offer an ideal setting for understanding generic features of the entanglement dynamics after a quantum quench [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35]. Indeed, for integrable systems the dynamics of entanglement-related quantities can be understood within the so-called quasiparticle picture [8]. Here we focus on the out-of-equilibrium dynamics of the entanglement entropy (von Neumann entropy), which is defined as [36, 37, 38, 39]

[TABLE]

with $\rho_{A}$ the reduced density matrix of a macroscopic subsystem $A$ (see Fig. 1 for a one-dimensional setup). In the quasiparticle picture, pairs of entangled quasiparticles are produced after the quench. As these pairs propagate ballistically, they entangle larger regions of the system (see Fig. 1 (a)). At a given time after the quench, the von Neumann entanglement entropy is the sum of the individual contributions coming from each pair that is shared between $A$ and its complement. This picture has been explicitly verified in free-fermion models [7]. It has been shown recently that it holds true also in the presence of interactions [32].

The quasiparticle prediction for the entanglement entropy of subsystem $A$ of length $\ell$ after a quench in generic integrable systems reads as [32]

[TABLE]

Here the index $\alpha$ labels the different types of quasiparticles present in the model, and $\lambda$ is the so-called rapidity, which distinguishes different modes of the same type of quasiparticles. The quasiparticle picture is built on two important ideas. First, at long times the entanglement entropy between $A$ and its complement coincides with the thermodynamic entropy of the statistical ensemble that describes local steady-state properties in $A$. For generic translationally invariant steady states, this ensemble is a Generalized Gibbs Ensemble [2, 3, 4, 5, 6, 40, 41, 42, 43] (GGE). Second, the velocity of the entangling quasiparticles are the group velocities of the particle-hole excitations [44] over the GGE steady state. A crucial observation is that, in contrast with free-fermion models, $v_{\alpha,\lambda}$ depends on the GGE describing the steady state after the quench (or equivalently on the pre-quench initial state), and it is “dressed” by the interactions. For generic lattice models with local interaction the velocities are finite, i.e, for any $\alpha,\lambda$, $|v_{\alpha,\lambda}|<\infty$. The function $s_{\alpha,\lambda}$, i.e., the entropy carried by quasiparticle pairs of type $\alpha$ in mode $\lambda$, is the contribution of these quasiparticles to the thermodynamic entropy of the GGE describing subsystem $A$.

The approach of Ref. [32] has been generalized in [45, 46, 47] to calculate the Rényi entropies in the steady-state, whereas calculating their full-time dynamics is a challenging open problem. Remarkably, the quasiparticle picture allows to obtain the dynamics of the logarithmic negativity [15, 48], which is a proper entanglement measure for mixed states.

Recently, there has been a growing interest in understanding the entanglement dynamics after quenches from piecewise homogeneous initial states. In the standard setup (see Fig. 1 (b)) two homogeneous chains in a different state ($L,R$ in Fig. 1) are joined together at $t=0$. One then studies the ensuing dynamics under an integrable globally homogeneous Hamiltonian. During the time evolution a lightcone spreads from the interface between the chains. For typical initial states, in the limit of long times and large distances $x$ from the origin, the system reaches at each fixed ray $\zeta\equiv x/t$ (see Fig. 1), a Local Quasi Stationary State [49] (LQSS) which is described by a GGE. These $\zeta$-dependent GGEs can be described analytically within the Generalized Hydrodynamics (GHD) formalism [50, 51]. Furthermore, by combining the quasiparticle picture with the GHD approach [52, 46, 53, 54] it is possible, in principle, to generalize the quasiparticle picture [32] to inhomogeneous settings. However, actually calculating the full-time entanglement dynamics in this case is a demanding task. The main difficulty is that, unlike in homogeneous quenches, the trajectories of the quasiparticles are not straight lines. Explicit results are available only in the short time limit [54] $t/\ell\to 0$ and the long time limit $t/\ell\to\infty$.

In order to overcome the above difficulties, we use a mapping between the GHD and the “molecular dynamics” of a system of classical particles called flea gas [55]. This mapping allows one to obtain the out-of-equilibrium dynamics in quantum integrable systems by performing classical simulations [56, 57, 58]. So far, the flea gas method has been employed only to integrable field theories, such as the Lieb-Liniger gas, but not to lattice models. In this paper, we discuss a generalization of the flea gas approach to the spin-$1/2$ anisotropic XXZ chain. An important remark is that, unlike for the Lieb-Liniger model [55], for the XXZ chain it is not straightforward to show analytically that the flea gas dynamics is fully equivalent to the GHD. In general for the XXZ chain the flea gas dynamics is expected to be different from the GHD. Our results suggest that that deviations are present, although much larger systems would be needed to ensure that the system is in the scaling limit. Still, we observe that deviations of the flea gas from the GHD are small. This results in a surprisingly good agreement between the flea gas and the GHD for the dynamics of local observables and the entanglement entropy.

Our main goal is to show that with moderate computational effort the flea gas framework allows one to compute accurate numerical predictions for the full-time entanglement dynamics after quenches from arbitrary inhomogeneous initial conditions, provided that the Generalized Hydrodynamics approach can be applied.

The flea gas algorithm requires to know the bare quasiparticle velocity and the scattering matrix between quasiparticles, which are easily obtained for a generic integrable model. The third and only non-trivial input to this method is the initial condition of the flea gas dynamics. For quenches from homogeneous initial states, the initial condition of the flea gas dynamics is given by the GGE describing the long time limit after the quench. The idea is that at comparatively short times after the quench the quasiparticles are described by the GGE density. For piecewise homogeneous setups (see, e.g., Fig. 1 (b)), the initial condition is given by the GGEs describing each homogeneous chain at $t=0$ in the hydrodynamic timescale.

We stress that it is also necessary to know the structure of quantum correlations in the initial states, i.e., how entanglement is shared among the quasiparticles. Here we restrict ourselves to the situation in which only entangled pairs are present. Other situations, for instance the case of entangled “triplets”, have been considered, at least in free models [26, 27].

As a benchmark of the method, we show that for quenches from homogeneous states our numerical results are in perfect agreement with (2). Moreover, for quenches from inhomogeneous initial states, in the limit $t,\ell\to\infty$ with $t/\ell\ll 1$ our results confirm the analytical predictions in Ref. [54]. Finally, to show the versatility of the method we provide results for the full-time dynamics of the entanglement entropy and of the mutual information after a quench from inhomogeneous settings, which are not easily accessible analytically [54].

The manuscript is organized as follows. In Section 2 we introduce the XXZ chain and the quenches considered in this study. In Section 3 we discuss the Bethe ansatz treatment of generic thermodynamic ensembles. The thermodynamic Bethe ansatz framework (TBA) is introduced in Section 3.1, which is followed by the description of steady states after homogeneous quenches in Section 3.2 and a summary of the GHD approach for inhomogeneous quenches in Section 3.3. In Section 4 we introduce the flea gas method, its implementation for the XXZ chain (see Section 4.1), and the calculation of entanglement-related quantities (see Section 4.2). Our numerical results are discussed in Section 5. In Section 5.1 we benchmark the method for homogeneous quenches. In Section 5.2 we provide results for entanglement entropy after quenches from piecewise homogeneous initial states. Finally, in Section 5.3 we discuss the mutual information. Section 6 concludes the article by mentioning some interesting future directions.

2 Model and quenches

The flea gas method that we propose for calculating the entanglement dynamics is expected to work for generic interacting integrable models, both on the lattice and in the continuum. However, here we provide results only for a prototypical lattice model, the spin-$1/2$ XXZ chain, which describes a system of interacting spins on a ring, and it is defined by the Hamiltonian

[TABLE]

Here $S_{i}^{+,-,z}$ are spin-$1/2$ operators, and $\Delta$ is the so-called anisotropy parameter. We restrict ourselves to the region with $\Delta>1$, where the system is gapped in the thermodynamic limit, although the method can be applied for $\Delta\leq 1$ as well. We impose periodic boundary conditions by identifying sites $1$ and $L+1$. We construct our initial states by joining two homogeneous blocks that are prepared in either the translationally invariant tilted Néel state or in the translationally invariant Majumdar-Ghosh (dimer) state. The method is applicable, in principle, to any low-entangled initial state.

The translationally invariant tilted Néel state is denoted as $|\mathrm{N},\theta\rangle$, and it is obtained by rotating the Néel state $\left|\uparrow\downarrow\uparrow\dots\right\rangle$ around the $\hat{z}$ axis and making it translationally invariant, i.e.,

[TABLE]

Here $\theta$ is the tilting angle and $\mathcal{T}$ is the one site translation operator to the right. The Néel state is recovered for $\theta=0$. Similarly, the translationally invariant dimer state $|\mathrm{D}\rangle$ is defined as

[TABLE]

In the homogeneous setup (Fig. 1 (a)), the chain is prepared in one of the states  (4) or (5) at $t=0$, and the system is let to evolve under (3). In the inhomogeneous case (Fig. 1 (b)) we consider quenches from the initial state $|\Psi_{0}\rangle=|\mathrm{N},\theta\rangle\otimes|\mathrm{D}\rangle$.

3 Bethe ansatz description of thermodynamic macrostates in the

XXZ chain

Here we introduce the thermodynamic Bethe ansatz (TBA) treatment of the XXZ chain [59], focusing on the features that are needed in the implementation of the flea gas method. First, we summarize the general TBA framework in Section 3.1. Then we report the TBA description of the steady states after the considered homogeneous quenches in Section 3.2. Finally, we summarize the generalized hydrodynamics (GHD) framework for quenches from inhomogeneous states in Section 3.3.

3.1 Thermodynamic Bethe Ansatz (TBA)

The XXZ chain is solved by the Bethe ansatz [59], which allows one to construct the eigenstates of (3). In the Bethe ansatz, the eigenstates are constructed with respect to the reference state with all spins up $|\uparrow\uparrow\cdots\uparrow\rangle$. Since the total magnetization $\sum_{j}S_{j}^{z}$ commutes with (3), the eigenstates are characterized by the total number $N$ of down spins, which is a good quantum number of the state. We refer to $N$ as the number of particles. In this study we focus on the thermodynamic limit, i.e., the limit $L,N\to\infty$, with particle density $N/L$ fixed.

A distinctive feature of integrable models is that their eigenstates can be interpreted as a collection of well-defined, i.e., having infinite lifetime, quasiparticles. For generic integrable models the quasiparticles are labelled by a set of real parameters $\{\lambda_{\alpha,j}\}_{\alpha,j}$, which are called rapidities. For the XXZ chain at $\Delta>1$ one has $\lambda_{\alpha,j}\in[-\pi,\pi]$. In general, there can be different species of quasiparticles. These are distinguished by the integer index $\alpha$. The total number of species depends on $\Delta$. For instance, at $\Delta\geq 1$ there is an infinite number of them. Quasiparticles with $\alpha=1$ correspond to magnon-like excitations, whereas for $\alpha>1$ they are bound states of $\alpha$ magnons ($\alpha$-strings [59]).

In the thermodynamic limit it is impossible to consider the individual rapidities $\lambda_{\alpha,j}$ of the quasiparticles. Instead, the standard TBA framework [59] uses the density of quasiparticles in rapidity space $\rho_{\alpha,\lambda}$, which are real functions of $\lambda$ for each species $\alpha$. One can also define the hole density $\rho_{\alpha,\lambda}^{\scriptscriptstyle(\mathrm{h})}$ as the density of unoccupied states in rapidity space. Another important quantity is the total density of states $\rho_{\alpha,\lambda}^{\scriptscriptstyle(\mathrm{t})}=\rho_{\alpha,\lambda}+\rho_{\alpha,\lambda}^{\scriptscriptstyle(\mathrm{h})}$. For a generic integrable model, this is a non-trivial function of $\lambda$, whereas for non-interacting systems the total density is a constant, reflecting that the rapidities are equally spaced. For the following, it is also convenient to define the filling functions $\vartheta_{\alpha,\lambda}$ and the functions $\eta_{\alpha,\lambda}$ as

[TABLE]

The densities $\rho_{\alpha,\lambda}$ and $\rho_{\alpha,\lambda}^{\mathrm{h}}$ are constrained by the Bethe equations arising from the periodic boundary conditions. In the thermodynamic limit, the Bethe equations become the Bethe–Gaudin–Takahashi (BGT) equations [59]

[TABLE]

where the functions $a_{\alpha,\lambda}$ are

[TABLE]

and $\eta\equiv\mathrm{arccosh}(\Delta)$. In (7), the scattering matrix $T_{\alpha,\beta}$ is defined as

[TABLE]

where $a_{\alpha,\lambda}$ is the same as in (8). The matrix $T_{\alpha,\beta}(\lambda-\mu)$ encodes all the information about the scattering between quasiparticles of type $(\alpha,\lambda)$ and $(\beta,\mu)$, and it will be crucial in the implementation of the flea gas algorithm (see section 4).

In the TBA framework, any set of particle and hole densities $\rho_{\alpha,\lambda},\rho_{\alpha,\lambda}^{\scriptscriptstyle(\mathrm{h})}$ identifies a thermodynamic macrostate. All the information about thermodynamic expectation values of local and quasi local operators is encoded in the functions $\rho_{\alpha,\lambda}$. These expectation values are determined by summing over the quasiparticles species and integrating over their rapidity. For instance, for the XXZ chain the energy of a macrostate identified by a set of densities $\rho_{\alpha,\lambda}$ reads [59]

[TABLE]

Besides the energy, integrable models have an infinite set of quasilocal charges that commute with the Hamiltonian. In the case of the XXZ model, these charges are obtained as the derivatives of the transfer matrix. In the thermodynamic limit, a conserved quasilocal charge ${\hat{Q}}$ is expressed as

[TABLE]

where ${q}_{\alpha,\lambda}$ is a known function, the density of the charge.

An important quantity that we will use is the bare velocity of the quasiparticles $v_{\alpha,\lambda}^{\mathrm{bare}}$. This is the group velocity defined from the bare quasiparticles dispersion as

[TABLE]

where $\epsilon_{\alpha,\lambda}$ is the bare energy of a quasiparticle defined in (10), and $p_{\alpha,\lambda}$ its bare momentum with $p^{\prime}_{\alpha,\lambda}=2\pi a_{\alpha,\lambda}$.

Interestingly, for generic integrable models the quasiparticle velocities depend on the thermodynamic macrostate [44], and they are “dressed” by the interactions. This happens because in interacting integrable models the addition or removal of a single quasiparticle causes a global rearrangement of the rapidities of the other quasiparticles. The net effect is a “dressing” of the bare quasiparticles properties, including the energies $\epsilon_{\alpha,\lambda}$, and hence the group velocities.

The correspondence between thermodynamic macrostates and microscopic eigenstates of (3) is not one-to-one. In fact, the densities $\rho_{\alpha,\lambda}$ and $\rho_{\alpha,\lambda}^{\scriptscriptstyle(\mathrm{h})}$ do not uniquely determine a microscopic eigenstate. In the thermodynamic limit, the number of microscopic eigenstates that give rise to the same set of macroscopic densities diverges exponentially with the system size. The number of these thermodynamically equivalent eigenstates is given in terms of the so-called Yang–Yang entropy [59] as

[TABLE]

with the entropy density function $s_{\alpha,\lambda}$ being

[TABLE]

The TBA formalism has been applied to describe thermal properties of the XXZ chain [59]. The corresponding thermodynamic ensemble is the Gibbs ensemble, and the Yang-Yang entropy (13) is the usual thermal entropy. However, the TBA framework can also be used to describe the thermodynamic macrostate arising after a quantum quench, i.e., the macrostate described by a generalized Gibbs ensemble (GGE) that takes into account the conservation of all the quasilocal charges. Then, the corresponding Yang-Yang entropy becomes the GGE thermodynamic entropy. Remarkably, this Yang-Yang entropy coincides with the von Neumann entanglement entropy of the post-quench steady state [60, 61, 62, 33, 63, 32], and it is one of the main ingredients to reconstruct the full-time dynamics of the entanglement entropy (as it is clear from (2)).

3.2 TBA treatment of the steady state after quenches from homogeneous states

Integrable models possess an extensive number of local and quasilocal conserved quantities. Their expectation value in the initial state is preserved during the dynamics. Thus, the post-quench dynamics is strongly constrained, implying that a Generalized Gibbs Ensemble, instead of the standard Gibbs one, has to be used to describe local properties of the steady state in the long time limit. The GGE can be thought of as emerging from a generalized microcanonical principle [64]. The eigenstates entering in the microcanonical average are the ones that have the correct expectation value of the local and quasilocal conserved quantities. In the thermodynamic limit the vast majority of these eigenstates give rise to the same set of densities $\rho_{\alpha,\lambda}$. This set of densities is called the representative state. The corresponding hole densities $\rho^{\scriptscriptstyle(h)}_{\alpha,\lambda}$ are obtained from the BGT equations (7). The representative state encodes all information about local properties of the steady state.

Within the TBA approach there are several techniques to determine the densities $\rho_{\alpha,\lambda}$ that describe the representative state. For instance, they can be determined from the overlaps between the eigenstates of the XXZ chain with the initial state, by using the so-called Quench Action method [5]. Alternatively, they can be calculated from the knowledge of the initial values of the local and quasilocal conserved quantities [65]. The latter method allows to deal, in principle, with any translationally invariant initial state.

It is customary to describe the representative state in terms of $\eta_{\alpha,\lambda}$. These functions satisfy the so-called Y-system [66], which leads to the set of recursive equations for the functions $\eta_{\alpha,\lambda}$ (cf. (6)) as

[TABLE]

with the convention that $\eta_{0,\lambda}\equiv 0$. Once $\eta_{1,\lambda}$ is known, then the functions $\eta_{\alpha,\lambda}$ for $\alpha>1$ can be computed using (15). The corresponding particle densities $\rho_{\alpha,\lambda}$ can be computed by substituting the $\eta_{\alpha,\lambda}$ in the BGT equations (7).

Clearly, Eq. (15) implies that to determine the steady-state properties after a homogeneous global quench, one needs to calculate only $\eta_{1,\lambda}$. For both the tilted Néel state and the Majumdar-Ghosh state, which are relevant for this work, the function $\eta_{1,\lambda}$ is exactly known. For the tilted Néel state $|\mathrm{N},\theta\rangle$ with tilting angle $\theta$, one has [67]

[TABLE]

where

[TABLE]

and

[TABLE]

For the dimer state, the funcion $\eta_{1,\lambda}$ reads  [65]

[TABLE]

3.3 Quenches from piecewise homogeneous initial states: Generalized Hydrodynamics

Here we consider quenches from piecewise homogeneous initial states (as described in Fig. 1 (b)). Two semi-infinite chains $L$ and $R$ are prepared in the translationally invariant states $|\Psi_{L}\rangle$ and $|\Psi_{R}\rangle$, and are suddenly joined together at $t=0$. The ensuing dynamics is governed by the globally translational invariant Hamiltonian (3). Recently, it has been shown that a Generalized Hydrodynamics (GHD) approach allows to study this quench [51, 50] in the long time limit at large spatial scales. Physically, during time evolution a light-cone spreads from the interface between the two chains. Outside of this lightcone, the properties of the system are the same as after the homogeneous quenches from the states $|\Psi_{L}\rangle$ and $|\Psi_{R}\rangle$ (see Fig. 1 (a)). Inside the lightcone and at late times, the expectation values of local and quasilocal observables become functions of $\zeta=x/t$. This reflects the propagation of stable quasiparticles between the two chains. This also suggests the emergence of a local quasi-stationary state for each $\zeta$. For integrable models, this corresponds to a $\zeta$-dependent GGE. Within the TBA framework, the GGE is represented by a set of TBA densities $\{\rho_{\alpha,\lambda}(\zeta)\}_{\alpha=1}^{\infty}$.

The key result of the GHD is that because of infinitely many conserved charges, the densities $\rho_{\alpha,\lambda}(\zeta)$ satisfy a simple continuity equation [51, 50] as

[TABLE]

In the bipartite quench the dependence of local observables on the space-time coordinates $x,t$ is only through $\zeta$. Here $v_{\alpha,\lambda}(\zeta)$ are the dressed group velocities (the same as in (2)), which are solutions of the system of integral equations [44]

[TABLE]

where $T_{\alpha\beta}$ is the scattering matrix (cf. (9) for the result for the XXZ chain) and $a_{\alpha,\lambda}$ is defined in (8). The functions $v_{\alpha,\lambda}^{\scriptscriptstyle\mathrm{bare}}$ are the bare velocities defined in (12).

Physically, Eq. (23) reflects that, due to integrability, the scattering between the quasiparticles is elastic, and the only effect of the interactions is to renormalize the quasiparticles velocities. Indeed, the term $\rho_{\beta,\mu}|v_{\alpha,\lambda}-v_{\beta,\mu}|$ in (23) is the number of quasiparticles with rapidity $\mu$ and of species $\beta$ that scatter in the unit time with the quasiparticle of species $\alpha$ and rapidity $\lambda$. The ratio $T_{\alpha\beta}(\lambda-\mu)/a_{\alpha,\lambda}$ can be interpreted as an effective shift of the trajectory of the quasiparticle with label $\alpha,\lambda$ due to the scatterings. This interpretation underlies the flea gas method (cf. section 4).

The GHD approach has been successfully applied to describe transport properties in spin systems, one-dimensional integrable field theories, both classical and quantum [68, 69, 70, 71, 58, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86]. Very recently, it has been shown that GHD provides a precise framework to describe experiments with trapped cold atoms [56]. A recent interesting direction is to generalize the approach to include diffusive corrections [87, 88, 89, 90, 91]. Finally, the GHD approach can be used to study the entanglement dynamics after quenches from piecewise-homogeneous initial states [52, 46, 54].

Unfortunately, calculating the full-time entanglement dynamics is in general a demanding task. The reason is that inside the lightcone (see Fig. 1 (b)) the trajectories of the quasiparticles are not straight lines. Explicit results are easily obtained only in some regimes. For instance, the steady-state value of the von Neumann entropy of a finite interval placed next to the interface between the two chains [52], as well as the growth rate of the entanglement entropy between two-semi-infinite systems [54], can be calculated in terms of the $\zeta=0$ macrostate only.

4 Flea gas approach for out-of-equilibrium integrable systems

The flea gas approach was introduced in Ref. [92] as an effective numerical method to simulate the GHD by employing classical “molecular dynamics” techniques. The method allows to simulate the dynamics of a quantum system starting from any thermodynamic macrostate, both homogeneous as well as inhomogeneous. So far, the approach has been implemented for the Lieb-Liniger gas but not for lattice systems such as the XXZ model.

The method was inspired by the correspondence between the continuity equation (22) and the hydrodynamic equations of a system of classical particles (hard-rod gas). Hard rods are classical one-dimensional objects undergoing elastic scattering. Here we denote their length as $d$. The hard rods dynamics is as follows. Hard rods move like free particles with bare velocity $v_{\mathrm{b}}$. When the distance between the centers of two hard rods equals $d$, they exchange their velocities. Following Ref. [92], here we adopt an alternative interpretation. One can think of hard rods as point-like objects. The scattering is then implemented as follows. When two hard rods are at the same point in space they scatter. The scattering consists of an instantaneous displacement by length $d$ of the positions of the two particles. Precisely, after assuming $d>0$ we impose that the particle on the left (right) is shifted by $d$ to the right (left).

Let us define the density of rods with “bare” velocity $v_{\mathrm{b}}$ as $\rho(v_{\mathrm{b}})$. The number of rods with velocity between $v_{\mathrm{b}}$ and $v_{\mathrm{b}}+dv$ and in the spatial interval $dx$ is $\rho(v_{b})dvdx$. During time evolution, many scatterings occur. The net effect is a renormalization of the velocity of the hard rods. Let us define this space-time dependent renormalized or “dressed” velocity as $v(v_{\mathrm{b}};x,t)$ The dressed velocity is a function of the bare velocity $v_{\mathrm{b}}$. The density $\rho(v_{\mathrm{b}})$ obeys the continuity equation [93]

[TABLE]

The renormalized velocity is given by the integral equation [93]

[TABLE]

Equation (25) has the same structure as (23), and it admits a simple interpretation. The expression $\rho(w)|v(v_{\mathrm{b}})-v(w)|$ is the number of scatterings per unit time between the hard rod with velocity $v$ and the ones with velocity $w$. The second term on the right hand side in (25) is the total shift that happens in the unit time to the trajectory of the hard rod with bare velocity $v$ due to the scatterings with other hard rods.

4.1 Flea gas for the XXZ chain and numerical implementation

We now discuss the application of the flea gas method for the XXZ chain. Before describing the method for the quenches from an inhomogeneous initial state, it is useful to consider the case of homogeneous ones. The TBA densities $\rho_{\alpha,\lambda}$ describing the steady state after the quench are stationary and homogeneous, i.e., they do not depend on $x,t$. The group velocity $v_{\alpha,\lambda}$ of the quasiparticles are obtained by solving the TBA system (23). The crucial observation is that Eq. (23) has the same structure as the equation for the hard rod gas (25). Equation (23) can be interpreted as the dressing equation for the velocities of a system of multi-species and point-like classical particles undergoing elastic scattering. Now each particle is identified by a double index $(\alpha,\lambda)$, and $T_{\alpha,\beta}(\lambda-\mu)/a_{\alpha,\lambda}$ is interpreted as a scattering length. Thus, we define $d_{\alpha,\beta}(\lambda,\mu)$ as

[TABLE]

Similarly to the hard rods, the particles move freely with bare velocities $v^{\scriptscriptstyle\mathrm{bare}}$ (now given by (12)). Scatterings occur when two particles are at the same point in space. If the particle coming from the left has labels $(\alpha,\lambda)$, and the particle coming from the right has labels $(\beta,\mu)$, then the particle coming from the left will jump $d_{\alpha,\beta}(\lambda,\mu)$, and the particle coming from the right will jump $-d_{\beta,\alpha}(\mu,\lambda)$.

A crucial remark is in order. Unlike the case of the Lieb-Liniger model [55] it is not straightforward to show that the dynamics outlined above reproduces the correct dressing for the bare velocities of the particles, i.e., Eq. (23). First, the displacement of the trajectory of a given particle is given as $\Delta x=v_{\alpha,\lambda}\Delta t$, which defines the dressed velocity $v_{\alpha,\lambda}$. The dressing of the velocity arises from the scattering with the other particles. The number of scatterings per unit time between a particle with label $(\alpha,\lambda)$ and particles with label $(\beta,\mu)$ is given as $\rho_{\beta,\mu}|v_{\alpha,\lambda}-v_{\beta,\mu}|\Delta t$.

The key issue is how to determine the direction of the jump. During the flea gas dynamics the particles move with their bare velocities. If a particle moving at $v_{\alpha,\lambda}^{\scriptscriptstyle\mathrm{bare}}$ scatters with another one with velocity $v_{\beta,\mu}^{\scriptscriptstyle\mathrm{bare}}$ its trajectory gets shifted by $\mathrm{sign}(v_{\alpha,\lambda}^{\scriptscriptstyle\mathrm{bare}}-v_{\beta,\mu}^{\scriptscriptstyle\mathrm{bare}})d_{\alpha,\beta}(\lambda,\mu)$. For the Lieb-Liniger gas one can show that the dressed velocities are monotonic functions of the bare ones, which implies that $\mathrm{sign}(v_{\alpha,\lambda}^{\scriptscriptstyle\mathrm{bare}}-v_{\beta,\mu}^{\scriptscriptstyle\mathrm{bare}})=\mathrm{sign}(v_{\alpha,\lambda}-v_{\beta,\mu})$. This ensures that the jumped length is $\mathrm{sign}(v_{\alpha,\lambda}-v_{\beta,\mu})d_{\alpha,\beta}(\lambda,\mu)$. By summing over $\beta$ and integrating over $\mu$, one obtains the term on the right-hand-side in (23). This shows that the flea gas dynamics gives the correct dressing for the group velocities of the particles. On the other hand, for the XXZ chain the dressed velocities are not monotonically increasing functions of the bare ones. An important consequence is that now $\mathrm{sign}(v_{\alpha,\lambda}^{\scriptscriptstyle\mathrm{bare}}-v_{\beta,\mu}^{\scriptscriptstyle\mathrm{bare}})\neq\mathrm{sign}(v_{\alpha,\lambda}-v_{\beta,\mu})$. This implies that for the XXZ chain one cannot conclude that the total jumped length is given as $(v_{\alpha,\lambda}-v_{\beta,\mu})d_{\alpha,\beta}(\lambda,\mu)$. To overcome this problem, our strategy here is to use the flea gas dynamics as outlined above, showing numerically that, at least in the quenches that we consider, it gives the correct dressing for the group velocities of the particles (see section 5.1).

We now discuss the details of the implementation of the flea gas method for the XXZ chain. The system is in the continuum, and it is of length $L$. Both space and time are treated as continuous variables. For a homogeneous quench, the initial state of the simulation is prepared as follows. First, we create a total number of particles $N_{p}$. The particles are described by the TBA densities $\rho_{\alpha,\lambda}$, which contain the full information about the post-quench GGE (see section 3.2 for the results for the quenches considered here). $N_{p}$ is chosen such that one has the correct value of the particle density, i.e.,

[TABLE]

Note that $N_{p}$ is not the total number of down spins $N$, which is given as $N=L\sum_{\alpha}\int d\lambda\alpha\rho_{\alpha,\lambda}$. This simply reflects that in the simulation multi-spin bound states are treated as individual point-like particles. The particles are labeled as $1,\dots,{N_{p}}$. Here we restrict ourselves to the situation in which only pairs of entangled quasiparticles with opposite rapidities [32] are present. For convenience, particles forming an entangled pair are labelled by consecutive integers as $(2\gamma,2\gamma+1)$ with $\gamma=1,\dots,N_{p}/2$. To each pair we assign a species label $\alpha$ with probability $r_{\alpha}$ given as

[TABLE]

Similarly, rapidities $\lambda_{2\gamma}=-\lambda_{2\gamma+1}$ are assigned to the pairs with probability $\rho_{\alpha,\lambda}=\rho_{\alpha,-\lambda}$. The position of each pair is random in the interval $[-L/2,L/2]$. Note that entangled particles are produced at the same point in space, implying $x_{2\gamma}=x_{2\gamma+1}$. However, to avoid spurious scatterings when the dynamics starts, we impose a tiny displacement between entangled particles. Finally, we assign to each pair their contribution to the Yang-Yang entropy, which is $s_{\alpha,\lambda}/\rho_{\alpha,\lambda}$ (cf. (14)).

During the time evolution, the particles move with the bare velocities $v^{\scriptscriptstyle\textrm{\tiny bare}}_{\alpha,\lambda}$, given as (cf (12))

[TABLE]

Here $a_{\alpha,\lambda}$ are defined in (8) and $a^{\prime}_{\alpha,\lambda}\equiv da_{\alpha,\lambda}/d\lambda$. During the simulation only the position of the particles are updated, whereas their labels, velocities, and entropies remain the same. Particles can collide, jumping backward and forward of distance $d_{\alpha,\beta}(\lambda-\mu)$ (cf. (26)). This happens as follows. Let us denote two colliding particles as $P_{1}$ and $P_{2}$, $P_{1}$ being the left particle and $P_{2}$ the right one, respectively. Let us assume that $P_{1}$ and $P_{2}$ have labels $(\alpha,\lambda)$ and $(\alpha^{\prime},\lambda^{\prime})$, respectively. Thus, $P_{1}$ jumps to the right of distance $d_{\alpha,\alpha^{\prime}}(\lambda,\lambda^{\prime})$, whereas particle $P_{2}$ jumps to the left of distance $d_{\alpha^{\prime},\alpha}(\lambda^{\prime},\lambda)=-d_{\alpha,\alpha^{\prime}}(\lambda,\lambda^{\prime})$. It is crucial to observe that while jumping, $P_{1}$ and $P_{2}$ can scatter with other particles that are within $d_{\alpha,\alpha^{\prime}}(\lambda,\lambda^{\prime})$. For example, if the trajectory of $P_{1}$, after scattering with $P_{2}$, crosses that of a third particle $P_{3}$ with label $(\alpha^{\prime\prime},\lambda^{\prime\prime})$, $P_{1}$ scatters with $P_{3}$, as well. This means that, in principle, there is a “cascade” of scatterings initiating when $P_{1}$ and $P_{2}$ collide.

The complete flea gas algorithm is illustrated in Fig. 2. The first step is to identify the pair of particles $P_{1}$ and $P_{2}$ that scatter first, and the corresponding scattering time $t_{\mathrm{\mathrm{coll}}}$. This is performed by the routine $\textsc{Find}(P_{1},P_{2},t_{\mathrm{coll}})$ in Fig. 2. This can be done efficiently by using standard methods in molecular dynamics simulations (see for instance Ref. [94]). Then, all the particles are evolved until $t_{\mathrm{coll}}$, when the scattering between $P_{1}$ and $P_{2}$ occurs. This is described by the procedure Collide in Fig. 2. $P_{1}$ and $P_{2}$ are instantaneously displaced by a distance $d_{1,2}$ and $d_{2,1}$ (cf. (26)). The displacement of the particles is implemented with the procedures JumpLeft and JumpRight, which are described in Fig. 3. Note that the two scattering particles are marked before starting the collision (see procedure Mark). This is to prevent that, while scattering with near particles, $P_{1}$ and $P_{2}$ scatter again with each other. Marked particles, instead of scattering, cross each other. After the scattering cascade starting with their first collision happened, $P_{1}$ and $P_{2}$ are unmarked.

4.2 Entanglement dynamics in flea gas simulations

The entanglement entropy at a given time is computed by counting the entangled pairs (weighted with their Yang-Yang entropy) that are shared between the subsystem of interest $A$ (cf. Fig. 1) and the rest, i.e., the number of pairs $(P_{\gamma},P_{\gamma+1})$, such that $x_{\gamma}$ and $x_{\gamma+1}$ are in different subsystems. The result for the entanglement entropy reads as

[TABLE]

Here the average $\langle\rangle_{t}$ is over different realizations of the flea gas dynamics up to time $t$. The sum is over the pairs that are shared between the two subsystems. Importantly, the factor $1/\rho_{\alpha,\lambda}$ takes into account that different types $(\alpha,\lambda)$ of particles appear in the sum (30) with a frequency $\ell\rho_{\alpha,\lambda}$.

5 Numerical results

We now provide numerical results showing the validity of the flea gas method to calculate the dynamics of the entanglement entropy after a generic quench in integrable systems. In section 5.1 we present some preliminary benchmarks of the approach. In section 5.2 we discuss the bipartite inhomogeneous quench depicted in Fig. 1 (b). Finally, in section 5.3 we discuss the mutual information between two intervals.

5.1 Preliminary benchmarks

A crucial feature of the flea gas dynamics is that it gives rise to the correct dressing of the group velocities of the quasiparticles given by (23). While this can be proven for the flea gas algorithm for the Lieb-Liniger model, this is not the case for the XXZ chain. Here we provide numerical evidence that, at least for the quenches that we consider, the flea gas dynamics presented in Section 4 gives rise to the correct dressing of the group velocities.

We first consider the quench from a homogeneous chain prepared in the Néel state. Our results are presented in Fig. 4 (a) and (b). The results are for the XXZ chain with $\Delta=1.5$. The figures show the dressed velocities of the first two strings $v_{\alpha,\lambda}$ with $\alpha=1,2$ plotted versus the rapidity $\lambda$. The full lines are the flea gas results. These are obtained as $\Delta x_{\gamma}/t$, where $\Delta x_{\gamma}$ is the displacement of the particles with respect to their initial position. The data are averaged over $10^{4}$ realizations of the flea gas dynamics. As it is clear from the Figure, there are large fluctuations in the central region around $\lambda=\pi/2$. This is because the density $\rho_{1,\lambda}$ is large at the edges of the interval, whereas it is suppressed at the center, for instance, for $\Delta=1.5$ by a factor of $\sim 50$. This effectively reduces the statistics for the central rapidities. For $\alpha=2$ the density $\rho_{2,\lambda}$ has a maximum around $\lambda=\pi/2$. However, it is in general much smaller than $\rho_{1,\lambda}$, again resulting in large fluctuations for the group velocities of the two strings. In Fig. 4 the dashed-dotted lines are the analytical results for the dressed velocities, which are obtained by solving numerically (23). Clearly, the agreement with the flea gas results is very good.

We also considered the dressed velocities in homogeneous thermal states. Panels (c) and (d) show results for the dressed velocities in the state described by the thermal density matrix

[TABLE]

where $\beta$ is the inverse temperature and $h$ a transverse magnetic field. The data are for $\beta=0.5$ and $\beta h=0.25$. The continuous lines are flea gas results for the dressed velocities of the first two strings, which perfectly match the analytical results of TBA (dashed-dotted lines).

We now move the quenches from piecewise homogeneous states. Here we consider the initial density matrix as

[TABLE]

where quantities with the subscript $L/R$ refer to the left and right chains (see Fig. 1 (b)). The quench from (32) was investigated in Ref. [76] using GHD. Here we consider $\beta_{\mathrm{L}}=0,\,(\beta h)_{\mathrm{L}}=1$, and $\beta_{\mathrm{R}}=0,\,(\beta h)_{\mathrm{R}}=2$, with $h$ being the magnetic field ($\Delta=2$). Due to the inhomogeneous initial condition now the dressed velocities depend on $\zeta\equiv x/t$ (see Fig. 1 (b)). To check the validity of the flea gas method, in principle, one has to check that the flea gas gives the correct result for $v_{\alpha,\lambda}(\zeta)$ for any $\zeta$. Here, instead, we consider the space-time dependence of the local energy density $E(\zeta)$ and magnetization $S^{z}(\zeta)$ plotted versus $\zeta=x/t$, with $t$ the time after the quench, and $x$ the distance from the origin of the lightcone. Both the quantities for $x,t\to\infty$ become functions of $\zeta$. In Figure 4 (e) and (f) the square symbols are the results of the flea gas simulation for a chain with $L=2000$ and $t=100$, whereas the full lines are the analytical results obtained in Ref. [76] by solving the GHD equations. The agreement between the flea gas and the GHD results is spectacular.

As a further check of the validity of the flea gas method we now discuss results for the dynamics of the von Neumann entanglement entropy after a quench from homogeneous initial states, for which analytical results (cf. Eq. (2)) are available. Our results are discussed in Fig. 7. The figure shows data for the XXZ chain with $\Delta=2$, quenching form the Néel state (see section 2). The rescaled entropy $S/\ell$ is plotted versus $t/\ell$, with $\ell$ the subsystem size. In the simulation we considered $\ell=100$ and a chain of length $L=2000$. The data are obtained by averaging over $\sim 10^{4}$ independent realizations of the flea gas dynamics. The continuous line is the flea gas result (30) up to $t/\ell\approx 1.5$, although results for larger times can be easily obtained. The dashed-dotted line is the analytical result (2) obtained in Ref. [32]. The agreement between the two is excellent.

Some remarks are in order. First, the flea gas picture is expected to capture correctly only the ballistic part of the entanglement dynamics, i.e., the leading behavior in $t/\ell$. Note that, however, subleading corrections, for instance diffusive corrections as ${\mathcal{O}}(\sqrt{t})$, are generically expected in the entanglement dynamics. In the flea gas framework diffusive corrections arise because of the average over the different realizations of the initial state, and are associated with the fluctuations of the particles trajectories. On the other hand, the diffusive corrections that are present in the flea gas are not expected to be the same as the quantum diffusive corrections of the XXZ chain. The origin of diffusion in interacting integrable models and in the flea gas have been under constant investigation in the last few years [88, 89, 90, 91, 87]. Finally, as it is clear from Fig. 7, subleading corrections are small. Only for very short times some deviations from (2) are present, which disappear in the scaling limit $t,\ell\to\infty$.

5.1.1 Effect of the wrong collisions

In section 4.1 we stressed that the fact that the group velocity of the quasiparticles are not monotonic increasing functions of the bare ones. This implies that $\mathrm{sign}(v_{\alpha,\lambda}-v_{\beta,\mu})\neq\mathrm{sign}(v_{\alpha,\lambda}^{\scriptscriptstyle\mathrm{bare}}-v_{\beta,\mu}^{\scriptscriptstyle\mathrm{bare}})$. As a consequence two quasiparticles colliding with velocities $v_{\alpha,\lambda}^{\scriptscriptstyle\mathrm{bare}}$ and $v_{\alpha,\lambda}^{\scriptscriptstyle\mathrm{bare}}$ are shifted by the “wrong” distance $\mathrm{sign}(v_{\alpha,\lambda}^{\scriptscriptstyle\mathrm{bare}}-v_{\beta,\mu}^{\scriptscriptstyle\mathrm{bare}})d_{\alpha,\beta}(\lambda,\mu)$.

Here we investigate this effect. We consider the initial state defined by the thermal density matrix (31) with $\beta=0.5$ and $\beta h=0.25$. We consider the XXZ chain with $\Delta=1.5$. The reason is that for $\Delta\to 1$ the dressing of the quasiparticles velocities is larger, which should enhance the effect of the wrong scatterings.

In Fig. 5 we compare the bare velocities and the dressed ones (dotted and continuous line, respectively). We show results only for $\alpha=1,2$. The effect of the dressing is clearly visible in the figure. Importantly, one consequence of the dressing is that the maximum of the velocities are shifted, as compared with the bare ones. This already implies that for some values of $\alpha,\lambda$ and $\beta,\mu$ one has that $\mathrm{sign}(v_{\alpha,\lambda}-v_{\beta,\mu})\neq\mathrm{sign}(v_{\alpha,\lambda}^{\scriptscriptstyle\mathrm{bare}}-v_{\beta,\mu}^{\scriptscriptstyle\mathrm{bare}})$, meaning that a priori the flea gas dynamics is not fully equivalent to the GHD. On the other hand, the behavior of the bare and dressed velocities is similar as a function of $\lambda$, suggesting that the effect of the wrong scatterings should be “small”.

In Fig. 6 we investigate the effect of the wrong collisions on the shift of the quasiparticles trajectories. In panel (a) we show the average total shift $\Delta x$ experienced by the quasiparticle with rapidity $\lambda$ due to wrong scatterings (continuous line) and the total number of scatterings (dotted line). The results are for the quasiparticlew with $\alpha=1$, i.e., the unbound quasiparticles. We observe that the wrong scatterings have a small effect on $\Delta x$, which is barely visible in the figure. Similar behavior is observed for the two-particle bound states. The results are shown in panel (b). In panel (c) we show the effect of the wrong scatterings on the trajectories of the quasiparticles with $\alpha=1$. The figure shows the same data as in (a). The average shift due to the wrong scatterings is $\approx 10^{-3}$.

Finally, two observations are in order. First, although panel (c) suggests a finite contribution of the wrong scattering to $\Delta x$, larger numerical simulations would be needed to ensure that the data are in the scaling limit $N,L\to\infty$ with $N$ the number of quasiparticles. Second, in principle, it should be possible to correct the effect of the wrong scatterings by imposing that upon colliding the displacement of the quasiparticle trajectories is the correct one $\mathrm{sign}(v_{\alpha,\lambda}-v_{\beta,\mu})$. Note that this requires knowing the dressed velocities $v_{\alpha,\lambda}$, which could be calculated during the simulation.

5.2 Entanglement dynamics after a quench from inhomogeneous initial conditions

Having established the validity of the flea gas method to simulate the entanglement dynamics after homogeneous quenches, we now consider the case of the inhomogeneous initial state in Fig. 1 (b). The calculation of the entanglement dynamics within the GHD framework is in general a complicated task. Explicit analytic results can be obtained only in few cases. For instance, the steady-state value of the von Neumann entanglement entropy for a finite subsystem placed next to the interface between the two chains (see Fig. 1 (b)) can be easily calculated. This corresponds to the limit $\ell/t\to 0$. In this limit, the entire subsystem is described by the GGE with $\zeta=0$ (see section 3.3). Following Ref. [32], the density of the steady-state von Neumann entanglement entropy coincides with that of the GGE entropy with $\zeta=0$. One has [52]

[TABLE]

Here $s_{\alpha,\lambda}(0)$ is the Yang-Yang entropy (cf. (14)) of the GGE with $\zeta=0$, which is obtained by using the GHD (see section 3.3). The result does not depend on which side of the system one places the interval, as long as $\ell$ is finite. Interestingly, one can show that the $\zeta=0$ macrostate describes the entanglement growth at short times, i.e., the limit $1\ll t\ll\ell$ as well [52, 54]. First, the entanglement entropy is expected to grow linearly at short times. Here we refer to the slope of the linear growth as the entanglement production rate [52, 46, 54]. The entanglement growth is due to the quasiparticles that cross the interface between the two chains. This suggests that the entanglement production rate is described by the $\zeta=0$ GGE. Indeed, if subsystem $A$ is the semi-infinite chain, the entanglement production rate $S/t$ is given as [54]

[TABLE]

Here $v_{\alpha,\lambda}(0)$ is the group velocities of the particle-hole excitations around the $\zeta=0$ GGE, which are obtained from (23). For a finite subsystem, the slope of the linear growth depends on the details of the bipartition. For simplicity we now consider the case of a finite interval of length $\ell$ placed in one of the two chains next to the interface (see Fig. 1 (b)). Clearly, the entanglement entropy gets contributions from both the edges of the subsystem. For short enough times but still in the linear regime, i.e, for large $t$ with $t/\ell\ll 1$, the contributions of the two edges decouple and can be summed independently. As in (34), one of the edges of $A$ is described by the GGE with $\zeta=0$. On the other hand, the other one is described by the GGE with $\zeta=\pm\infty$, depending on which side subsystem $A$ is placed in. The entanglement production rate is given as

[TABLE]

where $\sigma=\pm$ identifies the side in which subsystem $A$ is placed.

To illustrate how these features emerge in the flea gas simulations, in Fig. 8 we present numerical results for the quench from the initial state obtained by joining the Néel state and the dimer state $|\mathrm{N}\rangle\otimes|\mathrm{D}\rangle$ (see section 2 for the definition of these states, and section 3.2 for their TBA treatment). The results are for the XXZ chain with $\Delta=2$. The full and dotted lines correspond to the bipartitions with interval $A$ being $[-\ell,0]$ (in the Néel region) and $[0,\ell]$ (in the dimer region) respectively. In both cases we consider $L=2000$ and $\ell=100$. The results are obtained by averaging over $10000$ realizations of the “flea” gas dynamics (see section 4), and using (30).

For the quench from $|\mathrm{N}\otimes\mathrm{dimer}\rangle$, we observe that at $\Delta=2$ one has $s_{\alpha,\lambda}(+\infty)\approx s_{\alpha,\lambda}(-\infty)$ and $v_{\alpha,\lambda}(+\infty)\approx v_{\alpha,\lambda}(-\infty)$. From (35) one obtains that the entanglement production rate depends very mildly on which region the subsystem is placed. The theory predictions (35) for the entanglement production rates are not distinguishable on the scale of the figure and are reported as dashed-dotted line.

At intermediate $\zeta=t/\ell$, the entanglement entropy depends on all the values of $\zeta$. This happens because the entangling quasiparticles explore macrostates with different $\zeta$ as they travel in subsystem $A$ (see Fig. 1). Although it is possible, in principle, to write an analytic formula [54] for the evolution of the entanglement entropy at any time, its numerical evaluation is a demanding task. In contrast, the flea gas method allows to access easily the full-time entanglement dynamics, as it is clear from Fig. 8.

In Fig. 9 we present further checks of the validity of the flea gas method for inhomogeneous quenches. We consider the initial state obtained by joining the tilted Néel state and the dimer state, i.e., $|\mathrm{N},\theta\rangle\otimes|\mathrm{dimer}\rangle$ (see section 2), where $\theta$ is the tilting angle. Panel (a) and (b) show results for $\theta=\pi/3$, whereas in (c) and (d) we consider $\theta=\pi/6$. In all the cases we choose $\ell=100$ and total system size $L=2000$. In (a) and (d) the subsystem is placed on the Néel side ($A=[-\ell,0]$), whereas in (b) and (c) is in the dimer side ($A=[0,\ell]$). The fact that the production rate depends on the position of the interval is now apparent. The dashed-dotted lines are the theory predictions (cf. (35)) for the entanglement production rates. In Fig. 9 (a) some deviations from (35) are visible. These, however, are due to finite-size and finite-time effects. In the inset of Fig. 9 we report results for $\ell=500$, which are now in perfect agreement with (35). We observe that in general very large subsystems are needed to provide a robust numerical check of the GHD prediction (35).

5.3 Mutual information after quenches from inhomogeneous initial conditions

It is interesting to investigate the dynamics of the mutual information between two intervals. To this purpose, we now consider a tripartite system. Subsystem $A$ is made of two intervals $A_{1}$ and $A_{2}$ at a distance $d$. Here we consider only $d=0$, although the method works also for $d>0$. The two subsystems are embedded in an infinite system. The mutual information $I_{A_{1}:A_{2}}$ is a measure of the correlation shared between $A_{1}$ and $A_{2}$, although it is not a proper measure of the entanglement between them. $I_{A_{1}:A_{2}}$ is defined as

[TABLE]

where $S_{A_{1}}$, $S_{A_{2}}$, and $S_{A_{1}\cup A_{2}}$ are the von Neumann entanglement entropies of $A_{1}$, $A_{2}$ and $A_{1}\cup A_{2}$ with the rest of the system.

In the quasiparticle picture, the mutual information is proportional to the entangled pairs that are shared only between $A_{1}$ and $A_{2}$. On the other hand, the contribution of the quasiparticles to the mutual information is, again, the GGE thermodynamic entropy. Thus, the flea gas formula for $I_{A_{1}:A_{2}}$ is the same as (30) where the sum is restricted to the pairs of quasiparticles shared between $A_{1}$ and $A_{2}$.

The qualitative behavior of the mutual information is as follows. For two disjoint intervals at a distance $d$, the mutual information is zero at short times. At $t\sim d/t$, $I_{A_{1}:A_{2}}$ exhibits a linear increase. This corresponds to entangled pairs starting to be shared between $A_{1}$ and $A_{2}$. The growth persists up to $t\sim(d+\ell)/t$, when the mutual information starts to decrease. In systems with only one quasiparticle with perfect linear dispersion (as in CFT systems), the decrease is linear. In generic integrable models a much slower decrease is observed [46]. This is due to the fact that quasiparticles have a nontrivial dispersion, and slow quasiparticles entangle the two subsystems at long times. The mutual information can, in principle, be used as a tool to reveal the quasiparticle content of an integrable model. Typically, different species have different maximum velocities $v_{\alpha,M}$. This implies that if the distance between the two intervals is large enough, the mutual information exhibits a multi-peak structure in time, each peak corresponding to a different species [20, 95].

The mutual information after quenches from inhomogeneous initial states has not been investigated yet. In contrast with homogeneous global quenches [46], deriving the quasiparticle picture for the mutual information in inhomogeneous settings is a formidable task. Again, the reason is that the quasiparticles trajectories are non trivial functions of time. We now show that the flea gas approach allows to simulate effectively the full-time dynamics of the mutual information. We restrict ourselves to the case of two adjacent intervals, although the method works for disjoint intervals as well.

We present our results in Fig. 10, focusing on the XXZ chain with $\Delta=2$. The initial state that we consider is $|\mathrm{N},\theta\rangle\otimes|\mathrm{D}\rangle$. Different panels show different values of $\theta$. The data are for two equal-length adjacent intervals $[-\ell,0]$ and $[0,\ell]$ with $\ell=100$. The total chain length is $L=2000$. As expected, the mutual information is initially zero, it grows linearly at intermediate times, and it eventually decays to zero at asymptotically long times. Importantly, the initial slope of the mutual information depends only on the GGE with $\zeta=0$ because the interface between $A_{1}$ and $A_{2}$ is at the origin. In particular, the slope of the inital growth of the mutual information coincides with the entanglement production rate for the two semi-infinite chains (see Eq. (34)). This initial growth is reported in Fig. 10 as dashed-dotted line, and it perfectly describes the behavior of the flea gas results.

6 Conclusions

In this work we showed that the so-called flea gas method put forward in Ref. [92] provides a versatile tool for simulating the entanglement dynamics after quenches from generic initial states in integrable systems. We benchmarked the method in the Heisenberg XXZ chain, although it can be applied, in principle, to any integrable model. The method works for arbitrary initial states, both globally homogeneous as well as piecewise homogeneous. For globally homogeneous quenches the approach requires only the GGE macrostate that describes the steady-state. For piecewise homogeneous states, the key ingredients are the GGE macrostates describing the steady-state in the bulk of the two systems. Although in this case the entanglement dynamics can be obtained, in principle, by combining the quasiparticle picture with the Generalized Hydrodynamics, obtaining explicit formulas [54] is a demanding task because the trajectories of the quasiparticles are non trivial functions of time. Indeed, results can be obtained only in some limits. In contrast, in this work we showed that the flea gas approach allows to obtain easily the full-time dynamics of the entanglement entropy and of the mutual information between two intervals. Thus, the method paves the way to the study of entanglement dynamics using “molecular dynamics” simulations.

Our results open several possible new research directions. First, it would be important to investigate whether it is possible to prove analytically that for the XXZ chain the flea gas dynamics as described in section 4.1 gives the correct dressing for the group velocities.

Also, it would be useful to apply the method to more complicated setups, such as multipartite systems, or different initial states. Also, it would be important to go beyond the ballistic regime, studying corrections to the linear entanglement growth. This requires first to understand the subleading diffusive corrections in the flea gas method. Second, it requires to modify the flea gas dynamics to correctly reproduce the diffusive corrections that arise from the quantum fluctuations [88, 89]. An interesting direction would be to generalize the flea gas approach to study the entanglement dynamics in the presence of defects or impurities. Finally, it would be enlightening to understand whether it is possible to treat the entanglement of operators in integrable spin chains by using the flea gas approach, generalizing the results of Ref. [96] for the Rule $54$ chain.

7 Acknowledgements

MM acknowledges support from ERC under Consolidator grant number 771536 (NEMO). VA acknowledges support from the D-ITP consortium, a program of the NWO. VA also acknowledges support from the European Research Council under ERC Advanced grant 743032 DYNAMINT.