Entanglement and thermalization in open fermion systems
aa r X i v : . [ c ond - m a t . m e s - h a ll ] D ec Entanglement and thermalization in open fermion systems
Cristian Zanoci ∗ and Brian G. Swingle † Stanford Institute for Theoretical Physics, Stanford University, Stanford CA 94305 Department of Physics, Harvard University, Cambridge MA 02138 Martin Fisher School of Physics, Brandeis University, Waltham MA 02453December 16, 2016
Abstract
We numerically study two non-interacting fermion models, a quantum wire model and a Cherninsulator model, governed by open system Linblad dynamics. The physical setup consists of a unitarilyevolving “bulk” coupled via its boundaries to two dissipative “leads”. The open system dynamics ischosen to drive the leads to thermal equilibrium, and by choosing different temperatures and chemicalpotentials for the two leads we may drive the bulk into a non-equilibrium current carrying steady state.We report two main results in this context. First, we show that for an appropriate choice of dynamicsof the leads, the bulk state is also driven to thermal equilibrium even though the open system dynamicsdoes not act directly on it. Second, we show that the steady state which emerges at late time, even in thepresence of currents, is lightly entangled in the sense of having small mutual information and conditionalmutual information for appropriate regions. We also report some results for the rate of approach to thesteady state. These results have bearing on recent attempts to formulate a numerically tractable methodto compute currents in strongly interacting models; specifically, they are relevant for the problem ofdesigning simple leads that can drive a target system into thermal equilibrium at low temperature.
One of the most challenging problems in quantum many-body physics is the calculation of electrical andthermal currents in strongly interacting systems [1–5]. The formalism of linear response provides a generalmethod to compute currents in the limit of weak bias, but the needed correlation functions involve realtime computations at finite temperature. When the coupling is strong these computations are typicallyintractable or rely on uncontrolled approximations. A variety of other methods can be brought to bear uponthe problem, ranging from numerical approaches to the use of AdS/CFT duality [6–16], but calculatingtransport properties of strongly interacting systems generally remains challenging. Given the wealth ofexperimental data on electrical and thermal currents in quantum many-body systems, it is important toaddress this ongoing challenge. We are particularly interested in low temperature dynamics, where collectivemodes can dominate the physics, e.g. at quantum critical points.While a general approach to the problem of transport is difficult to imagine, recent ideas from the theoryof quantum entanglement have offered some new hope in this direction. The starting point of this approachis to directly consider the physical properties of the state ρ NESS corresponding to the non-equilibrium steadystate (NESS) of the system of interest (this is a starting point of many approaches, e.g. [1, 2, 11, 17–24]).The steady state is imagined to carry a current at some finite temperature. Then while a general densitymatrix ρ NESS can be exponentially complex, i.e. is a 2 n × n matrix for n spin-1/2s, it might be that theNESS steady state has relatively little entanglement, similar to a thermal state [25,26], and can be effectivelycompressed to a much smaller set of physically meaningful data. Working in the context of one-dimensional ∗ [email protected] † [email protected] L such that ρ NESS isequal to the time-independent steady state of L , L ( ρ NESS ) = 0. The parameters defining the NESS, e.g. thetemperature gradient, would be built into the dynamics L . One could then find the NESS by simulatingthe time evolution ∂ t ρ = L ( ρ ) until it converges to the steady state. Of course, simulating the dynamics of L may be challenging in general, but simulating the dynamics within a suitable class of low entanglementtensor network states could be feasible [27, 30, 63–65].Ref. [27] found that a small bath consisting of just a few sites was sufficient to obtain interesting physicsat high temperatures (compared to microscopic scales). It was also verified that the details of the bath didnot strongly effect the results, again at high temperature. However, since we are particularly interested inlow temperature physics, the problem of designing an open system dynamics whose steady state is the desiredNESS is potentially non-trivial. To form a useful component of any computational method, such an opensystem dynamics must have three properties: (1) it must be able to drive the system to thermal equilibriumat low temperatures, (2) it must drive the system into the steady state in a reasonable (non-exponential)amount of time, and (3) it must not be excessively complex, e.g. it should not use the detailed properties ofmany-body energy eigenstates. In this work we investigate the problem of designing a suitable dynamics L using a non-interacting fermion model.The dynamical system has two components, illustrated in Figure 1: one is the system of interest (the“target”) and the other one is the designer leads which thermalize the target (or more generally drive itout of equilibrium). In this work both the target and the lead are non-interacting fermion systems. Wepropose that when studying an interacting system, the target should include the interactions of the systemof interest, but the lead can still be taken to be non-interacting. Based on the general expectation thatinteracting systems which thermalize do so independent of the details of the bath, a non-interacting leadshould be adequate to induce thermalization so long as it has the correct gross features. Taking the lead tobe non-interacting immediately answers point (3) above, since such a lead, as we review below, is relativelysimple. In the remainder of this work, we study points (1) and (2) in both one- and two-dimensional non-interacting fermion models. The basic parameter of our lead model is the size of the lead; such a lead withmany sites we call an extended lead.Our results are as follows. First, thermal equilibrium at relatively low temperatures can be reached usingsuch an extended lead, but the required lead size typically grows with decreasing temperature. This addressespoint (1) above. Second, the time to reach the steady state in a metallic state is an inverse polynomial inthe system-plus-lead size. This addresses point (2) above. Third, the above two conclusions are modifiedwhen the system is insulating, either due to an energy gap or due to disorder. In these cases we find that thetime to reach the steady state is significantly longer, scaling exponentially with the system size. Fourth, thesteady states in our extended lead model are low entanglement states and exhibit approximate conditional2 L , μ L β R , μ R L W R
Figure 1: Schematic representation of the setup in which the total system is divided into two parts, the“target” wire in the middle and two “leads” on the left and right. The wire and the leads are alwaysarranged in this quasi-one-dimensional geometry. In the quantum wire model the wire and leads are strictlyone-dimensional, but the Chern insulator model the wire and leads are quasi-one-dimensional strips of atwo-dimensional system.independence for appropriate regions. To the best of our knowledge, we study lower temperatures and largersystems than previously considered in the literature. Compared to the results of Ref. [59], we study in moredetail the structure of the steady state, both in the unbiased and biased case, and we study a two-dimensionalmodel to investigate the validity of our conclusions in more than one dimension.In the remainder of this introductory section we introduce our models, define our observables, and discussour methods. We then present two sets of results, one set for the one-dimensional quantum wire model andone set for the two-dimensional Chern insulator model [66]. We conclude with a brief discussion of futurework and open questions.
We study the transport and entanglement properties of two models. The setup for each model consists ofthree disjoint regions: the left lead (L), the middle part representing the wire (W) whose transport propertieswe investigate, and the right lead (R) (see Figure 1). The three parts consist of N L , N W , and N R fermionsites respectively. The leads are held at inverse temperatures β L = 1 /T L , β R = 1 /T R and chemical potential µ L , µ R . This is accomplished using open system dynamics which drives the decoupled leads to thermalequilibrium with the indicated parameters. We first describe the Hamiltonian dynamics and then discussthe implementation of dissipation. The first model is a one-dimensional quantum wire of fermions. Each of thesegments described above is characterized by a hopping Hamiltonian which is quadratic in fermion creationand annihilation operators H a = − w N a − X x =1 ( c † x c x +1 + c † x +1 c x ) , (1.1)where a = L, W, R , and c x is the annihilation operator at site x . The creation and annihilation operatorsobey the standard anti-commutation relation { c x , c † y } = δ x,y . The leads couple to the middle part of thewire through similar hopping terms: − w ′ ( c † N L c N L +1 + h.c.) and − w ′ ( c † N L + N W c N L + N W +1 + h.c.). One canalso add an on-site potential term V x c † x c x for every fermion in the middle region W. Chern insulator model:
The second model is a two-band tight-binding Hamiltonian in two dimensionswhich exhibits the physics of the quantum anomalous Hall effect. This model is called the Chern insulatormodel and was first studied in [66]. We consider two states of fermions at each lattice site, which can beinterpreted as s and p orbitals. The Hamiltonian in momentum space k = ( k x , k y ) is given by H W = X k [ ǫ ( k ) + V h ( k )] , (1.2)3here ǫ ( k ) and h ( k ) are 2 × ǫ ( k ) = (cid:20) − t cos k x − t cos k y − t cos k x − t cos k y (cid:21) (1.3) h ( k ) = (cid:20) c (2 − cos k x − cos k y − e s ) sin k y + i sin k x sin k y − i sin k x − c (2 − cos k x − cos k y − e s ) (cid:21) (1.4)We can perform a Fourier transform and write the Hamiltonian in position space H W = N W,x X x =1 N W,y − X y =1 c † x,y − (cid:18) V c t (cid:19) iV iV (cid:18) V c − t (cid:19) c x,y +1 + h.c.+ N W,x − X x =1 N W,y X y =1 c † x,y − (cid:18) V c t (cid:19) − V V (cid:18) V c − t (cid:19) c x +1 ,y + h.c.+ N W,x X x =1 N W,y X y =1 c † x,y (cid:20) V c (2 − e s ) 00 − V c (2 − e s ) (cid:21) c x,y , (1.5)where N W,x , N W,y are the dimensions of the lattice, c x,y = [ c x,y, , c x,y, ] T is a two-component columnvector containing the annihilation operators of the two fermion states at lattice position ( x, y ). We treat α = ( x, y, i ) as a composite index labeling the fermions, which obey the anti-commutation relation { c α , c † β } = δ α,β = δ x,x ′ δ y,y ′ δ i,i ′ . Note that the Hamiltonian contains only nearest-neighbor couplings and there are noperiodic boundary conditions imposed.An important feature of the above Hamiltonian is the existence of edge states. If we rewrite H W in the( x, k y ) space, with periodic boundary conditions in y-direction and open boundary conditions in x-direction,we obtain a two-band energy spectrum with two edge states.The lead Hamiltonians are also given by Eq. (1.2), where we set V = 0 and t = w . The resulting leadHamiltonian is similar to the 1 D lead Hamiltonian because it involves only simple nearest neighbor hopping.The coupling between the left lead and the Chern insulator is given by H LW = N L,y X y =1 c † N L,x ,y (cid:20) − w ′ − w ′ − w ′ − w ′ (cid:21) c N L,x +1 ,y + h.c. (1.6)An analogous term H W R can be written for the coupling between the Chern insulator and the right lead.
We describe the interaction of our leads with an environment using Markovian open system dynamics. Thetime evolution of the density matrix ρ of the system is given by Lindblad’s equation dρdt = L ( ρ ) = − i [ H, ρ ] + X j L j ρL † j − X j { L † j L j , ρ } , (1.7)where H is the full system Hamiltonian and L j are jump operators describing the coupling to theenvironment. The operator L is called the Liouvillean operator. The jump operators are taken to belinear in c α and c † α and therefore can be written as L j = u † j c, or (1.8) L j = c † v j , (1.9)4here u j , v j , c = [ c , . . . , c n ] T are column vectors. In the case of a non-interacting fermion system withbath operators linear in the fermion modes, the Lindblad equation reduces to a single particle equation forthe Green’s function [67], as reviewed in Appendix A of [59].We construct the jump operators so that, absent a coupling to the wire, each lead would be driven tothe thermal equilibrium state appropriate to its decoupled lead Hamiltonian. The precise construction isdiscussed in Appendix B of [59]; we review it here. The right and left jump operators are computed in thesame way, so in what follows we focus only on the jump operators for the left bath. Let H L be the singleparticle Hamiltonian of the left bath, and let ǫ j and ψ j be its eigenvalues and eigenvectors. The eigenvectorsare written in the { c α } basis and can be viewed as a function of lattice site. Then we will have N L in rates v L,j = √ γ in ,j ψ j and N L out rates u L,j = √ γ out ,j ψ j , where ψ j is treated as a column vector. The onlyconstraint on the rates is that γ in ,j e β L ( ǫ j − µ L ) = γ out ,j . (1.10)Since only the ratio of the two matters, we can set γ in ,j = γ for all j = 1 , , . . . , N L . Therefore we getthe jump operators L out L,j = q γe β L ( ǫ j − µ L ) ψ † j c, (1.11) L in L,j = √ γc † ψ j . (1.12)Similar formulas are obtained for the right bath. For each of the two models defined above, we compute several observables that reveal the transport andentanglement properties of non-equilibrium steady states. We study electrical currents to probe out-of-equilibrium physics, energy occupation numbers to probe thermalization, decay rates of the open systemdynamics to probe timescales, and mutual information and conditional mutual information to probe theentanglement structure of the state.
Currents:
For the one-dimensional wire, the current operator through the link ( j, j + 1) is I = − i ( c † j c j +1 − c † j +1 c j ) . (1.13)For the two-dimensional lattice, we define the current flowing in the x -direction across the link (( x, y ) , ( x +1 , y )) I x = ic † x,y − (cid:18) V c t (cid:19) − V V (cid:18) V c − t (cid:19) c x +1 ,y + h.c. (1.14)and the current flowing in the y -direction across the link (( x, y ) , ( x, y + 1)) I y = ic † x,y − (cid:18) V c t (cid:19) iV iV (cid:18) V c − t (cid:19) c x,y +1 + h.c. (1.15) Occupation Numbers:
Next we are interested in computing the occupation numbers of the energyeigenstates of the bulk in NESS and comparing them to the thermal equilibrium distribution. Write thebulk Hamiltonian H W as a N W × N W matrix in { c j } basis. Let ǫ k and ψ k be the eigenvalues (single-particleenergies) and eigenvectors of H W . Then, for each energy mode k , we can define an annihilation operator c ǫ k = N W X j =1 ψ ∗ k ( j ) c j , (1.16)5 L , μ L β R , μ R L W RA B C
Figure 2: The partition of the system into three disjoint regions A , B , and C to compute the mutualinformation M I ( A : C ) and the conditional mutual information CM I ( A : C | B ).where the sum is over all the fermion modes in the bulk wire. The number operator for energy mode k is given by h c † ǫ k c ǫ k i = N W X i =1 N W X j =1 ψ k ( i ) ψ ∗ k ( j ) h c † i c j i . (1.17)We will compare these expectation values in two different states, one is the exact thermal state and theother one is the steady state of the open system dynamics. Decay Rates:
We will also compute the rate of relaxation ∆ to the steady state. This is obtained fromthe spectrum of the Liouvillean operator L as discussed in Section 1.3 and Appendix A. The inverse of thisrate determines the time needed to come exponentially close to the steady state. The rate typically decreaseswith increasing system size n , either as an inverse polynomial ∆ ∼ /n a or exponentially ∆ ∼ e − bn . Entanglement:
Finally, to study the entanglement structure of the current-carrying states we computethe mutual information and conditional mutual information. If we consider two regions, A and B, of oursystem, then the mutual information is defined as
M I ( A : B ) = S ( A ) + S ( B ) − S ( AB ) , (1.18)where S ( X ) denotes the von Neumann entropy of region X . The conditional mutual information betweenthree regions, A, B, and C (see Figure 2), of a system is given by CM I ( A : C | B ) = S ( AB ) + S ( BC ) − S ( B ) − S ( ABC ) . (1.19)Notice that the quantities defined above are expressed in terms of von Neumann entropy only. To computethe entropy of a region X , we first define the correlation matrix G Xαβ = h c † α c β i , (1.20)where α and β denote orbitals within X . It is worth mentioning that G X has eigenvalues between 0 and1. We can then define the entropy as S ( X ) = − Tr( G X ln G X + (1 − G X ) ln(1 − G X )) . (1.21)The interpretation of vanishing mutual information is simple. If M I ( A : C ) = 0 then ρ AC = ρ A ⊗ ρ C andthe AC is uncorrelated. The interpretation of vanishing conditional mutual information is more subtle. If CM I ( A : C | B ) = 0 then A → B → C forms a “quantum Markov chain” [54] meaning that C is independentof A given the state of B . As discussed at length in Refs. [53–60], a quantum Markov chain generalizesthe classical notion of a Markov chain in the sense that the total state ρ ABC can be recovered from themarginals ρ AB and ρ BC . This does not imply that A is uncorrelated with C , but it does imply that theyare unentangled. Hence vanishing conditional mutual information indicates a certain kind of short-range6ntanglement in the state. It is this physics of approximate conditional independence which was used inRefs. [59, 60] to show that certain NESS have tensor network representations. In order to compute our observables, we use a technique developed by Prosen [68], which allows us to solvethe Lindblad master equation using a canonical quantization in the Fock space of operators. The method isapplicable under the condition that the Hamiltonian is quadratic and the jump operators L k are linear infermionic operators. The key idea is to write the Liouvillean (see Appendix A) in terms of adjoint Majoranamaps and diagonalize it in term of normal master modes (anti-commuting operators acting on the Fock spaceof density operators). Then the physical NESS is given by the zero-mode eigenvector. Below we review themain results of Ref. [68]. A more detailed discussion of the method is included in Appendix A.First, rewrite the canonical fermion operators in terms of Majorana operators w j − = c j + c † j , w j = i ( c j − c † j ) (1.22)satisfying the anti-commutation relation { w j , w k } = 2 δ j,k . Throughout this section the labels i, j, ... runover all fermion modes in the problem, spatial and otherwise. Then the Hamiltonians of our systems can bewritten as quadratic forms in Majorana fermions H = n X j,k =1 w j H jk w k (1.23)while the Lindblad operators can be written as linear combinations of w j L k = n X j =1 l k,j w j . (1.24)Ref. [69] shows that the properties of the NESS can be derived from the Liouvillean shape matrix A , whichis an antisymmetric 4 n × n matrix incorporating the parameters of the bath operators and Hamiltonian A j − , k − = − iH jk − M kj + M jk A j − , k = 2 iM jk A j, k − = − iM kj A j, k = − iH jk − M jk + M kj (1.25)where M is a Hermitian matrix with entries M jk = 12 X i l i,j l ∗ i,k . (1.26)It is worth mentioning that our formula for A is different from the one in Ref. [68] in that we swap theindices of M . We believe this merely reflects a typo in Ref. [68].Recall that the eigenvalues of a complex antisymmetric matrix of even dimension always come in pairs ± β . Therefore we can order the eigenvalues of A as β , − β , β , − β , . . . , β n , − β n , with Re( β ) ≥ Re( β ) ≥ . . . ≥ Re( β n ) ≥
0. Let v , v , . . . , v n be the corresponding eigenvectors, written as column vectors.Ref. [69] proves three key results that we use to compute our observables. If the eigenvalues of A havestrictly positive real parts, Re( β j ) >
0, then1. The non-equilibrium steady state is unique.2. The rate of exponential relaxation to NESS is given by ∆ = 2 Re( β n ). There is also an extra factor of in the definition of M which comes from a rescaling of Lindblad operators L k → L k √ relative to their definition in Ref. [68].
7. The expectation value of any quadratic observable w j w k in NESS is given by h w j w k i = 2 n X m =1 v m, j − v m − , k − . (1.27)These results may be understood intuitively by noting that, roughly speaking, the β j represent the“energies of excitations”. Hence the NESS is unique when all β j have non-zero real part because all otherstates decay. Similarly, the rate of relaxation is determined by the slowest decaying exctiation, correspondingto β n in the ordering we have chosen.Notice that the results above provide all the necessary information to compute our observables, sincethey are all expressed in terms of the two-point correlation functions h c † j c k i , which of course is quadratic inMajorana operators. One can further use Wick’s theorem to compute expectation values of any higher orderobservable with an even number of fermion operators.The implementation of these techniques involves computing the eigenvalue decomposition of very largematrices. We used the Multiprecision Computing Toolbox for MATLAB [70] to solve for the eigenstates andeigenvalues of A with high precision (32 digits). The extra precision was required to ensure that we cancorrectly group the eigenvalues into ( β, − β ) pairs and that the eigenvectors are orthonormal. In this subsection we describe our analysis of the steady state physics of the quantum wire model. Throughoutthis discussion we set the nearest neighbor hopping w = 1, so that we effectively measure all energies inunits of w and all times in units of 1 /w . We first study the physics of thermalization when the leads areunbiased. Then we consider the physics of NESS when the leads are biased. Next we discuss the structureof entanglement in the steady state. Finally, we discuss the effects of disorder. We begin our analysis of the quantum wire model by investigating thermalization with unbiased leads. Tothis end we take β L = β R = β and µ L = µ R = 0. Recall that the open system dynamics is designed todrive each decoupled lead (L and R) into its decoupled thermal state. However, the leads are coupled to themiddle wire which, except for the lead coupling at its boundaries, enjoys unitary dynamics. The Hamiltoniancouplings are taken to be w = w ′ = 1. The question is how well the wire is thermalized by its boundarycouplings to the leads.We answer this question by comparing the exact thermal occupation numbers of the decoupled wire (with w ′ = 0) with the expectation values of the corresponding number operators in the steady state. We expectagreement if (1) the leads are effectively thermalizing the wire and (2) the effects of the lead-wire coupling w ′ is small. To be clear, throughout the analysis we use the decoupled ( w ′ = 0) wire energy eigenstates andwe compare the expectation value of the number operator in thermal equilibrium and in the steady state ofthe open system dynamics. Data for β = . , , γ (the strength of the dissipative terms in L ) is shown in Panels (a), (b), and (c) of Figure 3.For β = . ,
1, we find that γ = .
01 leads to a steady state distribution of occupation numbers that agreeswell with the thermal result. Increasing the lead size does not dramatically effect the final steady state. Forthe the lowest temperatures reached, β = 5, a lead of size N L = N R = 40 is not sufficient to thermalize awire of size N W = 120. As demonstrated in Panel (d) of Figure 3, by taking larger leads we are able toachieve thermal equilibrium even for temperatures as low as β = 5.What lead size is required to have approximate thermalization of the wire as a function of inversetemperature β ? Figure 4 shows the minimal size N L = N R of the leads needed in order for the occupationnumbers to be close (within a few percent) to the thermal distribution for different temperatures. Note againthat the size of the contacts is only important at low temperatures. Similar results are also obtained whenstudying energy eigenstates of the entire coupled wire and lead system. We conclude that a sufficiently largelead is able to thermalize a wire even at low temperature via only boundary couplings.8 ε F e r m i on O cc upa t i on N u m be r γ = 1 γ = 0.1 γ = 0.01Fermi distribution ( a ) −2 −1.5 −1 −0.5 0 0.5 1 1.5 20.10.20.30.40.50.60.70.80.9 Single Particle Energy ε F e r m i on O cc upa t i on N u m be r γ = 1 γ = 0.1 γ = 0.01Fermi distribution ( b ) −2 −1.5 −1 −0.5 0 0.5 1 1.5 200.10.20.30.40.50.60.70.80.91 Single Particle Energy ε F e r m i on O cc upa t i on N u m be r γ = 1 γ = 0.1 γ = 0.01Fermi distribution ( c ) −2 −1.5 −1 −0.5 0 0.5 1 1.5 200.10.20.30.40.50.60.70.80.91 Single Particle Energy ε F e r m i on O cc upa t i on N u m be r γ = 1 γ = 0.1 γ = 0.01Fermi distribution ( d ) Figure 3: Comparison of energy eigenstate occupation numbers between thermal equilibrium and the steadystate of L . (a) Occupation numbers for β = 0 . N L = N R = 40 and N W = 120, (b) Occupation numbersfor β = 1 . N L = N R = 40 and N W = 120, (c) Occupation numbers for β = 5 . N L = N R = 40 and N W = 120, and (d) Occupations number for β = 5 . N L = N R = N W = 120. β R equ i r ed ba t h s i z e Figure 4: The minimal size of the leads required to approximately thermalize the wire as a function of β .Approximate thermalization is roughly defined as having energy occupation numbers within a few percentof the thermal result. Recall that we do not expect exact agreement with the decoupled ( w ′ = 0) thermalresult anyway since w ′ = 0 in the Liouvillean L . 9 β L i nea r r e s pon s e c ondu c t an c e G Numerically computedBallistic transport
Figure 5: Landauer formula for ballistic conductance (red) compared to the numerically evaluated conduc-tance in the steady state of L (blue) as a function of β . System parameters are N L = N W = N R = 80, w = w ′ = 1, and γ = 0 . We continue our analysis by probing the extent to which biased leads define a steady state with the expectedtransport physics. To this end we take β L = β R = β and µ L = − µ R . We compare the value of the currentflowing through the wire in steady state with the corresponding Landauer formula for conductance [71].Recall that this formula is obtained from an infinite wire model where the left and right moving particlesare emitted from separate reservoirs with potentially different temperatures and chemical potentials.The linear response conductance is defined as G = I/ ( µ L − µ R ) for small µ L − µ R . We study theconductance as a function of temperature β L = β R = β by computing the current in the middle of the chainin the steady state with weakly biased leads. As reviewed in Appendix D of [59], the current predicted bythe ballistic Landauer formula [71] is I = Z π dk π v k f ( ǫ k − µ L , β L ) + Z − π dk π v k f ( ǫ k − µ R , β R ) , (2.1)where ǫ k = − w cos k is the energy, v k = ∂ǫ k ∂k = 2 w sin k is the group velocity, and f ( ǫ, β ) = ( e βǫ + 1) − is the Fermi distribution. This can be explicitly integrated to give I = 12 π (cid:20) β L ln (cid:18) e β L ( µ L +2 w ) + 1 e β L ( µ L − w ) + 1 (cid:19) + 1 β R ln (cid:18) e β R ( µ R − w ) + 1 e β R ( µ R +2 w ) + 1 (cid:19)(cid:21) (2.2)Next, in Figure 5 we plot the ballistic transport conductance and the numerically computed conductancein the steady state as a function of β = β L = β R , for a system with 240 sites ( N L = N W = N R = 80)with w = w ′ = 1 and γ = 0 .
05. Although there is an interesting feature near β = 2, the steady statecurrent generally matches the Landauer formula quite well. We also obtained numerical evidence that thelow temperature Landauer result is well approximated by the steady state value in the limit of large leads.Together we take these data as evidence that, just as the unbiased steady state captures well the physicsof thermal equilibrium, the biased steady state captures well the physics of current carrying states in thenon-interacting fermion wire.We can also show that relaxation to the steady state takes place in a reasonable time frame, inversepolynomial in the system size. Recall that the relaxation rate is determined from the real part of thespectrum of the Liouvillean as described in Section 1.3. Figure 6 shows the relaxation rate ∆ as a functionof total system size n for different inverse temperatures β . Other parameters are w = w ′ = 1, γ = 0 .
01, and µ R = − µ L = 0 . It doesn’t matter where in the wire we measure the current because in a stationary state the current is conserved. .5 3 3.5 4 4.5 5 5.5 6 6.5−14−13−12−11−10−9−8−7−6−5 Log of System Size n Log o f R e l a x a t i on R a t e ∆ β = 0.1 β = 1.0 β = 5.0 Figure 6: Log of the relaxation rates ∆ as a function of system size n = N L + N W + N R . Log o f M u t ua l I n f o r m a t i on M I ( A : C ) β = 0.1 β = 0.2 β = 0.5 β = 1.0 β = 2.0 β = 3.0 β = 5.0 ( a ) Log o f C ond i t i ona l M u t ua l I n f o r m a t i on C M I ( A : C | B ) β = 0.1 β = 0.2 β = 0.5 β = 1.0 β = 2.0 β = 3.0 β = 5.0 ( b ) Figure 7: (a)
Log of
M I ( A : C ) as a function of separation between regions A and C for different inversetemperatures β , (b) Log of
CM I ( A : C | B ) as a function of size of region B for different inverse temperatures β . The region geometries are defined in Figure 2. We continue our analysis of the wire model by studying the entanglement structure of the steady state. Theregion configurations we study were defined in Ref. [59] and are motivated by the physics of approximatecondition independence [60], [53–58]. Consider a system with w = w ′ = 1, µ R = − µ L = 0 .
1, sufficiently largeleads N L = N R = N W = 120, and sufficiently small γ = 0 .
05 (these choices are based on the parametersthat give the best thermalization results). The mutual information
M I ( A : C ) and conditional mutualinformation CM I ( A : C | B ) for the partition in Figure 2, as a function of the size of the middle region B are presented in Figure 7. Different curves correspond to different values of β = β L = β R . Notice that themutual information is basically equal to the conditional mutual information, and both decrease exponentiallyfast. When the logarithm of the MI/CMI is less than −
30, this should be interpreted as the system having 0MI/CMI. The reason those quantities are not exactly 0 is due to the numerical precision of the computation.We may also study the slope s of the linear regions as a function of temperature. A plot of 1 /s vs β isshown in Figure 8 where we find an approximately linear scaling of 1 /s with β . As we briefly mentioned in Section 1.3 we can add a disorder term V x c † x c x at every site x in the region W ofthe wire. We choose V x from a continuous uniform distribution over the interval [ − V , V ]. We look at smalldisorder values, where V = 0 . w , although the results of this section hold true for larger V . For a systemwith w = w ′ = 1, µ R = − µ L = 0 . N L = N R = N W = 120, and γ = 0 .
05, we plot the mutual informationand conditional mutual information (averaged over 500 disorder realizations) in Figure 9. We see that thetwo quantities decay roughly the same way as in the absence of disorder.Similarly, we study the effect of disorder on relaxation rates. Consider the 1D system with w = w ′ = 1,11 β I n v e r s e s l ope s Figure 8: The inverse of the slope of the linear region of CMI as a function of β . The linear dependenceon β is expected, at least at large β , from the low temperature thermal physics in which there is a thermalcorrelation length proportional to β . Log o f M u t ua l I n f o r m a t i on M I ( A : C ) β = 0.1 β = 1.0 β = 5.0 ( a ) Log o f C ond i t i ona l M u t ua l I n f o r m a t i on C M I ( A : C | B ) β = 0.1 β = 1.0 β = 5.0 ( b ) Figure 9: (a)
Log of
M I ( A : C ) as a function of separation between regions A and C for different inversetemperatures β in the presence of disorder V = 0 . w , (b) Log of
CM I ( A : C | B ) as a function of size ofregion B for different inverse temperatures β in the presence of disorder V = 0 . w . Averaging is performedover 500 disorder realizations. 12
50 100 150 200 250 300−40−35−30−25−20−15−10−5 System Size n
Log o f R e l a x a t i on R a t e ∆ V = 0.1wV = 0.5wV = 1.0w Figure 10: Log of relaxation rates ∆ as a function of system size n for different disorder strengths V .Averaging is performed over 1000 disorder realizations. Note that unlike the translation invariant case, herethe decay rate ∆ scales exponentially with the total system size n , ∆ ∼ e − bn . γ = 0 . µ R = − µ L = 0 .
1, and β = 1 .
0. We compute the relaxation rates for three disorder strengths, V = 0 . w, . w, w , and average the results over 1000 disorder realizations. The results are shown in Figure10. Even though we present the plots only for β = 1 .
0, a similar behavior is observed at both high ( β = 0 . β = 5 . In this subsection we describe our analysis of the steady state physics of the Chern insulator model. Through-out this discussion we set the nearest neighbor hopping w = 1, so that we effectively measure all energiesin units of w and all times in units of 1 /w . We first study the physics of thermalization when the leadsare unbiased. Then we consider the physics of the NESS when the leads are biased. Finally we discuss thestructure of entanglement in the steady state. A related study of a bosonic symmetry protected state subjectto open system dynamics may be found in Ref. [72] Here we follow closely our analysis of the one-dimensional quantum wire in the case of the two-dimensionalChern insulator. In what follows, the parameters of the model are set to V = 3, c = 1, e s = 0 . t = 1, w = 1, w ′ = 0 .
1, so that the system has an energy gap with two edge states. Consider the unbiased casewhere β L = β R = β and the chemical potential is inside the gap µ L = µ R = − .
6. We find that the leadsare able to drive the system to thermal equilibrium in the high temperature limit. The energy occupationnumbers for a system with N L,x = N L,y = N W,x = N W,y = N R,x = N R,y = 10 and β = 0 . V →
0, the energy gap closes, and the system behaves like an ordinary metal. In thisregime, we find that sufficiently large leads can thermalize the system even at low temperatures. The energyoccupation numbers for β = 1 . V = 0 at the lowest temperatures. As weshow below, the rate of decay to the steady state is also slow for the Chern insulator, and these observationsmay be related. The structure of the current carrying state is generally more complex in the Chern insulator model relativeto the quantum wire model because there exist many more current carrying modes in the two-dimensionalmodel. However, when the chemical potential is chosen to sit within the bulk energy gap, then the edge13 ε F e r m i on O cc upa t i on N u m be r γ = 0.01 γ = 0.1Fermi distribution ( a ) −4 −3 −2 −1 0 1 2 3 400.10.20.30.40.50.60.70.80.91 Single Particle Energy ε F e r m i on O cc upa t i on N u m be r Fermi distribution γ = 0.1 ( b ) Figure 11: (a)
Occupation numbers for β = 0 . V = 3, (b) Occupation numbers for β = 1 . V → C u rr en t Figure 12: The magnitude of the current flowing through every link of the Chern insulator model. Thephysical parameters are N L,y = N W,y = N R,y = N W,x = 20, N L,x = N R,x = 10, µ L = − . µ R = − .
5, and β L = β R = 1. As expected for this parameter regime, most of the current is carried by the edge states withvery little bulk conduction.
50 100 150 200 250 300 350 400−26−24−22−20−18−16−14−12−10−8−6 System Size n
Log o f R e l a x a t i on R a t e ∆ β = 0.1 β = 1.0 β = 5.0 ( a ) Log o f R e l a x a t i on R a t e ∆ β = 0.1 β = 1.0 β = 5.0 ( b ) Figure 13: (a)
Log of relaxation rate ∆ as a function of system size n for the insulator state, (b) Log ofrelaxation rate ∆ as a function of log of system size n for the metallic state.14
10 20 30 40 50 60−18−16−14−12−10−8−6−4−20 Size of width of region B
Log o f M u t ua l I n f o r m a t i on M I ( A : C ) β = 0.1 β = 1.0 β = 5.0 ( a ) Log o f C ond i t i ona l M u t ua l I n f o r m a t i on C M I ( A : C | B ) β = 0.1 β = 1.0 β =5.0 ( b ) Figure 14: (a)
Log of
M I ( A : C ) as a function of x-coordinate separation between regions A and C fordifferent inverse temperatures β , (b) Log of
CM I ( A : C | B ) as a function of the width of region B fordifferent inverse temperatures β .states of the Chern insulator are the primary carriers of current at low temperature. The current for a systemwith N L,y = N W,y = N R,y = N W,x = 20 and N L,x = N R,x = 10 biased by µ L = − . µ R = − . β L = β R = 1 is shown in Figure 12. We see that, as expected, the current is carried primarilyby the edge states. We conclude that in addition to driving the system to thermal equilibrium, the leadscan also drive the system into a NESS with the expected properties.We also study the time needed to reach the NESS in the Chern insulator model. Fix N L,y = N W,y = N R,y = 8 and vary the length of the insulator and leads N L,x = N W,x = N R,x . Also fix the chemicalpotential inside the gap µ L = − . µ R = − . γ = 0 .
05. Figure 13 shows the relaxation rate ∆as a function of total system size for different values of inverse temperature β . Panel (a) shows that thedecay rate depends exponentially on system size in the Chern insulating phase. This indicates a rather slowapproach to the NESS and is similar to the decay rate obtained for the one-dimensional disordered insulator.Panel (b) shows that the decay rate for the metallic phase depends approximately inverse polynomially onsystem size. This is the same qualitative structure as obtained from the one-dimensional clean metallic wire. Finally, we turn again to an analysis of the entanglement structure of the NESS. Consider a 2D system with N L,x = N R,x = 10, N W,x = 40, and N L,y = N W,y = N R,y = 8. We set γ = 0 .
05 and place the chemicalpotential in the energy gap µ L = − . µ R = − .
5. The plots for MI and CMI with different values of β areshown in Figure 14. Note that the discontinuity in the graphs occurs exactly when region B becomes thewhole insulator, and regions A and C become the left and right lead respectively. As in all other cases, theMI and CMI decay exponentially with separation or with the size of B respectively. Hence the physics ofapproximate conditional independence survives in the NESS of a two-dimensional model. In this paper we have given detailed evidence that one can design a relatively simple lead that thermalizestarget non-interacting fermion systems in one and two dimensions. By biasing the leads, the systems canalso be driven out of equilibirum to give expected results. Furthermore, we found that not every lead isable to thermalize a given target, at least in the non-interacting limit considered here. In particular, lowtemperatures seem to require a larger lead. It is reasonable to conjecture that adding interactions to thetarget system will help render thermalization more universal, e.g. independent of the details of the bath,but much remains to be understood in this context. In particular, it would be desirable to have a betterunderstanding of what parameters, e.g. lead size and γ , lead to the best thermalization. More generally,this is just one small step towards a general framework for constructing leads that are able to effectivelythermalize a variety of models.We also showed that the time to reach a steady state is reasonable and scales with an inverse power of thesystem-plus-lead size, at least for the one-dimensional quantum wire and two-dimensional Chern insulator15n its metallic phase. We did find that the environment-lead interaction rate, γ , needed to be taken smallto achieve thermalization. Such a small γ further slows the rate of approach to the steady state, but we didnot find evidence that γ had to scale with system size. In the presence of disorder in one dimension or inthe insulating phase in two dimensions, we found that the time to reach a steady state was longer, scalingexponentially with the system size. This presents a potential challenge to the open system dynamics methodfor finding NESS, although at least in one dimension we know this slowdown is physical and is due to thephysics of localization [73]. Given that larger baths better approximate thermal equilibrium but have slowerrelaxation times, it would be interesting to investigate optimizing the bath size with respect to these twocompeting demands. It would also be desirable to investigate the effects of interactions on the relaxationtimes of the open system dynamics.Finally, we showed that the systems we investigated have little entanglement in their steady states. Thusalthough we used free fermion technology in our calculations, the low degree of entanglement implies thatthe states in question could also have been represented in a tensor network form. Such a low entanglementstructure has been argued to generalize beyond the non-interacting limit [27, 59]. It might be illuminat-ing to go through the exercise of writing the non-interacting fermion steady states using tensor networks,particularly as a first step towards including the effects of interactions.In summary, our results provide further support for the previously outlined entanglement-based approachto calculating electrical and thermal currents in strongly interacting systems. The main idea is to representthe state of the system using tensor network methods and then to dynamically evolve, using an appropriateopen system dynamics, to a current carrying steady state. Here we have shown how to design a non-interacting fermion lead which has three key properties: (1) it effectively thermalizes the system down to lowtemperatures, (2) it does so in a reasonable time frame, and (3) it is not overly complex. Pioneering matrixproduct based computations have already been carried out in one dimension at high temperature [27], andwhile there are still algorithmic and conceptual barriers to implementing the general program outlined inthe introduction, our results have shed light on the design of lead systems. We have also raised several issuesin the physics of thermalization of open non-interacting fermion models. Our future work will be concernedwith testing these ideas for interacting systems at low temperatures in a variety of dimensions. Acknowledgements:
CZ is supported by a Stanford UAR Major Grant and by the Clifford C. F. WongUndergraduate Scholarship Fund. BGS is supported by the Simons Foundation as part of the It From QubitCollaboration. We thank D. Freeman, R. Mahajan, N. Tubman, and P. Hayden for discussions on relatedissues. We also thank J. Andress for feedback on the manuscript.
A Review of Prosen’s method
In this appendix we discuss in more detail the method for solving the Lindblad master equation for a systemwith Hamiltonian and bath operators given by Equations (1.23) and (1.24) respectively. Begin with thedefinition of Liouvillean as the operator acting on ρ on the right hand side of Lindblad’s equation dρdt = − i [ H, ρ ] + X k L k ρL † k − X k { L † k L k , ρ } ≡ ˆ L ρ. (A.1)Consider a 4 n dimensional Liouville space of operators K , with an inner product defined as follows h x | y i = 2 − n tr( x † y ) , x, y ∈ K . (A.2)A complete orthonormal basis for this space is given by operator-products P α ,α ,...,α n = w α w α · · · w α n n , α j ∈ { , } . (A.3)It turns out that | P ~α i is a fermionic Fock basis and we can define creation and annihilation linear mapsover K ˆ c j | P ~α i = α j | w j P ~α i , ˆ c † j | P ~α i = (1 − α j ) | w j P ~α i , (A.4)which satisfy the canonical anti-commutation relations16 ˆ c j , ˆ c k } = 0 , { ˆ c j , ˆ c † k } = δ j,k . (A.5)The Liouvillean can be written as a quadratic form in terms of these mapsˆ L = − i n X j,k =1 ˆ c † j H jk ˆ c k + 12 X i n X j,k =1 l i,j l ∗ i,k ˆ L j,k , (A.6)where ˆ L j,k = ( + exp( iπ ˆ N ))(2ˆ c † j ˆ c † k − ˆ c † j ˆ c k − ˆ c † k ˆ c j )+ ( − exp( iπ ˆ N ))(2ˆ c j ˆ c k − ˆ c j ˆ c † k − ˆ c k ˆ c † j ) (A.7)and N = P j ˆ c † j ˆ c j is the number of adjoint fermions. Notice that the Liouvillean commutes with theparity operator ˆ P = exp( iπ ˆ N ) and hence the operator space can be decomposed into even an odd subspacesvia an orthogonal projection K ± = ( ± ˆ P ) K . Since we are interested in expectation values of observablesquadratic in fermion operators, we restrict ourselves to the subspace K + where ˆ L j,k has the formˆ L j,k |K + = 4ˆ c † j ˆ c † k − c † j ˆ c k − c † k ˆ c j . (A.8)Combining Equations (A.6) and (A.8) yields a compact representation of the Liouvillean on the evensubspace ˆ L + = − c † (2 iH + M + M T )ˆ c + 2ˆ c † ( M − M T )ˆ c † , (A.9)where ˆ c = [ˆ c , ˆ c , . . . , ˆ c n ] T is a column vector, M is the matrix containing information about the baths,defined in Equation (1.26), and H is the Hamiltonian matrix with entries H jk . We can further simplify thisrepresentation by introducing 4 n Hermitian Majorana maps ˆ a k ˆ a j − = 1 √ c j + ˆ c † j ) , ˆ a j = i √ c j − ˆ c † j ) . (A.10)In terms of these maps, Equation (A.9) becomesˆ L + = ˆ aA ˆ a − A , (A.11)where ˆ a = [ˆ a , ˆ a , . . . , ˆ a n ] T is a column vector, A is the antisymmetric matrix introduced in 1.25, and A = 2tr( M ).Assuming that A is diagonalizable with eigenvalues β , − β , β , − β , . . . , β n , − β n and eigenvectors v , v , . . . , v n , the Liouvillean can be written in normal formˆ L + = − n X j =1 β j ˆ b ′ j ˆ b j , (A.12)where ˆ b j and ˆ b ′ j are normal master mode maps defined byˆ b j = v j − · ˆ a, ˆ b ′ j = v j · ˆ a, (A.13)and obeying the anti-commutation relations { ˆ b j , ˆ b k } = 0 , { ˆ b j , ˆ b ′ k } = δ j,k , { ˆ b ′ j , ˆ b ′ k } = 0 . (A.14)This normal form representation of the Liouvillean leads to the results stated in Section 1.3. We encouragethe reader to refer to [68] for a more detailed derivation.17 eferences [1] Supriyo Datta. Electronic transport in mesoscopic systems . Cambridge University Press, 1997.[2] A. Blandin, A. Nourtier, and D.W. Hone. Localized time-dependent perturbations in metals : formalism andsimple examples.
J. Phys. France , 37(4):369–378, 1976.[3] A. J. Leggett, S. Chakravarty, A. T. Dorsey, Matthew P. A. Fisher, Anupam Garg, and W. Zwerger. Dynamicsof the dissipative two-state system.
Rev. Mod. Phys. , 59:1–85, Jan 1987.[4] Jong E. Han. Solution of electric-field-driven tight-binding lattice coupled to fermion reservoirs.
Phys. Rev. B ,87:085119, Feb 2013.[5] Herbert Spohn. Kinetic equations from hamiltonian dynamics: Markovian limits.
Rev. Mod. Phys. , 52:569–615,Jul 1980.[6] R¨udiger Schack and Todd A. Brun. A c++ library using quantum trajectories to solve quantum master equations.
Computer Physics Communications , 102(1–3):210 – 228, 1997.[7] Sze M Tan. A computational toolbox for quantum and atomic optics.
Journal of Optics B: Quantum andSemiclassical Optics , 1(4):424, 1999.[8] A. Vukics and H. Ritsch. C++qed: an object-oriented framework for wave-function simulations of cavity qedsystems.
The European Physical Journal D , 44(3):585–599, 2007.[9] J.R. Johansson, P.D. Nation, and Franco Nori. Qutip: An open-source python framework for the dynamics ofopen quantum systems.
Computer Physics Communications , 183(8):1760 – 1772, 2012.[10] Guifr´e Vidal. Efficient simulation of one-dimensional quantum many-body systems.
Phys. Rev. Lett. , 93:040502,Jul 2004.[11] Lothar M¨uhlbacher and Eran Rabani. Real-time path integral approach to nonequilibrium many-body quantumsystems.
Phys. Rev. Lett. , 100:176403, May 2008.[12] Sean A. Hartnoll and Diego M. Hofman. Locally Critical Resistivities from Umklapp Scattering.
Phys. Rev.Lett. , 108:241601, 2012.[13] Raghu Mahajan, Maissam Barkeshli, and Sean A. Hartnoll. Non-Fermi liquids and the Wiedemann-Franz law.
Phys. Rev. , B88:125107, 2013.[14] Simone Giombi, Shiraz Minwalla, Shiroman Prakash, Sandip P. Trivedi, Spenta R. Wadia, and Xi Yin. Chern-Simons Theory with Vector Fermion Matter.
Eur. Phys. J. , C72:2112, 2012.[15] Guy Gur-Ari and Ran Yacoby. Correlators of Large N Fermionic Chern-Simons Vector Models.
JHEP , 02:150,2013.[16] S. Bhattacharyya, S. Minwalla, V. E. Hubeny, and M. Rangamani. Nonlinear fluid dynamics from gravity.
Journal of High Energy Physics , 2:045, February 2008.[17] Mikhail Pletyukhov, Dirk Schuricht, and Herbert Schoeller. Relaxation versus decoherence: Spin and currentdynamics in the anisotropic kondo model at finite bias and magnetic field.
Phys. Rev. Lett. , 104:106801, Mar2010.[18] R. Egger, L. M¨uhlbacher, and C. H. Mak. Path-integral monte carlo simulations without the sign problem:Multilevel blocking approach for effective actions.
Phys. Rev. E , 61:5961–5966, May 2000.[19] Klaus Mølmer, Yvan Castin, and Jean Dalibard. Monte carlo wave-function method in quantum optics.
J. Opt.Soc. Am. B , 10(3):524–538, Mar 1993.[20] Jean Dalibard, Yvan Castin, and Klaus Mølmer. Wave-function approach to dissipative processes in quantumoptics.
Phys. Rev. Lett. , 68:580–583, Feb 1992.[21] C. W. Gardiner, A. S. Parkins, and P. Zoller. Wave-function quantum stochastic differential equations andquantum-jump simulation methods.
Phys. Rev. A , 46:4363–4381, Oct 1992.[22] Tomaˇz Prosen. Third quantization: a general method to solve master equations for quadratic open fermi systems.
New Journal of Physics , 10(4):043026, 2008.[23] R. Alicki. Invitation to Quantum Dynamical Semigroups. In P. Garbaczewski and R. Olkiewicz, editors,
Dy-namics of Dissipation , volume 597 of
Lecture Notes in Physics, Berlin Springer Verlag , pages 239–264, 2002.[24] A. A. Gangat, T. I, and Y.-J. Kao. Steady States of Infinite-Size Dissipative Quantum Chains via ImaginaryTime Evolution.
ArXiv e-prints , August 2016.
25] M. M. Wolf, F. Verstraete, M. B. Hastings, and J. I. Cirac. Area Laws in Quantum Systems: Mutual Informationand Correlations.
Physical Review Letters , 100(7):070502, February 2008.[26] F. Verstraete, J. J. Garc´ıa-Ripoll, and J. I. Cirac. Matrix product density operators: Simulation of finite-temperature and dissipative systems.
Phys. Rev. Lett. , 93:207204, Nov 2004.[27] Tomaˇz Prosen and M. ˇZnidariˇc. Matrix product simulations of non-equilibrium steady states of quantum spinchains.
Journal of Statistical Mechanics: Theory and Experiment , 2009(02):P02035, 2009.[28] Toby S. Cubitt, Angelo Lucia, Spyridon Michalakis, and David Perez-Garcia. Stability of local quantum dissi-pative systems.
Communications in Mathematical Physics , 337(3):1275–1315, 2015.[29] Fernando G. S. L. Brand˜ao, Toby S. Cubitt, Angelo Lucia, Spyridon Michalakis, and David Perez-Garcia.Area law for fixed points of rapidly mixing dissipative quantum systems.
Journal of Mathematical Physics ,56(10):102202, 2015.[30] Michael Zwolak and Guifr´e Vidal. Mixed-state dynamics in one-dimensional quantum lattice systems: A time-dependent superoperator renormalization algorithm.
Phys. Rev. Lett. , 93:207205, Nov 2004.[31] R A Blythe and M R Evans. Nonequilibrium steady states of matrix-product form: a solver’s guide.
Journal ofPhysics A: Mathematical and Theoretical , 40(46):R333, 2007.[32] Steven R. White. Density matrix formulation for quantum renormalization groups.
Phys. Rev. Lett. , 69:2863–2866, Nov 1992.[33] U. Schollw¨ock. The density-matrix renormalization group.
Rev. Mod. Phys. , 77:259–315, Apr 2005.[34] Ulrich Schollw¨ock. The density-matrix renormalization group in the age of matrix product states.
Annals ofPhysics , 326(1):96 – 192, 2011.[35] Michael J. Hartmann, Javier Prior, Stephen R. Clark, and Martin B. Plenio. Density matrix renormalizationgroup in the heisenberg picture.
Phys. Rev. Lett. , 102:057202, Feb 2009.[36] Jian Cui, J. Ignacio Cirac, and Mari Carmen Ba˜nuls. Variational matrix product operators for the steady stateof dissipative quantum systems.
Phys. Rev. Lett. , 114:220601, Jun 2015.[37] E. Mascarenhas, H. Flayac, and V. Savona. A Matrix-Product-Operator Approach to the Nonequilibrium SteadyState of Driven-Dissipative Quantum Arrays.
ArXiv e-prints , April 2015.[38] Jiasen Jin, Davide Rossini, Rosario Fazio, Martin Leib, and Michael J. Hartmann. Photon solid phases in drivenarrays of nonlinearly coupled cavities.
Phys. Rev. Lett. , 110:163605, Apr 2013.[39] Hendrik Weimer. Variational principle for steady states of dissipative quantum many-body systems.
Phys. Rev.Lett. , 114:040402, Jan 2015.[40] Zi Cai and Thomas Barthel. Algebraic versus exponential decoherence in dissipative many-particle systems.
Phys. Rev. Lett. , 111:150403, Oct 2013.[41] M. Kliesch, D. Gross, and J. Eisert. Matrix-product operators and states: Np-hardness and undecidability.
Phys.Rev. Lett. , 113:160503, Oct 2014.[42] Tomaˇz Prosen. Spectral theorem for the lindblad equation for quadratic open fermionic systems.
Journal ofStatistical Mechanics: Theory and Experiment , 2010(07):P07020, 2010.[43] D. Karevski, V. Popkov, and G. M. Sch¨utz. Exact matrix product solution for the boundary-driven lindblad xxz chain.
Phys. Rev. Lett. , 110:047201, Jan 2013.[44] S R Clark, J Prior, M J Hartmann, D Jaksch, and M B Plenio. Exact matrix product solutions in the heisenbergpicture of an open quantum spin chain.
New Journal of Physics , 12(2):025005, 2010.[45] T. Prosen and M. ˇZnidariˇc. Diffusive high-temperature transport in the one-dimensional hubbard model.
Phys.Rev. B , 86:125118, Sep 2012.[46] T. Prosen. Matrix product solutions of boundary driven quantum chains.
ArXiv e-prints , April 2015.[47] T. Prosen. Exact nonequilibrium steady state of a strongly driven open xxz chain.
Phys. Rev. Lett. , 107:137201,Sep 2011.[48] T. Prosen. Open xxz spin chain: Nonequilibrium steady state and a strict bound on ballistic transport.
Phys.Rev. Lett. , 106:217206, May 2011.[49] T. Prosen. Exact nonequilibrium steady state of an open hubbard chain.
Phys. Rev. Lett. , 112:030603, Jan 2014.[50] Tomaˇz Prosen. Comments on a boundary-driven open xxz chain: asymmetric driving and uniqueness of steadystates.
Physica Scripta , 86(5):058511, 2012.
51] Vladislav Popkov and T. Prosen. Infinitely dimensional lax structure for the one-dimensional hubbard model.
Phys. Rev. Lett. , 114:127201, Mar 2015.[52] Viktor Eisler and Zolt´an Zimbor´as. Area-law violation for the mutual information in a nonequilibrium steadystate.
Phys. Rev. A , 89:032321, Mar 2014.[53] D´enes Petz. Sufficient subalgebras and the relative entropy of states of a von neumann algebra.
Comm. Math.Phys. , 105(1):123–131, 1986.[54] P. Hayden, R. Jozsa, D. Petz, and A. Winter. Structure of States Which Satisfy Strong Subadditivity of QuantumEntropy with Equality.
Communications in Mathematical Physics , 246:359–374, 2004.[55] O. Fawzi and R. Renner. Quantum conditional mutual information and approximate Markov chains.
ArXive-prints , October 2014.[56] D. Sutter, O. Fawzi, and R. Renner. Universal recovery map for approximate Markov chains.
ArXiv e-prints ,April 2015.[57] M. M. Wilde. Recoverability in quantum information theory.
Proceedings of the Royal Society of London SeriesA , 471:20150338, October 2015.[58] M. Junge, R. Renner, D. Sutter, M. M. Wilde, and A. Winter. Universal recovery from a decrease of quantumrelative entropy.
ArXiv e-prints , September 2015.[59] Raghu Mahajan, C Daniel Freeman, Sam Mumford, Norm Tubman, and Brian Swingle. Entanglement structureof non-equilibrium steady states. arXiv preprint arXiv:1608.05074 , 2016.[60] B. Swingle and J. McGreevy. Mixed s-sourcery: Building many-body states using bubbles of Nothing.
ArXive-prints , July 2016.[61] J. Maldacena. The Large-N Limit of Superconformal Field Theories and Supergravity.
International Journal ofTheoretical Physics , 38:1113–1133, 1999.[62] S. Ryu and T. Takayanagi. Holographic Derivation of Entanglement Entropy from the anti de SitterSpace/Conformal Field Theory Correspondence.
Physical Review Letters , 96(18):181602, May 2006.[63] Simeng Yan, David A. Huse, and Steven R. White. Spin-liquid ground state of the s = 1/2 kagome heisenbergantiferromagnet.
Science , 332(6034):1173–1176, 2011.[64] E.M. Stoudenmire and Steven R. White. Studying two-dimensional systems with the density matrix renormal-ization group.
Annual Review of Condensed Matter Physics , 3(1):111–128, 2012.[65] A. Kshetrimayum, H. Weimer, and R. Orus. A simple tensor network algorithm for 2d steady states.
ArXive-prints , December 2016.[66] Xiao-Liang Qi, Yong-Shi Wu, and Shou-Cheng Zhang. Topological quantization of the spin hall effect in two-dimensional paramagnetic semiconductors.
Physical Review B , 74(8):085308, 2006.[67] C.-E. Bardyn, M. A. Baranov, C. V. Kraus, E. Rico, A. ˙Imamo˘glu, P. Zoller, and S. Diehl. Topology bydissipation.
New Journal of Physics , 15(8):085001, August 2013.[68] Tomaˇz Prosen. Third quantization: a general method to solve master equations for quadratic open fermi systems.
New Journal of Physics , 10(4):043026, 2008.[69] Tomaˇz Prosen and Bojan ˇZunkoviˇc. Exact solution of markovian master equations for quadratic fermi systems:thermal baths, open xy spin chains and non-equilibrium phase transition.
New Journal of Physics , 12(2):025016,2010.[70] LLC Advanpix. Multiprecision computing toolbox for matlab, 2006.[71] R. Landauer. Spatial variation of currents and fields due to localized scatterers in metallic conduction.
IBMJournal of Research and Development , 1(3):223–231, July 1957.[72] ´A. Rivas and M. A. Martin-Delgado. Topological Heat Transport and Symmetry-Protected Boson Currents.
ArXiv e-prints , June 2016.[73] P. W. Anderson. Absence of diffusion in certain random lattices.
Phys. Rev. , 109:1492–1505, Mar 1958., 109:1492–1505, Mar 1958.