New Journal of Physics

The open access journal at the forefront of physics

Deutsche PhysikalischeGeseUschaft DPG IOP Institute Of PhySjCS

Entanglement negativity in the harmonic chain out of equilibrium

Viktor Eisler and Zoltan Zimboras

1 Institute for Theoretical Physics, Eotvos Lorand University, Pazmany setany 1/a, H-1117 Budapest, Hungary

2 Department of Computer Science, University College London, Gower Street, WC1 E 6BT London, UK

3 Department of Theoretical Physics, University of the Basque Country UPV/EHU, PO Box 644, E-48080 Bilbao, Spain

E-mail: zimboras@gmail.com

Received 18 August 2014

Accepted for publication 14 October 2014

Published 8 December 2014

New Journal of Physics 16 (2014) 123020

doi:10.1088/1367-2630/16/12/123020

Abstract

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.

Keywords: entanglement in many-body systems, quantum statistical mechanics, non-equilibrum statistical physics, conformal field theory

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

New Journal of Physics 16 (2014) 123020

1367-2630/14/123020+18$33.00 © 2014 IOP Publishing Ltd and Deutsche Physikalische Gesellschaft

1. Introduction

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

t=o t>0

Figure 1. Geometry of the oscillator chain before and after the quench.

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.

2. Model and setting

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

H = 12 [p2 + K(Xn+1 - Xn)2 + m0x2]' (1)

where xn and pn denote the position and momentum operators of the nth oscillator with canonical commutation relations [xm, pn ] = i5mn and [xm, xn ] = [pm, pn ] = 0. The parameters K

and w0 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 x0 = xN+1 = 0 and p0 = pN +1 = 0. The initial state is given by a simple product of Gibbs states

p (0) = — Hl ® Hr, (2)

Z i Z r

with inverse temperatures fa and fir for the left and right half-chains, respectively. The Hamiltonian Hl (Hr) 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 p (t) = e-lHtp (0)elHt under the action of Hamiltonian (1). The change in the geometry of the chain is depicted in figure 1.

Since both the initial state as well as the time evolution operator is Gaussian, p (t) can be fully characterized by its 2N X 2N covariance matrix C (t). This can be written in a blockmatrix form, composed of symmetric 2x2 matrices

Cm, n (t) = (Rm (t) Rl (t) + Rn (t) Rl (t)) , Rn (t) =

' xn (t) N pn (t)

where Rn(t) is a vector of the position and momentum operators in the Heisenberg picture at site n, and {O (t)) = Tr[p (0)O (t)] denotes the average of observable O (t) taken with the initialstate density operator.

To obtain the time evolution of Rn(t), we first introduce new operators by the canonical transformation

£ (n)(n)

nkn N + 1

1,..., N,

and similarly for pk. Note, that (n) are just the eigenvectors of the dynamical matrix in the Hamiltonian, satisfying Dirichlet boundary conditions 4>k (0) = 4>k (N + 1) = 0. The Hamilto-nian is then transformed into

H = 2 £ ( + " k2 Xk2),

with a dispersion relation

Mk = J m0 + 4 K sin2(qk/2) , qk =

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 xk (t) = pk (t) and pk (t) = —mk xk (t), which yields R (t) = S (t) R (0) with a symplectic matrix S(t) of block form

Sm,n (t ) = £ 4>* (m) (n)

cos mkt m- 1 sin mkt —mk sin mkt cos mkt

The time-evolution of the covariance matrix then reads

' Cl 0 0 Cr

C (t) = S (t) C (0) S (t)T, C (0) =

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

Cm, n = £4>t,k(m) k (n)

ma,k 0

0 m a,k

with a = l, r and 1 ^ m, n ^ N/2. Note, that $ak and wa k are the eigenvectors and eigenvalues of the dynamical matrix of the half-chains and are obtained by exchanging N — N/2 in equations (4) and (6), respectively.

3. The steady state

Our primary goal is to calculate the entanglement in the steady state, i.e., in the asymptotic limit t —> to 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 — to. The eigenvectors of the

dynamical matrix are then Fourier modes, ( (n) ~ e iqn, and the corresponding limit for the elements of the symplectic matrix in equation (7) is given by

dq -n 2n

giq (m - n)

cos mqt mq sin mq t

mq 1 sin mq t cos mqt

where m, n G Z and the dispersion rnq has the same form as in equation (6), but with continuous quasi-momenta q G (—n, n).

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. 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 and lmI, In I » 1.

The formal derivation of C(t) in the limit t ^ to 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

. —e ^ ' i cq '—n 2n V q

with 2x2 matrices

limCm,n(t) = in dqeiq(m-n)(c+ + i sgn(q) Cq"),

t^ro J-n 2n v '

C++ = -

Cq-= 1

' ßr mqx

+ coth

- coth

ßl mq

ßl mq

(0 -01)

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 ^ /?, 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 am and creation am operators, defined as the Fourier transforms of the modes

Xq + —Pq

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

lim Tr (p (t) am an ) = f " ^eiq (m-n) nq,

t^ro v ' J-n 2n

where the bosonic mode-occupation number is given by

—-1- q e (—n, 0),

g Pi Wq — 1

6 j 1 (16) - q e (0, n).

g Pl Wq — 1 *

The form of nq has a simple physical interpretation. Namely, all the normal modes with positive (negative) group velocities — originate from the left (right) reservoir and thus their

occupation numbers are given by the corresponding Bose-Einstein statistics with inverse temperature Pl (pp). 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 Pl = pp. 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]

limp (t) = -e-Hgff, Hgff = £ ((Q+ + v— Q—), (17)

t z „

where the integrals of motion can be constructed in the N ^ o limit of a periodic chain and read [27]

Qn = 1 ^PmPm+n + (xm — xm+ 1)(xm+n — xm+n +1) + W0 xmxm+n(18)

m=—o 1 o

Qn = 1 £ (pmxm+n — xmpm+n — pm+nXm + xm+npm ). (19)

m=—o

The corresponding sets of Lagrange multipliers v± can be determined by taking the Fourier

transform of Hgff, rewriting it in terms of the bosonic operators aq and a J, and requiring that the resulting expression reproduces the mode-occupations in equation (16). This leads to the equations

£ (Wq cos nq + v— sin nq) = WqPl, (20)

£ ( (Wq cos nq — V— sin nq) = Wqpr, (21)

which can be solved and yield the Lagrange multipliers

Pl + PL — /„ „ X, _4

5 n,0, V— = ((i — Pr)(—1)n "71—7. (22)

2 'n ^ n 4n — 1

Note, that the only reflection-symmetric conserved charge appearing in the GGE is Q+ = H the original Hamiltonian itself and the corresponding multiplier is simply the average inverse temperature. On the other hand, all the Q— charges have nonvanishing multipliers, decaying slowly in ~ 1/n for n » 1 and alternating in sign. Therefore, the operator Heff 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

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, |i~ = 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.

4. Partial transpose and logarithmic negativity

Our aim is to characterize the amount of entanglement in the time-evolved state p (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 A1 and A2, with A = A1 U A2 and B denoting the rest of the chain. The logarithmic negativity is then defined through the partial transpose of the reduced density matrix pA (t) = TrB p (t) as

S = lnTr

pAT (t) , (23)

where the superscript T2 indicates a partial transposition with respect to the indices in subsystem A2. The logarithmic negativity thus detects only the negative eigenvalues of pj2(t). In

particular, if there is no entanglement between A1 and A2, then all the eigenvalues of pj2 (t) are positive [34] and S 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 pA (t) are related to the symplectic eigenvalues of the reduced covariance matrix CA (t). Moreover, the partial transpose can also be implemented on the level of the covariances, with Cj2 (t) denoting the matrix where the signs of all the momenta pn with n G A2 are reversed. In turn, the logarithmic negativity is obtained as [9]

S = — ^ ln min (vj,1) (24)

from the symplectic eigenvalues Vj of Cj2 (t) and IAI denotes the number of sites in A. Note, that only the eigenvalues Vj < 1 contribute.

The symplectic spectrum of Cj2 (t) can be obtained through ordinary diagonalization, by multiplying with the symplectic matrix

= © (»J). (25)

which leads to the spectrum

Spect^C,2 (t)) = { ±1*1..... ±iv|A|}. (26)

5. Steady-state entanglement

We are interested in the entanglement of two neighbouring subsets of oscillators A1 and A2, each of size I. Before presenting results for the logarithmic negativity S, 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, rn0 ^ 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 S is insensitive to the zero-mode contribution and converges as m0 ^ 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 m0 = 0.005//, i.e., reducing the gap with the size of the subsystems.

5.1. Equal temperatures

We first consider the simplest case fcl = fc = fc. 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 p ~ exp(-fcH). 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, fc » 1, as a function of I and fc. 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 l1 and l2 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]

c l1 l2

S = — ln--+ constant. (27)

4 h + l2 K '

where c is the central charge of the CFT. For finite temperatures, the CFT is defined on an infinite cylinder of circumference fc which, however, can be mapped into the plane by a simple conformal transformation. The overall effect of the mapping is to replace the lengths li ^ fc /n sinh(/;n/fc). For intervals of equal lengths, l1 = t2 = l, this leads to

1 , fc , ln .„_.

S = — ln — tanh--+ constant. (28)

4 n fc

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

The CFT formula (28) has the correct limiting behaviour S ~ 1/4 ln l for l ^ fc whereas in the opposite limit of large segment sizes, l » fc, it predicts the saturation value S ~ 1/4 ln fc.

5=100 -

B=100 *

20 40 60 80 100 120

/J/7rtanh(&r//?)

Figure 2. Left: thermal-state logarithmic negativity for two adjacent intervals of size I in an infinite chain for various p. Right: scaled data according to CFT prediction. The dotted line has slope 1/4.

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).

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 p ~ exp (_faxH+ —faH_) with H± 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 p, the action thus decouples in the chiral fields and that live on infinite cylinders of circumferences fax and fa, 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 £ fa, fa) 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± 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 Vj fa, fa) < 1 of the partial transposed covariance matrix which contribute to

Figure 3. The five smallest symplectic eigenvalues Vj ft ft) < 1 of the partial transposed steady-state covariance matrix for two adjacent intervals of size I = 100 and various pairs of ftl and ft. The inset shows the deviation from the geometric mean of Gibbs-state symplectic eigenvalues vjft) and vj(ft).

£ ft ft), see equation (24). The small deviations from the geometric mean ^Vj ft)Vj (ft) 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].

6. Time evolution of entanglement

We now consider the complete time evolution of £ after the wall separating the two sides of the chain is removed, see figure 1. The subsystems A1 and A2 are chosen to be adjacent intervals with their common boundary located at the initial cut. Clearly, at t = 0 one has £ = 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 = 0, the result can directly be inferred from previous CFT calculations. The simplest choice is to consider the evolution of £ for two halves of an infinite system, A1 = (-to, 0] and A2 = [1, to). Since the state p(t) is pure, the logarithmic negativity is just the Renyi entropy with index 1/2, and inserting this into the CFT formula of [42] one obtains

£ = — ln t + constant. (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 = 21, with a bipartition A1 = [1, 1 ] and A2 = [I + 1, 21 ] with the initial cut located between

2.5 2 1.5 1

- 1 a..-' ¿ = 25 ° I = 50 » i = 75 ° t = 100 -

2£/7rsin(7rt/2£)

Figure 4. Left: time evolution of the half-chain logarithmic negativity after a local quench in the finite geometry for various I. Right: scaled data according to CFT prediction for times t < 21. The dotted line has slope 1/2.

sites I and I + 1. The corresponding result for £ can also be found from a CFT calculation based on [43] and reads

£ = 1 In 121/n sin (ntl21 )| + constant.

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 I and shows a quasi-periodic structure with period 21, 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 I = 25 which is supposedly due to lattice effects, caused by the slower normal-mode excitations with velocity -q < 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.

Finally, one could consider £ 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 £. The ground-state covariance matrix of the semi-infinite chain, with matrix elements given by the N —> to limit of equation (9), can be evaluated explicitly for m0 = 0 and yields [15]

Cm n = Cr(m + n) — Cr(m — n),

C r(x) = 1

W (1/2 + x) 0

4x2 — 1

for the right-hand side of the chain, m, n ^ 1, with w (z) being the digamma function. The left-hand side covariance matrix C^ n for m, n ^ 0 is obtained by a reflection of the indices m — 1 — m and n — 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

2.5 2 1.5 1

J^r □ □ A 0

- ¿ = 25 °

¿ = 50 *

¿ = 75 °

£ = 100 -

2.5 2 1.5 1

¿ = 25 ¿ = 50 ¿ = 75 ¿=100

ta{t - tf{t + t)c/0+c

Figure 5. Left: time evolution of entanglement after a local quench in the infinite geometry for various I. Horizontal lines indicate the ground-state values of £ for I = 25, 50. Right: scaled data according to equation (33) with parameters a = 0.5, b a 0.15, c a 0.13 and d = —(b + c).

can make use of the light-cone structure of the matrix elements, implying that for \m — n I » t the entries Sm,n (t) are exponentially small. Thus, the sums involved in Cj2 (t) can be truncated and evaluated numerically to a very high precision.

The result for £ is shown in figure 5. For times t < l, the logarithmic negativity develops a plateau, followed by a sharp drop at t a l, where the propagating front leaves the subsystem A. The data then converges slowly to the ground-state value of £, 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

£ = a ln t + b ln (l — t) + c ln (l + t) + d ln l + constant, (33)

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

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 fa = fa = p, with results on the time evolution of £ 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 a l and after a sudden decrease £ relaxes slowly towards its thermal-state value (28).

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, £ (fal) and £ (far), is indistinguishable from the data £ (fal, fa) for unequal temperatures. There are, however, again some small deviations.

Figure 6. Time evolution of entanglement for adjacent intervals of size I = 50 in the infinite geometry. Left: equal temperatures ftx = ft = ft. Horizontal lines indicate the steady-state entanglement for ft = 5, 10. Right: entanglement evolution for unequal temperatures (symbols) compared to averages of equal-temperature curves (lines).

7. 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 nq, 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 ftx, ft » 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].

Acknowledgments

The work of VE was realized in the framework of TAMOP 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.

Appendix A. Steady-state covariance matrix

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 Cm,n (0), which survives in the limit t ^ œ of the time evolution is located asymptotically far away from the initial cut. There the covariances are translationally invariant and given by

lim Cmn (0) = a,

= fn dq

J-n 2n

dq giq (m—n)

ftteq 2 '

where the + (-) sign in the limit corresponds to the right (left) hand side with a = r (a = l). The Fourier transform of a^, n is denoted by a% and we introduce the notation a± = (a£ ± axq )/2. Now we split up the initial covariance matrix in three terms as

Cm,n (0) = cm, n (0) + cm, n (0) + cm, n (0) where

Cm, n (0)

am, n + am, n

Cm, n (0) = sgn (m)-

and Cm 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 equation (8), which have to be compensated by Cm n (0). The time evolved matrices are given by Cl (t) = S (t) Cl (0)S (t)T such that C (t ) = C1(t ) + C 2(t ) + C 3(t ).

First, we consider C1(t) which is a product of three Toeplitz matrices and thus the multiplication can be performed in Fourier space C^(t) = Sq(t)a+ST(t). The symbol of the matrix in equation (10) can be rewritten as a sum of diagonal and offdiagonal contributions

Sq (t) = COS (teq t) 1 + Sin (teq t) Tq,

/ 1 \ 0 te—

—teq 0

Multiplying out Cq(t) and taking t ^ to, the rapidly oscillating terms can be substituted by their average and we arrive at

lim C,

1, n(t) = £ f eiq(m—n)2(a+ + ra+rT),

which gives the first term in the integral of equation (11).

The next step is to obtain the asymptotics of C2(t). Working again in Fourier space, one has now to include the Fourier transform of the sign function in Cm2 n (0) which leads to

*—'tv

dq- ei

(t) = f dqeiq (m-n) P f —e-i(P - q) J-n 2n J-n 2n p

2i T -Sq(t)a~STp (t),

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

j —j

t^œ J - n p

-g(p)eif (p)t = ing(q)eif (q)t sgnf (q)],

where f and g denote some smooth functions. In turn, the principal value integral has the effect of interchanging the oscillatory terms cos(m pt) — -sin(ft>qt) and sin(^pt) — cos(mqt) in

Sp (t). In the remaining integral over q, one can again take the averages of the time dependent terms which yields

limClt n(t) = i' dqeiq(m-n)i sgn (®q)Uo-rTq - rqa-). (A.7)

t —>to J-n 2n 2 v '

Noting that in our case sgn(a>q) = sgn(q) and calculating the matrix products leads to the second term in equation (11).

We finally show that limt—toC3 (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) the block-matrix form

C 3(0)=-( 5 i,

—q„ i

rn dq J-n 2n

iq (m+n) aa °q ■

The aa in the diagonal are Hankel matrices (with symbol a^) that arise due to the Dirichlet boundary condition in the center, while the offdiagonal Toeplitz matrices aa compensate the extra contributions from C1 (0) + C2 (0), as mentioned before. The time-evolved matrix C3 (t) has thus four different contributions from the various Hankel/Toeplitz blocks which can be written as

a=l,r s=±

with the definition

Cm, n (t) = ££ fn —qP Jn —p fn ei(pm-p'n ra (q, p, p') Sp (t) aa ST (t), . J-n 2k J-n 2n J-n 2n

'-n 2n J-n 2n

r±(q, p, p') =

2n5 (p - q)

2nô (p' ± q) ±

(A.10)

and T± (q, p, p') = [ T± (q, p', p) J*. Note, that T± is just the Fourier transform of the product of step functions [1 + sgn(m)]/2 x [1 ± sgn(n )]/2 which singles out the lower right/left block. Carrying out the integrals over p and p' in equation (A.9) using formula (A.6), one arrives at

limP f * f * ei(pm-pnT±(q, p, p')Sp(t)a^S^(t)

t — to j -n 2k j -n 2k

eiq (m± n) p t

■[ Sq (t) aq SqT (t) - Sq (t) ö^ St (t)

(q )(S

sq (t) aq st (t ) + sq (t) aq s: (t)

(A.11)

z_, Z„ z.

Figure B1. Conformal map from the quench geometry to the half-plane. Red dots indicate the locations of twist-field insertions and their images under w(z).

where we defined Sq(t) = —sin(ft>qt) 1 + cos(a>qt)rq and used the symmetry properties

S-q (t) = Sq (t), S—q (t) = Sq (t) and sgn (ft>—q) = — sgn(ftq). Calculating the matrix products, one finds that the expression in the brackets vanishes identically. The same holds for the integrals with T± which concludes our proof.

Appendix B. CFT results for local quench

Here we briefly summarize the method used to obtain S 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

S = limlnTr([2(t)), (B.1)

where ne 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 pA (t) = TrB p (t) in terms of a path integral. This is done by considering

an ne-sheeted Riemann surface and sewing together the replicas of pA2 (t). Each copy can be represented by a 2D path integral with slits along the intervals A1 and A2 on the real axis, and partial transposition T2 corresponds to interchanging the edges of the corresponding slit A2.

Instead of working on a replicated world-sheet, one can introduce replicated fields and the so-called twist fields, TUe and Tn, 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 = 0) into two semi-infinite lines A1 and A2. One has then a single contact point and thus

Tr(pl2(t)f = (r 2e (zo)), (B.2)

where the expectation value of the twist fields has to be calculated on a world-sheet representing the time-evolved density operator p (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 p (o). The real time evolution takes place between the endpoints z = ±ie of the cut, and the twist field has to be inserted at z0 = it. Note, that the

parameter e is needed to regularize the path-integral and the analytical continuation to real times t = 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 [42]

w = —+

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

(r I (zo))

/ dw \ 1

V dz z0 Re(w0),

where w0

w (z0) and the scaling dimension of the operator rUe is given by [15]

/ \ e 2

An. = -

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

Tr (p>))Ue a f.

Continuing the result to real time t — it, taking the limit t » e 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 Re (z) = ±1, 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 A1 = [—I, 0] and A2 = [0, 1 ] is more involved. There, one has to consider the three-point function — 2

{T (z—1) TUe (zo) T (z1)) of the twist fields with z0 = iT and z±1 = ±1 + iT. 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.

References

[1] Amico L, Fazio R, Osterloh A and Vedral V 2008 Rev. Mod. Phys. 80 517

[2] Calabrese P, Cardy J and Doyon B 2009 J. Phys. A: Math. Theor. 42 500301

[3] Vidal G, Latorre J I, Rico E and Kitaev A 2003 Phys. Rev. Lett. 90 227902

[4] Calabrese P and Cardy J 2009 J. Phys. A: Math. Theor. 42 504005

[5] Eisert J, Cramer M and Plenio M B 2010 Rev. Mod. Phys. 82 277

[6] Plenio M and Virmani S 2007 Quantum Inf. Comput. 7 1

[7] Vidal G and Werner R F 2002 Phys. Rev. A 65 032314

[8] Plenio M B 2005 Phys. Rev. Lett. 95 090503

[9] Audenaert K, Eisert J, Plenio M B and Werner R F 2002 Phys. Rev. A 66 042327

[10] Adesso G 2007 and Illuminati F J. Phys. A: Math. Theor. 40 7821

[11] Marcovitch S, Retzker A, Plenio M B and Reznik B 2009 Phys. Rev. A 80 012325

[12] Wichterich H, Molina-Vilaplana J and Bose S 2009 Phys. Rev. A 80 010304(R)

[13] Wichterich H, Vidal J and Bose S 2010 Phys. Rev. A 81 032311

[14] Calabrese P, Cardy J and Tonni E 2012 Phys. Rev. Lett. 109 130502

[15] Calabrese P, Cardy J and Tonni E 2013 J. Stat. Mech. P02008

[16] Calabrese P, Tagliacozzo L and Tonni E 2013 J. Stat. Mech. P05002

[17] Alba V 2013 J. Stat. Mech. P05013

[18] Chung C, Alba V, Bonnes L, Chen P and Läuchli A M 2014 Phys. Rev. B 90 064401

[19] Lee Y A and Vidal G 2013 Phys. Rev. A 88 042318

[20] Castelnovo C 2013 Phys. Rev. A 88 042319

[21] Bernard D and Doyon B 2012 J. Phys. A: Math. Theor. 45 362001

[22] Bernard D and Doyon B 2014 Ann. Inst. Henri Poincare 1 49

[23] Eisler V and Zimboräs Z 2014 Phys. Rev. A 89 032321

[24] Ajisaka S, Barra F and Zunkovic B 2014 New J. Phys. 16 033028

[25] Wolf M M, Verstraete F, Hastings M B and Cirac J I 2008 Phys. Rev. Lett. 100 070502

[26] Spohn H and Lebowitz J L 1977 Commun. Math. Phys. 54 97

[27] Boldrighini C, Pellegrinotti A and Triolo L 1983 J. Stat. Phys. 30 123

[28] Ho T G and Araki H 2000 Proc. Steklov Inst. Math. 228 191

[29] Ogata Y 2002 Phys. Rev. E 66 016135

[30] Aschbacher W H and Pillet C A 2002 J. Stat. Phys. 112 1153

[31] De Luca A, Viti J, Bernard D and Doyon B 2013 Phys. Rev. B 88 134301

[32] Collura M and Karevski D 2014 Phys. Rev. B 89 214308

[33] Collura M and Martelloni G 2014 J. Stat. Mech. P08006

[34] Peres A 1996 Phys. Rev. Lett. 77 1413

[35] Toth G, Knapp C, Gühne O and Briegel H 2007 Phys. Rev. Lett. 99 250405

[36] Anders J and Winter A 2008 Quantum Inf. Comput. 8 0245

[37] Anders J 2008 Phys. Rev. A 77 062102

[38] Ferraro A, Cavalcanti D, Garcia-Saez A and Acin A 2008 Phys. Rev. Lett. 100 080502

[39] Cavalcanti D, Ferraro A, Garcia-Saez A and Acin A 2008 Phys. Rev. A 78 012335

[40] Eisler V and Peschel I 2007 J. Stat. Mech. P06005

[41] Eisler V, Karevski D, Platini T and Peschel I 2008 J. Stat. Mech. P01023

[42] Calabrese P and Cardy J L 2007 J. Stat. Mech. P10004

[43] Stephan J M and Dubail J 2011 J. Stat. Mech. P08019

[44] Stephan J M and Dubail J 2013 J. Stat. Mech. P09002

[45] Ares F, Esteve J, Falceto F and Sänchez-Burillo E 2014 J. Phys. A: Math. Theor. 47 245301

[46] Karrasch C, Ilan R and Moore J E 2013 Phys. Rev. B 88 195129

[47] Castro-Alvaredo O, Chen Y, Doyon B and Hoogeveen M 2014 J. Stat. Mech. P03011

[48] Coser A, Tagliacozzo L and Tonni E 2014 J. Stat. Mech. P01008

Copyright of New Journal of Physics is the property of IOP Publishing and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.