Entanglement evolution and generalised hydrodynamics: interacting integrable systems

We investigate the dynamics of bipartite entanglement after the sudden junction of two leads in interacting integrable models. By combining the quasiparticle picture for the entanglement spreading with Generalised Hydrodynamics we derive an analytical prediction for the dynamics of the entanglement entropy between a finite subsystem and the rest. We find that the entanglement rate between the two leads depends only on the physics at the interface and differs from the rate of exchange of thermodynamic entropy. This contrasts with the behaviour in free or homogeneous interacting integrable systems, where the two rates coincide.

. Here x is measured from the interface between the two chains. In this setting we are interested in the entanglement entropy of a region A of length . A semiclassical description in terms of free quasiparticles applies to the entanglement entropy. In the quasiparticle pictures pairs of entangled quasiparticles are created in the bulk of the two chains. Quasiparticles forming an entangled pair have rapidity of opposite sign and initially travel with opposite velocities. different effective descriptions have been put forward in the extreme cases of integrable [8] and chaotic [32] systems, which allow one to recover quantitatively the dynamics of entanglement entropies in the limit of large times and subsystem sizes, also called "space-time scaling limit". These descriptions are respectively known as "semiclassical quasiparticle picture" and "minimal membrane picture". In this paper we are interested in the entanglement spreading in integrable models, thus we focus on the former. In the semiclassical quasiparticle picture one interprets the initial state as a source of entangled pairs of quasiparticles, generated uniformly in space at the initial time and propagating as free classical objects with opposite velocities. Quasiparticles produced at the same point in space are mutually entangled, whereas quasiparticles created further apart are not. The entanglement between a subsystem A and the rest is proportional to the total number of quasiparticles created at the same point and shared between A and its complement at time t. Assuming that there are N s species of particles, whose dispersion relation is characterised by a real "rapidity" λ, we find where χ A (x) is the characteristic function of the interval A, the function f α,λ (x) denotes the contribution to the entanglement of the pair (α, λ), (α, −λ), with rapidity ±λ and species α, originated at x, and X α,λ (x, t) gives the position at time t of the quasiparticle that was created at position x at time 0. This picture can be extended to the case where the initial state consists of more complicated multiplets of correlated particles [25,26]. For the sake of simplicity, however, here we stick to cases where only pairs are produced and, moreover, we assume that the pairs are composed by particles of the same species.
The semiclassical quasiparticle picture, as described above, only gives qualitative predic-tions. To make it quantitatively accurate one needs to specify what are the semiclassical particles responsible for the entanglement spreading, what are their trajectories, and what is the contribution of a single entangled pair to the total entanglement. In translationally invariant interacting integrable systems, this task has been accomplished in Ref. [31]. This has been done by exploiting the exact knowledge of the stationary state describing finite subsystems at large times [3,4] to complement the semiclassical quasiparticle picture. In particular, Ref. [31] assumed that the entangling quasiparticles are the elementary excitations on the stationary macrostate. This gives where v α,λ is the group velocity of the excitation (α, λ) and S Y Y α,λ is its contribution to the thermodynamic entropy. The quantitative prediction of Ref. [31] is conjectured to apply for all integrable models treatable via thermodynamic Bethe ansatz and to become exact in the space time scaling limit (at least for a class of "integrable" initial states [49]).
In this paper we derive an analogous quantitative prediction for cases where the initial state is not homogeneous. The idea is to use the recently developed theory of Generalised Hydrodynamics (GHD) [50,51] (see Sec. 3) to determine the state of the system at large times and use it to complement the semiclassical quasiparticle picture. In particular, we focus on the paradigmatic case of initial states formed by the junction of two macroscopically different homogeneous states (leads) and determine the time evolution of the entanglement entropy of A = [x 1 , x 1 + ] (see Fig. 1) in the space-time scaling limit. Our prediction applies to generic integrable models and is tested in the concrete case of the anisotropic spin-1/2 Heisenberg chain.
Note that the entanglement spreading in inhomogeneous settings has been already considered in the recent literature, see e.g. Refs. [52,53,54,55,57,58]. In particular, Refs. [54] and [57] considered exactly the problem studied here. Both these references, however, considered special cases. Ref. [57] investigated free systems, while Ref. [54] focussed on x 1 = 0 and considered the limits /t → 0 and /t → ∞. Our work represents a non-trivial extension of these studies.
Specifically, our prediction displays a remarkable novel effect due to the combination of inhomogeneity and interactions. This is most easily explained by considering the entropy of one of the leads, namely A = [0, ∞[, and comparing it with the homogeneous and the free cases. In particular, in the homogeneous case one has In words: the entanglement entropy increases with the same rate at which the two leads exchange thermodynamic entropy. Ref. [57] confirmed that, if the system is free, this feature remains true also in the inhomogeneous case. This is because, in both these cases, the trajectories of quasiparticles are straight lines; consequently, only one of the two entangled particles in a given pair can cross the junction, and one can forget about the pair structure. Here, instead, inhomogeneity and interactions cause the trajectories to curve, making it possible to both particles of a given entangled pair to cross the junction, and, in turn, slow down the rate of entanglement growth. In particular, this implies that the conjecture for the entanglement production rate put forward in Ref. [54] is generically incorrect (although the correction is typically small).
The rest of the manuscript is organised as follows. In Section 2 we review the Thermodynamic Bethe Ansatz treatment of generic integrable models. In Section 3 we summarise the GHD formalism for quenches from piecewise homogeneous initial states. In Section 4 we identify the entangling quasiparticles and present a detailed analysis of their trajectories. In Section 5 we derive and simplify the quasiparticle picture prediction for the full-time dynamics of the entanglement entropy. In Section 6 we analyse the entanglement production rate for the semi-infinite chain. In Section 7 we provide numerical checks on the validity of our results, by presenting tDMRG data for several quenches in the XXZ chain. Finally, in Section 8 we draw our conclusions. Three appendices contain various technical details of our derivations.

Thermodynamic Bethe Ansatz treatment
In this work we consider interacting integrable quantum many-body systems describable by Thermodynamic Bethe Ansatz [59] (TBA). This description applies to a variety of integrable models both in the continuum and on the lattice. In most of this paper we will keep the discussion at a general level, without specifying any concrete model. In Sec. 7 our results will be tested against numerical simulations in the paradigmatic example of the "gapped" XXZ spin-1/2 chain where we denoted by σ β j (β = x, y, z) the Pauli matrices at position j, by L the volume of the system, and by ∆ > 1 the "anisotropy".
Let us now briefly summarise the main aspects of the TBA description that are needed in the rest of the paper. When a TBA description applies, in the thermodynamic limit L → ∞ the eigenstates of the Hamiltonian are characterised by a set of functions where the real variable λ ∈ C ⊂ R is customarily called "rapidity" and the number N s is the "number of species". The number N s , and the domain C depend on the details of the specific model considered. For instance, in the XXZ chain with ∆ > 1 one has N s → ∞ and can choose C = [−π/2, π/2]. The functions {ρ α,λ } are known as "root densities", and can be interpreted as rapidity distributions of the system's stable quasiparticles. The rapidity λ parametrises the quasiparticles' dispersion relation while the index α labels different families. For instance, in the XXZ chain, α = 1 corresponds to magnon-like excitations, whereas quasiparticles with α > 1 can be thought of as bound states of α magnons. A quasiparticle of the species α with rapidity λ has energy e α (λ) and quasimomentum p α (λ) given by [59] p α (λ) = 2 arctan tan(λ) tanh(αη/2) , e α (λ) = −J sinh η sinh (αη) cosh(αη) − cos(2λ) + 2hα .
In other words, the root densities generalise to interacting integrable models the notion of momentum occupation numbers in free systems.
Expectation values of local operators can be expressed as functional of the root densities {ρ α,λ }. For instance, the densities q x of local conserved charges, i.e., operators with local densities commuting with the Hamiltonian, are expressed as simple linear functionals as follows Here the "bare charges" q β,µ are functions specifying the charge density, and we denoted by |{ρ α,λ } a generic eigenstate characterised by the set of root densities {ρ α,λ }. The correspondence between eigenstates and root densities is generically not one-to-one: a large number of eigenstates corresponds to the same set of root densities. This fact is usually referred to by saying that the root densities specify a "thermodynamic macrostate" of the system, while the eigenstates of the Hamiltonian correspond to its "microstates". Quantitatively, for a finite system of size L there are ∼ exp[L α dλ S Y Y α,λ ] eigenstates that in the limit L → ∞ are described by the same set of densities {ρ α,λ }. Any of these eigenstates can be regarded as a finite-size representative eigenstate of the thermodynamic macrostate. Here we introduced the Yang-Yang entropy density where the auxiliary functionals {ρ t α,λ } called "total root densities" are defined as [59] ρ t α,λ = a α,λ + Ns β=1 dµ T α,λ;β,µ ρ β,µ .
The precise form of a α,λ and T α,λ;β,µ depends, again, on the specific model considered. In general, however, the kernel T α,λ;β,µ encodes all the information about the two-particle scattering matrix (the only non trivial one in integrable models). For instance, for the gapped XXZ spin-1/2 chain we have [59] a α,λ = 1 π where we introduced η = cosh −1 ∆ > 0. The total root densities are interpreted as the densities of all possible values that the rapidities can take (the constraint originates from the fact that in finite volume the rapidities obey a set of non-trivial quantisation conditions [59]). These densities are generically not uniform and, as a consequence of the non-trivial interactions, depend on the root densities.
To conclude, we note that the TBA description can be used also for some mixed states. This is true every time a generalised microcanonical representation applies [3,5]. In other words, if the expectation values of local observables in the mixed state can be reproduced, in the thermodynamic limit, by expectation values on a single (carefully chosen) eigenstate of the Hamiltonian. For example, the TBA description can be used to describe systems in Gibbs and Generalised Gibbs states [3,5].

GHD description of the local quasi-stationary state
After being initialised in an inhomogeneous stateρ 0 (see Fig. 1 for an example), the system performs a non-trivial time evolution, during which the expectation values of local observables display fast oscillations in both position x and time t. At large times, however, these fast oscillations dephase away, and the expectation values become slow functions of x and t. In this regime, it is reasonable to expect that the expectation values can be described by a quasistationary stateρ s (x, t) retaining some slow dependence on position and time. Namely, we expect tr where H is the Hamiltonian of the system, O x is a generic observable localised around the point x, andρ s is the density matrix describing the quasi-stationary state. In the context of quantum non-equilibrium dynamics the emergence of such a state was first proposed in Ref. [60], where it was called locally quasi-stationary state (LQSS). Specifically, based on the intuition developed for homogeneous quenches [2,3,4,5], it was argued that, at fixed (x, t), the stateρ s (x, t) is a GGE constructed with the charges of the time evolving Hamiltonian. This means thatρ s (x, t) is homogeneous, stationary, and admits a "microcanonical" representation in terms of a TBA representative eigenstate, or, equivalently, of a set of root densities {ρ α,λ (x, t)}. Note that the densities depend on space and time. Determiningρ s (x, t) without solving the full non-equilibrium dynamics (in the presence of integrable interactions) is the key result of the theory of generalised hydrodynamics (GHD) introduced in Refs. [50,51]. Specifically it was shown that, at the leading order in x and t, the position-dependent root densities fulfil the following continuity equation The quantity v α,λ (x, t) appearing in (13) is the velocity of the elementary excitations on the state described by {ρ α,λ (x, t)}, see Ref. [61], and it is defined through the following integral equation where The model-dependent function v b α,λ is the velocity of the excitations on the "vacuum state" (the state with ρ α,λ = 0) and is known as "bare velocity". For the gapped XXZ spin-1/2 chain we have where a α,λ is given in Eq. (11). Note that for non-trivially interacting models, i.e. when T α,λ;β,µ = 0, the velocity v α,λ (x, t) depends on the densities {ρ α,λ (x, t)}. This makes the equation (13) highly non-trivial.
The simplification introduced by Eq. (13) is remarkable. To determine the late-time properties of an integrable quantum many-body system one needs to solve a system of differential equations whose number is proportional to the system size, instead of solving the Schrödinger equation, which has instead the dimension of the Hilbert space. There is, however, a remaining non-trivial step to make before a solution can be obtained: one has to impose the initial conditions for ρ α,λ (x, t). This has been successfully done in a number of cases [50,51,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80], including the bipartite quench protocol considered here (see below), but in general the problem is still open. We note that, very recently, it has been shown that GHD provides the precise framework to describe experiments with trapped cold atoms [62].
Equation (13) admits a very simple interpretation in terms of a "kinetic theory" of free classical particles moving in an inhomogeneous background. One regards ρ α,λ (x, t) as the distribution function for classical particles of the species α = 1, . . . , N s with momentum λ, at position x, and at time t. Equation (13) describes the evolution of the distribution functions ρ α,λ due to the motion of the particles, which are nothing but a "coarse grained version" of the stable interacting quasiparticles characterising integrable models. Indeed, at the leading order in (x, t), the only effect of the interaction is a renormalisation of the group velocity v α,λ (x, t). Note that, at sub-leading orders, the effect of interactions might spoil this interpretation [79,80,81,82].
Finally, it is convenient to observe that, given a set of quantities {g α,λ (x, t)} fulfilling (13), we have Equation (17) has a material-derivative form and it is generically easier to solve than (13), for instance by using the method of characteristics [65,66]. This implies that, instead of solving (13), it is convenient to define the "so-called" filling functions For the upcoming analysis it also important to note that the (x, t)-dependent Yang-Yang entropy densities fulfil a continuity equation of the form (13) as it readily follows from (13) and (17).

Bipartite quench
Let us now specialise the GHD formalism of the previous section to what will be referred to as a "bipartite quench" (see Fig. 1). This is the time evolution of a state that, up to irrelevant 1 corrections localised around the junction, has the form withρ L(R) two macroscopically different homogeneous states. In our numerical tests in the XXZ chain, we considerρ where we defined the "Dimer state" |D , the "tilted Néel" state |N, θ , and the "tilted ferromagnetic state" |F, θ as follows In this setting, ifρ L(R) have cluster decomposition properties, at long enough times all the quantities generally become functions of the "ray" ζ = x/t [60,50,51]; in the limit t → ∞ the LQSS on each ray becomes exactly stationary, and all the sub-leading corrections to (13) vanish. Rewriting (17) In this case, whenever information propagates with a bounded velocity, it is possible to impose the initial conditions in (26) at ζ → ±∞ and solve it [60,50,51]. Indeed, at infinite distances from the interface between the two leads there are regions where no information on the inhomogeneity can arrive, and local observables evolve as if the system were homogeneous. This means that their expectation values are described by stationary states that can be computed using the standard techniques developed for homogeneous quenches [2,3,4,5]. These stationary states provide the boundary conditions ϑ for (26). We then find where θ H (x) is the step function. Equation (27) is only an implicit solution because it depends on the velocity v α,λ (ζ) that in turn depends on ϑ α,λ (ζ). To find ϑ α,λ (ζ), it is convenient to adopt an iterative approach, combining (27) with the infinite time limit of (14) and (15). Namely Finally, for the upcoming discussion it is useful to mention that the continuity equation for the Yang-Yang entropy in terms of ζ reads as Figure 2: Velocity field v α,λ after a quench in the XXZ chain for ∆ = 10, plotted against ζ ≡ x/t. The pre-quench initial state is obtained by joining the Néel state |N, 0 (24) (left) and the ferromagnet |F, 0 (25) (right). Here we only show results for α = 1. Different panels correspond to different values of rapidity λ. Note that for ζ < v min ≈ −2 and ζ > v max ≈ 1 the velocity field does not depend on ζ, as expected.

Entangling quasiparticles in inhomogeneous backgrounds
The first step to turn quantitative the quasiparticle picture for the entanglement dynamics is to identify the entangling quasiparticles and determine their trajectories. Here we perform this task in the case of interacting integrable models after bipartite quenches. The basic assumption of this work is that the entangling quasiparticles are the excitations on the locally quasi-stationary state. The logic of this assumption is that, since the stable quasiparticle excitations in integrable models are responsible for the spreading of any kind of information, it is natural to argue that they also spread the entanglement. An analogous assumption has been formulated in Ref. [31] in the homogeneous case and in Ref. [57] for free systems.
To be specific, we denote by X α,λ (x, t) the position at time t of a particle of species α = 1, . . . , N s with rapidity λ that started at position x. The position is measured from the interface between the two chains (see Fig. 1). We then assume that X α,λ (x, t) is determined by the following classical equation of motion where v α,λ (x, t) are the same velocities as in (13). Note that the trajectories generated by (31) are generically not straight, in contrast to the homogeneous [31] and free [57] cases. For bipartite quenches (see Fig. 1), v α,λ (x, t) are obtained by solving the GHD equations (27)- (29) and are functions of The function v α,λ (ζ) is taken to be fairly generic; we only make one assumption on its ζ dependence: Assumption 1. For fixed λ and α, the equation ζ − v α,λ (ζ) = 0 has a unique solution, which we call ζ α,λ .
This is the standard assumption of GHD. It has been verified in all the examples examined up to now, but it has not yet been proven rigorously. A number of interesting properties of the velocity field follow from (27)- (29). For example one can show (see Appendix A) The first property stems from the bounds on the velocity at which information propagates from the interface to the bulk of the two semi-infinite chains; in space-time regions outside the lightcone spreading from the origin, the system is described by the macrostates ρ (L/R) α,λ . A concrete example of a ζ-dependent velocity field obtained in the XXZ model is reported in Fig. 2.
We remark that the quasiparticle picture is based on coarse grain in space, momentum, and, in turn, time, so the coordinates x and t in (31) have some intrinsic indetermination. Just to keep track of potential effects of that, we naively capture such corrections with a single parameter t 0 , which is regarded as the time when the initial conditions for the classical problem (31) Keeping a finite t 0 > 0 regularises the initial value problem (31)(33). Indeed, properties (i) and (ii) ensure that v α,λ (x/t) is a Lipschitz continuous function of x for all t > 0. This implies that Cauchy's Theorem applies and a unique solution exists for any given initial condition x. Eventually we will take the limit t 0 /t → 0. An immediate consequence of the uniqueness of the solution of (31)(33) is that, for fixed (α, λ), trajectories with different initial conditions cannot cross. In particular, it is useful to observe that quasiparticles with "initial velocity" ζ α,λ (cf. Assumption 1), follow linear trajectories, i.e., s α,λ (t) = ζ α,λ t, solves (31) with initial value ζ α,λ t 0 . This implies Equivalently, this means that a trajectory starting on the left or on the right of s α,λ (t) remains as such at any time.
To gain information on the qualitative form of a generic trajectory X α,λ (x, t), it is useful to explicitly integrate Eq. (31). This is done by means of the following convenient rewriting Integrating the expression (35) from t 0 to t we obtain After taking the limit t → ∞, the right-hand-side diverges logarithmically. The only way for the left-hand-side to match this behaviour is by having Figure 3: Pictorial representation of the functions Φ α,λ (ζ) (cf. (40)), its inverse Z α,λ (φ) (cf. (41)), and of J α,−λ (ζ) (cf. (63)). In (a) a quasiparticle is produced at time 0 at position ζt.
The function Φ α,λ (ζ) gives its ray at time t. In (b) Z α,λ (φ)t is the position at time 0 of the particle that at time t is on the ray φ. In (c) two quasiparticles with opposite rapidities λ and −λ are created. The function J α,−λ (ζ) gives the ray at time t of the quasiparticle whose partner with rapidity −λ is at ray ζ.
In words, this means that the trajectories of the quasiparticles become linear at asymptotically long times.
Crucially, Equation (36) allows us to explicitly determine the initial condition x in terms of the position of the quasiparticle at time t. After a straightforward derivation (see Appendix A), we find Here the first and the second terms account for quasiparticles created on the left and right lead (see Fig. 1), respectively. To treat the large time limit it is useful to define the scaling function such that Φ α,λ (ζ, t, t 0 )t is the position at time t of the particle that was originated at position ζt at time t 0 . From (38) it is straightforward to find Furthermore, it is useful to introduce the inverse function Z α,λ (φ, t, t 0 ) = Φ −1 α,λ (φ, t, t 0 ). Z α,λ (φ, t, t 0 )t gives the position at time t 0 of the particle (α, λ) that at time t has position φt. Using (40) we obtain From now on we take the limit t 0 /t → 0, dropping the dependence on t 0 in Φ α,λ and Z α,λ . Again, Φ α,λ (ζ)t is the position at time t of the quasiparticle that started at position ζt at time 0, while Z α,λ (φ)t is the position at time 0 of the quasiparticle that is at position φt at time t. The functions Φ α,λ and of Z α,λ are pictorially illustrated in Fig. 3. Using these definitions, one immediately has Z α,λ (Φ α,λ (ζ)) = ζ .
We mention that the function Z α,λ (φ)t has been already introduced in Refs. [65,66], where it is called "characteristics". An important property of Z α,λ (φ) and Φ α,λ (ζ) is that they are monotonic functions of their arguments. This is a direct consequence of the non-crossing condition (34) of the trajectories. To show the monotonicity of Z α,λ (φ) it is useful to observe that (41) implies Moreover, a trivial consequence of (34) is Equation (43) and (44) imply that Z α,λ and, in turn, its inverse Φ α,λ (ζ) are monotonic.

Flow of conserved quantities
In this section we show how the information about the quasiparticles' trajectories can be used to determine the flow of conserved quantities. The latter are represented by functions of x and t written in terms of a single-particle contribution g α,λ (x, t) that satisfies the continuity equation (13), namely where This definition includes the expectation values of all local and quasilocal conserved-charge densities, for which (cf. (7)) but also the Yang-Yang entropy, as demonstrated by Eq. (20).
Let us now specialise this relation to our setting and assume that all the functions of x and t become functions only of the "ray" x/t. In this case (48) becomes where we took X α,λ (a, t) = ζ α,λ t, x = ζt (51) and used (39). The l.h.s. of (50) can be directly computed using that g α,λ (ζ) fulfils the continuity equation (30). In particular we find Plugging this into (50) gives All this has a very simple physical interpretation. This equation states that, since (48) is conserved, g α,λ (Φ α,λ (ζ)) is renormalised by the factor to account for changes in the "volume" X α,λ (x, t) − X α,λ (a, t). For instance, if the trajectories are straight lines the volume is constant and, accordingly, we have Finally, note that (53) can be equivalently expressed as

Prediction for the entropy dynamics
In this section we derive the quasiparticle picture's prediction for the entanglement dynamics after a quench from a piecewise homogeneous initial condition (21). Specifically, we assume thatρ L(R) are both pure, low-entangled, states producing pairs of entangled quasiparticles. Namely, we assume that we can use (1). As noted in the previous section, for large enough times after a bipartite quench, we can make the following replacement (cf. the discussion after (41)) where the function Φ α,λ (ζ) is defined in (40). Plugging this expression into (1) and changing variables to the ray ζ = x/t we have where we introduced ζ 1,2 ≡ x 1,2 /t. The identification (57), however, is not sufficient. To make Eq. (58) predictive, we also have to determine the weight function f α,λ (ζ), i.e., the contribution of the pair of quasiparticles {(α, λ), (α, −λ)} to the entanglement entropy.
Here we conjecture that the entanglement entropy is given by the Yang-Yang entropy S Y Y α,λ (±∞) of the quasiparticles created initially in the two leads. This means As discussed in Sec. 3.1, S Y Y α,λ (±∞) is the density of Yang-Yang entropy of the stationary states reached after the quench from the homogeneous initial conditions on left and right leads. The ansatz (59) is motivated by two main observations. First, Eq. (59) has been verified for quenches from piecewise homogeneous initial states in free systems [57]. Second, for homogeneous quenches Eq. (59) holds true even in the presence of interactions [31].
Equation (58), supplemented with (59), gives a complete quasiparticle prediction for the dynamics of the entanglement entropy in the interacting case. It is, however, very complicated to evaluate, as it involves the determination of both v α,λ (ζ) and Φ α,λ (ζ) for all values of ζ. In the following we show that it is possible to drastically simplify Eq. (58). We begin by changing variables from ζ to φ = Φ α,λ (ζ) where the function Z α,λ (φ) is defined in (41) and Φ α,λ (ζ) = ∂ ζ Φ α,λ (ζ). Then we note where in the first step we used (43) and in the second (53). Finally, we use the identity (see Appendix B for the proof) In (62) we introduced the function J α,λ (ζ), which characterises the trajectory of entangled quasiparticles pairs. Specifically, J α,λ (ζ)t gives the position at time t of the quasiparticle (α, λ) starting at the same point as the particle (α, −λ) that at time t is at position ζt (see Fig. 3 for a pictorial representation). In formulae we have Note that one trivially has J α,λ (J α,−λ (ζ)) = ζ .
Putting all together we find the following expression for the entanglement entropy This expression can be simplified further by using the continuity equation (30) for the Yang-Yang entropy density. This allows one to perform the integral over φ explicitly The last step is achieved by means of the following identity (see Appendix B for the proof) Using (67) we can finally rewrite (66) as follows where we assumed ζ 1 < ζ 2 . Contrary to (58), (68) does not require to integrate over ζ; eq. (68) is one of the main results of this paper. The terms on the two lines of (68) correspond to processes where the quasiparticle (α, λ) crosses the boundary at ζ 1 and ζ 2 respectively. Let us look more closely at one of them, say the one on the first line. First we note that the contribution is positive only if the final configuration has only one of the two quasiparticles in the system, namely if the particle with rapidity λ enters the system and its companion is outside or if the particle of rapidity λ exits the system and its companion is inside. The contributions with final configurations featuring both quasiparticles inside or outside are instead weighted with a negative sign. This is in perfect agreement with our intuition: the entanglement increases when a pair is shared. Second, we note that the numerical value of the contribution is (cf. Eq. (61)) This is nothing but the number of pairs (α, λ), (α, −λ) for which (α, λ) crossed the border at x 1 before time t times f α,λ (Z α,λ (ζ 1 )), the contribution of a single pair. This is, again, in agreement with our expectations. Analogous considerations hold for the term on the second line.
In Section 6 we discuss further qualitative features of the result (68), focusing on the dynamics of the entanglement of the semi-infinite chain. Before specialising (68) to that situation, however, we provide some consistency checks of its validity.

Check I: Infinite time limit on a fixed ray
As a first check of (68), we consider the evolution of the entropy when the size of the subsystem A does not scale with t; in this case the entire subsystem is described by a single ray, say ζ 1 . At infinite time the entanglement entropy is expected to coincide with the entropy of the LQSS that describes the steady state at ray ζ 1 , in complete analogy with what happens after homogeneous quenches [35,36,37,38,31]. Let us then verify that (68) is consistent with this picture. Fixing A = [ζ 1 t, ζ 2 t], with in Eq. (68) and taking the infinite time limit with fixed subsystem size , we find which is the expected result. In the second step we used that S Y Y α,λ (ζ) fulfils the continuity equation (20).

Check II: A lead with sub-extensive entropy
Another interesting case is when one of the homogeneous states joined in the bipartite quench protocol, say that on the right, has sub-extensive entropy (S Y Y α,λ (∞) = 0). This setup has been investigated in Ref. [54] for quenches in the XXZ chain.
In this case we expect the entanglement entropy to be independent of the position of the system's right boundary, as long as it cannot be reached by quasiparticles coming from the left lead. In the following we show that this expectation is confirmed by Eq. (68). We first note that, as a special case of (17), the Yang-Yang entropy density satisfies where ρ t α,λ (φ) is the total root density defined in (9). In the specific case considered we have S Y Y α,λ (∞) = 0, so we find Considering an interval A = [ζ 1 t, ζ 2 t] with ζ 2 > v max (cf. Point (i) in Sec. 4), and plugging (73) into the semiclassical expression (68), we find Here we used sgn(J α,−λ (ζ)−ζ) = sgn(ζ −J α,λ (ζ)), which is a consequence of the monotonicity of the functions J α,λ (ζ) and of (64). We now show that the expression in (74) does not depend on ζ 2 as long as ζ 2 > v max . First, from the definition (63) of J α,λ (ζ) follows where we used that for ζ 2 > v max the quasiparticles velocities do not depend on ζ. Second, we observe Here we used that both Z α,λ and Φ α,λ are monotonic functions of their arguments, and can then be applied to both members of a difference in the argument of the sgn function. We also used that Φ α,λ = Z −1 α,λ , and Z α,λ (ζ α,λ ) = 0 (see Section 4). Putting all together we find J α,λ (ζ 2 ) > ζ α,λ . Combining this with the step function θ H (ζ α,λ − ζ 1 ) appearing in the integrand of (74), we then conclude which implies that (74) does not depend on ζ 2 .
6 Entanglement entropy of a semi-infinite interval: Entanglement production rate Here we focus on the entanglement entropy of a semi-infinite interval. There are several reasons why this quantity is interesting. First, it is much easier to calculate via direct numerical methods, such as tDMRG (cf. Sec. 7). Moreover, in contrast with Formula (68), it is much simpler to evaluate. Specifically, (68) requires determining the functions J α,λ (ζ) (see Sec. 4), which relate the trajectories of the quasiparticles forming an entangled pair. As we now show, this is not the case for the entanglement between two semi-infinite chains.
In the first step of (79) we used that, since the particles never cross, the sign of J α,−λ (ζ) − ζ is the sign of the difference of the initial velocities (σ α,λ (ζ) discriminates whether these velocities are computed in the left or in the right state). In the second step we used that the velocities in the left and right macrostates are odd functions of λ. Finally, the last equality in (80) is due to which is a consequence of Assumption 1 and of the continuity in ζ of v α,λ (ζ). We remark that the conditions 1.-3. are not very restrictive. For example, they hold for all the quenches considered so far in the XXZ chain. In particular, they are verified in all the bipartite quenches in the XXZ chain considered in this work.
Case (c) is complementary to (b): both quasiparticles emitted from the right reach the boundary and those emitted from the left are deflected back. The contribution reads as Finally, we point out that, once crossed the junction, a particle obeying the equation of motion (31) cannot come back. The trajectories of the entangled pairs corresponding to the three different scenario are schematically represented in Figure 4. Note that (a) is the standard case: it is realised in quenches from piecewise homogeneous initial states in free systems [57] and in quenches in interacting integrable systems evolving from homogeneous states [31]. Cases (b) and (c), instead, are a landmark of the simultaneous presence of inhomogeneity and interactions. In these cases the entangled pair contributes only for a finite time, i.e., only when the two quasiparticles are in different leads. After both quasiparticles crossed the boundary they do not contribute anymore to the entropy dynamics. The different scenarios (a)-(c) are illustrated in Figure 5 for some concrete bipartite quenches in the XXZ chain.

Entanglement versus thermodynamic entropy production rate
Let us now focus on the the entanglement production rate. This quantity is defined as the slope of the linear growth of the entanglement entropy of a semi-infinite interval, or, in other words, as the time derivative of (78) at fixed ray ζ. As discussed in Section 6, it reads as (cf. (82)) It is particularly instructive to compare it with the rate at which the two semi-infinite intervals exchange thermodynamic entropy. This is defined as First of all, we observe that the exchange rate of thermodynamic entropy is an upper bound for the entanglement production rate, i.e., This is a direct consequence of At the specific ray ζ = 0, after quenches from homogeneous initial states [31] and after bipartite quenches in free models [57], the bound is saturated, S ent,0 = S th,0 . Indeed, in these cases the group velocity does not depend on ζ, implying sgn(v α,λ (σ α,λ (0)∞))v α,λ (0) = sgn(v α,λ (0))v α,λ (0) = |v α,λ (0)| .  Figure 6: Entanglement (full simbols) versus thermodynamic (empty symbols) entropy production rate after the quench from a piecewise homogeneous initial state in the XXZ chain. The x-axis shows that chain anisotropy. Different symbols are used for different initial states obtained by joining the tilted Néel state (24), the dimer state (23), and the tilted ferromagnetic state (25). Full symbols are obtained using Eq. (88) while empty symbols are obtained using Eq. (89).
After bipartite quenches in interacting integrable systems, instead, S ent,0 and S th,0 are generically different. Using the simplifying assumption (83) we see that only scenario (a) in Figure 4 ensures S ent,0 = S th,0 . On the contrary, scenarios (b) and (c) imply S ent,0 = S th,0 . An illustration of the generic behaviour is given in Figure 6, which reports results for three bipartite quenches in the XXZ chain.
We point out that, if one of the scenarios (b) and (c) is allowed, the equality could hold only in the special case when the pairs reaching the junction originate on a lead with subextensive entropy (see Section 5.2). This is the case, for instance, after bipartite quenches in the XXZ chain if one of the two leads is prepared in the ferromagnetic state [58].

tDMRG benchmark
Here we present numerical checks of the results of Sections 5 and 6, focussing on the evolution of the half-chain entanglement entropy after a bipartite quench in the XXZ spin-1/2 chain with ∆ > 1 (cf. (4)).
We first focus on the quench from the state |N, 0 ⊗ |F, θ (cf. (24) and (25)). The results are presented in Figure 7. The different continuous curves are tDMRG results for different values of the tilting angle θ and of the chain anisotropy ∆. In the simulations we considered maximal bond dimensions χ max ≈ 400, which allowed us to reach times t ≈ 25 and half-chain entropy S ≈ 3. The dashed-dotted lines in the figure are the prediction (82). The agreement between the numerics and the analytical result is satisfactory for all quenches.
Let us now consider the difference between the entanglement production rate (88) and the production rate (89) of thermodynamic entropy. As demonstrated in Figure 6, this difference is generically quite small. This makes its numerical observation a highly non-trivial task. In practice, we verified that for all the quenches discussed in Fig. 7, the difference S ent − S th is not visible within the times and system sizes accessible with tDMRG.
To accentuate the discrepancy between S ent and S th , we consider the case that in Fig. 6 shows the largest difference |S ent −S th |, namely, the quench from the state |N ⊗|N, θ (cf. (24)). In this case, however, the entanglement growth is much faster, posing a severe limitation to the timescales accessible by tDMRG. The time evolution of the entanglement entropy after the quench for different values of the tilting angle θ and of the anisotropy ∆ is reported in Figure 8. For short times the tDMRG data exhibit large finite-time effects and are not described by (82). On the other hand, for t 6 the numerical data become compatible with the slope S ent . Still, much larger timescales are needed to provide a robust verification of (82).

Conclusions
We investigated the dynamics of the entanglement entropy after quenches from a piecewise homogeneous initial states in interacting integrable systems. By combining the quasiparticle picture for the entanglement spreading with the GHD approach, we derived an analytic prediction for the entropy evolution after the quench. Remarkably, the entanglement production  88)). The slope of the dotted lines is S th , the exchange rate of thermodynamic entropy at ζ = 0 (cf. (89)). rate, i.e., the growth rate of the entanglement between two half-infinite chains is described by a simple formula that we provided. This depends only on the thermodynamic macrostate (GGE) that describes local properties near the interface between the two chains at infinite time, as it was pointed out in Ref. [54]. We showed, however, that the entanglement production rate is different from the rate of exchange of thermodynamic entropy between the two half-infinite chains. This is in contrast with quenches in free-fermion models [57] and in homogeneous systems [31] and it is a genuine effect of the combination of inhomogeneity and interactions. Our work calls attention to several interesting directions for future research. An immediate one is to provide a more robust independent numerical check, going beyond the tDMRG time scales that we accessed in this work. Moreover, our analytic formula for the entanglement dynamics of a finite interval (cf. (68)) requires the quasiparticle trajectories, which have to be determined numerically. A promising alternative route is to apply the so-called "flea gas" approach [67]. There, the dynamics of out-of-equilibrium quantum systems is simulated by a gas of point-like particles travelling ballistically and scattering elastically.
Another interesting direction is to extend our framework to describe the dynamics of Rényi entropies. Indeed, as recently shown in Ref. [83], from the the dynamics of Rényi entropies one can extract that of the logarithmic negativity [84,85,86,87,88,89,90,91], which is a good entanglement measure for mixed states. Unfortunately, even if the steady-state value of the Rényi entropies is known [92,93,58], computing their full dynamics remains a highly challenging task. A severe complication is that the thermodynamic macrostate describing the Rényi entropies does not coincide with that describing local operators and depends nontrivially on the Rényi index.
Finally, a local breaking of integrability around the interface between the two chains is expected to have dramatic effects on the entanglement production rate along any ray. We wonder, however, whether the ray ζ = 0 is somehow exceptional, displaying a completely different qualitative behaviour. We leave this question to future investigations.