CClosed hierarchies and non-equilibrium steady states of drivensystems
Israel Klich
Department of Physics, University of Virginia, Charlottesville, 22904, VA, USA
We present a class of tractable non-equilibrium dynamical quantum systems whichincludes combinations of injection, detection and extraction of particles interspersedby unitary evolution. We show how such operations generate a hierarchy of equationstying lower correlation functions with higher order ones. The hierarchy closes forparticular choices of measurements and leads to a rich class of evolutions whoselong time behavior can be simulated efficiently. In particular, we use the method todescribe the dynamics of current generation through a generalized quantum exclusionprocess, and exhibit an explicit formula for the long time energy distribution in thelimit of weak driving.
Contents
I. Introduction II. General framework G III. Non-Equilibrium Steady State Equation
IV. Examples of dynamics and dependence on initial condition V. Summary References a r X i v : . [ qu a n t - ph ] F e b I. INTRODUCTION
Significant activity has been devoted to the study of quantum systems out of equilibrium,with a rapid increase in interest due to the relevance to experiments with ultra-cold atomicgases, whose coherent evolution may be effectively controlled and decoupled from dissipationto a heat bath [1–3]. Non equilibrium dynamics is typically studied in processes such asexternal driving, repeated quantum measurements and quantum quenches. The fundamentalquestion that arises in such cases is what is the long term behavior of the system: does iteventually reach a non-equilibrium steady state? What is the nature of such a state?In studying the aforementioned non-equilibrium situations, some highly successful tools ofequilibrium statistical physics, such as linear response theory, may easily fail. Thus, thereis a need to develop new methods to deal with some of these problems. Here we focus onone such idea - that of establishing closed hierarchies in order to get tractable equationsfor correlation functions. Specifically, in many statistical mechanics problems, it is possibleto make a systematic connection between the evolution of n body density functions with n + 1 density functions. A prime example for such a set of relations is the Bogoliubov-Born-Green-Kirkwood-Yvon (BBKGY) hierarchy, which is the essential structure leading to theBoltzmann equation. In the Boltzmann equation, single particle densities are tied to higherorder correlation functions represented in the collision integral (see, e.g. [4]). In this letter,we describe the requirements on obtaining a hierarchy under general quantum operationson fermions. We then show how the hierarchy may be closed for a quantum system that isperiodically evolved, detected, and injected with current. Finally, we use the idea to describedynamics of current buildup, and the energy distribution in the long term non-equilibriumsteady state.To begin the discussion, consider the most general evolution of a density matrix, describingunitary evolution, measurements and interaction with the environment. Written as ρ → L ( ρ ) = Σ ν A ν ρA † ν ; Σ ν A † ν A ν = 1 (1)This form ensures ρ remains a non-negative matrix, and the normalization condition on theKrauss operators A i ensures that Tr ρ = 1 is preserved under the evolution.In general, there is no simple relation between correlation functions computed in state ρ before and after the evolution (1), which necessitates working in an exponentially largeHilbert space and is therefore often un-tractable.Hierarchy structures have been used before in the context of Kossakowski-Lindblad evolu-tion, which is a particular limit of (1). For example, the steady state of a dissipative XXspin chain in the presence of driving and dissipation has been studied extensively [5–8].Also, conditions for a closed hierarchy in the continuous time frame work where also statedin [9–11]. Here we concentrate on a discrete time framework, but also supply correspondingKossakowski-Lindblad results as a special limit. In other processes, the possibility of gettinga closed equation for Kossakowski-Lindblad evolution of noise averaged expectation valueswas studied in [12], to explore the stability of fractional charges to noisy hopping processes.We utilize the power of this approach to study a non-equilibirum process of current gener-ation, as schematically depicted in Fig. 1 (a). In this process, we connect site a to a lead,where a current is injected, and particles are allowed to go out at site b (two choices for b are shown). The process is explicitly described by ρ −→ U ((1 − r ) ρ + rα [ (cid:15) a (2 − (cid:15) a ) a † a ρa a + (1 − (cid:15) a (1 − n a )) ρ (1 − (cid:15) a (1 − n a ))] + (2) r (1 − α )[ (cid:15) b (2 − (cid:15) b ) a b ρa † b + (1 − (cid:15) b n b ) ρ (1 − (cid:15) b n b )]) U † , where n a/b = a † a/b a a/b checks for the presence of a fermion on the injection/extraction site,and U = e − iτ (cid:80) h nm a † n a m describes evolution between attempts during a time interval τ .Here r is the overall attempt rate, α is the relative probability of injecting vs extractingattempts, and (cid:15) a,b are related to the efficiency of the injection/extraction attempts: when (cid:15) a,b = 1, particle injection or removal happens with probability 1 if an attempt is made. Weshow below that this process leads to a closed equation (36) for the two point function ofthe system, which can be then computed numerically. It is important to emphasize that thelong time steady state reached by the system is not a thermal equilibrium state, in that theenergy occupation is very different from a Fermi-Dirac distribution governed by the singleparticle Hamiltonian h governing the evolution U .For small r , we find a remarkable asymptotic formula for the steady state distribution Φ k ≡(cid:104) steady | a † k a k | steady (cid:105) . Here k labels the eigenstates | k (cid:105) of the single particle hamiltonian h nm , h | k (cid:105) = E k | k (cid:105) . Let p a,k = |(cid:104) a | k (cid:105)| , p b,k = |(cid:104) a | k (cid:105)| be overlaps of these states with thesites a, b . Then Φ k is a function of the ratio p a,k /p b,k :Φ k = A + B p a,k p b,k (1 − α ) (cid:15) b + α(cid:15) a p a,k p b,k (3)Note the appearance of the relative injection rates/extraction rates: α(cid:15) a and (1 − α ) (cid:15) b .The coefficients A , B are given below in Eq. (39). We emphasize that this expression isvalid for any system obeying the form (2), and is non perturbative.In the limit of low tunneling probability, (cid:15) a , (cid:15) b →
0, the result depends only on the ration ofinjection to removal rates and simplifies to:Φ k ∼ α(cid:15) a (1 − α ) (cid:15) b p a,k p b,k α(cid:15) a (1 − α ) (cid:15) b p a,k p b,k (4)This last expression has a simple interpretation: the probability of occupying a given mode k is determined by the ratio between the effective tunneling probability into energy k fromsite a compared to the effective tunneling rate of the state k through site b . The limitof r, (cid:15) a , (cid:15) b →
0, also corresponds to the limit where a Kossakowski-Lindblad equation canbe used to describe (2). Indeed, as we show below, one can obtain (4) from Kossakowski-Lindblad treatment of the process (2).We stress that in the low tunneling limit, the steady state Φ k does not depend on systemdetails except the tunneling rates and the probabilities p a/b,k . However, going back to theformula (3), the details of the distribution depend of sensitively on the choice of parameters.In particular, we note that even if p a,k = 0, i.e. there is no overlap between a given energymode and the insertion site (or mode), Φ k can be non vanishing, due to higher order pro-cesses, a feature which is absent in the simpler Kossakowski-Lindblad limit expression (4).This feature illustrates the non-perturbative dependence of Φ k on the system parameters(and on (cid:15) a , (cid:15) b , α ).For illustration, we consider hopping on a chain of length N , with the standard Hamiltonian H hop = (cid:80) N − i =1 a † i a i +1 + h.c. corresponding to Dirichlet boundary conditions. In this case p a,k /p b,k = sin ( πakN +1 ) / sin ( πbkN +1 ). In Fig. 1 we illustrate the result with N = 100, andinjection at a = 1. We evolve the system from an initial vacuum state at t = 0. The resultsfor extraction at the final and penultimate sites b = 100 ,
99 respectively, show sensitivityto the choice of operation sites. The energy distribution is computed numerically at longtimes and is clearly seen to approach Φ in the long time limit. We stress that once drivinghas stopped, the energy distribution Φ will remain the stationary distribution under thesubsequent free evolution. Fig 2 shows the actual evolution of the density as we inject theparticles into the system. (a)(b) Φ ( k ) Fermi - Dirac n k - � - � � � � � � ��������������� � � � ( � ) � α = ��� ���� � = � ���� � = ��� � � = ���� Φ ( k ) Fermi - Dirac n k - � - � � � � � � ��������������� � � � ( � ) � α = ��� ���� � = � ���� � = �� � � = ���� Figure 1: (a) The fermion hopping model. (b) Approach to Φ( k ) for r = 0 . , α = 0 . , τ = 0 . b = 100 (upper panel) and b = 99 (lower panel). 300 iterations betweensuccessive curves. For reference a Fermi-Dirac distribution is shown.Figure 2: The fermion hopping model: Evolution of local density, (cid:104) a † i a i (cid:105) , in space and time(red/blue corresponds to high/low density, same parameters). II. GENERAL FRAMEWORK
We now turn to establishing the framework for our processes. We consider a system offermions on a lattice of N sites. In (1) we take Krauss operators of the form A ν = m ν U ν ,where U ν is an evolution under a non-interacting hamiltonian, and m ν is a polynomial oforder r ν in fermion operators a † , a . The evolution under L of a general correlation function, (cid:104) a † i ..a † i l a i ( l ..a i ( l l (cid:105) ≡ Tr ρ a † i ..a † i l a i ( l ..a i ( l l (5)is given by (cid:104) a † i ..a † i l a i ( l ..a i ( l l (cid:105) −→ (cid:104) a † i ..a † i l a i ( l ..a i ( l l (cid:105) + (6) (cid:80) ν T r ρU † ν m † ν [ a † i ..a † i l a i ( l ..a i ( l l , m ν U ν ]where the normalization relation in (1) was used.The assumption that the U ν are non interacting, means that U † ν a i U ν = u ν ; ij a j for someunitary matrix u ν ∈ U ( N ). As a consequence the evolution of the l + l correlation func-tion (5), is related in (6) to correlation functions of an order at most l + l + 2 max ν ( r ν ),establishing a hierarchy of equations.We emphasize that the resulting state may be arbitrarily complex. Indeed, even whenstarting with a non-interacting thermal state, ρ ∼ exp ( − h ij a † i a j ) and taking each A ν a noninteracting unitary, ρ evolves into a sum of exponentials of fermion bi-linears. Such a statecan be used to approximate any interacting state whose determinant quantum Monte Carlodescription does not suffer from a sign problem [13].Below, we list several fundamental operations under which the hierarchy closes at the twopoint function level, for G ij ≡ (cid:104) a † i a j (cid:105) , inducing a map G → K ( G ). We start with the obviousone:(I) The non-interacting evolution L u ( ρ ) = U ρU † , as described above, induces a map G ij → K u ( G ) ij ≡ ( u † Gu ) ij (7)We augment the free evolution with the following types of operations acting on a single par-ticle mode: particle detection, injection and extraction. Below, for simplicity of presentationwe will associate the operation with the mode associated with site i .Denote P i the matrix ( P i ) mn = δ im δ in the projection on site i , and P ⊥ i = 1 − P i , we introduce:(II) Particle detection at site i : L D,i ( ρ ) = n i ρn i + (1 − n i ) ρ (1 − n i ) (8)where n i = a † i a i . The induced map on G is: K D,i ( G ) = P ⊥ i GP ⊥ i + P i GP i . (9)The process (II) may be viewed as a “decoherence” of the correlations G in between site i and the rest of the lattice. As a linear super-operator on matrices, the measurement K D,i has a simple spectrum. It acts as identity on matrices which do not mix site i with the rest,hence the non-zero subspace of matrices has a dimension 1 + ( dimP ⊥ i ) . The complementaryzero subspace is spanned by the off diagonal blocks, of dimensionality 2( dimP ⊥ i ).(III) Removal of a particle from site i is described by L out,i ( ρ ) = a i ρa † i + (1 − n i ) ρ (1 − n i ) (10)with the induced map on G : K out,i ( G ) = P ⊥ i GP ⊥ i . (11)As a super operator this simple map may be viewed as a projection on the space of matricesthat do not have an ( i, j ) or ( j, i ) element for any j .(IV) Finally, this operation injects a particle at site i : L in,i ( ρ ) = a † i ρa i + n i ρn i (12)and induces the map K in,i ( G ) = P i + P ⊥ i GP ⊥ i . (13)We note that in contrast with ( I − III ), the injection K in,i is an in-homogenuous transfor-mation on matrices, a property which we use below to compute steady states.It is also possible to add another two operations which correspond to ”softer” particle motioninto and out of the system, without performing a direct measurement on the system. Theseare described by:( ˜III) Soft removal site i is described by L out,i,(cid:15) ( ρ ) = (cid:15) (2 − (cid:15) ) a i ρa † i + (1 − (cid:15)n i ) ρ (1 − (cid:15)n i ) (14)with the induced map on G : K out,i,(cid:15) ( G ) = P ⊥ i GP ⊥ i + (1 − (cid:15) ) P i GP ⊥ i + (1 − (cid:15) ) P ⊥ i GP i + (1 − (cid:15) ) P i GP i . (15)Here 0 ≤ (cid:15) ≤
1, with (cid:15) = 1 corresponding to the operation (III). Similarly, we have:( ˜IV) Soft injection at site i : L in,i,(cid:15) ( ρ ) = (cid:15) (2 − (cid:15) ) a † i ρa i + (1 − (cid:15) (1 − n i )) ρ (1 − (cid:15) (1 − n i )) (16)and induces the map K in,i,(cid:15) ( G ) = P ⊥ i GP ⊥ i + (1 − (cid:15) ) P i GP ⊥ i + (1 − (cid:15) ) P ⊥ i GP i + (1 − (cid:15) ) P i GP i + (cid:15) (2 − (cid:15) ) P i . (17)Below, unless remarked differently, we will refer to both soft and hard process together,ommiting the ˜ notation. We can combine any of the site operations (II-IV) with the uni-tary evolutions (I) mixing the the addressed site i with the rest of the sites. When noparticle injection is present, the particle extraction map will generically drive G to 0, i.e.( K u K Out,i ) n → K u K In,i ) n , with no extractionpresent, will result in G ij → δ ij , when n → ∞ , which is the state where all sites are occupied.On the other hand the unitary evolution ( I ) and the detection process ( II ) preserve theaverage particle number, i.e. (cid:104) (cid:80) i a † i a i (cid:105) = Tr G remains constant under K u , K M . A. Universality of the transformations (I,II,III,IV) on G The set of transformations (I,II,III,IV) generate all possible transformations on the two pointfunction G , keeping G a valid to point function by construction. In other words, given twovalid correlation matrices G and G , there is a set of operations of the form (I,II,III,IV)that will take us from G to G . Proof:
We have already seen that it is possible to get G = 0 by emptying the system. It istherefore enough to show that we can get any G starting from the zero matrix.To do so, let u be unitary matrices that diagonalize G , i.e.: u † Gu = K u ( G ) = diag ( λ , .., λ N ) (18)Observing the operation (17), and noting that λ = (2 − (cid:15) ) (cid:15) for (cid:15) = 1 − √ λ we have: diag ( λ , , ..,
0) = K in, ,(cid:15) (0) (cid:15) = (1 − (cid:112) λ ) (19)similarly: diag ( λ , λ , ..,
0) = K in, ,(cid:15) ( K in, ,(cid:15) (0)) (cid:15) = (1 − (cid:112) λ ) (20)We can continue this way to populate the diagonal and get diag ( λ , .., λ N ). Finally, we undothe unitary u and have G = K u † ( K in,N,(cid:15) N ( K in,N − ,(cid:15) N − ( .... ))) (21)with (cid:15) i = (1 − √ λ i ) at step i . B. Soft extraction by tunneling and removal from auxiliary site.
We note that it is possible to induce the Kraus operators corresponding to the transformation L out, ( ρ ) = (cid:15) (2 − (cid:15) ) a ρa † + (1 − (cid:15)n ) ρ (1 − (cid:15)n ) (22)with the induced map on G : K out, ( G ) = P ⊥ GP ⊥ + (1 − (cid:15) ) P GP ⊥ + (1 − (cid:15) ) P ⊥ GP + (1 − (cid:15) ) P GP . (23)without carrying out any direct measurement on the system, instead the measurements arecarried out on outside the system. We can represent the operation of removing a particlefrom site 0 by coupling the site by a tunneling Hamiltonian to an auxiliary site e , and makingthe ”hard” removal on the site e .To make the derivation clear, let us denote by ρ S the density matrix of our system of n fermionic sites. And the density matrix including the extra site e is ρ S + e . We first performoperations on the larger system ρ S + e , and compute the change in ρ S = Tr e ρ S + e followingthe process.The protocol is as follows.(1) Site e is decoupled from our system, and an operation of particle removal from e is done.Thus ρ S + e −→ (1 − n e ) ρ S + e (1 − n e ) + a e ρ S + e a † e ρ S .(2) We apply the evolution with a tunneling between site e and 0, using the Hamiltonian H t ∝ i ( a † e a − a † a e ). i.e. we evolve ρ S + e with: U θ = e θ ( a † e a − a † a e ) (24)Following these operations, we have to compute how ρ S transformed ρ S → ρ S, → ρ S, . Thiscan be done explicitly by choosing a basis for the Fock space. With Fermions we have to fixan ordering, and we take: | m, σ, → k (cid:105) = (cid:0) a † e (cid:1) m (cid:16) a † (cid:17) σ (cid:16) a † (cid:17) k .. (cid:16) a † N − (cid:17) k N − (cid:16) a † N (cid:17) k N | Ω (cid:105) (25)where m, σ, k i ∈ { , } . The reduced density matrix is computed as: (cid:104) σ, → k | ρ S | σ (cid:48) , → k (cid:48) (cid:105) = (cid:104) , σ, → k | ρ S + e | , σ (cid:48) , → k (cid:48) (cid:105) + (cid:104) , σ, → k | ρ S + e | , σ (cid:48) , → k (cid:48) (cid:105) = (26) (cid:104) , σ, → k | ρ S + e + a e ρ S + e a † e | , σ (cid:48) , → k (cid:48) (cid:105) We now follow the steps outlined above.(1) After step 1, the total density matrix after particle removal from site e is: ρ S + e, = (1 − n e ) ρ S + e (1 − n e ) + a e ρ S + e a † e (27)and the system density matrix is: ρ S, = ρ S . We can also see this is to explicitly writing (cid:104) σ, → k | ρ S, | σ (cid:48) , → k (cid:48) (cid:105) = (cid:104) , σ, → k | ρ S + e, + a e ρ S + e, a † e | , σ (cid:48) , → k (cid:48) (cid:105) = (28)= (cid:104) , σ, → k | ρ S + e, | , σ (cid:48) , → k (cid:48) (cid:105) = (cid:104) σ, → k | ρ S | σ (cid:48) , → k (cid:48) (cid:105) ⇒ ρ S, = ρ S (2) We now apply the evolution U θ . We have ρ S + e, = U θ ρ S + e, U † θ , and therefore: (cid:104) σ, → k | ρ S, | σ (cid:48) , → k (cid:48) (cid:105) = (cid:104) , σ, → k | U θ ρ S + e, U † θ + a e U θ ρ S + e, U † θ a † e | , σ (cid:48) , → k (cid:48) (cid:105) = (29) (cid:104) , , → k | a σ U θ ρ S + e, U † θ (cid:16) a † (cid:17) σ (cid:48) + a σ a e U θ ρ S + e, U † θ a † e (cid:16) a † (cid:17) σ (cid:48) | , , → k (cid:48) (cid:105) To compute the matrix elements, we use the following properties of U θ : U † θ | , , → k (cid:48) (cid:105) = | , , → k (cid:48) (cid:105) ; U † θ | , , → k (cid:48) (cid:105) = | , , → k (cid:48) (cid:105) (30)1and the transformation: U † θ a U θ = cos( θ ) a − sin( θ ) a e (31) U † θ a e U θ = cos( θ ) a e + sin( θ ) a (32)By commuting the U θ operators through the a e , a operators we can now express the newmatrix elements as function of θ . We find that: (cid:104) σ, → k | ρ S, | σ (cid:48) , → k (cid:48) (cid:105) = δ σ δ σ (cid:48) (cid:104) , , → k | ρ S + e, | , , → k (cid:48) (cid:105) + sin ( θ ) δ σ δ σ (cid:48) (cid:104) , , → k | ρ S + e, | , , → k (cid:48) (cid:105) +(33) δ σ δ σ (cid:48) cos ( θ ) (cid:104) , , → k | ρ S + e, | , , → k (cid:48) (cid:105) +cos( θ ) δ σ δ σ (cid:48) (cid:104) , , → k | ρ S + e, | , , → k (cid:48) (cid:105) + cos( θ ) δ σ δ σ (cid:48) (cid:104) , , → k | ρ S + e, | , , → k (cid:48) (cid:105) We can identify the transformation on ρ S as: ρ S, = (1 − n ) ρ θS (1 − n ) + sin ( θ ) a ρ θS a † + cos ( θ ) n ρ S n + (34)cos( θ ) n ρ S (1 − n ) + cos( θ )(1 − n ) ρ θS n Rearranging the terms we finally have: ρ S → ρ S, = (1 − (1 − cos( θ ) n )) ρ S (1 − (1 − cos( θ ) n )) + sin ( θ ) a ρ S a † (35)Identifying (cid:15) = 1 − cos θ , and noting that sin θ = (cid:15) (2 − (cid:15) ), we have recovered the map (22). III. NON-EQUILIBRIUM STEADY STATE EQUATION
There are a myriad possible processes described by combinations of the operations ( I − IV ). Here we concentrate on current generation processes as described by Eq. (2), involvesoperations I, III , IV resulting in the map: G → (1 − r ) u † Gu + ru † { α ((1 − (cid:15) a P a ) G (1 − (cid:15) a P a ) + (2 (cid:15) a − (cid:15) a ) P a ) + (36)(1 − α )((1 − (cid:15) b P b ) G (1 − (cid:15) b P b )) } u. This simple model allows for a substantial reduction of complexity from the full quantumproblem of describing the evolution of ρ into an evolution equation for the two point function G ij , which can be tractable by either analytical or numerical methods. It is clear at thisstage that we can access very interesting situations.2To compute the eventual non-equilibrium steady state for (36) it is convenient to view thetransformation on G from a point of view of a super-operator. Here the N × N matrix G is viewed as an N dimensional vector, and the action of the evolution L on ρ translates in(36) into: G → Λ G + g, (37)where Λ is an N × N matrix, and g is the inhomogeneous contribution due to the particleinjection processes (13), and corresponding to the term rα (2 (cid:15) a − (cid:15) a ) u † P a u in (36).In general, whenever g = 0, the long time behavior will be determined as usual by the largesteigenvectors of Λ. However when g (cid:54) = 0, the situation is somewhat different: Indeed, fromEq. (37), we see that when (1 − Λ) is invertible, there exists a unique stationarity G , thatmay be written in the form: G steady = (1 − Λ) − g (38)If Λ − G r = G r , it means that the evolution u has an invariant subspace which does not include the sites a, b . In this case one has towork with a generalized inverse of (Λ − G steady ∼ G r + (1 − Λ) − g . While inhomogenous equations are acommon occurrence in the study of steady states in classical driven systems, they are usedless in quantum processes, where evolution is unitary. A recent example of such a non-homogenous equation in a quantum context is the calculation of the expectation values ofspin components in the steady state of a spin undergoing periodic laser pulses [14, 15].We now apply these ideas to our current injection process described by (2) and (36). Per-forming the inversion in superoperator space as in (38) in general is a daunting task. Inthe limit of r (cid:28)
1, we were able to solve exactly for the degenerate perturbation theory tolowest order in r , obtaining for the energy distribution Φ the result (3). The derivation issomewhat lengthy and given in the next section.The A , B coefficients in (3) are given below. Define: A = α (2 − (cid:15) a ) (1 − α ) (cid:15) b (cid:15) a Q ab ((2 − (cid:15) a )(2 − (cid:15) b )+2 Q ab (cid:15) a (cid:15) b ( α (2 − (cid:15) a )+(1 − α )(2 − (cid:15) b ))) (39) B = α (2 − (cid:15) a ) (cid:15) a (2 − (cid:15) b +2 αQ ab (cid:15) a (cid:15) b )((2 − (cid:15) a )(2 − (cid:15) b )+2 Q ab (cid:15) a (cid:15) b ( α (2 − (cid:15) a )+(1 − α )(2 − (cid:15) b ))) where: µ k = 2 ( α(cid:15) a p a,k + (1 − α ) (cid:15) b p b,k ) ; Q ab = Σ k p a,k p b,k µ k r . How canwe understand this? Note that at r = 0, there are infinitely many steady states (any G such that [ G, h ] = 0). However, when r (cid:54) = 0, Λ stops being degenerate and it singles out aparticular direction of breaking the degenerate space of matrices. A. Steady state distribution: Derivation
Here we derive the formulas (3),(39) for the non-equilibrium steady state energy distributionΦ. We will study the steady state equation associated with the process (36), taking (cid:15) a , (cid:15) b = 1for simplicity, however the derivation with (cid:15) a , (cid:15) b (cid:54) = 1 follows along exactly the same lines. G steady = (1 − r ) u † G steady u + (40) u † rα ( P a + P a ⊥ G steady P a ⊥ ) u + u † r (1 − α )( P b ⊥ G steady P b ⊥ ) u where u = e − iτh .Below we label the eigenstates of h by n , h | n (cid:105) = E n | n (cid:105) , and would like to find theprobability to find a state with energy E n occupied in the steady state. This probability isgiven by Φ n ≡ T r ( ρa † n a n ) = (cid:104) n | G | n (cid:105) .For r = 0, all states where [ G, h ] = 0, are immediately invariant under time evolution.Therefore, in the limit of r (cid:28) G steady which isapproximately diagonal. Let us write, in the energy basis, the ansatz: G steady = diag( { Φ , ... } ) + rD, (41)where Φ n = (cid:104) n | G steady | n (cid:105) are the steady states occupations, and D is an off-diagonal matrixin energy space. Eq. (40) becomes:Φ + rD = (1 − r )Φ + (1 − r ) ru † Du + rαu † P a u + (42) rαu † ( P a ⊥ Φ P a ⊥ ) u + u † r (1 − α )( P b ⊥ Φ P b ⊥ ) u + O ( r )We note that the zeroth order is eliminated and we wind up with: D = − Φ + u † Du + αu † P a u + αu † ( P a ⊥ Φ P a ⊥ ) u + u † (1 − α )( P b ⊥ Φ P b ⊥ ) u (43)4Furthermore, note that both D, u † Du are off-diagonal in energy. Therefore we have a closedequation for the diagonal elements:0 = − Φ n + αp a,n + α ( P a ⊥ Φ P a ⊥ ) nn + (1 − α )( P b ⊥ Φ P b ⊥ ) nn . (44)Explicitly,( P a ⊥ Φ P a ⊥ ) nn = (Φ − P a Φ − Φ P a + P a Φ P a ) nn = Φ n − p a,n Φ n + Σ l P a, nl Φ l P a, ln (45)where we have denoted p a,n = (cid:104) n | P a | n (cid:105) (and similarly p b,n = (cid:104) n | P b | n (cid:105) ) and P a,ln = (cid:104) l | P a | n (cid:105) .Note that using Σ l P a, nl P a, ln = p a,n we can write Eq. (44) as:0 = αp a,n − Φ n ( αp a,n + (1 − α ) p b,n ) + Σ l (Φ l − Φ n )( αP a, nl P a, ln + (1 − α ) P b, nl P b, ln ) . (46)At this point it is possible to argue that on the right, | Σ l (Φ l − Φ n )( αP a, nl P a, ln + (1 − α ) P b, nl P b, ln ) | is small, giving us a first guess for the answer:Φ n ∼ αp a,n α p a,n +(1 − α ) p b,n (47)However, as we see below, it is possible to do better and solve equation (44) exactly withoutthis condition. To do so notice that: P a, nl P a, ln = |(cid:104) n, a (cid:105)| |(cid:104) a, l (cid:105)| ≡ p a,n p a,l (48)Going back to (44) we write it as:0 = αp a,n − αp a,n + (1 − α ) p b,n )Φ n + Σ l ( αp a,n p a,l + (1 − α ) p b,n p b,l )Φ l (49)We rewrite the equation as an in-homogenous linear equation: Q → Φ = αZ F → F + V → Φ . (50)Here → F is a unit vector defined by: → F = p a,n Z F ; Z F = (cid:112) Σ n p a,n , (51) Q is a diagonal matrix Q nm = δ nm √ µ n ; µ n = 2( α p a,n + (1 − α ) p b,n ) , (52)5and V can be written in the form V nm = αp a,n p a,m + (1 − α ) p b,n | g m | = αZ F | F (cid:105)(cid:104) F | + (1 − α ) Z G | G (cid:105)(cid:104) G | . (53)The solution is given formally by:( Q − V ) → Φ = αZ F → F = ⇒ → Φ = Q − V αZ F → F = αZ F Q − −Q − V Q − Q − → F . (54)Next, we define the unit vector | F Q (cid:105) as | F Q (cid:105) = Z − Q − | F (cid:105) , ; Z = Σ n p a,n µ n Z F . (55)Note the normalization (cid:107) F Q (cid:107) = 1. Similarly we define | G Q (cid:105) = Z − Q − | G (cid:105) , ; Z = Σ n p b,n µ n Z G . (56)Using these, (54) is expressed as: → Φ = α Q − ( Z F Z FQ − αZ F Z | F Q (cid:105)(cid:104) F Q |− (1 − α ) Z G Z | G Q (cid:105)(cid:104) G Q | ) | F Q (cid:105) (57)In the next step we use the following relation: a | v (cid:105)(cid:104) v | + b | u (cid:105)(cid:104) u | | v (cid:105) = a + b + ab (1 −|(cid:104) v,u (cid:105)| ) { (1 + b ) | v (cid:105) − b (cid:104) u, v (cid:105)| u (cid:105)} , (58)which holds for normalized vectors || u || = || v || = 1. We are not aware if the expression(58) appears in the literature, but it can be verified explicitly by multiplying both sides by(1 + a | v (cid:105)(cid:104) v | + b | u (cid:105)(cid:104) u | ).We will use (58) on (57), with | F Q (cid:105) , | G Q (cid:105) playing the role of | u (cid:105) , | v (cid:105) . Thus, we take in (58): a → − αZ F Z ; b → − (1 − α ) Z G Z , (59)and c ≡ (cid:104) F Q | G Q (cid:105) = Σ n Z G Z F Z FQ Z GQ p b,n p a,n µ n = Z G Z F Z FQ Z GQ Σ n p b,n p a,n µ n (60)noting Z F Z FQ = (cid:113) Σ n p a,n µ n ; Z G Z GQ = (cid:113) Σ n p b,n µ n (61)6we have c = (cid:114) (Σ l p a,lµl )(Σ l p b,lµl ) Σ n p b,n p a,n µ n (62)Using these expressions with (58) and (57) we find:Φ n = αZ F Z FQ √ µ n (( − αZ F Z | F Q (cid:105)(cid:104) F Q |− (1 − α ) Z G Z | G Q (cid:105)(cid:104) G Q | ) | F Q (cid:105) ) n = αZ F Z FQ √ µ n − αZ F Z − (1 − α ) Z G Z + αZ F Z (1 − α ) Z G Z (1 −| c | ) ×{ (1 − (1 − α ) Z G Z ) (cid:104) n | F Q (cid:105) + (1 − α ) Z G Z c ∗ (cid:104) n | G Q (cid:105)} = αµ n (1 − (1 − α ) Z G Z ) p a,n +(1 − α ) Z F Z FQ Z G Z GQ c ∗ p b,n − αZ F Z − (1 − α ) Z G Z + αZ F Z (1 − α ) Z G Z (1 −| c | ) . Denoting Q aa = Z F Z = Σ l p a,l µ l ; Q bb = Z G Z = Σ l p b,l µ l ; Q ba = Σ l p a,l p a,l µ l , (63)we find that: Φ n = αµ n (1 − (1 − α ) Q bb ) p a,n +(1 − α ) Q ba p b,n − αQ aa − (1 − α ) Q bb + α (1 − α )( Q aa Q bb − Q ba2 ) . (64)As a final simplification we note that:2 αQ ba + 2(1 − α ) Q bb = Σ l p b,l (2 αp a,l + 2(1 − α ) p b,l )2 ( α p a,l + (1 − α ) p b,l ) = Σ l p b,l = 1 , (65)and similarly we have: 2 αQ ba + 2(1 − α ) Q bb = 1. Using these relations in Eq. (64), we canexpress the final result in terms of Q ab alone, finding:Φ k = A + B p a,k p b,k (1 − α ) + α p a,k p b,k ; A = 2(1 − α ) αQ ba Q ba , B = α (1 + 2 αQ ba )1 + 2 Q ba . (66)as mentioned, the derivation with (cid:15) a , (cid:15) b (cid:54) = 1 follows exactly the same line, giving the coeffi-cients (39). B. Kossakowski-Lindblad limit
In the Kossakowski-Lindblad limit, the treatment is considerably simpler. Starting with:˙ ρ = 1 i (cid:126) [ H, ρ ] + γ a ( a † a ρa a − { ρ, a a a † a } γ b ( a b ρa † b − { ρ, a † b a b } G is: ˙ G = − i (cid:126) [ G, h t ] − { γ a P a + γ b P b , G } + γ a P a . (67)The steady state obeys ˙ G steady = 0, we again set G steady = Φ + rD where Φ is diagonal and D is strictly off diagonal in energy, and assume that r → γ a and γ b are approachingzero. We take a diagonal matrix element of the equation to find, in lowest order in r thatΦ k ( γ a P a + γ b P b ) kk = γ a ( P a ) kk ⇒ Φ k ≡ γ a p a,k γ a p a,k + γ b p b,k . (68)Setting γ a = rα(cid:15) a , γ b = r (1 − α ) (cid:15) b , as representing the appropriate rates in the processdescribed in (2) we recover (3). IV. EXAMPLES OF DYNAMICS AND DEPENDENCE ON INITIALCONDITION
The dependence of the dynamics on the initial condition is of interest by itself. While inFig. 1, we started the evolution from the vacuum state, in Fig. 3, we describe such a processwhere the system is started off as the ground state of H hop . The evolution happens in stages.In the initial stage of evolution we observe two shock wave fronts: one propagating with aregion of reduced density from the right, collides with a front of enhanced density propagatedfrom the left. It is interesting to note that the evolution is on a faster time scale than thespeed of propagation of a wave-packet localized at a point by free evolution. In the contextof classical non equilibrium processes, shock waves have been described for the asymmetricexclusion process in e.g.[16] (It is possible to use the present system also to describe suchsituations, however this will be done elsewhere).As the fronts collide the imbalance between the left and right sides of the chain startsto decrease. Finally, soliton like density packets of different velocities, are observed atlonger time scales, and may be related to the soliton described in [17] in the context of theorthogonality catastrophe. It is interesting to note the injected particles traveling from theleft travel with faster velocities compared to their partners from the other side.In Fig. 3 we show the average particle density ¯ n ≡ N − T rG . One of the interesting featuresobserved is a qualitative change in the slope of ¯ n ( t ) around 350 iterations. This changeseems to correspond to the annihilation of the high density front coming from the left. To8 n n r Density Depletion
Figure 3: Density depletion in a system where particles are extracted from the right at higherrate than injected on the left, here r = 1 and α = 0 .
3, initial state is the half filled ground state of H hop . Left: Real space evolution of local density. Right: Evolution of space averaged density. check this behavior, we consider, in Fig. 4 the evolution when the initial stage is asymmetricitself: Here in the initial stage all sites i on the left, i < i >
100 are occupied. This state evolves through four fronts that collide andeventually annihilate. Note that for coherent evolution from such an initial state, it hasbeen shown that the front propagation has a scaling 1 /t [18]. In the context of evolution9 n r Density Depletion
Figure 4: Density evolution under the same dynamics as Fig. 3, however with a domain wall as aninitial state. Here the depletion happens in two steps, but eventually reaches the same asymptoticvalue (¯ n ∼ .
38) as in Fig. 3. of magnetization in a spin chain the evolution of initial domain wall was studied in [19].Comparing the density evolution in Fig. 3 and Fig. 4, we see that there is a transient be-havior associated with the different nature of the initial states, and their stages of evolution.In Fig. 4, there is a noticeable change in depletion rate around 100 and 300 iterations, thefirst kink corresponds the initial high density region on the right hitting the left side: at thatpoint injection of particles becomes harder for a while and | ∂ t ¯ n | decreases until the densitygoes down enough on the left. The second kink is observed when the high density region isreflected back to the right: extracting particles on the right is then easier and | ∂ t ¯ n | grows.At long times the density seems to decay asymptotically as 1 /t towards the non-equilibriumsteady state density.0 V. SUMMARY
We presented a class of non equilibrium quantum processes that correspond to closed hier-archies of evolution equations, and can thus be studied numerically efficiently. We have usedthis idea to explore non-equllibrium generation of currents and approach to steady states.We remark that the resulting states may also be viewed as Floquet states, and we have thussupplied a particular way of engineering such states, that may be of interest in the contextof topological Floquet states[20–22] and generation of topological states via dissipation [23].Moreover, the energy distribution Φ k should be studied further: one can hope to test theresulting highly excited current carrying steady states in a variety of settings from coldatoms to mesoscopic systems and spin chains. We emphasize that our result does not relyon integrability in the sense of Bethe Anzats that is useful in one dimension and has beenused in studies of dissipative spin chains. Thus, our treatmentis available for periodicallydriven fermion systems that do not correspond to spin chains, and most importantly, tohigher dimensional systems. Acknowledgement