Entanglement dynamics of a many-body localized system coupled to a bath
EEntanglement dynamics of a many-body localized system coupled to a bath
Elisabeth Wybo,
1, 2
Michael Knap,
1, 3, 2 and Frank Pollmann
1, 2 Department of Physics, Technical University of Munich, 85748 Garching, Germany Munich Center for Quantum Science and Technology (MCQST), Schellingstr. 4, D-80799 M¨unchen Institute for Advanced Study, Technical University of Munich, 85748 Garching, Germany (Dated: April 29, 2020)The combination of strong disorder and interactions in closed quantum systems can lead to many-body localization (MBL). However this quantum phase is not stable when the system is coupled to athermal environment. We investigate how MBL is destroyed in systems that are weakly coupled toa dephasive Markovian environment by focusing on their entanglement dynamics. We numericallystudy the third R´enyi negativity R , a recently proposed entanglement proxy based on the negativitythat captures the unbounded logarithmic growth in the closed case and that can be computedefficiently with tensor networks. We also show that the decay of R follows a stretched exponentiallaw, similarly to the imbalance, with however a smaller stretching exponent. I. INTRODUCTION
Interacting quantum systems subjected to strong dis-order can realize an exotic many-body localized (MBL)phase of matter [1–4]. Similarly to a non-interacting An-derson insulator [5], the MBL phase is characterized bythe absence of conventional transport and by spatial cor-relations that decay exponentially with distance. How-ever, there are also important differences for example inthe frequency-dependent response [6, 7] and in the en-tanglement dynamics [8–10]. Most prominently, the en-tanglement entropy grows logarithmically in time [8–10],due to effective interactions between the localized orbitals(localized integrals of motion, LIOMs) [11, 12]. Evidencefor an MBL phase has also been found experimentallyin systems of ultracold atoms, trapped ions or super-conducting qubits by the observation of a non-thermalsaturation value of local densities [13–16] and by entan-glement dynamics [17–20].An important question is on which time scales signa-tures of MBL are observable in real systems, which arenever truly isolated. In general, we expect that a cou-pling of the system to a bath leads to delocalization astransport is restored [21–26]. When considering Marko-vian dephasing noise described by the Lindblad equation,the interference between the LIOMs is destroyed and theMBL state is driven into a featureless infinite tempera-ture state [22–24]. It has been argued that local densities(e.g. the imbalance [13]) show a universal slower thanexponential (specifically, a stretched exponential) decaythat can be explained in terms of the LIOMs [22, 23].The stretched exponential decay has also been observedexperimentally in a cold atom setup [15]. These worksexplain the dynamics of the imbalance in dephasive MBLsystems by means of purely classical rate equations, thatthus only consider the hopping of diagonal states in thedensity matrix. Given the recent experimental focus onentanglement dynamics in MBL systems, it is an openquestion of how Markovian noise affects pure quantumcorrelations.In this work, we investigate how the MBL phase is dy- . . . W ∆Γ FIG. 1. A sketch of our setup: a spin chain with disorderedz-directed fields h i ∈ [ − W, W ], spin exchange ∆ and couplingΓ to a bath. namically destroyed by a dephasive coupling to a bath,and focus on the decay of quantities that are sensitiveto quantum correlations over bipartitions of the system.This is motivated by the fact that one of the most strik-ing dynamical signatures of MBL is a generic logarithmicgrowth of entanglement under a quench which is com-pletely absent in case of Anderson localization. Our goalis to investigate a recently introduced entanglement mea-sure for open quantum systems based on the negativity,the third R´enyi negativity [27–30], that can dynamicallycapture the difference between Anderson localization andMBL provided the dissipation is sufficiently weak. Wemotivate our choice of this quantity, investigate how itscales, and discuss its experimental relevance.The remainder of this paper is organized as follows:We start by introducing our model and setup in Sec. II.In Sec. III we provide some background about entan-glement measures for open quantum systems, and moti-vate our choice to consider the third R´enyi negativity. InSec. IV we present our computations of this quantity andits distribution over disorder realizations, and concludein Sec. V. a r X i v : . [ c ond - m a t . s t r- e l ] A p r II. MODEL AND SETUP
The random-field XXZ Hamiltonian on a chain withopen boundary conditions is given by H = J (cid:34) L − (cid:88) i =1 (cid:0) S xi S xi +1 + S yi S yi +1 + ∆ S zi S zi +1 (cid:1)(cid:35) + L (cid:88) i =1 h i S zi , (1)where S x,y,zi are the spin- operators, and the h i ’sare randomly and uniformly distributed in the interval[ − W, W ]. Here, we will fix the disorder to W = 5 J suchthat our systems are in the MBL phase [31].The time dependence of the density matrix is given bythe Lindblad master equation [32]˙ ρ = − i [ H, ρ ( t )]+ (cid:88) i (cid:18) L i ρ ( t ) L † i − (cid:110) L i L † i , ρ ( t ) (cid:111)(cid:19) ≡ L ( ρ ) , (2)which models the coupling of the system to a Markovian,i.e. memoryless, bath. We consider the Lindblad oper-ators for dephasing noise L i = √ Γ S zi , with this choicethe decoherent part of the Lindblad equation also con-serves the total magnetization (cid:80) i S zi in the system. TheLindblad equation then takes the form˙ ρ = − i [ H, ρ ( t )] + Γ (cid:88) i (cid:18) S zi ρ ( t ) S zi − ρ ( t ) (cid:19) . (3)We sketch our system in Fig. 1. In the limit of a purelydephasive coupling H = 0 the off-diagonal matrix ele-ments of the density matrix just decay exponentially. Sothis specific choice of Lindblad operators removes entan-glement over time.To simulate the time evolution accoring to the Lind-blad equation, we use exact diagonalization and the timeevolving block decimation (TEBD) algorithm on densitymatrices [33, 34], where the time-evolution operator isgiven by U ( t ) = exp( L t ) with the superoperator L = − iH ⊗ + i ⊗ H +Γ (cid:88) i (cid:18) S zi ⊗ S zi − ⊗ (cid:19) . (4)As initial state of the quench we use the N´eel productstate ρ = | ψ (cid:105) (cid:104) ψ | with | ψ (cid:105) = | . . . (cid:105) , so we work inthe sector with total magnetization M = (cid:80) i (cid:104) S zi (cid:105) = 0.The time-evolution operator acts on a vectorized ver-sion of the density matrix in which the spin indices arecombined as | ρ ( t + d t (cid:105) = exp( L d t ) | ρ ( t ) (cid:105) . We note thatthe efficiency of our TEBD simulation of the density ma-trix critically depends on the entropy of | ρ ( t ) (cid:105) viewedas a pure state in operator space. This operator-spaceentropy cannot be easily related to the quantum entan-glement of the density matrix as it also contains clas-sical correlations [35]. In our setup, it has been shownthat this quantity grows logarithmically which allows fora simulation over long times [8, 24]. At late times theoperator-space entropy converges to a value set by the steady state, which in our case is the identity restrictedto the M = 0 sector, and scales as log L [24]. We providesome further numerical details in Appendix A. III. ENTANGLEMENT IN OPEN QUANTUMSYSTEMS
In order to motivate the entanglement quantities weare interested in, let us review some general terminologyabout bipartite mixed state entanglement [36]. Considerthe density matrix ρ of a generic quantum system. If itis possible to decompose the density matrix as a convexsum over a bipartition ρ = (cid:88) i p i ρ Ai ⊗ ρ Bi (5)with (cid:80) i p i = 1 and p i (cid:62)
0, the density matrix is calledseparable, meaning that there are only classical correla-tions between part A and B . If it is not possible to writesuch a decomposition the density matrix is quantum en-tangled as it is impossible to produce it locally. Now thequestion is to find an entanglement measure that is onlysensitive to these quantum correlations. Typical entan-glement monotones in the pure state formalism such asthe von Neumann entanglement entropy and the R´enyientropies S q ( | ψ (cid:105) ) = 11 − q log Tr ρ qB (6)with reduced density matrix ρ B = Tr A | ψ (cid:105) (cid:104) ψ | are notgood measures in the mixed case since they are also sen-sitive to classical correlations.It is a hard problem to determine whether or not adensity matrix is separable over a bipartition as in (5);in fact it has been proven to be NP-hard and no generalsolution is known [37]. However there are some easilycomputable criteria that determine if the density matrixis not separable. They can in some cases distinguish en-tangled and separable states. The most important ofsuch criteria are based on positive maps, or on entangle-ment witnesses [36].The well known criterion of separability of Peres isbased on a positive map (transposition) and investigateswhether the density matrix remains positive under par-tial transposition of a subsystem: if ρ is separable ⇒ ρ T B is positive [38]. The negativity is an entanglement mono-tone that measures the violation of this criterion [39] N ( ρ ) = (cid:107) ρ T B (cid:107) − , (7)in which the trace norm is defined as (cid:107) O (cid:107) = Tr √ OO † ,hence a summation over the absolute values of the eigen-values. The trace is invariant under partial transposition,i.e. 1 = Tr ρ = Tr ρ T B , in which we assumed a normalizeddensity matrix. Therefore N ( ρ ) = (cid:88) λ i < | λ i | , (8)is effectively just a summation over the negative eigenval-ues introduced in the density matrix under partial trans-position. We can also define the logarithmic negativity E ( ρ ) = log (cid:107) ρ T B (cid:107) . The negativity has been previouslystudied in the context of MBL in closed quantum systemsin Ref. [40, 41] and has been experimentally measuredbetween two quibts in Ref. [20]. Note that the nega-tivity dynamics in open quantum systems can be non-asymptotic, however in our setup we avoid this ‘sudden-death dynamics’ by explicit spin conservation as we ex-plain in Appendix B. The negativity is difficult to com-pute using tensor-network approaches. A more accessiblequantity is the R´enyi negativity, much in the spirit of theR´enyi entropy in the pure-state formalism, using a replicaconstruction E q ( ρ ) = log Tr (cid:0) ρ T B (cid:1) q . (9)This replica construction for entanglement negativitieshas been proposed in the context of field theories inRefs. [27, 28, 42]. Unlike the R´enyi entropies for purestates, the R´enyi negativites for mixed states are no en-tanglement monotones. However the moments of thepartially transposed density matrix can be used to es-timate the negativity as shown in Ref. [29]. It isalso easy to see that the analytic continuations of (9)are different for even ( q e ) and odd ( q o ) powers. Wehave that lim q e → E q e ( ρ ) = E ( ρ ) while lim q o → E q o ( ρ ) =log Tr ρ T B = 0 due to the the normalization of the den-sity matrix. For pure states, we can work out the powersof the partially transposed density matrix in terms of thereduced density matrix [28]Tr (cid:0) ρ T B (cid:1) q = (cid:40) Tr ρ q o B , q o odd,(Tr ρ q e / B ) , q e even. (10)By taking the limit q e →
1, we see that the logarithmicnegativity can be linked to the R´enyi entropy of order1 / E ( | ψ (cid:105) ) = 2 log Tr ρ / B = S / . (11)As entanglement proxy we will consider R q ( ρ ) = − log (cid:32) Tr (cid:0) ρ T B (cid:1) q Tr ρ q (cid:33) = log Tr ρ q − E q ( ρ ) , (12)as this quantity remains zero for diagonal density ma-trices. For q = 1 , R q vanishes, such that the first non-trivial quantity is R . This R´enyi negativity gives alwayszero for product states, but is not necessarily zero for allseparable (classically correlated) states, and hence is noentanglement monotone. We work out R for the 2-qubitWerner state as an example in Appendix C. It is howevera computable, and potentially measurable, entanglementestimator. Note that for a pure state | ψ (cid:105) R ( | ψ (cid:105) ) = − log (cid:0) ρ T B (cid:1) = − log Tr ρ B = 2 S . (13) M l − M l M l M l M M M M M M M l − M l − M l − M l M l M l − log ········· ········· M M M M M M A B M l − M l − FIG. 2. A sketch of the R entanglement measure, see Eq.(12), that we compute using tensor-network techniques. Wemake a bipartition of the system into two subsystems A and B , and partially transpose the degrees of freedom of subsys-tem B before taking the trace. R has been previously studied in Ref. [30] in the con-text of finite-temperature phase transitions. In our con-text, the N´eel state at t = 0 and the maximally mixedstate at t = + ∞ have a value of R = 0 as their densitymatrices do not have any negative values. At intermedi-ate times, when there is some entanglement, the trace ofthe partially transposed density matrix is reduced mean-ing that tr ( ρ TB ) tr ρ <
1, and thus R >
0. We sketchhow R is computed using tensor-network techniques inFig. 2. So basically we need to compute 1 . A is used in the same way inthe numerator and denominator.Another possibility to quantify entanglement in openquantum systems is given by entanglement witnesses.They have the advantage that they can be experimen-tally relevant, because often they rely on simple expecta-tion values. Their disadvantage is that the most optimalwitness in general requires an optimization over the fullHilbert space, we discuss this further in Appendix D. IV. RESULTSA. Closed system
As we have shown in the previous paragraph R re-duces to the third R´enyi entropy in the closed quantumsystem. To check that it indeed shows the same char-acteristic features as the von Neumann entropy in thethermodynamic limit, we have plotted its behavior inFig. 3. We observe the typical logarithmic growth in R for MBL systems, and a fast saturation for Andersonlocalized systems.As can be seen in Fig. 4, the distributions of R arevery broad, and it is our goal to understand the shapeof the distribution. Therefore we start from a simplifiedLIOM Hamiltonian neglecting couplings between threeor more LIOMs H = (cid:88) i h i τ zi + (cid:88) i>j J ij τ zi τ zj (14)in which the J ij ’s decay exponentially with the distancebetween the spins J ij = J e − r/ξ with r = i − j .Assuming an initial product state of two spins | ψ (0) (cid:105) = √ ( | (cid:105) + | (cid:105) ) ⊗ √ ( | (cid:105) + | (cid:105) ), which are generated for theLIOMs because our initial state is prepared in a productstate for the physical spins. We obtain for the entan-glement generated under time-evolution with Hamilto-nian (14), that R ( t ; J r ) = − log (cid:18) tJ r )8 (cid:19) , (15)hence the maximum R that can be generated betweenthe spins is 2 log 2 as expected.The couplings J r have been shown to be distributedaccording to a log-normal distribution [43] P J ( J ; r, ξ , ξ ) = (cid:114) ξ πr J exp (cid:32) − (log J + 2 r/ξ ) r/ξ (cid:33) . (16)The parameters ξ and ξ characterize respectively thegrowth of the mean and variance with distance betweenthe spins. Then we can estimate the distribution of R for a bipartition of size L by summing over L values of R that are calculated by sampling the J ’s from L dis-tribution functions (16) R ( t ) = L (cid:88) r =1 R ( t ; J r ) . (17)The average and some histograms given by this model arecompared to the numerics in Fig. 4 for L = 16 and ∆ = 1.We have taken the parameters ξ and ξ of the distribu-tion such that the growth of our model has the sameslope ( ξ log t ) as our data, and such that ξ /ξ ≈ R = 2 log 2 both in the model as in the numerics cor-responding to a singlet bound over the bipartition [44].Secondly we note that there are two simplifications inour model (i) the difference between τ -spins and physicalspins and (ii) the fact that we only took into account 2-spin couplings. The first simplification is reflected in theshort-time dynamics. The second simplification inducestoo long tails in the model towards low entanglement, aswe did not take into account multi-spin couplings whichcan also provide significant contributions to the entan-glement.The distribution of R for Anderson localized systemswould decay quickly for values higher than the singletbond R = 2 log 2 as entanglement cannot propagatethrough the system, in contrast to what we observe forthe MBL system. − tJ . . . . . . R ∆ = 0 .
0∆ = 0 .
1∆ = 0 .
5∆ = 1 . tJ . . . R ∆ = 1 . FIG. 3. Quench dynamics of R = 2 S in the closed quantumsystem at fixed system size L = 14 for various interactionstrengths. In the inset we show the finite size scaling of R at fixed interaction strength ∆ = 1, from bottom to top L =10 , , ,
16. The errorbars show the standard error of themean. tJ . . . R (a) modelnumerics − − R )0 . . . . tJ = 103 . − − R )0 . . . . tJ = 547 . − − R )0 . . . . tJ = 1000 . FIG. 4. (a) Comparison of the average of R obtained by thenumerics over 2500 ensembles with system size L = 16 and∆ = 1 and by the model described in the text with ξ = 0 . ξ = 0 . B. Open system
1. Stretched exponential decay of R When we switch on the dephasing, R stays a good en-tanglement measure unlike the von Neumann and R´enyientropies. In the open system, R undergoes a char-acteristic stretched exponential decay starting at timescales (cid:38) / Γ as shown in Fig. 5. Such a decay canbe understood as a superposition of local exponentialdecays, and has also been observed in the imbalancein Refs. [22, 23, 25] and is also experimentally con-firmed [15]. We observed such a decay as well in ourexact simulations for other entanglement measures likethe negativity and the Fisher information as we show inAppendix D.Next, we quantitatively extract the stretching expo-nent b of R ∼ e − (Γ t/a ) b by considering different systemsizes. To this end, we use the TEBD algorithm on densitymatrices, and compute R for the rather large couplingstrength Γ = 0 . J . From Fig. 5 we see that the exponentis around b ≈ .
25, and our data shows that interactionsin the Hamiltonian do not influence this exponent much.This is expected to hold true as long as interactions aresmall compared to the disorder in the system ∆ < W [25].Ref. [22, 23] respectively reported a stretching exponent b ≈ .
38 and b ≈ .
42 for the imbalance. We observean exponent in R that is significantly smaller indicatingthat entanglement is more robust than transport underdephasing.
2. MBL at intermediate time scales
As we have seen in the previous section, interactions inthe system are only a subleading effect in the stretchedexponential tails, and hence these tails do not allow oneto distinguish MBL from Anderson localization. Theentanglement quantity R is however able to distin-guish MBL from Anderson by determining the maximallyreached value of R at intermediate time scales, on thecondition that the dephasing is sufficiently weak and theinteractions are sufficiently strong, i.e. Γ /J < ∆ seeFig. 6. This is because the dephasing typically set in ata time t ∼ / Γ and the logarithmic growth at a time t ∼ (∆ J ) − . We also compare the entanglement dynam-ics to the relaxation dynamics of the imprinted densitypattern, as measured by the imbalance, I = (cid:104) S ze − S zo (cid:105)(cid:104) S ze + S zo (cid:105) ,where S ze /S zo sums S zi over even/odd sites.We also compute the interquartile range (IQR) of ourdata, shown in the insets of Fig. 6, which is a measure forthe spread of a distribution that is less sensitive to thetails than the variance. We choose this measure becauseof the limited number of ensembles that are numericallyfeasible for open quantum systems, which implies thatwe only have limited access to the tails. The dephasingis driving the state into a trivial steady state which im-plies that the distribution of R will converge to a δ -peakat zero entanglement. We indeed see a clear dip in thespread of the distribution as the dephasing sets in. Thefull distributions of R , shown in Fig. 7, possess strongtails even in the presence of dephasing noise, however,their width is decreasing over time.The fact that R carries only traces from MBL at inter-mediate times, when the dephasing is not yet completely − − t Γ10 × × × × − l og ( R ) L = 20 L = 32 L = 36 .
03 0 .
04 0 . /L . . e x p o n e n t b − − t Γ10 × × × × − l og ( R ) L = 20 L = 32 L = 36 .
03 0 .
04 0 . /L . . e x p o n e n t b FIG. 5. Quench dynamics of the R´enyi negativity R in a sys-tem with disorder W = 5 J and coupling Γ = 0 . J . Top: MBLwith interaction strength ∆ = 1. Bottom: Non-interactingsystem ∆ = 0. In the inset we show the best fitting parame-ter for the stretching exponent b with 3 σ errorbars obtainedby the least-square method. The blue line shows one of thefitting functions with b ≈ . dominating the dynamics can also be seen from the dis-tributions: if we want to detect traces of MBL, we needto have some larger entanglement clusters remaining overthe biparition in some ensembles. Thus the distributionof R must contain some part that has an entanglementthat is higher than the singlet entanglement R = 2 log 2,for MBL to be detectable. This criterion is more sensi-tive than just looking at the mean of R , since it focuseson the upper part of its full distribution. V. CONCLUSION
We have discussed a novel entanglement measure foropen quantum systems in the context of many-body lo-calization. We have seen that the third R´enyi negativityforms a promising measure to study the entanglementdynamics of an MBL system that is slightly coupled toa dephasing environment. R can distinguish MBL fromAnderson up to intermediate time scales as it reproducesthe logarithmic growth of entanglement in the clean MBLsystem. In addition, we conclude that all quantities, en- − tJ . . . . . R ∆ = 0Γ = 0 . J ∆ = 1Γ = 0 . J Γ = 0 . J − tJ I Q R ( l og R ) − tJ . . . . . I − tJ . . . I Q R ( I ) FIG. 6. Dynamics of the R´enyi negativity (top) and the im-balance (bottom) under a quench in the open or closed system( L = 20). The disorder is fixed to W = 5 J and the interac-tion strength is ∆ = 1 for the MBL system, and ∆ = 0 for thenon-interacting Anderson localized system. At intermediatetime scales and for sufficiently weak dephasing strength, R distinguishes MBL from single particle localization. A fea-ture that is absent in the imbalance. In the insets we showthe interquartile range (IQR) which forms a measure for thespread of the distribution. tanglement and transport, decay according to a stretchedexponential. However the stretching exponents are foundto be smaller for the entanglement quantities, meaningthat the late time entanglement dynamics is slower thanfor instance the dynamics of the imbalance under dephas-ing.The quantities Tr ρ and Tr (cid:0) ρ T B (cid:1) are measurablewithout the need of full state tomography by performingjoint measurements on n = 3 copies of the state [29, 45–48]. Alternatively, one could link Tr ρ n and Tr (cid:0) ρ T B (cid:1) n tothe statistical correlations of random measurements ona single copy of the state [49], by further developing themeasurement protocols proposed for the R´enyi entangle-ment entropies [17, 50–54].For future work it would be interesting to investigatewhether these novel protocols to measure entanglementin open quantum systems could be potentially experi-mentally as relevant as the protocols to measure R´enyientropies in closed quantum systems. From the theo- − . − . − . . . . . . . tJ = 61 . Γ = 0 J Γ = 0 . J Γ = 0 . J − . − . − . . . . . . . tJ = 394 . − . − . − . . R )0 . . . . . Γ = 0 J Γ = 0 . J − . − . − . . R )0 . . . . . FIG. 7. Histograms at two different times (columns) corre-sponding to the averages shown in Fig. 6. (a)-(b) Interactingsystem with ∆ = 1. (c)-(d) Non-interacting system ∆ = 0.The green line indicates R = 2 log 2 which corresponds tothe maximal entanglement between 2-spins. retical perspective it would be interesting to investigatehow R behaves under different forms of dissipation. Inparticular for non-hermitean types of Lindblad operatorsit would be interesting to investigate which signaturesthe entanglement structure of a non-trivial steady statecontains. VI. ACKNOWLEDGMENTS
We thank A. Elben, S. Gopalakrishnan and T. Groverfor useful discussions. Our tensor-network calcula-tions were performed using the TeNPy Library [55].We acknowledge support from the Technical Univer-sity of Munich - Institute for Advanced Study, fundedby the German Excellence Initiative and the EuropeanUnion FP7 under grant agreement 291763, the DeutscheForschungsgemeinschaft (DFG, German Research Foun-dation) under Germanys Excellence Strategy-EXC-2111-390814868, Research Unit FOR 1807 through grantsno. PO 1370/2-1, DFG TRR80 and DFG grant No.KN1254/1-1, and from the European Research Coun-cil (ERC) under the European Unions Horizon2020 re-search and innovation programme (grant agreements No.771537 and 851161).
Appendix A: Numerical details
The density matrix ρ is represented as a matrix prod-uct operator (MPO) ρ = d (cid:88) i ,...,i L j ,...,j L M i ,j [1] M i ,j [2] . . . M i L ,j L [ L ] | i , . . . , i L (cid:105)(cid:104) j , . . . , j L | = d (cid:88) i ,...,i L M i [1] M i [2] . . . M i L [ L ] ( σ i ⊗ σ i ⊗ . . . ⊗ σ i L L ) , where d is the local dimension of the Hilbert space, andwhere σ ,...,d l is a set of d × d Hermitean and tracelessmatrices (Gell-Mann matrices), and where σ l = . Inour spin-1 / d = 2, these are just the Paulimatrices σ , , l = ( σ xl , σ yl , σ zl ). The matrix M i l l has di-mension χ l − × χ l with max( χ l − , χ l ) ≤ χ , where χ isthe maximal bond dimension. We can impose a canoni-cal form on the MPO in the same way as for the MPS. Weapply the TEBD algorithm on this MPO [33, 34, 39, 56],which means that we make a Suzuki-Trotter decom-position of the time-evolution superoperator exp( L d t )in terms of the two-site time-evolution superoperators U i,i +1 ( δt ) = exp( L i,i +1 δt ) with L i,i +1 = − iH i,i +1 ⊗ + i ⊗ H i,i +1 − (cid:18) Γ i + Γ i +1 (cid:19) ⊗ + Γ i ( S zi ⊗ S zi ) + Γ i +1 (cid:0) S zi +1 ⊗ S zi +1 (cid:1) . (A1)In our simulations we use the usual fourth order Trotterdecomposition scheme which does in principle destroy thecanonical form while performing the updates on all evenor odd bonds, in addition because dissipation makes thetime evolution non-unitairy. However this effect is verysmall for Γ (cid:46) J , therefore we can still use this schemewith very good accuracies. In case the dissipation takesthe form of dephasing noise, the complexity of the statecan be reduced. (Note that the infinite temperature state ∼ can be represented by a MPO with bond dimension χ = 1.) We truncate the singular values after acting with U i,i +1 ( δt ) on a bond by only keeping the χ largest onesor by only keeping the ones that are larger than a certain (cid:15) trunc . Note that, unlike for the pure state case, this doesnot fully correspond to truncating in the entanglementof the density matrix, but rather in its complexity, orso-called operator-space entropy [35].In the remainder of this section we show a compari-son between various parameters of the TEBD on MPOalgorithm. We show results for one particular disorderensemble such that we can compare the errors caused bythe algorithm, without the statistical errors from the av-eraging. In Fig. 8 we compare the fourth order TEBDalgorithm with the exact results, by showing the rela-tive error in the R´enyi negativity. As expected this errorincreases in time and with decreasing bond dimension.The main source of error that declares the small devia-tions from the exact result at maximal bond dimension − − − t Γ10 − − − − − − | ∆ ( R ) | / R χ = 256 χ = 200 χ = 100 10 − t Γ10 − − − | ∆ ( R ) | / R χ = 256 χ = 200 χ = 100 χ = 5010 − − − t Γ10 − − − − − − | ∆ ( R ) | / R dt = 0 . dt = 0 . dt = 0 . − t Γ10 − − − − − − | ∆ ( R ) | / R dt = 0 . dt = 0 . dt = 0 . FIG. 8. Relative error of the fourth order TEBD scheme forvarious simulation parameters for a small system of L = 8spins. In the left column the coupling is Γ = 0 . J , in theright one Γ = 1 J . In the top line we vary the maximal bonddimension. In the bottom line we vary the time step (in unitsof J − ) at maximal bond dimension χ = 256. − − − t Γ10 − − − − − − − | R χ − R χ m a x | / R χ m a x χ = 200 χ = 300 χ = 400 10 − − − t Γ10 − − − − | R d t = . − R d t = . | / R d t = . χ = 200 χ = 300 χ = 400 χ = 500 FIG. 9. Relative error of the fourth order TEBD scheme for L = 40 spins with a coupling strength of Γ = 0 . J . On theleft we show the relative difference between simulations withbond dimension χ and bond dimension χ max = 500 with atime step of d t = 0 .
05. On the right we show the differencebetween simulation with time step is d t = 0 . t = 0 . χ . is the Trotter error because of the splitting of the time-evolution operator. However this error can be controlledby choosing a small enough time step, as can be deducedfrom Fig. 8 where we also plotted the performance of theTEBD scheme at various time steps at the exact bonddimension.In Fig. 9 we show the relative error with respect tothe largest bond dimension that was easily computablefor a system size of L = 40, as well as a comparisonbetween different time steps. From this we conclude thatwe maximally need a bond dimension around χ = 400,and time step d t = 0 . Appendix B: Sudden death dynamics of thenegativity
Entanglement quantities in open quantum systemsmay decay non-asymptotically, unlike transport quanti-ties. This so-called sudden death dynamics is a knownphenomenon, that imposes challenges on the stability of − − − t Γ10 − − − E Γ = 0 . J Γ = 0 . J Γ = 0 . J − − − t Γ10 − − − E Γ = 0 . J Γ = 0 . J Γ = 0 . J FIG. 10. Negativity dynamics in a small system of L = 8spins. When spin-conservation is broken by adding a term g (cid:80) i S xi to the Hamiltonian, the negativity dynamics stopsabruptly at a finite time. In the left plot g = 0 . J , and in theright one g = 0 . J . quantum memories [57, 58]. In our setup this specific dy-namics only occurs in the negativity when we explicitlybreak the spin-conservation symmetry as illustrated inFig. 10. In this section we investigate why the negativitydecay is always asymptotic when the evolution conservesthe total spin. The U (1)-symmetry leads to entries of thedensity matrix that are always zero, only one sub-blockof the density matrix that corresponds to the consideredspin-sector is occupied. In the M = (cid:80) i (cid:104) S zi (cid:105) = 0 sector,for a chain with an even number of spins, the dimensionis m = C LL/ . Partial transposition maps at least partof the off-diagonal elements of the occupied sub-block toother spin sectors. Consider for instance the two qubitmatrix element in the M = 0 sector c | (cid:105) (cid:104) | , after par-tially transposing the second qubit index this becomes c | (cid:105) (cid:104) | , which is a matrix element outside the M = 0sector. Clearly under spin-conserving dynamics this ma-trix element would have remained zero.As the diagonal elements remain of course invariant un-der partial transposition, we can split ρ T B into two blocks B out and B in , corresponding to occupied elements insideor outside the original spin sector ρ T B = B out ⊕ B in (B1)with the blocks of the generic form, B out = (cid:18) AA † (cid:19) ∈ C n, n and B in = ( B in ) † ∈ C m,m (B2)with m + 2 n = dim( H ), because of the Hermiticity of theoriginal density matrix. From this simple argument wecan make no a priori assumptions about the structure ofthe eigenvalues of B in . However, it is easy to see that fora density matrix of the form B out the eigenvalues comein pairs with opposite signs ± λ , ± λ , · · · ± λ n . The factthat there are always negative eigenvalues present due toinherent structure of the partially transposed density ma-trix of a system with spin conservation, prevents suddendeath dynamics in the negativity. . . . λ . . . . E . . . λ . . . R FIG. 11. The logarithmic negativity E and the third R´enyinegativity R for the two-qubit Werner state defined by λ . R is not an entanglement monotone as it takes a non-zerovalue in the separable regime λ < / Appendix C: R is not an entanglement monotone By the partial transposition criterion of Peres, it fol-lows that each separable state has a positive partialtranspose. Therefore each separable state has negativ-ity E ( ρ sep ) = 0, however this is not true for the R´enyinegativity. Consider for instance the two-qubit Wernerstate ρ ( λ ) = λ | φ (cid:105) (cid:104) φ | + (1 − λ ) I with λ ∈ [0 ,
1] and | φ (cid:105) aBell pair, in the positive partial transpose (PPT) regime λ < /
3. In this regime ρ T B has only positive eigenvalues,and as the PPT criterion is a sufficient for separability inthe two qubit case, ρ ( λ < /
3) is separable. However theeigenvalues of ρ T B and ρ are not the same, which implies anon-trivial value of R . Explicitly the eigenvalues of ρ ( λ )are { × (1 − λ ) , × (3 λ + 1) } , while the eigenvalues of ρ T B ( λ ) are { × (1+ λ ) , × (1 − λ ) } . So in this case R takes a non-trivial value, while E = 0. This is illustratedin Fig. 11. Entanglement monotones satisfy invarianceunder LOCC (local operations and classical communica-tion). However a separable state is transformable intoany other separable by means of LOCC. (Note that lo-cal unitairy transformations fall into the class of LOCCtransformations.) Therefore an entanglement monotonemust remain constant over the set of separable states,which is clearly not the case for the R´enyi negativity inour example. Appendix D: Quantum Fisher information asentanglement witness
Here we discuss the quantum Fisher information (QFI)which quantifies the sensitivity of a state to a unitairytransformation e iθO generated by a linear Hermitean op-erator of the form O = (cid:80) i n i · S i , where n i is a unitvector and S i is the vector of spin matrices ( S xi , S yi , S zi ).Therefore it measures the spread of quantum correlationsvia the operator O . The QFI witnesses entanglement ina state if its value is larger than the system size F Q > L ,and by other conditions it can also witness multipartiteentanglement [59]. For pure states, the QFI is given by − − − t Γ10 − − l og f Q a: 0.10 ± ± . J Γ = 0 . J Γ = 0 . J − t Γ05 f Q − − − t Γ10 − l og E a: 0.17 ± ± . J Γ = 0 . J Γ = 0 . J − t Γ0 . . . E − − − t Γ10 − l og R a: 0.02 ± ± . J Γ = 0 . J Γ = 0 . J − t Γ0 . . . R FIG. 12. Lindblad quench dynamics in the XXZ model atdisorder W = 5 J , interaction strength ∆ = 1 and system size L = 8. Stretched exponential fits e − ( Γ ta ) b for the decay thatstart approximately at time ∼ fit the data very well. Weshow from top the bottom: the Fisher information density f Q = F Q /L , the negativity and the third R´enyi negativity.The insets show the same data in a different scale. The er-rorbars show the standard error of the mean upon averagingover around 1000 disorder ensembles. the variance of OF Q ( | ψ (cid:105) , O ) = 4 (cid:16) (cid:104) ψ | OO | ψ (cid:105) − | (cid:104) ψ | O | ψ (cid:105)| (cid:17) . (D1)For mixed states the QFI cannot be related to simple ex-pectation values, instead the full spectral decompositionof the density matrix ρ = (cid:80) i p i | s i (cid:105) (cid:104) s i | is necessary [60] F Q ( ρ, O ) = 2 (cid:88) i,jp i + p j > ( p i − p j ) p i + p j | (cid:104) s j | O | s i (cid:105)| , (D2)and can only be computed using exact diagonalization,unless ρ takes the form of a thermal state [61]. The QFIrelies on the choice of generator O , and for simplicity wewill choose the staggered magnetization O = (cid:80) i ( − i S zi which seems a natural choice to consider the quenchdynamics from an initial N´eel state. Note that thechoice O = (cid:80) i S zi would imply a vanishing Fisherinformation due to spin conservation, while O = (cid:80) i S xi would imply that the Fisher information is equal to thesystem size for the N´eel state F ( t =0) Q = L and underdephasing dynamics again converges to the systemsize F ( t → + ∞ ) Q = L . The QFI has been experimentallymeasured in the context of MBL in Ref. [14].In Fig. 12 we have computed the QFI, the negativityand R by exact calculations. From this we see that theQFI also decays according to a stretched exponential, andthat the stretching exponents of the negativity and R are indeed approximately equal. [1] D. Basko, I. Aleiner, and B. Altshuler, Metalinsulatortransition in a weakly interacting many-electron systemwith localized single-particle states, Annals of Physics , 1126 (2006).[2] R. Vosk and E. Altman, Many-body localization in onedimension as a dynamical renormalization group fixedpoint, Phys. Rev. Lett. , 067204 (2013).[3] R. Nandkishore and D. A. Huse, Many-body lo-calization and thermalization in quantum statisti-cal mechanics, Annual Review of Condensed MatterPhysics , 15 (2015), https://doi.org/10.1146/annurev-conmatphys-031214-014726.[4] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Col-loquium: Many-body localization, thermalization, andentanglement, Rev. Mod. Phys. , 021001 (2019).[5] P. W. Anderson, Absence of diffusion in certain random lattices, Phys. Rev. , 1492 (1958).[6] S. Gopalakrishnan, M. M¨uller, V. Khemani, M. Knap,E. Demler, and D. A. Huse, Low-frequency conductivityin many-body localized systems, Phys. Rev. B , 104202(2015).[7] K. Agarwal, E. Altman, E. Demler, S. Gopalakrishnan,D. A. Huse, and M. Knap, Rare-region effects and dy-namics near the many-body localization transition, An-nalen der Physik , 1600326 (2017).[8] M. ˇZnidariˇc, T. c. v. Prosen, and P. Prelovˇsek, Many-body localization in the heisenberg xxz magnet in a ran-dom field, Phys. Rev. B , 064426 (2008).[9] J. H. Bardarson, F. Pollmann, and J. E. Moore, Un-bounded growth of entanglement in models of many-bodylocalization, Phys. Rev. Lett. , 017202 (2012).[10] M. Serbyn, Z. Papi´c, and D. A. Abanin, Universal slow growth of entanglement in interacting strongly disorderedsystems, Phys. Rev. Lett. , 260601 (2013).[11] M. Serbyn, Z. Papic, and D. Abanin, Local conservationlaws and the structure of the many-body localized states,Phys. Rev. Lett. , 127201 (2013).[12] D. A. Huse, R. Nandkishore, and V. Oganesyan, Phe-nomenology of fully many-body-localized systems, Phys.Rev. B , 174202 (2014).[13] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. L¨uschen,M. H. Fischer, R. Vosk, E. Altman, U. Schneider, andI. Bloch, Observation of many-body localization of inter-acting fermions in a quasirandom optical lattice, Science , 842 (2015).[14] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. Hess,P. Hauke, M. Heyl, D. Huse, and C. Monroe, Many-bodylocalization in a quantum simulator with programmablerandom disorder, Nature Physics , 10.1038/nphys3783(2015).[15] H. P. L¨uschen, P. Bordia, S. S. Hodgman, M. Schreiber,S. Sarkar, A. J. Daley, M. H. Fischer, E. Altman, I. Bloch,and U. Schneider, Signatures of many-body localizationin a controlled open quantum system, Phys. Rev. X ,011034 (2017).[16] P. Bordia, H. L¨uschen, S. Scherg, S. Gopalakrishnan,M. Knap, U. Schneider, and I. Bloch, Probing slow re-laxation and many-body localization in two-dimensionalquasiperiodic systems, Phys. Rev. X , 041047 (2017).[17] T. Brydges, A. Elben, P. Jurcevic, B. Vermersch,C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F.Roos, Probing r´enyi entanglement entropy via random-ized measurements, Science , 260 (2019).[18] M. Rispoli, A. Lukin, R. Schittko, S. Kim, M. E. Tai,J. L´eonard, and M. Greiner, Quantum critical behaviourat the many-body localization transition, Nature ,385 (2019).[19] A. Lukin, M. Rispoli, R. Schittko, M. E. Tai, A. M. Kauf-man, S. Choi, V. Khemani, J. L´eonard, and M. Greiner,Probing entanglement in a many-body–localized system,Science , 256 (2019).[20] B. Chiaro, C. Neill, A. Bohrdt, M. Filippone, F. Arute,K. Arya, R. Babbush, D. Bacon, J. Bardin, R. Barends,S. Boixo, D. Buell, B. Burkett, Y. Chen, Z. Chen,R. Collins, A. Dunsworth, E. Farhi, A. Fowler, B. Foxen,C. Gidney, M. Giustina, M. Harrigan, T. Huang,S. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi,J. Kelly, P. Klimov, A. Korotkov, F. Kostritsa, D. Land-huis, E. Lucero, J. McClean, X. Mi, A. Megrant,M. Mohseni, J. Mutus, M. McEwen, O. Naaman, M. Nee-ley, M. Niu, A. Petukhov, C. Quintana, N. Rubin,D. Sank, K. Satzinger, A. Vainsencher, T. White,Z. Yao, P. Yeh, A. Zalcman, V. Smelyanskiy, H. Neven,S. Gopalakrishnan, D. Abanin, M. Knap, J. Martinis, andP. Roushan, Growth and preservation of entanglement ina many-body localized system, 1910.06024v1.[21] A. Carmele, M. Heyl, C. Kraus, and M. Dalmonte,Stretched exponential decay of majorana edge modesin many-body localized kitaev chains under dissipation,Phys. Rev. B , 195107 (2015).[22] M. H. Fischer, M. Maksymenko, and E. Altman, Dynam-ics of a many-body-localized system coupled to a bath,Phys. Rev. Lett. , 160401 (2016).[23] E. Levi, M. Heyl, I. Lesanovsky, and J. P. Garrahan,Robustness of many-body localization in the presence ofdissipation, Phys. Rev. Lett. , 237203 (2016). [24] M. V. Medvedyeva, T. c. v. Prosen, and M. ˇZnidariˇc,Influence of dephasing on many-body localization, Phys.Rev. B , 094205 (2016).[25] B. Everest, I. Lesanovsky, J. P. Garrahan, and E. Levi,Role of interactions in a dissipative many-body localizedsystem, Phys. Rev. B , 024310 (2017).[26] I. Vakulchyk, I. Yusipov, M. Ivanchenko, S. Flach, andS. Denisov, Signatures of many-body localization insteady states of open quantum systems, Phys. Rev. B , 020202 (2018).[27] P. Calabrese, J. Cardy, and E. Tonni, Entanglement neg-ativity in quantum field theory, Phys. Rev. Lett. ,130502 (2012).[28] P. Calabrese, J. Cardy, and E. Tonni, Entanglement neg-ativity in extended systems: a field theoretical approach,Journal of Statistical Mechanics: Theory and Experi-ment , P02008 (2013).[29] J. Gray, L. Banchi, A. Bayat, and S. Bose, Machine-learning-assisted many-body entanglement measure-ment, Phys. Rev. Lett. , 150503 (2018).[30] K.-H. Wu, T.-C. Lu, C.-M. Chung, Y.-J. Kao, andT. Grover, Entanglement renyi negativity across a finitetemperature transition: a monte carlo study, (2019),1912.03313.[31] A. Pal and D. A. Huse, Many-body localization phasetransition, Phys. Rev. B , 174411 (2010).[32] H.-P. Breuer and F. Petruccione, The Theory of OpenQuantum Systems (Oxford University Press, 2007).[33] F. Verstraete, J. J. Garc´ıa-Ripoll, and J. I. Cirac, Ma-trix product density operators: Simulation of finite-temperature and dissipative systems, Phys. Rev. Lett. , 207204 (2004).[34] M. Zwolak and G. Vidal, Mixed-state dynamics in one-dimensional quantum lattice systems: A time-dependentsuperoperator renormalization algorithm, Phys. Rev.Lett. , 207205 (2004).[35] T. Prosen and I. Piˇzorn, Operator space entanglemententropy in a transverse ising chain, Phys. Rev. A ,032316 (2007).[36] O. G¨uhne and G. Toth, Entanglement detection10.1016/j.physrep.2009.02.004 (2008), 0811.2803v3.[37] L. Gurvits, Classical deterministic complexity of ed-monds’ problem and quantum entanglement, quant-ph/0303055v1.[38] A. Peres, Separability criterion for density matrices,Phys. Rev. Lett. , 1413 (1996).[39] G. Vidal and R. F. Werner, Computable measure of en-tanglement, Phys. Rev. A , 032314 (2002).[40] J. Gray, A. Bayat, A. Pal, and S. Bose, Scale invariantentanglement negativity at the many-body localizationtransition, (2019), 1908.02761v1.[41] C. G. West and T.-C. Wei, Global and short-range entan-glement properties in excited, many-body localized spinchains, (2018), 1809.04689v1.[42] M. Rangamani and M. Rota, Comments on en-tanglement negativity in holographic field theories10.1007/JHEP10(2014)060, 1406.6989v2.[43] V. K. Varma, A. Raj, S. Gopalakrishnan, V. Oganesyan,and D. Pekker, Length scales in the many-body localizedphase and their spectral signatures, Phys. Rev. B ,115136 (2019).[44] R. Singh, J. H. Bardarson, and F. Pollmann, Signaturesof the many-body localization transition in the dynamicsof entanglement and bipartite fluctuations, New Journal of Physics , 023046 (2016).[45] J. Cai and W. Song, Novel schemes for directly measur-ing entanglement of general states, Phys. Rev. Lett. ,190503 (2008).[46] H. A. Carteret, Noiseless quantum circuits for the peresseparability criterion, Phys. Rev. Lett. , 040502 (2005).[47] F. Mintert, M. Ku´s, and A. Buchleitner, Concurrence ofmixed multipartite quantum states, Phys. Rev. Lett. ,260502 (2005).[48] A. J. Daley, H. Pichler, J. Schachenmayer, and P. Zoller,Measuring entanglement growth in quench dynamics ofbosons in an optical lattice, Phys. Rev. Lett. , 020505(2012).[49] Y. Zhou, P. Zeng, and Z. Liu, Single-copies estimation ofentanglement negativity, (2020), 2004.11360.[50] S. J. van Enk and C. W. J. Beenakker, Measuring Tr ρ n on single copies of ρ using random measurements, Phys.Rev. Lett. , 110503 (2012).[51] Y. Nakata, C. Hirche, M. Koashi, and A. Winter, Ef-ficient quantum pseudorandomness with nearly time-independent hamiltonian dynamics, Phys. Rev. X ,021006 (2017).[52] A. Elben, B. Vermersch, M. Dalmonte, J. I. Cirac,and P. Zoller, R´enyi entropies from random quenches inatomic hubbard and spin models, Phys. Rev. Lett. ,050406 (2018).[53] B. Vermersch, A. Elben, M. Dalmonte, J. I. Cirac, andP. Zoller, Unitary n -designs via random quenches inatomic hubbard and spin models: Application to themeasurement of r´enyi entropies, Phys. Rev. A , 023604(2018).[54] A. Elben, B. Vermersch, C. F. Roos, and P. Zoller, Sta-tistical correlations between locally randomized measure-ments: A toolbox for probing entanglement in many-body quantum states, Phys. Rev. A , 052323 (2019).[55] J. Hauschild and F. Pollmann, Efficient numericalsimulations with Tensor Networks: Tensor NetworkPython (TeNPy), SciPost Phys. Lect. Notes , 5 (2018),code available from https://github.com/tenpy/tenpy ,arXiv:1805.00055.[56] F. Verstraete, J. I. Cirac, and V. Murg, Matrix prod-uct states, projected entangled pair states, and varia-tional renormalization group methods for quantum spinsystems 10.1080/14789940801912366, 0907.2796v1.[57] M. P. Almeida, F. de Melo, M. Hor-Meyll, A. Salles,S. P. Walborn, P. H. S. Ribeiro, and L. Davidovich,Environment-induced sudden death of entanglement, Sci-ence , 579 (2007).[58] T. Yu and J. H. Eberly, Sudden death of entanglement,Science , 598 (2009).[59] P. Hyllus, W. Laskowski, R. Krischek, C. Schwem-mer, W. Wieczorek, H. Weinfurter, L. Pezz´e, andA. Smerzi, Fisher information and multiparticle entan-glement, Phys. Rev. A , 022321 (2012).[60] S. L. Braunstein and C. M. Caves, Statistical distanceand the geometry of quantum states, Phys. Rev. Lett. , 3439 (1994).[61] P. Hauke, M. Heyl, L. Tagliacozzo, and P. Zoller, Mea-suring multipartite entanglement through dynamic sus-ceptibilities, Nature Physics12