New Neumerical Method to Calculate Time-Dependent Quantum Properties in Finite Temperature Based on the Heisenberg Equation of Motion

For the purpose of computer calculation to evaluate time-dependent quantum properties in finite temperature, we propose new numerical method expressed in the forms of simultaneous differential equations. At first we derive the equation of motion in finite temperature, which is found to be same expression as Heisenberg equation of motion except for the c-number. Based on this equation, we construct numerical method to estimate time-dependent physical properties in finite temperature precisely without using analytical procedures such as Keldysh formalism. Since our approach is so simple and is based on the simultaneous differential equations including no terms related to self-energies, computer programming can be easily performed. It is possible to estimate exact time-dependent physical properties, providing that Hamiltonian of the system is taken to be a one-electron picture. Furthermore, we refer to the application to the many body problem and it is numerically possible to calculate physical properties using Hartree Fock approximation. Our numerical method can be applied to the case even when perturbative Hamiltonians are newly introduced or Hamiltonian shows complex time-dependent behavior. In this article, at first, we derive the equation of motion in finite temperature. Secondly, for the purpose of verification and of exhibiting the usefulness, we show the derivation of gap equation of superconductivity and of sum rule of electrical conductivity and the application to the many body problem. Finally we apply this method to these two cases: the first case is most simplified resonance charge transfer neutralization of an ion and the second is the same process but impurity potential is newly introduced as perturbative Hamiltonian. Through both cases, it is found that neutralization process is not so sensitive to temperature, however, impurity potential as small as 10 meV strongly influences the neutralization of ion.


Introduction
In the previous work [1], we proposed a new and numerical method so as to evaluate various time-dependent properties and to extend the theoretical forms easily even when Hamiltonian shows complex time-dependence or perturbative Hamiltonian is newly introduced.Our proposed method is based on the Heisenberg equations of motion as shown below where operator i is in Heisenberg representation.By using Equation (1.1), the differentiation of can be expressed as By taking the expectation value of each term, we obtain the following differential equation: However, as shown above, theoretical treatments are restricted to the ground state, i.e., T = 0 K; thus various time-dependent physical properties in finite temperature S.-I.KONDO 1538 have remained unresolved.
As to the theoretical evaluation of time-dependent properties in finite temperature, linear response theory [2] is an excellent method and can be applied to many cases such as calculations of electrical conductivity and magnetic susceptibility under the oscillating field and it is proved that this method gives exact solutions within the range of linearity.Usually time-dependent perturbative is given in the form of separation of variables in the linear response theory: , can be easily deduced as where BA g t

   
BA is a retarded thermal Green function and defined by and As shown in Equations (1.5)- (1.8), if a perturbative Hamiltonian can not be decomposed into the form of separation of variables of Equation (1.4) or can show complex time-dependent behaviours, it seems to be a complicated task to evaluate time-dependent physical properties av B t   R because of analytical and/or numerical difficulty in the estimation of retarded thermal Green function BA g t and of analytical corrections concerning to Equations (1.5)- (1.8).In addition, linear response theory can not be applied unless linearity is observed.
Although Green function methods are excellent and smart in analyzing the various time dependent phenomena, integral schemes are essential in diagram methods (Dyson equation), which easily leads to numerical and/or analytical difficulties and awkwardness for the explanation of experimental data.For example, let us consider the simplest case of atom-surface collision where the electron transfer matrix between surface and an ion and energy level of the ion shows strong time-dependence due to the change in atom-surface distance, we can write down the time-dependent Newns-Anderson Hamiltonian as below  k k (1.11) where z is an atom-surface distance with showing timedependence ( v: velocity of ion); k is the energy of a conduction electron with momentum k; k and k are creation and annihilation operators of a conduction electron with momentum k; a denotes the energy level of an ion and usually depends on the z; a and a C are creation and annihilation operators of a state of the ion, respectively; and ak is the electron transfer matrix element from the conduction electron k to a state of the ion and can be expressed as a function of z.Keldysh formalism is usually adopted in treating such a non-equilibirum problem and Dyson equation is given in the following form: ,  is retaded Green function of ion and is expressed as  are retarded Green function of non-perturative state of ion and self-energy, respectively.

 
, , The number of electron on the ion at t , i.e., , which we should seek and calculate, is In above equation, n k denotes the number of electron occupying the band state k at initial state and is given As shown in the series of Equations (1.12) to (1.21), many calculations and procedures are necessary to evaluate , which suggest that analytical efforts to evaluate are usually failed except for a few cases.Since energy state of ion , which is experimenttally determined, is considered complex, we can easily guess analytical form of can not be obtainable and numerical expression of is only possible.
Additionally analytical procedures to treat the term   Usually the term of 1 2 of Equation (1.23) can not be evaluated in the analytical form except for the specific cases such as or , therefore, the procedure on the basis of (1.12) seems too complicated even if perturbative method is used.Consequently awkward and complex analytical and/or numerical schems for calculating integral parts and a lot of approximations ignoring experimental condi-tions are essentially required to evaluate .Actually, theoretical solutions are obatined in the limited conditions such as wide band limit.
Furthermore, if there exists impurity atoms or dislocations on surface, we should take account of the presence of impurity potential  as shown below When such a perturbative Hamiltonian is intoroduced into the system, it should be noted that the methods on the basis of Keldysh formalism seems almost impossible for evaluating   n  precisely because of a lots of awkward analytical and/or numerical tasks and approximations ignoring experimental conditions.
Concerning the other approaches, we can mention the works performed by Brako and Newns, [3,4].Starting from the equation of motion method proposed by Bloss and Hone [5], they performed calculations by solving the simultaneous differential equations while regarding operators (Q-number) as c-number.Based on these approaches, their method are expressed in the integral forms and exact solutions can be attainable in the case of wide band limit with assuming .Actually experimental data on the polarisation of the light emitted in the electronic transition of H atoms scattered on Ni(111) surface [6] were theoretically interpreted and examined by their method [7].
However, their method also seems to meet the numerical and analytical difficulties and nuisances when perturbative Hamiltonian is introduced or the system shows complicated time-dependence because of integral schemes being involved.Actually, their method can not be applied to the case when impurity potential as shown in Equation (1.24) is introduced into the system.
Concerning these theoretical methods as stated above, the main reason for complicating calculations is that theoretical expressions are given not in the differential forms but in the integral ones.In this article, therefore, for the purpose of constructing the method which can be applicable to experimental analysis, we propose a new numerical methods on the basis of differential forms.At first we derive the extension of Heisenberg equation of motion to finite temperature.Then, based on the such an extension, we construct simultaneous differential equations which can evaluate time-dependent physical propav B t erties in finite temperature.Our method can be applied to the case when a perturbative Hamiltonian can show so complex time-dependent behaviours that analysis on the basis of linear response theory and/or of Green function methods seem to meet analytical difficulties and nuisance and even when perturbative Hamiltonian is newly introduced.
In Section 2, we show the derivation of theoretical formalism and verify our proposed method for two cases: superconductivity and electrical conductivity.Furthermore we refer to the application to the many body problem using Hatree Fock approximation.After confirming our calculations, in Section 3, we apply our method to the two cases of resonant charge transfer processes, in which time-dependent properties play important roles on deciding dynamics.The first case is the most simlified resonant charge transfer process that a single charged ion is neutralzed in front of metal.During calculating, to simplify the discussion we assume electron transfer matrix remains constant.However it is sufficiently possible to estimate physical properties under the condition that electron transfer matrix shows k dependence.In the second case, we numerically treat the same process when impurity potential is introduced as perturbative Hamiltonian.The second case is more complicated than the first one and is usually considered as a difficult and awkward task to analyze on the basis of conventional ways such as Keldysh formalism because perturbative Hamiltonian often complicates the scheme of the evaluation of self energy.Furthermore, temperature dependence of resonant charge transfer process under the presence of impurity potential remains unresolved.Finally, in Section 4 we conclude this paper and discuss the remaining problems and further developments.

Derivation of Theoretical Framework
Firstly, let us start from the Neumann equation, which is given by where operator  ˆt  is a time-dependent density operator.Additionally we set 0  (chemical potential) = 0 hereafter for simplifying the discussion.Then time-dependent physical properties By taking the the differentiation of the above equation, we obtain Substituting the Neumann equation of Equation (2.1) into the above equation, we get , we obtain the following formula with respect to the the dif- It should be easily noticed that Equation (2.5) takes the same form as Heisenberg equation of motion in Equation (1.1) except for the c-number.The above equation corresponds to the extension of Heisenberg equation of motion to finite temperature.Equation (2.5) directly leads to the fact that time-dependent physical properties in finite temperature are expressed not in the integral forms as Dyson equations but in the differential forms (simultaneous differetial equations); this, therefore, indicates analytical and/or numerical nuisances will be expected to greatly decrease even if Hamiltonian includes complex terms and/or perturbative Hamiltonian is introduced.

Derivation of Gap Equation in Superconductivity
Next, in order to verify the above theoretical results, let us apply the above equation to the superconductivity.we consider the following Hamiltonian: (2.6) where In the above equations, k is the energy of a conduction electron with momentum k.
Using mean field approximation and Equation (2.10), Hamiltonian Equation (2.6) can be reduced to oneelectron picture and given by Next, we take operator as Accordingly, is Then, based on the Equation (2.5), the differential In the limit of , Consequently from Equation (2.14) where (2.17) By combining the above relation with Equation (2.10), we obtain Next, using Bogoliubov transformation to the Hamiltonian of Equation (2.11)  Then Hamiltonian can be easily diagonalized and is given by .Then   . Consequently, the term of By substituting the above result into Equation (2.18), we finally obtain gap equation.

Calculation of Electrical Conductivity
Next, we apply the Equation (2.5) to the electrical conduction.At first, we assume the following Hamiltonian: x y z j j j j j j j E E E x y z We obtain j , the time differential form of , from the relation where Since we focus on the electrical conduction, we take operator as an electric current density along  direction, 1 1 , where is a volume of the system.Let us consider the case of magnetic field , i.e., , then we obtain the following from Equations (2.27) and (2.29), , .
By applying the above result to Equation (2.5), we obtain the following differential equation with respect to the electric current density along  direction: It should be noted that the first term in the right hand expresses the force related to the potential while the second expresses term the force related to the applied electrical field.For example, let us consider the simplest case where  is a velocity of electron along the where In the above equation, the term of  indicates the force due to the potential V; thus Equation (2.37) corresponds to the classical equation of motion.
Next, let us back to the Equation (2.34).We solve this equation under the initial condition of The solution of Equation (2.34) can be easily evaluated as below If we focus on the term of in the right hand of the above equation, the exact expression of this term is:

 
The density operator     f of Equation (2.39) can be extended into the series of the power of  .
Ĉonsequently,using the above extension of   obtain the following result when     Considering the case where there is no pertu H bative amiltonian of 1 Ĥ , we can easily conclude no electrical current, which m s From the above result, even if pertubative Hamiltonian Since l in the Equation (2.44) co er rresponds to the powof Hamiltonian 1 Ĥ , it is sufficient to take 1 l  within the range of linearity, i.e., the first order of Hailtonian 1 ˆm H . Accordingly, on the basis of linear response theory, is Consequently, within the range of linearity, Finally we obtain the sum rule (2.50)

Application to Many Body Problem
, let us lim i Concerning the application to many body problem consider the follwing Anderson Hamiltonian: It should be noted that our method is quite sim cause these non-linear simulataneous differential equations do not have any terms related to comlicated selfenergies or vertex corrections, which should be essentially involved in considering the electrons correlation.
Next, let us discuss the above approxima view points of mean field approximation (H approximation).Firstly, we define ˆd n ple betion in the atree Fock as where ˆd n  denotes average value of ˆd n  .Then, ig- , we obtain the following Hamiltonian by substituting the above equation into Coulomb repulsive term of Equation (2.53):   [1], a resonant electron charge transfer (RCT) between a metal surface and an atom is interesting and important event , which Newns first proposed and discussed to explain the chemical adsorption in 1969 [8] on the basis of Anderson-Newns Model.In considering the various surface phenomena, RCT have played important roles to determ lecular dissociation and chemisorption.In the theoreti mework to trea the analysis on time-dependent Anderson-Newns model (TDAN) have been reported and discussed extensively to account for the various experimental results.The evaluation of time-dependent properties, therefore, are essential in analyzing various results related to such surface dynamical processes.Concerning to these analysis, it is well known that the Keldysh formalism [9] has been often used, because this method was developed for the purpose of being applicable to non-equilibrium states such as TDAN.Consequently analysis on the basis of the above formalism have been reported and discussed in the field of surface science [10][11][12][13][14][15][16].
In this section, firstly, let us apply our theoretical results to the most simplified RCT neutralization process, where a singly charged ion approaches a metal surface and moves away from it after the ion-surface collision.To describe the electronic interaction between a metal surface and an ion, we consider the following TDAN while ignoring spin orientations: In the above equations, all variables and operator are same as alleady defined in the section of introduction.Since z can be expressed as a function of time (z = v|t| v: By combining Equat (2.5), we obtain ions (3.1)-(3.4)with equation where we define   av ji n t is obtained from the relation .
Consequently, we obtain the simultaneous differential equations corresponding to the Hamiltonian of Equation (3.1).It should be noted that evaluated differential equations are the same expressions as alread previous article [1].The differences between finite temperature case and ground state (0 K) are initial conditions.The merit of this method is that one can obtain numerical s (3.5)-(3.7)under various initial conditions by using computer codes introduction, since Equations (3.5)- (3.rential f rd the solubasic s as integral schemes ascribed to son equation.Consequently our pr y stated in the solutions of Equation .As stated in the 7) are given in the diffe orms, numerical procedures towa tions are ally free from such nuisance Dy oposed method seems to be applicable to the cases even when Hamiltonian includes complicated time dependent terms.
When calculating these simultaneous differential equations, we follow the same procedures as stated in the previous article [1].Concerning the surface-particle (ion) interaction, we adopt the Gaussian form of occupying momentum k and the complete ionic state at t   , respectively.We assume that Fermi level μ 0 remains in the middle of the band at bandwidth D and energy interval E h  are 0.8 eV and 0.01 eV, respectively.The range of energy  neutralization rate, it is found that the RCT is numerically confirmed at the range of 1.4 eV a exp V a 0 −2 ) and ion ve    2 eV.As illustrated in Figure 1, neutralization rate is not so sensitive to temperature, but rather slow even in the range of 600 K, which seems consistent with the report by Sulston et al. [17].
Secondly, let us apply our method to the same RCT process when impurity potential is newly introduced.Actually, such a theoretical calculations as evaluating the temperture dependent properties of RCT under the preence of impurity potential hav not been reported and discussed.Usually the presence of perturvative Hamiltocedures, bec ly reobtain exact s rres s e nian will easily cause the difficulty in the analytical proause the complicated and/or awkward schemes for the evaluation of self energy is usual quired.It seems, therefore, to be a very difficult task to olution in this case, using conventinal methods such as Green functions based on Keldysh fo malism.However, our proposed methods straightforwardly gives exact solution and procedu of analysis are very simple.Let us consider the case where impurity potential exists on crystal surface.Hamiltonian is (3.12)Concerning   a E z , we take the same expression as Equation (3.9).Based on the above Hamiltonian, we obtain the following equations;

 
Fermi Level is As shown in Equations (3.13)- (3.15), one can easily see that our proposed method is enough simple to evaluate numerical solutions using computer codes and is able to gives exact solutions because impurity potenti regarded as one body potential. for V 0 = +10 meV at T = 0 K, 300 K an 600 K together with the data of V 0 = 0 m gure 2, the shape for d eV at T = 0 K.As shown in Fi of neutralization rate drastically changes when impurity potential as small as 10 meV is newly introduced into the system: single peak for V 0 = 0 meV while about three peaks V 0 = +10 meV.From the values of neutralization rate, neutralization is numerically confirmed at the range of about 1.4 eV    2.4 eV for V 0 = +10 meV, showing the extension in comparison with the case of V 0 = 0 meV.This extension of neutralization is also n erically confirmed in the previous article [ uld be remarkable that the all as 10 meV greatly extends t um 1]; thus it sho presence of impurity potential as sm he neutralization range.H mplified RCT and RCT w our proposed owever, neutralization rate is not so sensitive to temperature, which shows same tendency as Figure 1.
In this subsection, for the purpose of verifying the usefulness of our proposed method, we apply this approach to the two cases of RCT: most si ith impurity potential newly induced.method straightforwardly offers exact solutions by such a simple procedure as numerically solving the simultaneous differential equations of (3.5)-(3.7)and (3.13)- (3.15).It should be easily noticed again that the only differences between present work and previous article [1] are initial conditions.Thus we can easily conclude that simultaneous differential equations for seeking temperature dependence of neutralization probability take the same forms as those in previous article [1] except for initial cconditions.

Conclusions
Starting from Neumann equation, we obtain equation of motion in finite temperature, i.e., Equation (2.5).As illustrated below, it should be obvious that our theoretical   at T = 0 K. Through the analysis of both cases, neutralization process is found to be not so sensitive to temperature up to 600 K while impurity potential as small as 10 meV greatly changes the neutralization process.
In addition, as already stated in the previous article [1], small energy intervals ΔE can give accurate and reliable calculation results, which means that the number of electrons N is large because of Δ N D E  (D: band width).Furthermore, assuming that the system consists of N electrons, then numerical solutions of 2 N simultaneous differential equations are required.For example, let us consider actual case.Usually D takes several eV.If we set energy interval Δ 1 E  meV, the value of N is in the range of several thousand, which indicates that we have to solve 10 6 -10 7 simutaneous differential equations.In this article, we adopt algorithm based on the Runge-Kutta-Fehlberg formula and solve numerical 10 3 simultaneous differential equations, ve a hug difequations a from points of mo ula ly about 6 × however a fast algorithm to sol e number of simultaneous ferential ccurately is strongly desirable the view re precise sim tions.
," Progress of Theoretical

4 )
By defining perturbative Hamiltonian in the form of separation of variables as shown in Equation (1.4), various time-dependent physical properties av B t B in finite temperature, which originates from the operator operators of a conduction electron with momentum k and spin orientations  .Con- constant, v: ion velocity) without an imaginary part and momentum k dependence.For the purpose of simplifying the discussion, we assume S.However, as shown in Equations (3.5)-(3.7), it is sufficiently possible to calculate numerically in the case of

Figure 1
Figure 1 illustrates neutralization rate

Figure 2
Figure 2. Temperature dependence of neutralization rate     aa as a function of n    en yz oud state (T = 0 K) and finite temperature are simple enough to change initial conditions: Ĉ  and h band, i.e., −0.4 eV ≤ E k ≤ 0 eV is filled with electrons at T = 0 K.Additionally, when the ion is close to the surface, the a E on the surface-ion distance z can be expressed as below in the atomic unit, It should be noted that our proposed method can be easily applied to the many body problem.For example, when there exists the Coulomb repulsive term between d  and d  electrons as given as below, † By introducing such an approximation, we obtain nonlinear simultaneous differential equations with showing closed form, which means that numerical solutions are obtainable.However, the above approximation does not include spin cross term such as proximations than Hartree Fock are considered to have the spin cross terms, we could expect such approximat s are ex r ssed as functions of For further research,it seems more importatnt to find and determine optimized function of F and G. physical properties even if the sysime-dependent behaviors.The ded as difficult to analyze be-cause of the presence of impurity potential.The presence of perturbative Hamiltonian such as impurity pot tial usually causes comlicated and awkward schemes for the evaluation of self energy.However, our propsoed method straightforwardly deduces the solutions without any assumution, which is considered a great merit in anal ing the case where pertubative Hamiltoniam is introduced.Based on our method, the only differences between gr