Entanglement negativity in the harmonic chain out of equilibrium

We study the entanglement in a chain of harmonic oscillators driven out of equilibrium by preparing the two sides of the system at different temperatures, and subsequently joining them together. The steady state is constructed explicitly and the logarithmic negativity is calculated between two adjacent segments of the chain. We find that, for low temperatures, the steady-state entanglement is a sum of contributions pertaining to left- and right-moving excitations emitted from the two reservoirs. In turn, the steady-state entanglement is a simple average of the Gibbs-state values and thus its scaling can be obtained from conformal field theory. A similar averaging behaviour is observed during the entire time evolution. As a particular case, we also discuss a local quench where both sides of the chain are initialized in their respective ground states.


Introduction
In the last decade it has been recognized that the entanglement content of manybody quantum states carries essential information about the underlying system [1,2]. In the simplest case, when the system is in a pure state, the entanglement between two complementary parts is measured by the entanglement entropy. One of the most remarkable results for the case of one-dimensional systems is that the entanglement entropy of ground states displays a universal behaviour. At critical points it grows logarithmically in the subsystem size [3], with a prefactor governed by the central charge of the underlying conformal field theory (CFT) [4]; while for noncritical chains it saturates to a finite value, i.e., an area law holds [5].
However, the use of the entropy as an entanglement measure is restricted to pure states and a bipartite setting. The entanglement between non-complementary subsystems, embedded in a larger system, thus needs a different characterization, since the reduced state is in general mixed. Among the numerous proposals to quantify mixed-state entanglement [6], the logarithmic negativity [7,8] turns out to be a particularly useful and easily computable measure. It is especially simple to evaluate for coupled harmonic oscillators [9] and in general for Gaussian states of continuous variable systems [10]. In particular, the logarithmic negativity was used to characterize tripartite entanglement in the ground state of the one-dimensional harmonic chain [11]. Furthermore, results obtained for various quantum spin chains [12,13] indicate universal features in the entanglement of disjoint intervals at criticality. These universal features have recently been understood via CFT methods [14,15]; the corresponding analytical predictions have been compared to numerical simulations in various 1D critical systems [16,17,18] and a very good agreement was found. In two dimensions, entanglement negativity has also been shown to detect topological order [19,20].
In spite of this renewed interest, the behaviour of the entanglement negativity has not yet been investigated out of equilibrium. Here we consider such a problem for a chain of harmonic oscillators that is released from an initial state where the two sides of the chain are kept at different temperatures. The system driven by this thermal gradient evolves, for long times, into a nonequilibrium steady state (NESS) carrying a constant flux of energy. In the gapless limit of the harmonic chain and for sufficiently low reservoir temperatures, the NESS is believed to be described by a nonequilibrium CFT with universal behaviour for e.g. the energy flow [21,22] and thus one expects that universal signatures may also appear in the entanglement negativity.
The choice of our nonequlibrium setup is further motivated by recent studies of a free-fermion chain where the same type of NESS was shown to violate the area law for the mutual information [23]. Similar logarithmic violations were found for the XY spin chain [24] which shows a marked contrast to the thermal-state behaviour where a strict area law can be proven to hold [25]. Since the mutual information measures only the total (classical + quantum) correlations between subsystems, it is natural to ask whether such a singular behaviour would be present for the steady-state entanglement negativity. The NESS of the harmonic chain is an ideal candidate to attack this question since, in contrast to free fermions or spin chains, the logarithmic negativity can easily be extracted from the covariance matrix [9].
Here we focus on the entanglement between adjacent segments of the chain and in a low-temperature limit of the NESS. Our main result is to show that the logarithmic negativity is, to a very good approximation, a sum of two contributions from left-and right-moving normal-mode excitations emitted from the reservoirs. They both carry onehalf of the corresponding thermal-state entanglement that can be found from a simple generalization of the ground-state CFT calculations [14,15]. Hence the steady-state entanglement of the harmonic chain obeys a strict area law.
In the next section we define the model and describe the initial and time evolved states, which is followed by a discussion of the steady state properties for the infinite chain in Sec. 3. Next, we present the covariance matrix formalism in Sec. 4 that is used to obtain the logarithmic negativity. The method is then applied to calculate the steady-state entanglement in Sec. 5 whereas the full time evolution starting from the initial state is presented in Sec. 6. The results are discussed in Sec. 7 and some details of the calculations are presented in two Appendices.

Model and setting
The Hamiltonian of the harmonic chain of length N (with units m = = 1) is given by where x n and p n denote the position and momentum operators of the n-th oscillator with canonical commutation relations [x m , p n ] = iδ m,n and [x m , x n ] = [p m , p n ] = 0. The parameters K and ω 2 0 set the strength of the nearest neighbour coupling and the external harmonic confining potential at each site, respectively. On the chain ends we impose Dirichlet boundary conditions (fixed walls) which imply x 0 = x N +1 = 0 and The initial state is given by a simple product of Gibbs states with inverse temperatures β l and β r for the left and right half-chains, respectively. The Hamiltonian H l (H r ) is defined by a similar expression as in Eq. (1) with the limits of the sum given by 0, N/2 (N/2, N) and Dirichlet boundary conditions imposed at sites 0 and N/2 + 1 (N/2 and N + 1). The initially disconnected halves are joined together at time t = 0 and the state evolves unitarily ρ(t) = e −iHt ρ(0)e iHt under the action of Hamiltonian (1). The change in the geometry of the chain is depicted in Fig. 1. Figure 1. Geometry of the oscillator chain before and after the quench.
Since both the initial state as well as the time evolution operator is Gaussian, ρ(t) can be fully characterized by its 2N × 2N covariance matrix C(t). This can be written in a block-matrix form, composed of symmetric 2 × 2 matrices where R n (t) is a vector of the position and momentum operators in the Heisenberg picture at site n, and O(t) = Tr [ρ(0)O(t)] denotes the average of observable O(t) taken with the initial-state density operator. To obtain the time evolution of R n (t), we first introduce new operators by the canonical transformation and similarly for p k . Note, that φ k (n) are just the eigenvectors of the dynamical matrix in the Hamiltonian, satisfying Dirichlet boundary conditions φ k (0) = φ k (N + 1) = 0. The Hamiltonian is then transformed into with a dispersion relation In the following we will choose K = 1, which sets the maximal velocity of the normalmode excitations to unity. The Heisenberg equations of motion are given byẋ The time-evolution of the covariance matrix then reads where C(0) has a block-diagonal form, composed of N × N covariance matrices of Gibbs states on the two disconnected half-chains. Their matrix elements are given by with α = l, r and 1 ≤ m, n ≤ N/2. Note, that φ α,k and ω 2 α,k are the eigenvectors and eigenvalues of the dynamical matrix of the half-chains and are obtained by exchanging N → N/2 in Eqs. (4) and (6), respectively.

The steady state
Our primary goal is to calculate the entanglement in the steady state, i.e., in the asymptotic limit t → ∞ of the time evolution. For this limit to be well defined, one should set 1 ≪ t ≪ N in order to avoid reflections of the induced wavefront from the fixed ends of the chain. This can be achieved by working directly in the thermodynamic limit N → ∞. The eigenvectors of the dynamical matrix are then Fourier modes, φ q (n) ∼ e −iqn , and the corresponding limit for the elements of the symplectic matrix in Eq. (7) is given by where m, n ∈ Z and the dispersion ω q has the same form as in Eq. (6), but with continuous quasi-momenta q ∈ (−π, π).
The time evolution on such an infinite chain was considered for the case of classical harmonic oscillators before, and the existence of a steady state was shown [26,27]. Note, however, that the time evolution operator S(t) is exactly the same for quantum oscillators and the only difference between the two problems is the form of the initial covariance matrix C(0). Moreover, it was shown in Ref. [27] that, for a large class of initial conditions, the covariance matrix converges locally under the time evolution S(t) in the limit t → ∞. This is true, in particular, if C(0) is translationally invariant asymptotically far away from the cut on both left-and right-hand sides, which is satisfied by the Gibbs-state covariances in Eq. (9) in the limit N → ∞ and |m|, |n| ≫ 1.
The formal derivation of C(t) in the limit t → ∞ was given in Ref. [27]. However, there is a small mistake in the calculation and thus we reiterate the main steps with correct formulas in Appendix A. In turn, the asymptotic limit of the covariance matrix is obtained as with 2 × 2 matrices Note, that the limit in (11) holds by fixed indices m, n and gives a steady-state covariance matrix which is locally translational invariant. In case β l = β r , the steady state supports a nonzero energy current which is encoded in the offdiagonal matrix elements of Eq. (13). The nature of the NESS is however best understood by rather considering the correlation functions of the bosonic annihilation a m and creation a † m operators, defined as the Fourier transforms of the modes which bring the Hamiltonian (5) into diagonal form. The asymptotic form of the bosonic correlation functions is obtained as where the bosonic mode-occupation number is given by The form of n q has a simple physical interpretation. Namely, all the normal modes with positive (negative) group velocities dωq dq originate from the left (right) reservoir and thus their occupation numbers are given by the corresponding Bose-Einstein statistics with inverse temperature β l (β r ). Note, that this behaviour is completely analogous to the one found for free-fermion related models [28,29,30,31,32], and also agrees with recent results obtained within CFT [21,22].

GGE form of the steady state
The steady state given by the bosonic mode-occupation numbers (16) does not correspond to a Gibbs ensemble, unless β l = β r . However, taking into account all the (infinite set of) local conservation laws for the harmonic chain, one can express the NESS density matrix as a generalized Gibbs ensemble (GGE) [27] lim t→∞ where the integrals of motion can be constructed in the N → ∞ limit of a periodic chain and read [27] The corresponding sets of Lagrange multipliers µ ± n can be determined by taking the Fourier transform of H eff , rewriting it in terms of the bosonic operators a q and a † q , and requiring that the resulting expression reproduces the mode-occupations in Eq. (16). This leads to the equations which can be solved and yield the Lagrange multipliers Note, that the only reflection-symmetric conserved charge appearing in the GGE is Q + 0 = H the original Hamiltonian itself and the corresponding multiplier is simply the average temperature. On the other hand, all the Q − n charges have nonvanishing multipliers, decaying slowly µ − n ∼ 1/n for n ≫ 1 and alternating in sign. Therefore, the operator H eff has long-range couplings. Interestingly, the steady state of a freefermion chain, starting from the same initial condition, has a completely analogous GGE description [23,29], the only difference being that multipliers with odd indices vanish identically there, and thus one has no sign alternation in H eff .

Partial transpose and logarithmic negativity
Our aim is to characterize the amount of entanglement in the time-evolved state ρ(t) of the harmonic chain. In the following, we focus on a particular measure of entanglement, the logarithmic negativity [7]. In the most general case, one is interested in a tripartite setting where entanglement is to be measured between subsystems A 1 and A 2 , with A = A 1 ∪ A 2 and B denoting the rest of the chain. The logarithmic negativity is then defined through the partial transpose of the reduced density matrix ρ A (t) = Tr B ρ(t) as where the superscript T 2 indicates a partial transposition with respect to the indices in subsystem A 2 . The logarithmic negativity thus detects only the negative eigenvalues of ρ T 2 A (t). In particular, if there is no entanglement between A 1 and A 2 , then all the eigenvalues of ρ T 2 A (t) are positive [33] and E vanishes due to normalization. Instead of working with density matrices, the logarithmic negativity of the harmonic chain is easier to obtain using the covariance matrix formalism [9]. Indeed, the eigenvalues of ρ A (t) are related to the symplectic eigenvalues of the reduced covariance matrix C A (t). Moreover, the partial transpose can also be implemented on the level of the covariances, with C T 2 A (t) denoting the matrix where the signs of all the momenta p n with n ∈ A 2 are reversed. In turn, the logarithmic negativity is obtained as [9] from the symplectic eigenvalues ν j of C T 2 A (t) and |A| denotes the number of sites in A. Note, that only the eigenvalues ν j < 1 contribute.
The symplectic spectrum of C T 2 A (t) can be obtained through ordinary diagonalization, by multiplying with the symplectic matrix which leads to the spectrum

Steady-state entanglement
We are interested in the entanglement of two neighbouring subsets of oscillators A 1 and A 2 , each of size ℓ. Before presenting results for the logarithmic negativity E, one should point out a subtlety of the numerical treatment. In all the following calculations, we are interested in the gapless limit of the harmonic chain, ω 0 → 0, where the Hamiltonian has an underlying CFT with central charge c = 1. Note, however, that both expressions (12) and (13) diverge for q → 0 in this limit due to the zero-mode. Nevertheless, we observe numerically that E is insensitive to the zero-mode contribution and converges as ω 0 → 0. To ensure that we always remain in the gapless regime, all the calculations are carried out by setting ω 0 = 0.005/ℓ, i.e., reducing the gap with the size of the subsystems.

Equal temperatures
We first consider the simplest case β l = β r = β. The local perturbation in the center, due to the change in the boundary condition, is propagated away and the NESS converges locally to the Gibbs equilibrium state ρ ∼ exp(−βH). The thermal-state entanglement in both spin models [34] and oscillator chains [35,36] has been studied previously with a focus on the critical temperature above which the state becomes separable.
Here, instead, we consider the scaling of the entanglement negativity of two adjacent intervals in the low-temperature regime, β ≫ 1, as a function of ℓ and β. It turns out that CFT methods, applied recently to calculate the ground-state entanglement in the oscillator chain [15], can easily be generalized to finite temperatures in this case. In particular, for two adjacent intervals of lengths ℓ 1 and ℓ 2 embedded in an infinite chain, the ground-state logarithmic negativity can be expressed via three-point functions of twist fields on the complex plane and yields [15] where c is the central charge of the CFT. For finite temperatures, the CFT is defined on an infinite cylinder of circumference β which, however, can be mapped into the plane by a simple conformal transformation. The overall effect of the mapping is to replace the lengths ℓ i → β/π sinh(ℓ i π/β). For intervals of equal lengths, ℓ 1 = ℓ 2 = ℓ, this leads to where we have set c = 1 corresponding to the harmonic chain. The CFT formula (28) has the correct limiting behaviour E ∼ 1/4 ln ℓ for ℓ ≪ β whereas in the opposite limit of large segment sizes, ℓ ≫ β, it predicts the saturation value E ∼ 1/4 ln β. This is indeed what we observe from the numerical data, obtained by the method in Section 4, and shown on the left of Fig. 2. When scaled according to the CFT variable, as shown on the right of Fig. 2, we observe an excellent data collapse and agreement with the prediction of Eq. (28).

Unequal temperatures
Turning to the nonequilibrium case, we now present some simple heuristic arguments how the steady-state entanglement can be related to the equilibrium value in Eq. (28). Following the results of Section 3, the NESS density operator can be written in the form ρ ∼ exp (−β l H + − β r H − ) with H ± describing the evolution of right-and leftmoving particles, respectively. In the CFT context [22], they correspond to mutually commuting chiral components of the Hamiltonian H = H + + H − . In the path-integral representation of ρ, the action thus decouples in the chiral fields φ + and φ − that live on infinite cylinders of circumferences β l and β r , respectively. Due to this separation, the partition functions involved in the calculation of the logarithmic negativity [15] are also supposed to factorise into chiral components and thus their contribution is additive. Finally, making the natural assumption that the entanglement content of the chiral theories is half of that of the full CFT, one expects the relation with the steady-state entanglement E(β l , β r ) being the average of the thermal-state values (28) corresponding to the left and right reservoirs. Our numerical calculations confirm the validity of Eq. (29) to an extremely good precision. The tiny deviations from the equality are supposed to be a consequence of the zero-mode which has been neglected in the above CFT reasoning. In fact, the presence of the zero-mode couples the two chiral branches H ± and thus the factorization of the NESS density matrix is not perfect. However, the effect of this coupling seems to be rather small for the range of temperatures we have considered. This is illustrated in Fig. 3 on the level of the symplectic eigenvalues ν j (β l , β r ) < 1 of the partial transposed covariance matrix which contribute to E(β l , β r ), see Eq. (24). The small deviations from the geometric mean ν j (β l )ν j (β r ) of the Gibbs-state symplectic eigenvalues are shown on the inset. The deviations seem to increase with increasing temperatures and the relation is expected to break down approaching the critical temperature where the entanglement vanishes [35,36].

Time evolution of entanglement
We now consider the complete time evolution of E after the wall separating the two sides of the chain is removed, see Fig. 1. The subsystems A 1 and A 2 are chosen to be adjacent intervals with their common boundary located at the initial cut. Clearly, at  Figure 3. The five smallest symplectic eigenvalues ν j (β l , β r ) < 1 of the partial transposed steady-state covariance matrix for two adjacent intervals of size ℓ = 100 and various pairs of β l and β r . The inset shows the deviation from the geometric mean of Gibbs-state symplectic eigenvalues ν j (β l ) and ν j (β r ). t = 0 one has E = 0 and the entanglement has to increase to reach its steady-state value, discussed in the previous section.

Local quench at zero temperature
The special case, when both sides of the chain are prepared in their respective ground states will be treated first. In fact, this is the same situation, also known as a local quench, which was studied before for free-fermion chains [37,38], as well as in the context of CFT [39,40], in various geometries and with a focus on the entanglement entropy.
In some simple bipartite settings, B = ∅, the result can directly be inferred from previous CFT calculations. The simplest choice is to consider the evolution of E for two halves of an infinite system, A 1 = (−∞, 0] and A 2 = [1, ∞). Since the state ρ(t) is pure, the logarithmic negativity is just the Rényi entropy with index 1/2, and inserting this into the CFT formula of Ref. [39] one obtains Alternatively, one can follow the line of Ref. [15] and work out the CFT representation of the partial-transpose density matrix. The calculation is sketched in Appendix B and leads to the same result. In order to test the result numerically, however, one has to choose a finite system of size N = 2ℓ, with a bipartition A 1 = [1, ℓ] and A 2 = [ℓ + 1, 2ℓ] with the initial cut located between sites ℓ and ℓ + 1. The corresponding result for E can also be found from a CFT calculation based on Ref. [40] and reads E = 1 2 ln |2ℓ/π sin(πt/2ℓ)| + const.
This can now be compared to numerical results obtained with the exact time-dependent covariance matrix, Eqs. (7)(8)(9), following again the steps in section 4, with the result shown in Fig. 4. On the left, the data is plotted for various half-chain lengths ℓ and shows a quasi-periodic structure with period 2ℓ, similar to the evolution of entanglement entropy in the same geometry. This is due to reflections of the front, induced by the quench, from the fixed ends of the chain. Note, however, the slight upwards shift of the curve for ℓ = 25 which is supposedly due to lattice effects. Nevertheless, when plotted against the CFT scaling variable, the data shows a very good collapse and is seen to reproduce Eq. (31) for large arguments, as shown on the right of Fig. 4. Finally, one could consider E for the infinite chain in the tripartite setting of section 5. As discussed in Appendix B, the CFT treatment is rather involved, requiring knowledge of higher order twist-field expectation values. Thus, we resort to the numerical evaluation of E. The ground-state covariance matrix of the semi-infinite chain, with matrix elements given by the N → ∞ limit of Eq. (9), can be evaluated explicitly for ω 0 = 0 and yields [15] C r m, for the right-hand side of the chain, m, n ≥ 1, with ψ(z) being the digamma function. The left-hand side covariance matrix C l m,n for m, n ≤ 0 is obtained by a reflection of the indices m → 1 − m and n → 1 − n in Eq. (32).
The time evolved covariance matrix, Eq. (8), requires to carry out the matrix product with the symplectic evolution matrix S(t) which, in principle, is infinitely large. However, we can make use of the light-cone structure of the matrix elements, implying that for |m−n| ≫ t the entries S m,n (t) are exponentially small. Thus, the sums involved in C T 2 A (t) can be truncated and evaluated numerically to a very high precision. The result for E is shown in Fig. 5. For times t < ℓ, the logarithmic negativity develops a plateau, followed by a sharp drop at t ≈ ℓ, where the propagating front leaves the subsystem A. The data then converges slowly to the ground-state value of E, shown by horizontal lines on the left of Fig. 5. The plateau region is again reminiscent of the behaviour of the entanglement entropy in the corresponding geometry [37,39]. We thus propose the ansatz where the coefficients are fitting parameters. In fact, some of them can be fixed by requiring that in the limit ℓ → ∞ we recover the result (30), implying a = 1/2 and d = −(b + c). For the remaining two we obtain b ≈ 0.15 and c ≈ 0.13 from a fit to the data with ℓ = 100. The data is then scaled together using these values in the right of Fig. 5 and shows a nice collapse.

Finite temperatures
We finally study the case where the initial states on left and right-hand side are prepared at finite temperatures. We consider again the infinite geometry, where the initial covariance matrices are given by Eq. (9) with the sum replaced by an integral. First, we consider the unbiased case β l = β r = β, with results on the time evolution of E shown on the left of Fig. 6. When compared to the local quench results in Fig. 5, one sees that the curves become flatter and eventually saturate in time for increasing temperatures. Nevertheless, one observes the same light-cone effect at t ≈ ℓ and after a sudden decrease E relaxes slowly towards its thermal-state value (28).
The case of unequal temperatures is shown on the right of Fig. 6. The result (29) for the steady state suggests, that the same relation might be true for the time evolution as well. Indeed, the average of the time evolved entanglement with equal temperature initial conditions, E(β l ) and E(β r ), is indistinguishable from the data E(β l , β r ) for unequal temperatures. There are, however, again some small deviations.  Figure 6. Time evolution of entanglement for adjacent intervals of size ℓ = 50 in the infinite geometry. Left: Equal temperatures β l = β r = β. Horizontal lines indicate the steady-state entanglement for β = 5, 10. Right: Entanglement evolution for unequal temperatures (symbols) compared to averages of equal-temperature curves (lines).

Discussion
In conclusion, we have studied the entanglement, measured by the logarithmic negativity, in a simple steady state of the harmonic oscillator chain driven by thermal reservoirs at different temperatures. The steady-state density matrix factorises into two Gibbs-like states, with Hamiltonians given by only the left-or right-moving particles, and the entanglement is found to be the sum of the chiral contributions. These are simply equal to half of the thermal-state entanglement corresponding to each reservoir, which can be found by CFT calculations. The above additivity property depends crucially on the assumption that the effect of the zero-mode, which couples the chiral branches, can be neglected. This seems not to be valid for the steady state of a free-fermion chain, prepared with the same initial condition as considered here. Indeed, in the latter case the mode-occupation n q , given by the Fermi-statistics analogue of Eq. (16), develops a jump singularity around q = 0. This leads to a contribution in the mutual information of two adjacent segments, which scales logarithmically in the segment size [23,41], and thus the additivity of the chiral contributions does not hold for this measure. Note, however, that in the CFT limit β l , β r ≫ 1, the prefactor of the logarithm is exponentially suppressed as a function of the inverse temperatures, and numerical results indicate that the additivity is practically recovered even for very large subsystem sizes.
In the opposite limit of high temperatures, the additivity property for the logarithmic negativity is weakly violated for the harmonic chain, however, the area law is still strictly obeyed. The question thus remains open, whether the logarithmic violation of the area law found for the mutual information in the fermionic NESS could persist for the entanglement negativity. As a further extension, one could test the additivity property of the entanglement for interacting models, such as the NESS of the XXZ chain [42] or of special integrable quantum field theories [43].
where the + (−) sign in the limit corresponds to the right (left) hand side with α = r (α = l). The Fourier transform of σ α m,n is denoted by σ α q and we introduce the notation σ ± q = (σ r q ± σ l q )/2. Now we split up the initial covariance matrix in three terms as C m,n (0) = C and C 3 m,n (0) describes the remaining terms. Note, that the sum of the first two terms gives the correct asymptotic behaviour in (A.1), however, they generate nonzero matrix elements in the offdiagonals of C(0) in Eq. (8), which have to be compensated by C 3 m,n (0). The time evolved matrices are given by C i (t) = S(t)C i (0)S(t) T such that First, we consider C 1 (t) which is a product of three Toeplitz matrices and thus the multiplication can be performed in Fourier space C 1 q (t) = S q (t)σ + q S T q (t). The symbol of the matrix in Eq. (10) can be rewritten as a sum of diagonal and offdiagonal contributions S q (t) = cos(ω q t) 1l + sin(ω q t) Γ q , Multiplying out C 1 q (t) and taking t → ∞, the rapidly oscillating terms can be substituted by their average and we arrive at which gives the first term in the integral of Eq. (11). The next step is to obtain the asymptotics of C 2 (t). Working again in Fourier space, one has now to include the Fourier transform of the sign function in C 2 m,n (0) which leads to where the P sign indicates that the Cauchy principal value has to be taken. It can be evaluated using the identity where f and g denote some smooth functions. In turn, the principal value integral has the effect of interchanging the oscillatory terms cos(ω p t) → − sin(ω q t) and sin(ω p t) → cos(ω q t) in S T p (t). In the remaining integral over q, one can again take the averages of the time dependent terms which yields Noting that in our case sgn(ω ′ q ) = sgn(q) and calculating the matrix products leads to the second term in Eq. (11).
We finally show that lim t→∞ C 3 (t) = 0. Using the form of C(0) in Eqs. (8) and (9) as well as the expression for the eigenvectors in (4), one has in the thermodynamic limit (N → ∞) the block-matrix form Theσ α in the diagonal are Hankel matrices (with symbol σ α q ) that arise due to the Dirichlet boundary condition in the center, while the offdiagonal Toeplitz matrices σ α compensate the extra contributions from C 1 (0) + C 2 (0), as mentioned before. The time-evolved matrix C 3 (t) has thus four different contributions from the various Hankel/Toeplitz blocks which can be written as with the definition and F l ± (q, p, p ′ ) = F r ± (q, p ′ , p) * . Note, that F r ± is just the Fourier transform of the product of step functions [1 + sgn(m)] /2 × [1 ± sgn(n)] /2 which singles out the lower right/left block. Carrying out the integrals over p and p ′ in Eq. (A.9) using formula (A.6), one arrives at lim t→∞ P π −π dp 2π π −π dp ′ 2π e i(pm−p ′ n) F r ± (q, p, p ′ )S p (t)σ r q S T p ′ (t) = where we definedŜ q (t) = − sin(ω q t)1l + cos(ω q t)Γ q and used the symmetry properties S −q (t) = S q (t),Ŝ −q (t) =Ŝ q (t) and sgn(ω ′ −q ) = −sgn(ω ′ q ). Calculating the matrix products, one finds that the expression in the brackets vanishes identically. The same holds for the integrals with F l ± which concludes our proof.

Appendix B. CFT results for local quench
Here we briefly summarize the method used to obtain E for a local quench. For a detailed discussion of ground-state entanglement negativity within CFT, we refer to Ref. [15]. The essential step is to rewrite Eq. (23) as where n e is an even integer which, at the end of the calculation, has to be analytically continued to one. Thus, one first needs to express the trace of an even power of the partial-transposed reduced density matrix ρ A (t) = Tr B ρ(t) in terms of a path integral. This is done by considering an n e -sheeted Riemann surface and sewing together the replicas of ρ T 2 A (t). Each copy can be represented by a 2D path integral with slits along the intervals A 1 and A 2 on the real axis, and partial transposition T 2 corresponds to interchanging the edges of the corresponding slit A 2 .
Instead of working on a replicated world-sheet, one can introduce replicated fields and the so-called twist fields, T ne and T ne , that cyclically permute replicas in one of two directions. Each time when replicas are sewn together, one has to calculate expectation values of products of the two twist fields, inserted at the endpoints of the slits. Whenever the slit corresponds to the partial transposed subsystem, the twist operators have to be interchanged. The simplest choice for the subsystems is a bipartition (B = ∅) into two semiinfinite lines A 1 and A 2 . One has then a single contact point and thus Tr(ρ T 2 A (t)) ne = T 2 ne (z 0 ) (B.2) where the expectation value of the twist fields has to be calculated on a world-sheet representing the time-evolved density operator ρ(t). This is depicted on the left of Fig. B1, where the cut in the middle corresponds to the imaginary time evolution of two decoupled CFTs with fixed boundary conditions, yielding the initial state ρ(0). The real time evolution takes place between the endpoints z = ±iǫ of the cut, and the twist field has to be inserted at z 0 = iτ . Note, that the parameter ǫ is needed to regularize the path-integral and the analytical continuation to real times τ = it must be carried out at the end of the calculation. Although the original world-sheet has a complicated geometry, one can apply a conformal mapping [39] which transforms it to the half-plane, as shown on the right of Fig. B1. On the halfplane the expectation value of a one-point function is known and one can then use the conformal transformation formula to obtain where w 0 = w(z 0 ) and the scaling dimension of the operator T 2 ne is given by [15] ∆ ne = c 6 Continuing the result to real time τ → it, taking the limit t ≫ ǫ and finally using Eq. (B.1), one arrives to the result in Eq. (30).
The local quench for the finite system can be treated in a similar way. There the world-sheet has a double-pants geometry, with fixed boundary conditions along Re(z) = ±ℓ, and the proper conformal mapping to the half-plane was given in Ref. [40]. Carrying out the analogous steps, one arrives to the formula in Eq. (31).
On the other hand, the tripartite case in the infinite system with line segments A 1 = [−ℓ, 0] and A 2 = [0, ℓ] is more involved. There, one has to consider the three-point function T (z −1 )T 2 ne (z 0 )T (z 1 ) of the twist fields with z 0 = iτ and z ±1 = ±ℓ + iτ . This is mapped under (B.3) into a three-point function on the half-plane which, however, has the complexity of a six-point function on the full plane. Although some recent progress has been made in the derivation of higher order twist-field correlators in Ref. [44], their structure is rather involved and we have not been able to tackle this case analytically.