Entanglement negativity in the harmonic chain out of equilibrium

In the last decade it has been recognized that the entanglement content of many-body 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 (1D) 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 1D 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–18] and a very good agreement was found. In two dimensions (2D), 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 one-half 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 section 3. Next, we present the covariance matrix formalism in section 4 that is used to obtain the logarithmic negativity. The method is then applied to calculate the steady-state entanglement in section 5 whereas the full time evolution starting from the initial state is presented in section 6. The results are discussed in section 7 and some details of the calculations are presented in two appendices.

The Hamiltonian of the harmonic chain of length N (with units $m=\hbar =1$) is given by

$$\begin{eqnarray}&&H=\frac{1}{2}\mathop{\sum }\limits_{n=0}^{N}\left[ p_{n}^{2}+K{{({{x}_{n+1}}-{{x}_{n}})}^{2}}+\omega _{0}^{2}x_{n}^{2} \right],\end{eqnarray} \tag{ 1 }$$

where x_n and p_n denote the position and momentum operators of the nth oscillator with canonical commutation relations $[{{x}_{m}},{{p}_{n}}]={\rm i}{{\delta }_{m,n}}$ and $[{{x}_{m}},{{x}_{n}}]=[{{p}_{m}},{{p}_{n}}]=0$. The parameters K and $\omega _{0}^{2}$ 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 ${{p}_{0}}={{p}_{N+1}}=0$.

The initial state is given by a simple product of Gibbs states

$$\begin{eqnarray}&&\rho (0)=\frac{1}{{{Z}_{{\rm l}}}}{{{\rm e}}^{-{{\beta }_{{\rm l}}}{{H}_{{\rm l}}}}}\otimes \frac{1}{{{Z}_{{\rm r}}}}{{{\rm e}}^{-{{\beta }_{{\rm r}}}{{H}_{{\rm r}}}}},\end{eqnarray} \tag{ 2 }$$

with inverse temperatures ${{\beta }_{{\rm l}}}$ and ${{\beta }_{{\rm r}}}$ for the left and right half-chains, respectively. The Hamiltonian H_l (H_r) is defined by a similar expression as in equation (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 $\rho (t)={{{\rm e}}^{-{\rm i}Ht}}\rho (0){{{\rm e}}^{{\rm i}Ht}}$ under the action of Hamiltonian (1). The change in the geometry of the chain is depicted in figure 1.

Zoom In Zoom Out Reset image size

Since both the initial state as well as the time evolution operator is Gaussian, $\rho (t)$ can be fully characterized by its $2N\times 2N$ covariance matrix $C(t)$. This can be written in a block-matrix form, composed of symmetric 2 × 2 matrices

$$\begin{eqnarray}&&{{C}_{m,n}}(t)=\langle {{R}_{m}}(t)R_{n}^{T}(t)+{{R}_{n}}(t)R_{m}^{T}(t)\rangle ,\quad {{R}_{n}}(t)=\left( \begin{array}{ccccccccccccccc} {{x}_{n}}(t) \\ {{p}_{n}}(t) \\ \end{array} \right),\end{eqnarray} \tag{ 3 }$$

where R_n(t) is a vector of the position and momentum operators in the Heisenberg picture at site n, and $\langle \mathcal{O}(t)\rangle ={\rm Tr}\left[ \rho (0)\mathcal{O}(t) \right]$ denotes the average of observable $\mathcal{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

$$\begin{eqnarray}&&{{x}_{k}}=\mathop{\sum }\limits_{n=1}^{N}{{\phi }_{k}}(n){{x}_{n}},\quad {{\phi }_{k}}(n)=\sqrt{\frac{2}{N+1}}{\rm sin} \frac{\pi kn}{N+1},\quad k=1,...,N,\end{eqnarray} \tag{ 4 }$$

and similarly for p_k. Note, that ${{\phi }_{k}}(n)$ are just the eigenvectors of the dynamical matrix in the Hamiltonian, satisfying Dirichlet boundary conditions ${{\phi }_{k}}(0)={{\phi }_{k}}(N+1)=0$. The Hamiltonian is then transformed into

$$\begin{eqnarray}&&H=\frac{1}{2}\mathop{\sum }\limits_{k=1}^{N}\left( p_{k}^{2}+\omega _{k}^{2}x_{k}^{2} \right),\end{eqnarray} \tag{ 5 }$$

with a dispersion relation

$$\begin{eqnarray}&&{{\omega }_{k}}=\sqrt{\omega _{0}^{2}+4\;K{{{\rm sin} }^{2}}({{q}_{k}}/2)},\qquad {{q}_{k}}=\frac{\pi k}{N+1}.\end{eqnarray} \tag{ 6 }$$

In the following we will choose K = 1, which sets the maximal velocity of the normal-mode excitations to unity.

The Heisenberg equations of motion are given by ${{\dot{x}}_{k}}(t)={{p}_{k}}(t)$ and ${{\dot{p}}_{k}}(t)=-\omega _{k}^{2}{{x}_{k}}(t)$, which yields $R(t)=S(t)R(0)$ with a symplectic matrix S(t) of block form

$$\begin{eqnarray}{{S}_{m,n}}(t)=\mathop{\sum }\limits_{k=1}^{N}{{\phi }_{k}}*(m){{\phi }_{k}}(n)\left( \begin{array}{ccccccccccccccc} {\rm cos} {{\omega }_{k}}t & \omega _{k}^{-1}{\rm sin} {{\omega }_{k}}t \\ -{{\omega }_{k}}{\rm sin} {{\omega }_{k}}t & {\rm cos} {{\omega }_{k}}t \\ \end{array} \right).\end{eqnarray} \tag{ 7 }$$

The time-evolution of the covariance matrix then reads

$$\begin{eqnarray}C(t)=S(t)C(0)S{{(t)}^{T}},\qquad C(0)=\left( \begin{array}{ccccccccccccccc} {{C}^{l}} & 0 \\ 0 & {{C}^{r}} \\ \end{array} \right),\end{eqnarray} \tag{ 8 }$$

where $C(0)$ has a block-diagonal form, composed of $N\times N$ covariance matrices of Gibbs states on the two disconnected half-chains. Their matrix elements are given by

$$\begin{eqnarray}C_{m,n}^{\alpha }=\mathop{\sum }\limits_{k=1}^{N/2}{{\phi }_{\alpha ,k}}*(m){{\phi }_{\alpha ,k}}(n)\left( \begin{array}{ccccccccccccccc} \omega _{\alpha ,k}^{-1} & 0 \\ 0 & {{\omega }_{\alpha ,k}} \\ \end{array} \right){\rm coth} \frac{{{\beta }_{\alpha }}{{\omega }_{\alpha ,k}}}{2},\end{eqnarray} \tag{ 9 }$$

with $\alpha =l,r$ and $1\leqslant m,n\leqslant N/2$. Note, that ${{\phi }_{\alpha ,k}}$ and $\omega _{\alpha ,k}^{2}$ are the eigenvectors and eigenvalues of the dynamical matrix of the half-chains and are obtained by exchanging $N\to N/2$ in equations (4) and (6), respectively.

Our primary goal is to calculate the entanglement in the steady state, i.e., in the asymptotic limit $t\to \infty$ of the time evolution. For this limit to be well defined, one should set $1\ll t\ll 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\to \infty$. The eigenvectors of the dynamical matrix are then Fourier modes, ${{\phi }_{q}}(n)\sim {{{\rm e}}^{-{\rm i}qn}}$, and the corresponding limit for the elements of the symplectic matrix in equation (7) is given by

$$\begin{eqnarray}{{S}_{m,n}}(t)=\int _{-\pi }^{\pi }\frac{{\rm d}q}{2\pi }{{{\rm e}}^{{\rm i}q(m-n)}}\left( \begin{array}{ccccccccccccccc} {\rm cos} {{\omega }_{q}}t & \omega _{q}^{-1}{\rm sin} {{\omega }_{q}}t \\ -{{\omega }_{q}}{\rm sin} {{\omega }_{q}}t & {\rm cos} {{\omega }_{q}}t \\ \end{array} \right),\end{eqnarray} \tag{ 10 }$$

where $m,n\in \mathbb{Z}$ and the dispersion ${{\omega }_{q}}$ has the same form as in equation (6), but with continuous quasi-momenta $q\in (-\pi ,\pi )$.

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 [27] that, for a large class of initial conditions, the covariance matrix converges locally under the time evolution S(t) in the limit $t\to \infty$. 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 equation (9) in the limit $N\to \infty$ and $|m|,|n|\gg 1$.

The formal derivation of C(t) in the limit $t\to \infty$ was given in [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

$$\begin{eqnarray}&&\mathop{{\rm lim} }\limits_{t\to \infty }{{C}_{m,n}}(t)=\int _{-\pi }^{\pi }\frac{{\rm d}q}{2\pi }{{{\rm e}}^{{\rm i}\,q(m-n)}}\left( C_{q}^{+}+{\rm i}\;{\rm sgn} (q)\;C_{q}^{-} \right),\end{eqnarray} \tag{ 11 }$$

with 2 × 2 matrices

$$\begin{eqnarray}C_{q}^{+}=\frac{1}{2}\left[ {\rm coth} \left( \frac{{{\beta }_{{\rm r}}}{{\omega }_{q}}}{2} \right)+{\rm coth} \left( \frac{{{\beta }_{{\rm l}}}{{\omega }_{q}}}{2} \right) \right]\left( \begin{array}{ccccccccccccccc} \omega _{q}^{-1} & 0 \\ 0 & {{\omega }_{q}} \\ \end{array} \right),\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}C_{q}^{-}=\frac{1}{2}\left[ {\rm coth} \left( \frac{{{\beta }_{{\rm r}}}{{\omega }_{q}}}{2} \right)-{\rm coth} \left( \frac{{{\beta }_{{\rm l}}}{{\omega }_{q}}}{2} \right) \right]\left( \begin{array}{ccccccccccccccc} 0 & -1 \\ 1 & 0 \\ \end{array} \right).\end{eqnarray} \tag{ 13 }$$

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 ${{\beta }_{{\rm l}}}\ne {{\beta }_{{\rm r}}}$, the steady state supports a nonzero energy current which is encoded in the offdiagonal matrix elements of equation (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}^{\dagger }$ operators, defined as the Fourier transforms of the modes

$$\begin{eqnarray}&&{{a}_{q}}=\sqrt{\frac{{{\omega }_{q}}}{2}}\left( {{x}_{q}}+\frac{{\rm i}}{{{\omega }_{q}}}{{p}_{q}} \right),\qquad a_{q}^{\dagger }=\sqrt{\frac{{{\omega }_{q}}}{2}}\left( {{x}_{-q}}-\frac{{\rm i}}{{{\omega }_{q}}}{{p}_{-q}} \right),\end{eqnarray} \tag{ 14 }$$

which bring the Hamiltonian (5) into diagonal form. The asymptotic form of the bosonic correlation functions is obtained as

$$\begin{eqnarray}&&\mathop{{\rm lim} }\limits_{t\to \infty }{\rm Tr}(\rho (t)\;a_{m}^{\dagger }{{a}_{n}})=\int _{-\pi }^{\pi }\frac{{\rm d}q}{2\pi }{{{\rm e}}^{{\rm i}q(m-n)}}\;{{n}_{q}},\end{eqnarray} \tag{ 15 }$$

where the bosonic mode-occupation number is given by

$$\begin{eqnarray}{{n}_{q}}=\left\{ \begin{array}{ccccccccccccccc} \frac{1}{{{{\rm e}}^{{{\beta }_{{\rm r}}}{{\omega }_{q}}}}-1} & q\in \left( -\pi ,0 \right), \\ \frac{1}{{{{\rm e}}^{{{\beta }_{{\rm l}}}{{\omega }_{q}}}}-1} & q\in \left( 0,\pi \right). \\ \end{array} \right.\end{eqnarray} \tag{ 16 }$$

The form of n_q has a simple physical interpretation. Namely, all the normal modes with positive (negative) group velocities $\frac{{\rm d}{{\omega }_{q}}}{{\rm d}q}$ originate from the left (right) reservoir and thus their occupation numbers are given by the corresponding Bose–Einstein statistics with inverse temperature ${{\beta }_{{\rm l}}}$ (${{\beta }_{{\rm r}}}$). Note, that this behaviour is completely analogous to the one found for free-fermion related models [28–33], and also agrees with recent results obtained within CFT [21, 22].

3.1. Generalized Gibbs ensemble (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 ${{\beta }_{{\rm l}}}={{\beta }_{{\rm 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 GGE [27]

$$\begin{eqnarray}&&\mathop{{\rm lim} }\limits_{t\to \infty }\rho (t)=\frac{1}{Z}{{{\rm e}}^{-{{H}_{{\rm eff}}}}},\qquad {{H}_{{\rm eff}}}=\mathop{\sum }\limits_{n=0}^{\infty }(\mu _{n}^{+}Q_{n}^{+}+\mu _{n}^{-}Q_{n}^{-}),\end{eqnarray} \tag{ 17 }$$

where the integrals of motion can be constructed in the $N\to \infty$ limit of a periodic chain and read [27]

$$\begin{eqnarray}&&Q_{n}^{+}=\frac{1}{2}\mathop{\sum }\limits_{m=-\infty }^{\infty }\left[ {{p}_{m}}{{p}_{m+n}}+({{x}_{m}}-{{x}_{m+1}})({{x}_{m+n}}-{{x}_{m+n+1}})+\omega _{0}^{2}{{x}_{m}}{{x}_{m+n}} \right],\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&Q_{n}^{-}=\frac{1}{4}\mathop{\sum }\limits_{m=-\infty }^{\infty }({{p}_{m}}{{x}_{m+n}}-{{x}_{m}}{{p}_{m+n}}-{{p}_{m+n}}{{x}_{m}}+{{x}_{m+n}}{{p}_{m}}).\end{eqnarray} \tag{ 19 }$$

The corresponding sets of Lagrange multipliers $\mu _{n}^{\pm }$ can be determined by taking the Fourier transform of ${{H}_{{\rm eff}}}$, rewriting it in terms of the bosonic operators a_q and $a_{q}^{\dagger }$, and requiring that the resulting expression reproduces the mode-occupations in equation (16). This leads to the equations

$$\begin{eqnarray}&&\mathop{\sum }\limits_{n=0}^{\infty }\left( \mu _{n}^{+}{{\omega }_{q}}{\rm cos} nq+\mu _{n}^{-}{\rm sin} nq \right)={{\omega }_{q}}{{\beta }_{l}},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&\mathop{\sum }\limits_{n=0}^{\infty }\left( \mu _{n}^{+}{{\omega }_{q}}{\rm cos} nq-\mu _{n}^{-}{\rm sin} nq \right)={{\omega }_{q}}{{\beta }_{r}},\end{eqnarray} \tag{ 21 }$$

which can be solved and yield the Lagrange multipliers

$$\begin{eqnarray}&&\mu _{n}^{+}=\frac{{{\beta }_{{\rm l}}}+{{\beta }_{{\rm r}}}}{2}{{\delta }_{n,0}},\qquad \mu _{n}^{-}=({{\beta }_{{\rm l}}}-{{\beta }_{{\rm r}}}){{(-1)}^{n}}\frac{4}{\pi }\frac{n}{4{{n}^{2}}-1}.\end{eqnarray} \tag{ 22 }$$

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 inverse temperature. On the other hand, all the $Q_{n}^{-}$ charges have nonvanishing multipliers, decaying slowly $\mu _{n}^{-}\sim 1/n$ for $n\gg 1$ and alternating in sign. Therefore, the operator ${{H}_{{\rm eff}}}$ has long-range couplings. Interestingly, the steady state of a free-fermion 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}_{{\rm eff}}}.$

To conclude this section, we remark that the non-local structure of the GGE is a direct consequence of the thermodynamic limit. Indeed, if boundary conditions are retained at the ends of the chain, one expects that the currents are suppressed for large times due to reflections, $\mu _{n}^{-}=0$ for all n, and thus the GGE becomes local. This was recently proven for a free-fermion field theory [32] and we believe that the same holds also in the continuum limit of the oscillator chain.

Our aim is to characterize the amount of entanglement in the time-evolved state $\rho (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₁ and A₂, with $A={{A}_{1}}\cup {{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 ${{\rho }_{A}}(t)={\rm T}{{{\rm r}}_{B}}\;\rho (t)$ as

$$\begin{eqnarray}&&\mathcal{E}={\rm lnTr}\left| \rho _{A}^{{{T}_{2}}}(t) \right|,\end{eqnarray} \tag{ 23 }$$

where the superscript T₂ indicates a partial transposition with respect to the indices in subsystem A₂. The logarithmic negativity thus detects only the negative eigenvalues of $\rho _{A}^{{{T}_{2}}}(t)$. In particular, if there is no entanglement between A₁ and A₂, then all the eigenvalues of $\rho _{A}^{{{T}_{2}}}(t)$ are positive [34] and $\mathcal{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 ${{\rho }_{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_{A}^{{{T}_{2}}}(t)$ denoting the matrix where the signs of all the momenta p_n with $n\in {{A}_{2}}$ are reversed. In turn, the logarithmic negativity is obtained as [9]

$$\begin{eqnarray}&&\mathcal{E}=-\mathop{\sum }\limits_{j=1}^{|A|}{\rm ln} \ {\rm min} ({{\nu }_{j}},1)\end{eqnarray} \tag{ 24 }$$

from the symplectic eigenvalues ${{\nu }_{j}}$ of $C_{A}^{{{T}_{2}}}(t)$ and $|A|$ denotes the number of sites in A. Note, that only the eigenvalues ${{\nu }_{j}}\lt 1$ contribute.

The symplectic spectrum of $C_{A}^{{{T}_{2}}}(t)$ can be obtained through ordinary diagonalization, by multiplying with the symplectic matrix

$$\begin{eqnarray}{{\Sigma }_{A}}=\mathop{\oplus }\limits_{j\in A}\left( \begin{array}{ccccccccccccccc} 0 & 1 \\ -1 & 0 \\ \end{array} \right),\end{eqnarray} \tag{ 25 }$$

which leads to the spectrum

$$\begin{eqnarray}&&{\rm Spect}\left( {{\Sigma }_{A}}C_{A}^{{{T}_{2}}}(t) \right)=\left\{ \pm {\rm i}{{\nu }_{1}},...,\pm {\rm i}{{\nu }_{|A|}} \right\}.\end{eqnarray} \tag{ 26 }$$

We are interested in the entanglement of two neighbouring subsets of oscillators A₁ and A₂, each of size ℓ. Before presenting results for the logarithmic negativity $\mathcal{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, ${{\omega }_{0}}\to 0$, where the Hamiltonian has an underlying CFT with central charge c = 1. Note, however, that both expressions (12) and (13) diverge for $q\to 0$ in this limit due to the zero-mode. Nevertheless, we observe numerically that $\mathcal{E}$ is insensitive to the zero-mode contribution and converges as ${{\omega }_{0}}\to 0$. Interestingly, this is in contrast with the behaviour of the entropy and mutual information, which both involve a divergent zero-mode contribution. To ensure that we always remain in the gapless regime, all the calculations are carried out by setting ${{\omega }_{0}}=0.005/\ell$, i.e., reducing the gap with the size of the subsystems.

5.1. Equal temperatures

We first consider the simplest case ${{\beta }_{{\rm l}}}={{\beta }_{{\rm r}}}=\beta$. 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 $\rho \sim {\rm exp} (-\beta H)$. The thermal-state entanglement in both spin models [35] and oscillator chains [36–39] 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, $\beta \gg 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 ℓ₁ and ℓ₂ 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]

$$\begin{eqnarray}&&\mathcal{E}=\frac{c}{4}{\rm ln} \frac{{{\ell }_{1}}{{\ell }_{2}}}{{{\ell }_{1}}+{{\ell }_{2}}}+{\rm constant},\end{eqnarray} \tag{ 27 }$$

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 ${{\ell }_{i}}\to \beta /\pi {\rm sinh} ({{\ell }_{i}}\pi /\beta )$. For intervals of equal lengths, ${{\ell }_{1}}={{\ell }_{2}}=\ell$, this leads to

$$\begin{eqnarray}&&\mathcal{E}=\frac{1}{4}{\rm ln} \frac{\beta }{\pi }{\rm tanh} \frac{\ell \pi }{\beta }+{\rm constant},\end{eqnarray} \tag{ 28 }$$

where we have set c = 1 corresponding to the harmonic chain.

The CFT formula (28) has the correct limiting behaviour $\mathcal{E}\sim 1/4{\rm ln} \ell$ for $\ell \ll \beta$ whereas in the opposite limit of large segment sizes, $\ell \gg \beta$, it predicts the saturation value $\mathcal{E}\sim 1/4{\rm ln} \beta$. This is indeed what we observe from the numerical data, obtained by the method in section 4, and shown on the left of figure 2. When scaled according to the CFT variable, as shown on the right of figure 2, we observe an excellent data collapse and agreement with the prediction of equation (28).

Zoom In Zoom Out Reset image size

5.2. 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 equation (28). Following the results of section 3, the NESS density operator can be written in the form $\rho \sim {\rm exp} (-{{\beta }_{{\rm l}}}{{H}_{+}}-{{\beta }_{{\rm r}}}{{H}_{-}})$ with ${{H}_{\pm }}$ describing the evolution of right- and left-moving 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 ${{\phi }_{+}}$ and ${{\phi }_{-}}$ that live on infinite cylinders of circumferences ${{\beta }_{{\rm l}}}$ and ${{\beta }_{{\rm 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

$$\begin{eqnarray}&&\mathcal{E}({{\beta }_{{\rm l}}},{{\beta }_{{\rm r}}})=\frac{\mathcal{E}({{\beta }_{{\rm l}}})+\mathcal{E}({{\beta }_{{\rm r}}})}{2},\end{eqnarray} \tag{ 29 }$$

with the steady-state entanglement $\mathcal{E}({{\beta }_{{\rm l}}},{{\beta }_{{\rm r}}})$ being the average of the thermal-state values (28) corresponding to the left and right reservoirs.

Our numerical calculations confirm the validity of equation (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}_{\pm }}$ 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 figure 3 on the level of the symplectic eigenvalues ${{\nu }_{j}}({{\beta }_{{\rm l}}},{{\beta }_{{\rm r}}})\lt 1$ of the partial transposed covariance matrix which contribute to $\mathcal{E}({{\beta }_{{\rm l}}},{{\beta }_{{\rm r}}})$, see equation (24). The small deviations from the geometric mean $\sqrt{{{\nu }_{j}}({{\beta }_{{\rm l}}}){{\nu }_{j}}({{\beta }_{{\rm 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 [36, 39].

Zoom In Zoom Out Reset image size

We now consider the complete time evolution of $\mathcal{E}$ after the wall separating the two sides of the chain is removed, see figure 1. The subsystems A₁ and A₂ are chosen to be adjacent intervals with their common boundary located at the initial cut. Clearly, at t = 0 one has $\mathcal{E}=0$ and the entanglement has to increase to reach its steady-state value, discussed in the previous section.

6.1. 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 [40, 41], as well as in the context of CFT [42, 43], in various geometries and with a focus on the entanglement entropy.

In some simple bipartite settings, $B=\varnothing$, the result can directly be inferred from previous CFT calculations. The simplest choice is to consider the evolution of $\mathcal{E}$ for two halves of an infinite system, ${{A}_{1}}=\left( -\infty ,0 \right]$ and ${{A}_{2}}=\left[ 1,\infty \right)$. Since the state $\rho (t)$ is pure, the logarithmic negativity is just the Rényi entropy with index $1/2$, and inserting this into the CFT formula of [42] one obtains

$$\begin{eqnarray}&&\mathcal{E}=\frac{1}{2}{\rm ln} t+{\rm constant}.\end{eqnarray} \tag{ 30 }$$

Alternatively, one can follow the line of [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\ell$, with a bipartition ${{A}_{1}}=\left[ 1,\ell \right]$ and ${{A}_{2}}=\left[ \ell +1,2\ell \right]$ with the initial cut located between sites ℓ and $\ell +1$. The corresponding result for $\mathcal{E}$ can also be found from a CFT calculation based on [43] and reads

$$\begin{eqnarray}&&\mathcal{E}=\frac{1}{2}{\rm ln} \left| 2\ell /\pi {\rm sin} (\pi t/2\ell ) \right|+{\rm constant}.\end{eqnarray} \tag{ 31 }$$

This can now be compared to numerical results obtained with the exact time-dependent covariance matrix, equations (7–9), following again the steps in section 4, with the result shown in figure 4. On the left, the data is plotted for various half-chain lengths ℓ and shows a quasi-periodic structure with period $2\ell$, 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 $\ell =25$ which is supposedly due to lattice effects, caused by the slower normal-mode excitations with velocity $\frac{{\rm d}{{\omega }_{q}}}{{\rm d}q}\lt 1$. Furthermore, there might be universal subleading contributions, originating from a more careful CFT treatment and breaking the periodicity [44], that are, however, difficult to identify in the numerics. Nevertheless, when plotted against the CFT scaling variable in equation (31), the data shows a very good collapse and is seen to reproduce the formula for large arguments, as shown on the right of figure 4.

Zoom In Zoom Out Reset image size

Finally, one could consider $\mathcal{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 $\mathcal{E}$. The ground-state covariance matrix of the semi-infinite chain, with matrix elements given by the $N\to \infty$ limit of equation (9), can be evaluated explicitly for ${{\omega }_{0}}=0$ and yields [15]

$$\begin{eqnarray}C_{m,n}^{{\rm r}}={{C}^{{\rm r}}}(m+n)-{{C}^{{\rm r}}}(m-n),\quad {{C}^{{\rm r}}}(x)=\frac{1}{\pi }\left( \begin{array}{ccccccccccccccc} \psi (1/2+x) & 0 \\ 0 & \frac{4}{4{{x}^{2}}-1} \\ \end{array} \right)\end{eqnarray} \tag{ 32 }$$

for the right-hand side of the chain, $m,n\geqslant 1$, with $\psi (z)$ being the digamma function. The left-hand side covariance matrix $C_{m,n}^{{\rm l}}$ for $m,n\leqslant 0$ is obtained by a reflection of the indices $m\to 1-m$ and $n\to 1-n$ in equation (32).

The time evolved covariance matrix, equation (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|\gg t$ the entries ${{S}_{m,n}}(t)$ are exponentially small. Thus, the sums involved in $C_{A}^{{{T}_{2}}}(t)$ can be truncated and evaluated numerically to a very high precision.

The result for $\mathcal{E}$ is shown in figure 5. For times $t\lt \ell$, the logarithmic negativity develops a plateau, followed by a sharp drop at $t\approx \ell$, where the propagating front leaves the subsystem A. The data then converges slowly to the ground-state value of $\mathcal{E}$, shown by horizontal lines on the left of figure 5. The plateau region is again reminiscent of the behaviour of the entanglement entropy in the corresponding geometry [40, 42]. We thus propose the ansatz

$$\begin{eqnarray}&&\mathcal{E}=a{\rm ln} t+b{\rm ln} (\ell -t)+c{\rm ln} (\ell +t)+d{\rm ln} \ell +{\rm constant},\end{eqnarray} \tag{ 33 }$$

where the coefficients are fitting parameters. In fact, some of them can be fixed by requiring that in the limit $\ell \to \infty$ we recover the result (30), implying $a=1/2$ and $d=-(b+c)$. For the remaining two we obtain $b\approx 0.15$ and $c\approx 0.13$ from a fit to the data with $\ell =100$. The data is then scaled together using these values in the right of figure 5 and shows a nice collapse.

Zoom In Zoom Out Reset image size

6.2. 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 equation (9) with the sum replaced by an integral. First, we consider the unbiased case ${{\beta }_{{\rm l}}}={{\beta }_{{\rm r}}}=\beta$, with results on the time evolution of $\mathcal{E}$ shown on the left of figure 6. When compared to the local quench results in figure 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\approx \ell$ and after a sudden decrease $\mathcal{E}$ relaxes slowly towards its thermal-state value (28).

Zoom In Zoom Out Reset image size

The case of unequal temperatures is shown on the right of figure 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, $\mathcal{E}({{\beta }_{{\rm l}}})$ and $\mathcal{E}({{\beta }_{{\rm r}}})$, is indistinguishable from the data $\mathcal{E}({{\beta }_{{\rm l}}},{{\beta }_{{\rm r}}})$ for unequal temperatures. There are, however, again some small deviations.

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 equation (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, 45], and thus the additivity of the chiral contributions does not hold for this measure. Note, however, that in the CFT limit ${{\beta }_{{\rm l}}},{{\beta }_{{\rm r}}}\gg 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 [46] or of special integrable quantum field theories [47].

The work of VE was realized in the framework of TÁMOP 4.2.4.A/1-11-1-2012-0001 'National Excellence Program'. The project was supported by the European Union and co-financed by the European Social Fund. ZZ acknowledges funding by the British Engineering and Physical Sciences Research Council (EPSRC), the Basque Government (Project No. IT4720-10), and by the European Union through the ERC Starting Grant GEDENTQOPT and the CHIST-ERA QUASAR project.

In this appendix we derive the NESS covariance matrix in equation (11), following the steps of [27] and correcting a factor of two error there. The main idea of the calculation is that the only relevant information about the initial state, stored in ${{C}_{m,n}}(0)$, which survives in the limit $t\to \infty$ of the time evolution is located asymptotically far away from the initial cut. There the covariances are translationally invariant and given by

$$\begin{eqnarray}\mathop{{\rm lim} }\limits_{m,n\to \pm \infty }{{C}_{m,n}}(0)=\sigma _{m,n}^{\alpha }=\int _{-\pi }^{\pi }\frac{{\rm d}q}{2\pi }{{{\rm e}}^{{\rm i}q(m-n)}}\left( \begin{array}{ccccccccccccccc} \omega _{q}^{-1} & 0 \\ 0 & {{\omega }_{q}} \\ \end{array} \right){\rm coth} \frac{{{\beta }_{\alpha }}{{\omega }_{q}}}{2},\end{eqnarray} \tag{ A.1 }$$

where the + (−) sign in the limit corresponds to the right (left) hand side with $\alpha ={\rm r}$ ($\alpha ={\rm l}$). The Fourier transform of $\sigma _{m,n}^{\alpha }$ is denoted by $\sigma _{q}^{\alpha }$ and we introduce the notation $\sigma _{q}^{\pm }=(\sigma _{q}^{{\rm r}}\pm \sigma _{q}^{{\rm l}})/2$.

Now we split up the initial covariance matrix in three terms as ${{C}_{m,n}}(0)=C_{m,n}^{1}(0)+C_{m,n}^{2}(0)+C_{m,n}^{3}(0)$ where

$$\begin{eqnarray}&&C_{m,n}^{1}(0)=\frac{\sigma _{m,n}^{{\rm r}}+\sigma _{m,n}^{{\rm l}}}{2},\quad C_{m,n}^{2}(0)={\rm sgn} (m)\frac{\sigma _{m,n}^{{\rm r}}-\sigma _{m,n}^{{\rm l}}}{2},\end{eqnarray} \tag{ A.2 }$$

and $C_{m,n}^{3}(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 equation (8), which have to be compensated by $C_{m,n}^{3}(0)$. The time evolved matrices are given by ${{C}^{i}}(t)=S(t){{C}^{i}}(0)S{{(t)}^{T}}$ such that $C(t)={{C}^{1}}(t)+{{C}^{2}}(t)+{{C}^{3}}(t)$.

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_{q}^{1}(t)={{S}_{q}}(t)\sigma _{q}^{+}S_{q}^{T}(t)$. The symbol of the matrix in equation (10) can be rewritten as a sum of diagonal and offdiagonal contributions

$$\begin{eqnarray}{{S}_{q}}(t)={\rm cos} ({{\omega }_{q}}t)\;{\mathbb{1}}+{\rm sin} ({{\omega }_{q}}t)\;{{\Gamma }_{q}},\qquad {{\Gamma }_{q}}=\left( \begin{array}{ccccccccccccccc} 0 & \omega _{q}^{-1} \\ -{{\omega }_{q}} & 0 \\ \end{array} \right).\end{eqnarray} \tag{ A.3 }$$

Multiplying out $C_{q}^{1}(t)$ and taking $t\to \infty$, the rapidly oscillating terms can be substituted by their average and we arrive at

$$\begin{eqnarray}&&\mathop{{\rm lim} }\limits_{t\to \infty }C_{m,n}^{1}(t)=\int _{-\pi }^{\pi }\frac{{\rm d}q}{2\pi }{{{\rm e}}^{{\rm i}q(m-n)}}\frac{1}{2}\left( \sigma _{q}^{+}+{{\Gamma }_{q}}\sigma _{q}^{+}\Gamma _{q}^{T} \right),\end{eqnarray} \tag{ A.4 }$$

which gives the first term in the integral of equation (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_{m,n}^{2}(0)$ which leads to

$$\begin{eqnarray}&&C_{m,n}^{2}(t)=\int _{-\pi }^{\pi }\frac{{\rm d}q}{2\pi }{{{\rm e}}^{{\rm i}q(m-n)}}P\int _{-\pi }^{\pi }\frac{{\rm d}p}{2\pi }{{{\rm e}}^{-{\rm i}(p-q)n}}\frac{2{\rm i}}{p-q}{{S}_{q}}(t)\sigma _{p}^{-}S_{p}^{T}(t),\end{eqnarray} \tag{ A.5 }$$

where the P sign indicates that the Cauchy principal value has to be taken. It can be evaluated using the identity

$$\begin{eqnarray}&&\mathop{{\rm lim} }\limits_{t\to \infty }P\int _{-\pi }^{\pi }\frac{{\rm d}p}{p-q}g(p){{{\rm e}}^{{\rm i}f(p)t}}={\rm i}\pi g(q){{{\rm e}}^{{\rm i}f(q)t}}\;{\rm sgn} \left[ f^{\prime} (q) \right],\end{eqnarray} \tag{ A.6 }$$

where f and g denote some smooth functions. In turn, the principal value integral has the effect of interchanging the oscillatory terms ${\rm cos} ({{\omega }_{p}}t)\to -{\rm sin} ({{\omega }_{q}}t)$ and ${\rm sin} ({{\omega }_{p}}t)\to {\rm cos} ({{\omega }_{q}}t)$ in $S_{p}^{T}(t)$. In the remaining integral over q, one can again take the averages of the time dependent terms which yields

$$\begin{eqnarray}&&\mathop{{\rm lim} }\limits_{t\to \infty }C_{m,n}^{2}(t)=\int _{-\pi }^{\pi }\frac{{\rm d}q}{2\pi }{{{\rm e}}^{{\rm i}q(m-n)}}{\rm i}\;{\rm sgn} (\omega _{q}^{\prime })\frac{1}{2}\left( \sigma _{q}^{-}\Gamma _{q}^{T}-{{\Gamma }_{q}}\sigma _{q}^{-} \right).\end{eqnarray} \tag{ A.7 }$$

Noting that in our case ${\rm sgn} (\omega _{q}^{\prime })={\rm sgn} (q)$ and calculating the matrix products leads to the second term in equation (11).

We finally show that ${{{\rm lim} }_{t\to \infty }}{{C}^{3}}(t)=0$. Using the form of $C(0)$ in equations (8) and (9) as well as the expression for the eigenvectors in (4), one has in the thermodynamic limit ($N\to \infty$) the block-matrix form

$$\begin{eqnarray}{{C}^{3}}(0)=-\left( \begin{array}{ccccccccccccccc} {{\tilde{\sigma }}^{{\rm l}}} & {{\sigma }^{{\rm l}}} \\ {{\sigma }^{{\rm r}}} & {{\tilde{\sigma }}^{{\rm r}}} \\ \end{array} \right),\qquad \tilde{\sigma }_{m,n}^{\alpha }=\int _{-\pi }^{\pi }\frac{{\rm d}q}{2\pi }{{{\rm e}}^{{\rm i}q(m+n)}}\sigma _{q}^{\alpha }.\end{eqnarray} \tag{ A.8 }$$

The ${{\tilde{\sigma }}^{\alpha }}$ in the diagonal are Hankel matrices (with symbol $\sigma _{q}^{\alpha }$) that arise due to the Dirichlet boundary condition in the center, while the offdiagonal Toeplitz matrices ${{\sigma }^{\alpha }}$ 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

$$\begin{eqnarray}&&C_{m,n}^{3}(t)=\mathop{\sum }\limits_{\alpha =l,r}\mathop{\sum }\limits_{s=\pm }\int _{-\pi }^{\pi }\frac{{\rm d}q}{2\pi }P\int _{-\pi }^{\pi }\frac{{\rm d}p}{2\pi }\int _{-\pi }^{\pi }\frac{{\rm d}p^{\prime} }{2\pi }{{{\rm e}}^{{\rm i}(pm-{{p}^{\prime }}n)}}\mathcal{F}_{s}^{\alpha }(q,p,p^{\prime} ){{S}_{p}}(t)\sigma _{q}^{\alpha }S_{p^{\prime} }^{T}(t),\end{eqnarray} \tag{ A.9 }$$

with the definition

$$\begin{eqnarray}&&\mathcal{F}_{\pm }^{{\rm r}}(q,p,p^{\prime} )=\frac{1}{2}\left[ 2\pi \delta (p-q)-\frac{2{\rm i}}{p-q} \right]\frac{1}{2}\left[ 2\pi \delta (p^{\prime} \pm q)\pm \frac{2{\rm i}}{p^{\prime} \pm q} \right],\end{eqnarray} \tag{ A.10 }$$

and $\mathcal{F}_{\pm }^{{\rm l}}(q,p,p^{\prime} )={{\left[ \mathcal{F}_{\pm }^{{\rm r}}(q,p^{\prime} ,p) \right]}^{*}}$. Note, that $\mathcal{F}_{\pm }^{{\rm r}}$ is just the Fourier transform of the product of step functions $\left[ 1+{\rm sgn} (m) \right]/2\times \left[ 1\pm {\rm sgn} (n) \right]/2$ which singles out the lower right/left block. Carrying out the integrals over p and $p^{\prime}$ in equation (A.9) using formula (A.6), one arrives at

$$\begin{eqnarray}&&\begin{array}{l} \mathop{{\rm lim} }\limits_{t\to \infty }P\int _{-\pi }^{\pi }\frac{{\rm d}p}{2\pi }\int _{-\pi }^{\pi }\frac{{\rm d}p^{\prime} }{2\pi }{{{\rm e}}^{{\rm i}(pm-{{p}^{\prime }}n)}}\mathcal{F}_{\pm }^{{\rm r}}(q,p,p^{\prime} ){{S}_{p}}(t)\sigma _{q}^{{\rm r}}S_{p^{\prime} }^{T}(t) \\ \quad =\frac{{{{\rm e}}^{{\rm i}q(m\pm n)}}}{4}\left[ {{S}_{q}}(t)\sigma _{q}^{{\rm r}}S_{q}^{T}(t)-{{{\hat{S}}}_{q}}(t)\sigma _{q}^{{\rm r}}\hat{S}_{q}^{T}(t) \right. \\ \qquad -\left. {\rm i}\;{\rm sgn} (\omega _{q}^{\prime })({{{\hat{S}}}_{q}}(t)\sigma _{q}^{{\rm r}}S_{q}^{T}(t)+{{S}_{q}}(t)\sigma _{q}^{{\rm r}}\hat{S}_{q}^{T}(t)) \right] \\ \end{array}\end{eqnarray} \tag{ A.11 }$$

where we defined ${{\hat{S}}_{q}}(t)=-{\rm sin} ({{\omega }_{q}}t){\mathbb{1}}+{\rm cos} ({{\omega }_{q}}t){{\Gamma }_{q}}$ and used the symmetry properties ${{S}_{-q}}(t)={{S}_{q}}(t)$, ${{\hat{S}}_{-q}}(t)={{\hat{S}}_{q}}(t)$ and ${\rm sgn} (\omega _{-q}^{\prime })=-{\rm sgn} (\omega _{q}^{\prime })$. Calculating the matrix products, one finds that the expression in the brackets vanishes identically. The same holds for the integrals with $\mathcal{F}_{\pm }^{{\rm l}}$ which concludes our proof.

Here we briefly summarize the method used to obtain $\mathcal{E}$ for a local quench. For a detailed discussion of ground-state entanglement negativity within CFT, we refer to [15].

The essential step is to rewrite equation (23) as

$$\begin{eqnarray}&&\mathcal{E}=\mathop{{\rm lim} }\limits_{{{n}_{e}}\to 1}{\rm lnTr}{{(\rho _{A}^{{{T}_{2}}}(t))}^{{{n}_{e}}}},\end{eqnarray} \tag{ B.1 }$$

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 ${{\rho }_{A}}(t)={\rm T}{{{\rm r}}_{B}}\;\rho (t)$ in terms of a path integral. This is done by considering an n_e-sheeted Riemann surface and sewing together the replicas of $\rho _{A}^{{{T}_{2}}}(t)$. Each copy can be represented by a 2D path integral with slits along the intervals A₁ and A₂ on the real axis, and partial transposition T₂ corresponds to interchanging the edges of the corresponding slit A₂.

Instead of working on a replicated world-sheet, one can introduce replicated fields and the so-called twist fields, ${{\mathcal{T}}_{{{n}_{e}}}}$ and ${{\bar{\mathcal{T}}}_{{{n}_{e}}}}$, 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=\varnothing$) into two semi-infinite lines A₁ and A₂. One has then a single contact point and thus

$$\begin{eqnarray}&&{\rm Tr}{{(\rho _{A}^{{{T}_{2}}}(t))}^{{{n}_{e}}}}=\langle \mathcal{T}_{{{n}_{e}}}^{2}({{z}_{0}})\rangle ,\end{eqnarray} \tag{ B.2 }$$

where the expectation value of the twist fields has to be calculated on a world-sheet representing the time-evolved density operator $\rho (t)$. This is depicted on the left of figure 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 $\rho (0)$. The real time evolution takes place between the endpoints $z=\pm {\rm i}\epsilon$ of the cut, and the twist field has to be inserted at ${{z}_{0}}={\rm i}\tau$. Note, that the parameter is needed to regularize the path-integral and the analytical continuation to real times $\tau ={\rm i}t$ must be carried out at the end of the calculation.

Zoom In Zoom Out Reset image size

Although the original world-sheet has a complicated geometry, one can apply a conformal mapping [42]

$$\begin{eqnarray}&&w=\frac{z}{\epsilon }+\sqrt{{{\left( \frac{z}{\epsilon } \right)}^{2}}+1}\end{eqnarray} \tag{ B.3 }$$

which transforms it to the half-plane, as shown on the right of figure B1. On the half-plane the expectation value of a one-point function is known and one can then use the conformal transformation formula to obtain

$$\begin{eqnarray}&&\langle \mathcal{T}_{{{n}_{e}}}^{2}({{z}_{0}})\rangle \propto {{\left( {{\left| \frac{{\rm d}w}{{\rm d}z} \right|}_{{{z}_{0}}}}\frac{1}{\operatorname{Re}({{w}_{0}})} \right)}^{{{\Delta }_{{{n}_{e}}}}}},\end{eqnarray} \tag{ B.4 }$$

where ${{w}_{0}}=w({{z}_{0}})$ and the scaling dimension of the operator $\mathcal{T}_{{{n}_{e}}}^{2}$ is given by [15]

$$\begin{eqnarray}&&{{\Delta }_{{{n}_{e}}}}=\frac{c}{6}\left( \frac{{{n}_{e}}}{2}-\frac{2}{{{n}_{e}}} \right),\end{eqnarray} \tag{ B.5 }$$

with the central charge c. Carrying out the calculations, one obtains

$$\begin{eqnarray}&&{\rm Tr}{{(\rho _{A}^{{{T}_{2}}}(\tau ))}^{{{n}_{e}}}}\propto {{\left( \frac{\epsilon }{{{\epsilon }^{2}}-{{\tau }^{2}}} \right)}^{{{\Delta }_{{{n}_{e}}}}}}.\end{eqnarray} \tag{ B.6 }$$

Continuing the result to real time $\tau \to {\rm i}t$, taking the limit $t\gg \epsilon$ and finally using equation (B.1), one arrives to the result in equation (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 $\operatorname{Re}(z)=\pm \ell$, and the proper conformal mapping to the half-plane was given in [43]. Carrying out the analogous steps, one arrives to the formula in equation (31).

On the other hand, the tripartite case in the infinite system with line segments ${{A}_{1}}=\left[ -\ell ,0 \right]$ and ${{A}_{2}}=\left[ 0,\ell \right]$ is more involved. There, one has to consider the three-point function $\langle \mathcal{T}({{z}_{-1}})\bar{\mathcal{T}}_{{{n}_{e}}}^{2}({{z}_{0}})\mathcal{T}({{z}_{1}})\rangle$ of the twist fields with ${{z}_{0}}={\rm i}\tau$ and ${{z}_{\pm }}1=\pm \ell +{\rm i}\tau$. 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 [48], their structure is rather involved and we have not been able to tackle this case analytically.