4 III. TIME EVOLUTION To calculate the retarded Green’s function G r $_{σ}$( t ) = − iθ ( t ) 〈 GS |{ c$_{σ}$ ( t ) , c † σ }| GS 〉 , one first applies the c † σ or c$_{σ}$ on the ground state from the previous VMPS calcula- tions: | φ (0) 〉 = c † σ | GS 〉 , | χ (0) 〉 = c$_{σ}$ | GS 〉 . Then let states evolute in real time to get | φ ( t ) 〉 = e − i ( H − E$_{G}$ ) $^{t}$c † σ | GS 〉 and | χ ( t ) 〉 = e i ( H − E$_{G}$ ) $^{t}$c$_{σ}$ | GS 〉 . Thus the retarded Green’s function can be calculated as The operation of the local operators on the ground state MPS is Where a string of JW signs and a d × d matrix s + act on the k -th site. After the multiplication, the new state is still represented as the the MPS of the same rank. The real-time evolution is performed by first split- ting the evolution operator into small pieces e − iHt = ( e − iH Δ $^{t}$)$^{M}$, where Δ t = t $_{M}$, then applying the Trotter decomposition. To the second order one has e − iH Δ t = e − iH$_{o}$ Δ t/ $^{2}$e − iH$_{e}$ Δ t/ $^{2}$e − iH$_{Z}$ Δ $^{t}$e − iH$_{e}$ Δ t/ $^{2}$e − iH$_{o}$ Δ t/ 2 + O ((Δ t ) $^{3}$) . Where H = H$_{o}$ + H$_{e}$ + H$_{Z}$ . H$_{e}$ and H$_{o}$ are two particle terms act on even and odd bonds and H$_{Z}$ includes the single particle terms. Since each term within H$_{α}$ , ( α = e, o, Z ) commutes one could apply the evolution operator U$_{α}$ (Δ t ) = e − iH$_{α}$ Δ t associate with each bond one after another without leading to further errors. In the particle-hole symmetric case, the onsite energies of the chain are zero and the chemical potential µ = U/ 2 can be put into the interaction term. Thus the single particle evolution operator U$_{Z}$ = 1 The operation of the single site evolution operator U$_{Z}$ is similar to the applying of the creation/annihilation op- erators, Eqn 4. The operation of the two site evolution operators U$_{e}$ or U$_{o}$ is shown schematically in Fig.4. The physical indices for the two sites are first merged and ex- changed, then the d 2 × d 2 operator is separated into two matrices U$_{1}$ and U$_{2}$ by the singular value decomposition (SVD). Each of them only carries the physical indices of one site. In both cases, the evolution operators are writ- ten as the matrix product operator (MPO) $^{35}$, which has four indices, two for the visual dimensions and two other for the physical indices, see Fig.2. Generally after apply- ing them onto the target states, the resulting states could not be written as the MPS with same visual dimension as before. Thus one needs to project the state into the orig- inal MPS space. It is done variationally by minimizing the norm: N = ‖| φ $^{′}$〉− U$_{α}$ (Δ t ) | φ ( t ) 〉‖ , and approximately one has | φ ( t + Δ t ) 〉 = | φ $^{′}$〉 . The minimization procedure can be carried out following the similar procedures as in Sec.II. It is performed every step after the applying of the evolution operator on the original state. In the course of evolution, we increase the bond dimension to Imaginary time Green’s function could also be calcu- lated similarly, in which we perform the evolution algo- rithm along the imaginary time axis. Imaginary time evo- lution technique could also be used to search the ground state. The imaginary time Green’s function G ( τ ) is a de- cayed function from n$_{σ}$ − 1 , and it is decaying faster for larger interacting strength U . For the particle-hole sym- metric case, the decaying form approaches single mode decay − 1 $_{2}$e − µτ in the atomic limit. See Figure 5 for the real time Green’s function for dif- ferent interacting strength. The retarded Green’s func- tion G $^{r}$( t ) evolutes from the Bessel function at the non- interacting limit to the cosine function at the atomic limit. In between one has the oscillated decaying curves, and the period of the oscillation gives the frequency of the Hubbard band. The normalization of the density of states (DOS) is automatically fulfilled since it implies that G $^{r}$( t = 0) = − 1 . There may be a concern on the orthogonality catastrophe (OC) 36,37 which states that in the thermodynamic limit local perturbation leads to com- plete reconstruction of the ground state of a fermionic system in such a way that the overlap of the "old" and "new" ground-state wave functions is proportional to N − $^{α}$. But since we are dealing with finite N here, the OC is irrelevant. $^{36,37}$. By comparing the real time Green’s function with different bath sites, we notice that even with small number of sites, one could reproduce the re- sulting spectral information for the thermal-dynamical limit, as long as the wavefront created by the adding (re- moving) an electron on the impurity site does not reach the boundary. maintain the smallest eigenvalue of the reduced density matrix larger than a threshold w$_{min}$ = 10 − $^{10}$. G r $_{σ}$( t ) = − iθ ( t )[ 〈 φ (0) | φ ( t ) 〉 + 〈 χ (0) | χ ( t ) 〉 ] (3) $_{c}$† k $_{↑}$| GS 〉 = ∏ l<k ( $^{p$_{l}$}$↑ )$_{σ}$$_{l}$$_{σ}$ ′ $_{l}$A [ σ ′ $^{$_{l}$]}$( s + $_{k}$)$_{σ}$$_{k}$$_{σ}$ ′ $_{k}$A [ σ ′ $_{k}$] ∏ m>k A [ σ$_{m}$ $^{]}$| ⃗ σ$_{l}$, σ$_{k}$, ⃗ σ$_{m}$ (4) 〉 Figure 4: (Color online) Writing the two particle evolution operator into a MPO. (a). The physical indices of the four indices tensor are merged, giveing a d 2 × d 2 matrix. (b). The order of the indices are exchanged. (c). SVD of the matrix separates it into two parts (d). Rearrangement of the physical indices gives the MPO representation of the two-sites evolution operator. U α ′ β ′ αβ U2 U1 β ′ β α ′ α U β ′ β α ′ α α U1 U2 α ′ β ′ β U β α α ′ β ′ SVD ( a ) ( b ) ( c ) ( d )