Solving real time evolution problems by constructing excitation operators

In this paper we study the time evolution of an observable in the interacting fermion systems driven out of equilibrium. We present a method for solving the Heisenberg equations of motion by constructing excitation operators which are defined as the operators A satisfying [H,A]=\lambda A. It is demonstrated how an excitation operator and its excitation energy \lambda can be calculated. By an appropriate supposition of the form of A we turn the problem into the one of diagonalizing a series of matrices whose dimension depends linearly on the size of the system. We perform this method to calculate the evolution of the creation operator in a toy model Hamiltonian which is inspired by the Hubbard model and the nonequilibrium current through the single impurity Anderson model. This method is beyond the traditional perturbation theory in Keldysh-Green's function formalism, because the excitation energy \lambda is modified by the interaction and it will appear in the exponent in the function of time.


I. INTRODUCTION
The interacting fermion systems have a central position in modern condensed matter physics. 1 Their properties have been intensively studied both in theory and experiment. Good agreements between the theoretical predictions and the experimental results have been obtained when these systems are in thermal equilibrium. Examples include the Coulomb blockade effect and the Kondo effect in quantum dots, [2][3][4][5][6][7][8][9] and the charge density wave in the low-dimensional materials, e.g. the carbon nanotubes and the edge states in fractional quantum Hall effects. [10][11][12][13][14][15] However, how to understand the real time evolution of an interacting system driven out of equilibrium remains a great challenge in spite of intense efforts in recent years. The standard tool for the quantum field theory of nonequilibrium states is the Keldysh-Green's function technique. But the Keldysh techniques suffer from several shortcomings. First it is based on the Wick's theorem, and then the initial states must be limited to the equilibrium states of the quadratic Hamiltonian. Secondly, the result from a perturbative expansion with respect to the interaction strength does not incorporate the physics of the problem when the Coulomb interaction is the largest energy scale. Finally, for some models, e.g., the Kondo model, the perturbation theory is plagued by the infrared divergence at low temperatures.
For these reasons, new theoretical tools have been developed to understand the interacting fermion systems in nonequilibrium. Numerical methods are developed like the real time Quantum Monte Carlo, [16][17][18][19] the time-dependent numerical renormalization group, 20-22 the scattering state numerical renormalization group, 23 the adaptive time-dependent density matrix renormalization group [24][25][26] and the non-equilibrium dynamical mean field theory. 27 The numerical method has the advantage that the result can be obtained in very high precision and the parameters of the model can be chosen arbitrarily. However, they also suffer the inconvenience that the large computer resources are required and the result does not incorporate the physics of the problem directly. The analytical methods are developed, which can supplement the shortcomings of the numerical methods, such as the scattering Bethe-ansatz, 28 the method based on integrability, 29, 30 the realtime renormalization group, 31-34 the nonequilibrium flow equation, 35-39 and various approximation  schemes building on the Green's function techniques. 40,41 However, no unique method is available which can cover all regimes of interest in different models. It is still necessary to study new methods.
The Hamiltonian of a typical interacting fermion system can be expressed aŝ whereĉ k (ĉ † k ) is the fermionic annihilation (creation) operator. And the normal ordering is with respect to the non-interacting ground state of the Fermi sea. In the traditional diagrammatic expansion the time evolution of a physical observable, denoted as Ô (t) whereÔ(t) satisfies the Heisenberg equation of motion is calculated by expanding the S-matrix and summing up a series of diagrams in a self-consistent way. In this paper we propose an alternate way to solve the Heisenberg equation of motion by decomposing the observable operatorÔ into the linear combination of the excitation operators which are defined as the operatorÂ satisfying the eigen equation This method circumvents the expansion of the S-matrix and then the contraction of field operators, and is applicable even the initial state is nontrivial. Furthermore, as we will show, this method is beyond the traditional perturbation theory because the excitation energy λ which appears in the exponent in the time function is modified by the interaction strength.
In Sec. II we give the definition of the excitation operators. In Sec. III we demonstrate how to construct the excitation operator and use it to calculate the time evolution of the creation operator in a toy model inspired by the Hubbard model. And the comparison between this model and the nonequilibrium flow equation approach which also contributes to solve the Heisenberg equation of motion is given. In Sec. IV we perform our method in single impurity Anderson model and calculate the evolution of the current operator. Sec. V is a short summary.

II. EXCITATION OPERATORS
The excitation operator 42 of a HamiltonianĤ is the operatorÂ which satisfies where λ is a real number denoting the excitation energy. We suppose that |ψ x and |ψ y are two arbitrary eigen-states ofĤ with the eigen-energies E x and E y respectively. Then the excitation operator can be expressed asÂ The corresponding excitation energy is λ x, y = E x − E y . An observable operatorÔ can be decomposed into the linear combination of the excitation operators, i.e., where O x,y = ψ x |Ô|ψ y . According to Eq. (4), we have e iĤtÂ x,y e −iĤt = e iλ x,y tÂ x,y . 012194-3 Pei Wang AIP Advances 2, 012194 (2012) The solution of the Heisenberg equation of motion can be expressed aŝ Once we know how to decompose an observable operator, we obtain the time evolution of it. In many-body problems it is generally difficult to find the eigen-states of the Hamiltonian and then calculate the excitation operator according to Eq. (5). An alternate way is to expandÂ into a series of products of field operators, substitute the expansion into Eq. (4) and decide the coefficients. For the Hamiltonian of Eq. (1), we could expressÂ aŝ where M k , M k k and M k 1 k 2 k 1 are the undetermined coefficients. The commutator of the Hamiltonian with the products of N field operators contains the products of (N + 2) operators, which indicate thatÂ should be an infinite series. We need to truncate the series to get an approximated expression ofÂ.
Since we employ the perturbative expansion in the expression ofÂ, our method is also a kind of perturbative approach. While the direct perturbative approach for solving the Heisenberg equations of motion will generate secular terms, which are uncontrolled at large timescales. 36,37 By first decomposing the observable into the linear combination of the excitation operators, we avoid the secular divergences and make the results about the long-time behavior more explicit.

A. Excitation operators
In this section, we use a toy model Hamiltonian to demonstrate how to construct excitation operators and use them to calculate the evolution of an observable operator. In Eq. (10) the first term denotes the single particle energy, and the second and third terms the interactions. Here σ = ↑, ↓ denotes the particle spin, andσ the opposite to σ . The normal ordering is with respect to the Fermi sea with chemical potential μ. The Fermi distribution function is denoted as n k = θ (μ − k ). It is not difficult to find that this model Hamiltonian is inspired by the Hubbard model. In fact it is the simplified Hubbard model in which most of the interaction terms are thrown off. Only the degenerate interaction term is left which is defined as the interaction keeping the sum of the single particle energies of the two annihilation operators exactly as same as that of the two creation operators. For example, the degenerate interaction in the Hamiltonian of Eq. (1) should be the interaction terms We are interested in calculating e iĤtĉ † kσ e −iĤt . We notice thatĉ † kσ increases the particle number of a state by one. Then it will be decomposed into the excitation operators in which arbitrary term contains one more creation operator than the annihilation operator. We use the symbolÂ kσ to represent the excitation operators. Again we could suppose that the first term ofÂ kσ isĉ † kσ in the simplest case. The commutator ofĤ withĉ † kσ contains the terms :ĉ † So we could suppose the expression of the excitation operator aŝ In Eq. (11) we neglect the high order terms consisting of more than three field operators. The eigen-equation is written as where λ k is the excitation energy. Substituting Eq. (10) and (11) in and comparing the coefficients beforeĉ † kσ between the left and right sides of the equation, we get But −n 2 k + n k = 0 whether n k = 0 or n k = 1. Then the excitation energy becomes Comparing the coefficients before :ĉ † kσĉ † k σĉ k σ :, :ĉ † k σĉ † kσĉ k σ : and :ĉ † kσĉ † k σĉ k σ : between the two sides of the eigen-equation for k = k , we get These equations are solved and the coefficients M are decided. When n k = n k we find that and M 3 k,k is an arbitrary number. What we want is a set of linearly-independent excitation operators. So we set M 3 k,k = 0, since :ĉ † kσĉ † k σĉ k σ : is itself an excitation operator as will be shown next. Eq. (16) has two linearly-independent solutions: M 1 k,k = 0, M 2 k,k = 1 2n k −1 ; and M 1 k,k = 1 2n k −1 , M 2 k,k = 0. One could choose arbitrary one forÂ kσ . We use an indicator function χ k,k = 0, 1 to denote our choice, and express the corresponding term inÂ kσ as k σĉ k σ : For n k = n k , the solution of Eq. (15) is In summary the excitation operator can be expressed aŝ k σĉ k σ : Different indicator functions will give differentÂ kσ . But we will show next that the choice of χ k,k does not affect the last result of e iĤtĉ † kσ e −iĤt . To decomposeĉ † kσ we need to cancel the terms : c † kσ c † k σ c k σ :, :ĉ † k σĉ † kσĉ k σ : and :ĉ † kσĉ † k σĉ k σ : inÂ kσ . So we construct the excitation operators with them as the leading terms. To avoid ambiguity we use the symbolB to denote these operators. For arbitrary k = k we suppose that where we again neglect the high order terms consisting of more than three field operators. This supposition is reasonable because the commutator [Ĥ ,B kk σ ] does not contain the terms like single creation operator. The eigen-equation is now written as By comparing the coefficients between the left and right sides of the equation we have In the matrix form of this equation we find that (N 1 k,k , N 2 k,k , N 3 k,k ) is in fact the eigenvector of a (3 × 3) matrix. When n k = n k Eq. (22) can be expressed in the matrix form as ⎛ This matrix has three linearly-independent eigenvectors. Correspondingly we get three linearly-independentB kk σ . When n k = n k Eq. (22) becomes ⎛ where we have used the identity n k + n k = 1. This matrix contributes another three eigenvectors.
Totally we obtain six types ofB kk σ , which are all listed in table I with the corresponding excitation energies.
In summary we get the excitation operatorsÂ kσ andB kk σ . We could also get the excitation operators of higher orders by similar analysis. Even we only keep the lowest order terms in the perturbative expansion ofB kk σ , we get a nontrivial result. The interaction strength U enters the expression of the excitation energy ofB kk σ . We conclude thatB kk σ reflects the characteristics of the collective excitations. When n k = 0 we find thatB 1 k σ )ĉ k σ and the corresponding excitation energy is k + 2U. This operator annihilates an electron and simultaneously creats a spin singlet. This procedure costs an extra energy of 2U.

B. The evolution of single creation operator
The operatorĉ † kσ can be decomposed into the linear combination of the operatorsÂ kσ andB kk σ . The second order terms inÂ kσ with n k = n k can be canceled byB 1 kk σ ,B 2 kk σ andB 3 kk σ , and those with n k = n k byB 4 kk σ ,B 5 kk σ andB 6 kk σ . It is easy to find The expressions of e iĤtÂ kσ e −iĤt and e iĤtB kk σ e −iĤt are got from Eq. (7) and the excitation energies listed in table I. Substituting back in the expression ofÂ kσ andB kk σ we get Note that the expression of e iĤtĉ † kσ e −iĤt is unique, independent to the indicator function. In Eq. (26) the function of time is U-dependent, which indicates that our method involves the renormalization of the excitation energies and is beyond the traditional perturbation theory. This equation is an operator identity. As the operatorĉ † kσ (t) acting on e iĤt |ψ where |ψ is an eigenstate, we get e iĤtĉ † kσ |ψ . So the physical meaning of Eq. (26) is that it shows how the system responses to a single particle excitation by producing the collective excitations.
It is worthwhile to mention that the nonequilibrium flow equation 36,37 is also a method designed for solving the Heisenberg equation of motion by expanding the operator into a power series. There are three main difference between our method and the flow equation approach. First in the flow equation approach a unitary transformation is performed to the Hamiltonian, while in our method the Hamiltonian keeps invariant. Secondly, we turns the problem into the diagonalization of the coefficient matrix instead of solving the differential equation in the flow equation approach. Finally we deal with the degenerate interaction in the same way as the non-degenerate ones, which suggests that this method is applicable in the model where the degenerate interaction is overwhelming, e.g. in the Tomonaga-Luttinger model. 15 While in the flow equation approach the degenerate interaction has a different scaling behavior in the flow.

IV. SINGLE IMPURITY ANDERSON MODEL
Next we consider the single impurity Anderson model, because it is important both in theory and experiment. Furthermore, many papers have contributed to study this model, so that we could compare our method with the others. We will use the excitation operators to calculate the time evolution of the current through the impurity.
The Hamiltonian of the Anderson impurity model is expressed aŝ where σ = ↑, ↓ denotes the electron spin and α = L, R the left and right leads. Hereĉ † kασ andd † σ are the single-electron creation operators in the leads and at the impurity site respectively. At initial time the leads are in thermal equilibrium with chemical potentials ± V sd 2 respectively. The coupling between leads and the impurity is switched on at time t = 0.
In the pre-diagonalization basis, 43 where B s 1 s 2 s 1 s 2 is the abbreviation of B s 1 B s 1 B s 2 B s 2 and the normal ordering is with respect to the initial state. The (anti)symmetric operator is defined asĉ k±σ = 1 √ 2 (ĉ kLσ ±ĉ k Rσ ) and the hybridization operator asĉ sσ = k V s − k B sĉk+σ + B sdσ . Here the coefficient B s = V √ 2 s + 2 and the linewidth = ρπV 2 , where ρ is the density of states of the lead. In the Anderson impurity model we will neglect the degenerate interaction by ordering B s 1 s 2 s 1 s 2 = 0 as s 1 + s 2 = s 1 + s 2 . This is reasonable in the Anderson impurity model, since the coefficient B s → 0 in the thermodynamic limit. The Fermi function in hybridization basis has sharp edges in the thermodynamic limit, i.e., c † s σ c sσ 0 = δ s,s n s , where n s = 1 2 (θ ( V sd 2 − s ) + θ (− V sd 2 − s )).