Nonequilibrium dynamics in the one-dimensional Fermi-Hubbard model: A comparison of the nonequilibrium Green functions approach and the density matrix renormalization group method
N. Schluenzen, J.-P. Joost, F. Heidrich-Meisner, and M. Bonitz
TL;DR
This paper compares the nonequilibrium Green functions approach with the density matrix renormalization group method for simulating the dynamics of strongly-correlated fermions in the one-dimensional Fermi-Hubbard model, highlighting their respective strengths.
Contribution
It provides an extensive comparison between NEGF and DMRG methods, demonstrating NEGF's reliability for weak to intermediate coupling and DMRG's efficiency at strong coupling.
Findings
NEGF is reliable for weak to intermediate coupling.
DMRG is efficient at strong coupling.
Both methods complement each other in simulating nonequilibrium dynamics.
Abstract
The nonequilibrium dynamics of strongly-correlated fermions in lattice systems have attracted considerable interest in the condensed matter and ultracold atomic-gas communities. While experiments have made remarkable progress in recent years, there remains a need for the further development of theoretical tools that can account for both the nonequilibrium conditions and strong correlations. For instance, time-dependent theoretical quantum approaches based on the density matrix renormalization group (DMRG) methods have been primarily applied to one-dimensional setups. Recently, two-dimensional quantum simulations of the expansion of fermions based on nonequilibrium Green functions (NEGF) have been presented [Schluenzen et al., Phys. Rev. B 93, 035107 (2016)] that showed excellent agreement with the experiments. Here we present an extensive comparison of the NEGF approach to numerically…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11
Figure 12
Figure 13
Figure 14
Figure 15Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Nonequilibrium dynamics in the one-dimensional Fermi-Hubbard model: A comparison of the nonequilibrium Green functions approach
and the density matrix renormalization group method
N. Schlünzen
Institut für Theoretische Physik und Astrophysik, Christian-Albrechts Universität Kiel, Leibnizstr. 15, Germany
J.-P. Joost
Institut für Theoretische Physik und Astrophysik, Christian-Albrechts Universität Kiel, Leibnizstr. 15, Germany
F. Heidrich-Meisner
Department of Physics and Arnold Sommerfeld Center for Theoretical Physics,Ludwig-Maximilians-Universität München, D-80333 München, Germany
M. Bonitz
Institut für Theoretische Physik und Astrophysik, Christian-Albrechts Universität Kiel, Leibnizstr. 15, Germany
Abstract
The nonequilibrium dynamics of strongly-correlated fermions in lattice systems have attracted considerable interest in the condensed matter and ultracold atomic-gas communities. While experiments have made remarkable progress in recent years, there remains a need for the further development of theoretical tools that can account for both the nonequilibrium conditions and strong correlations. For instance, time-dependent theoretical quantum approaches based on the density matrix renormalization group (DMRG) methods have been primarily applied to one-dimensional setups. Recently, two-dimensional quantum simulations of the expansion of fermions based on nonequilibrium Green functions (NEGF) have been presented [Schlünzen et al., Phys. Rev. B 93, 035107 (2016)] that showed excellent agreement with the experiments. Here we present an extensive comparison of the NEGF approach to numerically accurate DMRG results. The results indicate that NEGF are a reliable theoretical tool for weak to intermediate coupling strengths in arbitrary dimensions and make long simulations possible. This is complementary to DMRG simulations which are particularly efficient at strong coupling.
pacs:
67.85.-d, 05.30.Jp, 37.10.Jk
I Introduction
Experiments addressing the nonequilibrium dynamics of quantum many-body systems have made remarkable progress in recent years, both probing ultrafast dynamics in strongly correlated materials Giannetti et al. (2016); Gandolfi et al. (unpublished) and quantum quenches in interacting quantum gases (see Refs. Langen et al., 2015; Eisert et al., 2015; Gogolin and Eisert, 2016 for a review). Among the many ultracold quantum-gas experiments with fermions we mention the study of the expansion dynamics of strongly-correlated fermions in a two-dimensional optical lattice Schneider et al. (2012), the collapse and revival dynamics of Fermi-Bose mixtures Will et al. (2015) and the real-time decay of a density wave in one-dimensional lattices Pertot et al. (2014). Very recently, several experimental groups reported the successful implementation of fermionic quantum-gas microscopes Cheuk et al. (2015); Parsons et al. (2015); Haller et al. (2015); Edge et al. (2015); Omran et al. (2015); Cheuk et al. (2016a, b); Greif et al. (2016); Boll et al. (2016), which will give unprecedented access to both equilibrium and nonequilibrium properties of the Fermi-Hubbard model. Given the tremendous success of the earlier bosonic quantum-gas microscopes in exploring the nonequilibrium realm Cheneau et al. (2012); Fukuhara et al. (2013a, b); Hild et al. (2014); Preiss et al. (2015), a considerable experimental activity in studying quantum-quench dynamics in the Fermi-Hubbard model can be expected in the near future. Quantum-gas microscopes operate with two-dimensional systems which will push the efforts into this most challenging regime (see also Ref. Cocchi et al., 2016) while also allowing to study one-dimensional systems Boll et al. (2016).
A large body of theoretical work has concentrated on one-dimensional systems, the reason being both experiments Kinoshita et al. (2006); Hofferberth et al. (2007); Langen et al. (2013); Gring et al. (2012); Trotzky et al. (2012); Cheneau et al. (2012); Ronzheimer et al. (2013) as well as the availability of powerful theoretical tools based on field theory Giamarchi (2004), integrability Essler et al. (2005) or numerical methods. While exact diagonalization (ED) is still an indispensable tool (see, e.g., Refs. Manmana et al., 2005; Rigol et al., 2008; Prelovsek and Bonca, 2013), it is limited to small systems. Nonetheless, for problems restricted to the dynamics of a single charge carrier coupled to spin or phonon degrees of freedom, there exist Krylov-space approaches that operate in a subspace of the full Hilbert space constructed by selecting only those states accessible by the Hamiltonian dynamics Bonča et al. (1999). Such an exact diagonalization in a limited functional space has been applied quite extensively to two-dimensional problems as well (see, e.g., Refs. Vidmar et al., 2011; Mierzejewski et al., 2011; Bonča et al., 2012; Golež et al., 2014).
Time-dependent density matrix renormalization group (DMRG) methods Vidal (2004); Daley et al. (2004); White and Feiguin (2004) have been very widely applied to nonequilibrium problems and yield numerically accurate results, but are limited by the accessible times scales and are primarily useful for one-dimensional systems. A recent variant of the method Zaletel et al. (2015) has been tailored for long-range interactions and is thus better suited for coupled one-dimensional and two-dimensional systems Zaletel et al. (2015); Hauschild et al. (2015) but cannot overcome the exponential scaling of a matrix-product states ansatz with the number of coupled chains. The application of time-dependent tensor-network approaches that are based on ansatz states suitable for two-dimensional systems such as the projected entangled pair states has been very little explored Dorando et al. (2009); Haegeman et al. (2011); Lubasch et al. (2011).
Apart from time-dependent DMRG methods, there are other many-body methods for the real-time evolution including continuous-time quantum Monte Carlo Gull et al. (2011) and time-dependent dynamical mean-field theory approaches Schmidt and Monien (unpublished); Freericks et al. (2006); Eckstein et al. (2010). The former, while being able to achieve essentially exact results for short evolution times, can suffer from a dynamical sign problem Gull et al. (2011). The latter method often utilizes continuous-time quantum Monte Carlo as an impurity solver, while in more recent developments, time-dependent DMRG has also been successfully used for this purpose Wolf et al. (2014); Ganahl et al. (2015). Time-dependent DMFT methods are not exact in two dimensions either but are argued to capture better the physics of strongly-correlated systems in higher dimensions, leading to a wide range of applications in the context of nonequilibrium dynamics in the Hubbard model (see, e.g., Ref. Eckstein et al., 2010). Finally, the iterative equation-of-motion method for operators provides an alternative approach Uhrig (2009), which has also been applied to quantum quench problems in the 2D Fermi-Hubbard model Hamerla and Uhrig (2014).
Despite all these efforts, there still is a significant gap between the rapidly progressing experiments in the field of ultracold atoms and accurate quantum dynamics simulations when it comes to correlated systems in two dimensions. To contribute towards closing this gap, two of us have applied an alternative approach to the quantum simulation of the nonequilibrium mass-transport of correlated fermions studied in the experiment of Ref. Schneider et al., 2012: nonequilibrium Green functions (NEGF). Previously, this theory has been successfully applied to a variety of many-particle systems including the correlated electron gas Kwong and Bonitz (2000), electron-hole plasmas Kwong et al. (1998), nuclear matter Rios et al. (2011) and electrons in quantum dots Balzer et al. (2009); Lorke et al. (2006), for a recent overview see Ref. png, 2016. Extensive applications to finite Hubbard clusters were presented in Ref. Hermanns et al., 2014 and first applications of NEGF to mass transport in small lattice systems of correlated fermions were shown in Ref. Bonitz et al., 2015. Finally, in Ref. Schlünzen et al., 2016 these simulations were extended to strong coupling by using -matrix selfenergies as well as to substantially larger systems. Applying an extrapolation to the thermodynamic limit the nonequilibrium correlated quantum mass transport in two-dimensional fermion ensembles could be directly compared to the experiments of Ref. Schneider et al., 2012, and excellent agreement was observed. For an overview on the NEGF approach and its application to inhomogeneous Hubbard clusters, see Ref. Schlünzen and Bonitz, 2016.
Even though NEGF simulations are computationally demanding, they have a number of remarkable advantages. First, they do not exhibit an exponential scaling with system size, as is the case for exact diagonalization, and they do not have a dynamical sign problem as continuous time QMC methods. Second, they are not limited with respect to the system dimensionality, as opposed to matrix-product state methods. At the same time, in contrast to ED, NEGF simulations are not a first-principle method since they involve a many-body approximation—the selfenergy—which determines the accuracy and the quality of the results, similar to the approximate exchange-correlation energy in density functional theory. DMRG, on the other hand, also involves an approximation but the numerical errors depend on a control parameter, the discarded weight, and whenever this can be made sufficiently small, the results can become essentially exact for system sizes larger than what is accessible to ED Schollwöck (2005, 2011).
The accuracy of NEGF simulations of spatially inhomogeneous fermion systems was tested before for few-electron atoms Balzer et al. (2010) and small Hubbard clusters von Friesen et al. (2010) where exact diagonalization results are available. This analysis was extended to larger Hubbard systems, on the order of 10 sites in Refs. Hermanns et al., 2014; Lacroix et al., 2014, revealing a high accuracy of simulations with second-order Born selfenergies, for weak coupling and moderate times (on the order of 20 inverse hopping times). However, the quality of the results for larger systems has remained open until now, due to the lack of reliable benchmark data. On the other hand, for small Hubbard clusters, also problems were reported: in the case of a strong excitation, two-time NEGF simulations were found to exhibit an unphysical damping of the dynamics von Friesen et al. (2009, 2010). The origin of this behavior has been traced back to the selfconsistent nature of the used approximations. These deficiencies could be removed to a large extent by making the transition to single-time dynamics with the help of the generalized Kadanoff–Baym ansatz Lipavský et al. (1986) with Hartree-Fock propagators (HF-GKBA) Hermanns et al. (2014).
Thus, there is a clear need to further study the question of accuracy and predictive capability of NEGF simulations, in particular, for systems larger than those studied so far, for longer simulation times and beyond the weak-coupling limit. The goal of this paper is to present such an analysis by benchmarking NEGF results using a variety of different selfenergy approximations, in a two-time as well as in a single-time formalism (i.e., using the GKBA), against DMRG results. Due to the inherent properties of matrix-product states Schollwöck (2011), these comparisons have to focus on 1D fermion systems. We choose a set of four non-trivial cases of nonequilibrium dynamics in the Fermi-Hubbard model for which correlations play a crucial role.
The Hamiltonian of the Fermi-Hubbard chain is
[TABLE]
where {\mathchoice{\makebox[5.55557pt][l]{\displaystyle\hat{c\mspace{3.5mu}}}}{\makebox[5.55557pt][l]{\textstyle\hat{c\mspace{3.5mu}}}}{\makebox[5.55557pt][l]{\scriptstyle\hat{c\mspace{3.5mu}}}}{\makebox[5.55557pt][l]{\scriptscriptstyle\hat{c\mspace{3.0mu}}}}}_{\bm{s},\sigma\vphantom{{\mathchoice{\makebox[5.00002pt][l]{\displaystyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\textstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptscriptstyle\bar{\bar{x\mspace{3.0mu}}}}}}}}^{\dagger\vphantom{\dagger}} creates a fermion with spin in site and n_{\bm{s}\vphantom{{\mathchoice{\makebox[5.00002pt][l]{\displaystyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\textstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptscriptstyle\bar{\bar{x\mspace{3.0mu}}}}}}}}^{\sigma\vphantom{\dagger}}={\mathchoice{\makebox[5.55557pt][l]{\displaystyle\hat{c\mspace{3.5mu}}}}{\makebox[5.55557pt][l]{\textstyle\hat{c\mspace{3.5mu}}}}{\makebox[5.55557pt][l]{\scriptstyle\hat{c\mspace{3.5mu}}}}{\makebox[5.55557pt][l]{\scriptscriptstyle\hat{c\mspace{3.0mu}}}}}_{\bm{s},\sigma\vphantom{{\mathchoice{\makebox[5.00002pt][l]{\displaystyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\textstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptscriptstyle\bar{\bar{x\mspace{3.0mu}}}}}}}}^{\dagger\vphantom{\dagger}}{\mathchoice{\makebox[5.55557pt][l]{\displaystyle\hat{c\mspace{3.5mu}}}}{\makebox[5.55557pt][l]{\textstyle\hat{c\mspace{3.5mu}}}}{\makebox[5.55557pt][l]{\scriptstyle\hat{c\mspace{3.5mu}}}}{\makebox[5.55557pt][l]{\scriptscriptstyle\hat{c\mspace{3.0mu}}}}}_{\bm{s},\sigma\vphantom{{\mathchoice{\makebox[5.00002pt][l]{\displaystyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\textstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptscriptstyle\bar{\bar{x\mspace{3.0mu}}}}}}}}^{\vphantom{\dagger}}. is the hopping matrix element (set to unity in our simulations), denotes the onsite interaction and is the number of sites (the lattice spacing and are set to unity). The cases studied include (i) the symmetric and (ii) asymmetric expansion from a band insulator into an empty lattice, (iii) the expansion from a band insulator in a quasi-1D situation on a two-leg ladder, and (iv) the decay of an ideal charge-density wave state. These four initial situations are sketched in Fig. 1.
As a result of this analysis, the applicability range of NEGF simulations and relevant approximations is being mapped out. Our main results are the following: NEGF simulations with the HF-GKBA are reliable also for moderate coupling, , if the proper selfenergies are being used. These are the -matrix selfenergy–for small or large filling–and the third-order selfenergy (including all diagrams of third order, cf. Sec. II.1)–close to half filling. In all cases, two-time simulations are less accurate (due to the unphysical damping mentioned above) but they can be used to estimate the deviations of the HF-GKBA from the exact result as typically the latter is enclosed between single-time and two-time NEGF results. Finally, NEGF simulations fill the gap left open by DMRG by being capable to treat large systems (of any dimensionality) and to achieve long simulation times, for weak and moderate coupling, whereas the DMRG is advantageous and more efficient for strong coupling.
The remainder of this paper is as follows: in Sec. II we give a brief introduction into NEGF and time-dependent DMRG simulations. This is followed in Sec. III by numerical results. There we study the four cases introduced above and depicted in Fig. 1: a symmetric and asymmetric sudden expansion (confinement quench) in 1D, Secs. III.1 and III.2, respectively, the sudden expansion in a two-leg ladder, Sec. III.3, and a charge-density-wave initial state, Sec. III.4.
II Methods
II.1 Nonequilibrium Green functions (NEGF)
The central quantity in the nonequilibrium Green functions theory is the (single-particle) Green function . It is defined on the Schwinger–Keldysh contour Schwinger (1961); Keldysh (1965) via the time-ordering operator ,
[TABLE]
where, denotes the ensemble average. The Green function can be understood as a generalization of the nonequilibrium single-particle density matrix, n_{\bm{s},{\bm{s}}^{\prime}\vphantom{{\mathchoice{\makebox[5.00002pt][l]{\displaystyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\textstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptscriptstyle\bar{\bar{x\mspace{3.0mu}}}}}}}}^{\sigma\vphantom{\dagger}}(t), onto the two-time plane. Therefore, provides easy access not only to the observables related to n_{\bm{s},{\bm{s}}^{\prime}\vphantom{{\mathchoice{\makebox[5.00002pt][l]{\displaystyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\textstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptscriptstyle\bar{\bar{x\mspace{3.0mu}}}}}}}}^{\sigma\vphantom{\dagger}} but, in addition, also to the spectral properties of the system. However, the full -particle information is not directly available from , although, for example, the pair correlation function can be reconstructed from Bonitz et al. (2013).
The equations of motion for the single-particle Green function are the Keldysh–Kadanoff–Baym equations Kadanoff and Baym (1962),
[TABLE]
together with the adjoint equation ( denotes the matrix element of the single-particle Hamiltonian). denotes the selfenergy which is the only unknown of the theory, and with an exact the method would be exact. In practice, the selfenergy has to be approximated for which systematic many-body schemes (e.g. Feynman diagrams) exist that are applicable in equilibrium as well as in nonequilibrium, via the use of the time contour.
In the following we list the selfenergies that are used in the present paper. The contribution of the first order in the interaction is given by the Hartree–Fock (mean field) selfenergy,
[TABLE]
which is contained in each of the approximations used below. Many-body approximations that go beyond the mean field level (that are of higher than first order in ) contain, in addition, a correlation selfenergy, i.e., {\Sigma_{\bm{s}{\bm{s}}^{\prime}\vphantom{{\mathchoice{\makebox[5.00002pt][l]{\displaystyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\textstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptscriptstyle\bar{\bar{x\mspace{3.0mu}}}}}}}}^{\sigma\vphantom{\dagger}}=:\Sigma_{\bm{s}{\bm{s}}^{\prime}\vphantom{{\mathchoice{\makebox[5.00002pt][l]{\displaystyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\textstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptscriptstyle\bar{\bar{x\mspace{3.0mu}}}}}}}}^{\textnormal{HF},\sigma\vphantom{\dagger}}+\Sigma_{\bm{s}{\bm{s}}^{\prime}\vphantom{{\mathchoice{\makebox[5.00002pt][l]{\displaystyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\textstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptscriptstyle\bar{\bar{x\mspace{3.0mu}}}}}}}}^{\textnormal{cor},\sigma\vphantom{\dagger}}}. The first correlation correction is of second order and works well for weak coupling, i.e. , for a discussion see Ref. Hermanns et al. (2014). Here we want to go beyond the weak coupling regime. Therefore, we focus on two higher order many-body approximations. The first is the -matrix approximation (TMA) in the particle-particle channel and yields a selfenergy \Sigma_{\bm{s}{\bm{s}}^{\prime}\vphantom{{\mathchoice{\makebox[5.00002pt][l]{\displaystyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\textstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptscriptstyle\bar{\bar{x\mspace{3.0mu}}}}}}}}^{\textnormal{cor},\sigma\vphantom{\dagger}} which accounts for scattering processes up to infinite order (see Ref. Schlünzen and Bonitz, 2016 for a detailed discussion). This is realized by the selfconsistent, recursive structure of the -matrix which can be understood as an effective interaction that obeys its own equation of motion (the Lippmann-Schwinger equation), Eq. (6),
[TABLE]
The corresponding Feynman diagrams are shown in Fig. 2(a). The TMA is known to perform best in the limit of small (or large) density Schlünzen and Bonitz (2016); von Friesen et al. (2009); Hermanns (2016), i.e., when the interaction in the system is dominated by electron-electron or hole-hole scattering events. Around half-filling, however, electron-hole scattering gains in importance which is not captured by the particle-particle TMA. Therefore, we introduce, in addition, the third-order approximation (TOA) which contains all selfenergy contributions up to . In this approximation the correlation selfenergy, \Sigma_{\bm{s}{\bm{s}}^{\prime}\vphantom{{\mathchoice{\makebox[5.00002pt][l]{\displaystyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\textstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptscriptstyle\bar{\bar{x\mspace{3.0mu}}}}}}}}^{\textnormal{cor},\sigma\vphantom{\dagger}}, attains the following form Hermanns (2016),
[TABLE]
The corresponding diagrams are shown in Fig. 2(b). In the TOA, particle-particle and electron-hole scattering processes are considered on equal footing, yet only to third order inclusively. Both the TMA and TOA approach have been found to perform well for weak to moderate coupling strengths as long as the respective density conditions are fulfilled Schlünzen and Bonitz (2016); von Friesen et al. (2009); Hermanns (2016); Schlünzen et al. (2016); Hermanns et al. (2014); Bonitz et al. (2015); von Friesen et al. (2009, 2010); Balzer et al. (2016).
Finally, we introduce the generalized Kadanoff-Baym ansatz (GKBA) which is an approximation that reduces the complexity of the time structure of NEGF theory by separating the time-diagonal Green functions from the off-diagonal ones. The full KBE of Eq. (3) is solved only for , while, for , the Green function is reconstructed from its time-diagonal values, i.e., from the single-particle density matrix. For the latter, the less and greater component of which originate from the mapping of the time contour onto the real time axis Balzer and Bonitz (2013) are reconstructed according to Lipavský et al. (1986),
[TABLE]
where {\mathfrak{n}^{<,\sigma}_{\bm{s}{\bm{s}}^{\prime}}(t)=n_{\bm{s},{\bm{s}}^{\prime}\vphantom{{\mathchoice{\makebox[5.00002pt][l]{\displaystyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\textstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptscriptstyle\bar{\bar{x\mspace{3.0mu}}}}}}}}^{\sigma\vphantom{\dagger}}(t)} and {\mathfrak{n}^{>,\sigma}_{\bm{s}{\bm{s}}^{\prime}}(t)=n_{\bm{s},{\bm{s}}^{\prime}\vphantom{{\mathchoice{\makebox[5.00002pt][l]{\displaystyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\textstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptscriptstyle\bar{\bar{x\mspace{3.0mu}}}}}}}}^{\sigma\vphantom{\dagger}}(t)-\delta^{\vphantom{\dagger}}_{{\bm{s}},{{\bm{s}}^{\prime}}\vphantom{{\mathchoice{\makebox[5.00002pt][l]{\displaystyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\textstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptscriptstyle\bar{\bar{x\mspace{3.0mu}}}}}}}}}. The GKBA does not violate the attractive properties of the NEGF method, as it retains density and energy conservation, as well as time reversibility Scharnke et al. (2017). When using the GKBA, still the question remains how the retarded and advanced propagators Latini et al. (2014) are approximated. Here, we concentrate on Hartree–Fock propagators—the resulting approximation will be called HF-GKBA Hermanns et al. (2014). This approximation has been shown to eliminate (or drastically reduce) the artificial damping properties of two-time simulations for strongly excited systems and, at the same time substantially improving the computational performance.
II.2 Time-dependent density matrix renormalization group method (DMRG)
The density matrix renormalization group method White (1992); Schollwöck (2005, 2011) relies on approximating many-body wave-functions via matrix-product states of a finite bond dimension . A matrix-product state can be written as
[TABLE]
where are the local degrees of freedom at site and are matrices of dimensions (for details and the role of boundary conditions, see Ref. Schollwöck, 2011). Any wave-function can be brought into the form Eq. (9) by a sequence of singular value decompositions where in general, the bond dimension of the matrices will scale exponentially in system size. To illustrate this procedure, consider a one-dimensional system that is cut into two parts and . By denoting complete basis sets in the parts and by and , we can express a many-body wave-function as
[TABLE]
By means of a singular value decomposition of the rectangular matrix , this can be reexpressed in terms of new basis sets in and with a single index
[TABLE]
where the are the singular values and is the Schmidt number, which in general scales exponentially with system size. At the heart of the approximation used in DMRG and matrix-product states methods in general is a truncation in the number of states used to represent by keeping only those states with the largest Schmidt coefficients , i.e.,
[TABLE]
This is equivalent to diagonalizing the reduced density matrix of part and truncating in its eigenspectrum, which was White’s original formulation White (1992)
[TABLE]
While actual algorithms are described comprehensively in Ref. Schollwöck, 2011, we here want to explain for which many-body states Eq. (12) provides a useful approximation in the sense that few states (order of ) suffice to obtain numerically accurate results for observables . This obviously depends on how quickly the eigenvalues of the reduced density matrix decay. A correct intuition can be gained from relating the decay of to the entanglement entropy
[TABLE]
A fast decay of translates into a weakly entangled wave-function and vice versa. The crucial question is the scaling of the entanglement entropy with the system size. For ground-states of gapped Hamiltonians with short-range interactions, an area law holds Eisert et al. (2010)
[TABLE]
where is the linear dimension of subsystem and is the dimension. Scaling in translates directly into scaling in , i.e., the liner dimension of the full system. In , we obtain for where is the correlation length and this implies that the numerical effort (i.e., the number of states used to approximate ) does not increase with system size since Schollwöck (2011). For critical systems in one dimension, the entanglement entropy acquires a logarithmic correction (see Refs. Schollwöck, 2011; Eisert et al., 2010 and references therein). This reasoning explains why matrix-product states based techniques work primarily for one-dimensional systems since in 2D, even if an area law holds, the scaling is exponential in the width of the system Stoudenmire and White (2012). The other important issue is whether an efficient algorithm can be formulated based on matrix-product states. It turns out that most matrix-product state methods including DMRG scale as and linearly in Schollwöck (2011).
For real-time evolutions , the application of the time-evolution operator can be efficiently implemented via a Trotter-Suzuki decomposition (using ) into operators which is just a local two-site gate affecting two -matrices in Eq. (9) ( is the time step) Vidal (2004); Daley et al. (2004); White and Feiguin (2004). In general, a time-propagated many-body state will develop volume-law like entanglement even if the initial state was a product state Schollwöck (2011). For global quantum quenches (such as the time-evolution from a product state such as our case (iv), see Fig. 1(d)), the entanglement grows linearly in time , implying that the number of states needed to maintain the same quality of approximation to the true will increase exponentially. This limits the accessible times in global quenches to about , while in local quenches, geometric quenches such as the sudden expansion considered in examples (i)-(iii) or for slow parameter changes, a milder entanglement increase occurs.
The two main parameters that control the accuracy of time-dependent DMRG simulations are the time step and the discarded weight Gobert et al. (2005); Schollwöck (2011). The latter is defined as
[TABLE]
which is a measure for the error made per truncation. The quality of DMRG data has to be analyzed as a function of both and , with the latter the dominant parameter since the dependence of the error on can be reduced by using higher-order Trotter-Suzuki decompositions Schollwöck (2011). In this work, we use a time-dependent DMRG implementation as introduced in Refs. White and Feiguin, 2004; Daley et al., 2004 and we varied the time step between and the discarded weight with a maximum number of states.
III Numerical Results
In our simulations we consider four different nonequilibrium setups. The corresponding initial states are depicted in Fig. 1. In all cases, the incipient configuration consists of both, doubly occupied and empty Hubbard sites. In the first setup, the occupied sites are arranged on a straight line to form a one-dimensional band insulator (BI). During time propagation, and in the absence of any further potential, the density starts to expand symmetrically towards the left and right edges of the Hubbard chain, cf. Fig. 1(a). Next, in order to also investigate an asymmetric expansion dynamics the initial BI is placed onto the leftmost sites of the chain allowing the density to escape only to the right, see Fig. 1(b). Further, to analyze the effect of the dimensionality of the system we extend the asymmetric setup to a two-leg Hubbard ladder the leftmost rungs of which are initially doubly occupied, Fig. 1(c). The dynamics on such ladders is often used to investigate the 1D-to-2D crossover. Finally, we consider a setup that generates a final state at a constant and large density where correlation effects are expected to manifest themselves even stronger, cf. Fig. 1(d). Here, the initial state consists of a one-dimensional Hubbard chain with alternating occupation . During the evolution from this charge-density wave (CDW), the particles quickly form an entangled many-body state in which correlations play a crucial role. These four setups will be analyzed in detail in Secs. III.1–III.4.
III.1 Sudden expansion in 1D: symmetric case
We start the numerical analysis by considering a confinement quench giving rise to a sudden symmetric expansion of an ensemble of fermions into an empty lattice (see Fig. 1(a)). This setup has been studied in many papers, including experimental studies Schneider et al. (2012); Ronzheimer et al. (2013); Xia et al. (2014); Vidmar et al. (2015), and theoretically using DMRG methods Heidrich-Meisner et al. (2008, 2009); Kajala et al. (2011); Langer et al. (2012); Bolech et al. (2012); Vidmar et al. (2013); Mei et al. (2016); Hauschild et al. (2015) and NEGF Schlünzen et al. (2016). We exclusively study an initial density of on the sites that are occupied at , which was previously considered in Refs. Kajala et al., 2011; Keßler et al., 2013; Schlünzen et al., 2016 such that the general properties are well understood.
Here, we focus on the quantitative details of the time dynamics and compare our NEGF results to DMRG. We consider a chain of length with particles for . The evolution of the respective density profiles is shown in Fig. 3(a) for six consecutive time steps .
The solid, black lines correspond to the DMRG results and the dashed green lines belong to the NEGF calculations using the -matrix approximation (TMA) while the orange lines are obtained by additionally invoking the HF-GKBA, cf. Sec. II.1. As expected, the general trend of the density is to propagate outwards resulting in a bell-shaped profile which can be seen from all considered descriptions. For times exceeding , the site occupations start to deviate slightly in the three simulations. In the full two-time NEGF calculation the fermion expansion is slightly faster than in the DMRG, while in the HF-GKBA simulation the particles stay closer to the center and are in very good agreement with the DMRG.
A quantity more sensitive to correlations is the double occupancy,
[TABLE]
the dynamics of which are displayed in Fig. 3(b). It is evident that it follows the trend of the density by which it is dominated. Again, in the full two-time NEGF calculation d_{\bm{s}\vphantom{{\mathchoice{\makebox[5.00002pt][l]{\displaystyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\textstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptstyle\bar{\bar{x\mspace{3.5mu}}}}}{\makebox[5.00002pt][l]{\scriptscriptstyle\bar{\bar{x\mspace{3.0mu}}}}}}}}^{\vphantom{\dagger}} expands faster than in the DMRG result, where the deviations are larger than for the density. In contrast, the HF-GKBA is again very close to the latter.
To better quantify the discrepancies between DMRG and the two NEGF approaches we introduce the total density deviation between the two methods in the following way,
[TABLE]
where I and II denote the respective method. This quantity allows to analyze the time dependent difference of the density profiles. It should be noted that the quantitative value of has no direct interpretation. Instead, by dividing by the total number of Hubbard sites , one gets the average deviation per site. Adopting DMRG as the reference method, we investigate the dependence of the deviation on the interaction strength and time by calculating the total deviations for the two-time TMA simulation and the HF-GKBA results which are displayed in Fig. 4(a) and (b), respectively.
As expected, the total deviation grows in time for all cases. Interestingly, however, the deviations saturate around . A closer look reveals that, during the early propagation period, the growth appears to be superlinear, followed by a receding phase after which the growth becomes more fluctuating. The lengths of these time periods strongly depends on the interaction strength, as they become shorter for larger . As a consequence, for times around the total deviation increases with the interaction strength while for times around it decreases with . The overall trend is common between the full TMA results and the HF-GKBA simulations. The only noticeable difference is that remains a little smaller for larger times and small in the HF-GKBA calculations.
To better understand how the density deviations vary with and , we replot these quantities in Fig. 5(a) for the two time points, and . In addition, we compute the total deviation of the double occupancy which is defined, by analogy to Eq. (18), as
[TABLE]
which is shown in Fig. 5(b).
While for all results confirm the trend that the deviations grow with increasing , for the dependence is more irregular. In the latter case, the deviations between HF-GKBA and DMRG, in particular, are nearly independent of the interaction strength. The two-time NEGF results for the double occupancy, however, show large deviations for small . The decrease of and with increasing for later times can be understood from the direct dynamics of the density profiles. Since for large the particles predominantly remain in the center of the system the growth of the deviations is limited due to the absence of moving particles.
From the presented results, it turns out that the DMRG result is typically enclosed between the HF-GKBA and the two-time NEGF result in -matrix approximation. At the same time the HF-GKBA data are slightly closer to the DMRG results.
To further analyze the expansion behavior following a 1D sudden confinement switch, it is instructive to analyze the expansion velocity of the fermion cloud which is defined according to Schlünzen et al. (2016)
[TABLE]
This quantity measures the temporal growth of the particle cloud which has a mean square radius from which the initial size is subtracted. This quantity was analyzed in detail for 1D, 2D and 3D systems and a broad range of system size in Ref. Schlünzen et al., 2016. Here we focus on the time evolution of for 1D systems and compare again DMRG, two-time NEGF simulations and HF-GKBA. The results are shown in Fig. 6 for .
As one can see, for all cases starts from the same value which is the largest expansion velocity in an empty lattice (cf. Refs. Schneider et al., 2012; Langer et al., 2012; Schlünzen et al., 2016). A noninteracting gas expands with a a constant Schneider et al. (2012). For , decreases from its initial value until it slowly converges to an approximately constant value once becomes comparable to the bandwidth.
This behavior is explained by the large effective mass of doublons in the limit , where perturbation theory results in an effective hopping matrix elements for (see, e.g., Ref. Winkler et al., 2006). As a consequence, the doublons become inert on the accessible time scales and the system remains largely in a weakly-correlated state that is essentially a product state in the core region Heidrich-Meisner et al. (2009); Kajala et al. (2011); Bolech et al. (2012) (see also Refs. Hofmann and Potthoff, 2012; Muth et al., 2012; Rausch and Potthoff, unpublished). Therefore, is dominated by the few atoms that expand after some doublons have dissolved into single particles Kajala et al. (2011) and at long times, the expansion velocity is dominated by these fast atoms while the slow doublons do not contribute Vidmar et al. (2013).
For the applicability of the NEGF methods, this dynamical freezing of a build-up of correlations as implies that the methods become more accurate again, since the wave-functions acquire a simpler structure than at weak . This explains the a priori counterintuitive observation that the numerical deviations of the NEGF methods compared to DMRG (see the discussion of Fig. 4 and Fig. 5) become smaller as increases, even though the NEGF techniques are by construction weak-coupling methods. The regime of can much easier be accessed by DMRG with longer times becoming accessible Heidrich-Meisner et al. (2009), demonstrating the usefulness of NEGF and DMRG as complementary approaches for weak and strong coupling, respectively.
Interestingly, for large interaction strength both DMRG and the HF-GKBA propagation show oscillations in the expansion velocity with similar frequency. In contrast, for the two-time TMA calculations an onset of oscillations is seen only for , in all other cases the expansion velocity quickly damps monotonically approaching an asymptotic value. This is consistent with earlier observations that two-time propagations of the KBEs for strongly excited small Hubbard systems can be accompanied by unphysical damping in the density evolution von Friesen et al. (2009, 2010); Hermanns et al. (2014), as was noted in the introduction Sec. I. Since the initial confinement quench in our simulations constitutes such a strong excitation it is very likely that the missing of the oscillations of the expansion velocity in the two-time simulations are associated with this artificial damping.
On the other hand, the HF-GKBA is known to remove the artificial damping in strongly excited small systems Hermanns et al. (2014). Therefore, it is not surprising, that in the present setup, the HF-GKBA simulations exhibit better agreement with the DMRG for intermediate times, including the reproduction of the oscillations of the expansion velocity. This is particularly the case for small and moderate couplings, . For larger couplings, the long-time asymptotics of the expansion velocity of the two-time simulations is closer to the DMRG than the HF-GKBA result. This behavior is also consistent with the earlier observations for the evolution of the density profile and the double occupancy. This complementarity of the performance of the two-time and the HF-GKBA NEGF simulations are a particularly attractive feature.
Therefore, having both NEGF results at hand, allows one to estimate, e.g., the value of the asymptotic expansion velocity. For all , the DMRG solution of this asymptotic value lies within the NEGF results. Utilizing this observation, one can extract the exact value of with a relative error of , for all .
The experiment Schneider et al. (2012) used a different measure for the expansion velocity derived from the time evolution of the half-width-at-half-maximum, called core expansion velocity. In Ref. Schlünzen et al., 2016, a direct comparison of numerical results for this core expansion velocity to experimental data of Ref. Schneider et al., 2012 was presented, with a very good agreement. Our analysis of the errors of densities as a function of and time in the different NEGF schemes further corroborates the validity of the NEGF data used in that comparison.
We close this discussion by noting that in principle, it should be possible to compute the asymptotic expansion velocities from the Bethe ansatz, along the lines of Refs. Bolech et al., 2012; Mei et al., 2016. For instance, was computed for with an excellent agreement with DMRG results Mei et al. (2016). The extension of Ref. Mei et al., 2016 to initial densities is left for future research.
III.2 Sudden expansion in 1D: asymmetric case
It is now interesting to further investigate whether the observed accuracy and the complementary behavior of single-time and two-time NEGF simulations is just a special case of the symmetric expansion. To this end we now consider a modified setup [cf. Fig. 1(b)] where the confinement quench gives rise to a density expansion in only one direction. The results are presented in Fig. 7(a) for fermions and . The respective evolution of the double occupancy is shown in Fig. 7(b). As one can see, the results obtained by all considered methods lie very close to each other. To better distinguish between the particular profiles, we show the deviations to the DMRG results in Fig. 8. The subfigures Fig. 8(a)-(d) correspond to different interaction strengths .
As in the symmetric 1D setup before we observe a complementary behavior of the two-time result and the HF-GKBA. First, it is striking that again both approximations exhibit opposite deviations from the DMRG: while the two-time results show a slightly too fast expansion, the HF-GKBA results are retarded. Correspondingly, the deviations of the local densities from the DMRG results have opposite signs: the HF-GKBA (two-time TMA) densities are above (below) the DMRG result, on the originally doubly occupied sites and, vice versa, for the unoccupied sites. Also, the deviations have a similar dependence on the coupling strength as in the symmetric case, Sec. III.1: For the HF-GKBA that exhibits smaller density deviations than the two-time result for all considered , the maximum deviation is found for intermediate coupling strengths (). In contrast, for , the two-time result shows large deviations at the edge of the occupied region.
After considering the densities we again compute the width of the expanding particle cloud, using Eq. (20) (see Fig. 9(a)). This quantity confirms the observations made before for the density: compared to the DMRG, the expansion of the particle cloud is slightly accelerated (decelerated) for the two-time (GKBA) simulations. With increasing the two-time result becomes more accurate than the HF-GKBA. For large couplings, , the deviations between two-time and single-time approximations vanish. In this limit, also the two-time result for the expansion is retarded, in comparison to the DMRG. This analysis indicates that for a combination of two-time and single-time simulations is able to reproduce the cloud size with a relative error not exceeding . However, for large couplings, the inaccuracies grow and appear to arise from the inadequacy of the underlying -matrix approximation in the particle-particle channel. This makes it necessary to study additional many-body approximations.
The considered system of sites and initially doubly occupied sites obeys a high symmetry between the electron- (hole-) density on site and the hole- (electron-) density on site . To generalize our findings, we also present results for a system of initially occupied sites, for which this symmetry is broken. The corresponding widths of the particle cloud are shown in Fig. 9(b). As one can see, all trends agree with the previous results and especially the enclosing behavior of the NEGF methods seems not to depend on the symmetry of the system.
III.3 Asymmetric sudden expansion on a two-leg ladder
As mentioned before, the generalization to higher system dimensions constitutes a challenging problem to DMRG due to the additional degrees of freedom in the correlation growth. Therefore, only simple 2D toy models have been simulated so far with time-dependent DMRG, including the expansion of strongly interacting bosons on few-leg ladders and in small 2D clusters Hauschild et al. (2015). Since the NEGF method is not restricted with respect to dimensionality it is very interesting to compare the performance of both methods on a two-leg ladder to see if the good agreement of the previous 1D analysis can be confirmed for higher dimension.
As in Sec. III.2, we consider an asymmetric expansion setup, now with a ladder of ten rungs, the leftmost five of which are initially doubly occupied (cf. Fig. 1(c)). The resulting density evolution is illustrated in Fig. 10(a) where the density distributions for several time steps are shown in a simulation using the HF-GKBA+TMA approach with . As in the 1D case, the particles tend to move to the right. To quantify the growth of the width of the particle cloud we again use the reduced radius of Eq. (20) the time evolution of which is shown in Fig. 10(b) for all considered methods and . As one can see, the behavior is very similar to the 1D case (cf. Fig. 9). The slowing down of the expansion for increasing interaction strength is well predicted by all considered methods for small interaction strengths, whereas for larger , the DMRG curve lies between the NEGF results.
It should be mentioned that the evolutions of the reduced radius for all share a common short-time phase (this is also present in 1D but becomes more apparent on the ladder), for which behaves like the ideal system. This phase shortens with increasing interaction strength, which is due to the build-up of correlations. The behavior is similar for the symmetric expansion setup (cf. Sec. III.1) for which the dependence of the early expansion phases and the connection the onset of correlations are analyzed in detail in Refs. Schlünzen and Bonitz, 2016; Schlünzen et al., 2016.
III.4 Relaxation of charge-density wave state of doublons
We now turn to the fourth setup that is depicted in Fig. 1(d), an alternating sequence of doubly occupied and empty states corresponding to a charge-density wave. There have been a number of experiments starting from similar initial states with both fermions Pertot et al. (2014) or bosons Trotzky et al. (2012), but mostly of the type. Theoretically, there is much interest in the decay of charge-density waves or initial states with perfect Néel order in the Fermi-Hubbard model, with previous work on both its 1D version Enss and Sirker (2012); Bauer et al. (2015) and for higher-dimensional systems (see, e.g., Ref. Balzer et al., 2015). The decay from the bosonic version of our initial state was studied in Ref. Carleo et al., 2012.
Here the dynamics is governed by a short-time process in which particles move into the empty sites, provided that , which is the band-width. After that, a spreading and build-up of correlations sets in, for which the relevant velocity is typically strongly dependent on Cheneau et al. (2012); Enss and Sirker (2012); Bauer et al. (2015).
A useful quantity for the dynamics is the occupation imbalance which is defined as the difference of the densities on all even and all odd sites,
[TABLE]
where () sums up all densities of the even (odd) sites. The imbalance starts from and is then expected to decay. The results for fermions and five different couplings are shown in Fig. 11(a).
In the figure we compare DMRG results to NEGF simulations using a third-order approximation for the selfenergy. We show results for the single-time limit, i.e. after applying the HF-GKBA111For NEGF, we leave out the full two-time propagation results due to the artificial damping (cf. Fig. 6) that completely cancels out all occurring oscillations in the CDW setup.. The agreement is excellent for small and moderate couplings. Only once the interaction strength becomes as large as , small deviations are visible which grow for . This is not surprising because the third-order approximation does not capture higher order corrections. It contains, however, the third-order electron-hole diagram (cf. Fig. 2) which becomes essentially important at half filling. Therefore, TOA simulations are superior to the -matrix calculations with respect to the description of CDW dynamics. Similar trends are seen in the total double occupancy which is displayed in Fig. 11(b). The dependence of errors of the NEGF methods on (i.e., an increase as becomes order of the bandwidth) in this example is the expected generic behavior since these methods are by construction weak-coupling approaches. DMRG works particularly well for in such problems Enss and Sirker (2012), illustrating the complementary strength of NEGF versus DMRG time evolutions.
IV Summary and outlook
IV.1 Summary of main results
In this paper we thoroughly investigated the accuracy and applicability range of NEGF-based approaches in the description of the complex and correlated electron dynamics in strongly excited large Hubbard chains. The basis for this benchmark analysis were DMRG simulations performed for the same setups. Based on this analysis for the selected four setups we may conclude that NEGF simulations are reliable and accurate, thereby fully confirming earlier comparisons to exact diagonalization results for small clusters. Thus, NEGF simulations have predictive power, far beyond the present systems and situations. More precisely, our conclusions can be summarized as follows:
The quality of the NEGF results crucially depends on the choice of the selfenergy, , which is clearly dictated by the physical situation. For weak coupling, (not studied here, cf. Hermanns et al. (2014)), the second Born approximation is adequate. For moderate coupling, , proper approximations are the particle-particle -matrix (TMA) and the third-order approximation (TOA), which was investigated here for the first time. 2. 2.
For , the choice of depends on the local densities (filling): for densities close to zero (or close to one), TMA is appropriate, confirming earlier results for small clusters von Friesen et al. (2009), whereas near half filling TOA is significantly more accurate, as it contains contributions neglected in TMA. 3. 3.
For the present system sizes the HF-GKBA (with the relevant selfenergy) yields more accurate results compared to the corresponding two-time simulations (due to the artificial damping observed in the latter). While the envelopes of global dynamical quantities (energies, cloud size, expansion velocity, density imbalance etc.) are captured very accurately, oscillations of these quantities are reproduced only qualitatively, for . 4. 4.
Full two-time NEGF simulations can be used as a support of the HF-GKBA data, as typically the exact result is enclosed between the single-time and two-time simulations. One half of the difference of the two yields a (conservative) estimate of the numerical error, at least for couplings .
Based on this analysis of the NEGF capabilities, the main outcome of this paper is that NEGF and DMRG have, to a large degree, complementary strengths and limitations, with respect to the interaction strength. If does not exceed the bandwidth of the system, the NEGF approach has predictive power even for long-time propagations, and it is directly applicable to 2D and 3D systems Schlünzen et al. (2016). In contrast, the exponential spreading of entanglement narrows the DMRG approach to very short 1D simulations (somewhat larger times can be reached than presented here by using more states and possibly also by using variants of the algorithm Enss and Sirker (2012); Haegeman et al. (2011)). On the other hand, If is larger than the bandwidth, the NEGF approach, in its present form, does not describe the dynamics properly, due to the built-in perturbative character of the approximations, whereas the DMRG method provides the exact dynamics for rather long times, although being limited to 1D and small 2D systems.
IV.2 Complementarity of NEGF and DMRG simulations: A case study
To illustrate this complementarity and the reach of the two methods, we have performed additional long-time simulations and investigated the system size-dependence of the simulations, for the CDW setup (cf. Sec. III.4). As a particularly sensitive quantity, we introduce the average double occupation
[TABLE]
The time evolution of is shown in Fig. 12 for different chain lengths , ranging from to , corresponding to particles, for (a) , (b) and (c) . The NEGF simulations use the HF-GKBA with -matrix selfenergies. TOA simulations had a stability problem and are included only for shorter times (see insets).
Starting from the case of [cf. Fig. 12(a)], one can see that the short-time dynamics () of all considered chains are very similar. After the transient oscillations have decayed a quasi-stationary regime is observed. However, at some time the double occupation abruptly increases again [cf. Fig. 12(a)]. These revivals occur periodically, becoming weaker with increasing system size , and their periods increase nearly linearly with . This indicates particles with a critical velocity that pass through the entire system. It should be noted that, after the revivals, starts to fluctuate inhomogeneously with an amplitude that decreases with . We note that similar revivals and system-size dependencies were reported in Ref. Kiendl and Marquardt, unpublished. While for the DMRG simulations are restricted to very short times, e.g., for ( for ), the NEGF simulations easily allow to reach and more, for all chain lengths. The excellent agreement with DMRG is striking, suggesting that also the long time results are reliable.
Consider now the case , Fig. 12(b). Here the complementary behavior of the two approaches becomes particularly obvious. While DMRG simulations show an improved performance compared to and reach times of the order of , for (, for ), NEGF simulations still reach times of the order of , however, it is more difficult to achieve convergence. Simulations with the most accurate TOA selfenergy are stable only for short times, on the order of , for , similar to DMRG, and are in good agreement with the latter [cf. inset of Fig. 12(b)]. Long-time simulations are presently possible only with -matrix selfenergies which, however, exhibit a small upshift, compared to DMRG. Interestingly, the -dependent revivals that were observed in the NEGF simulations for are confirmed here as well by the NEGF results and, even more clearly in the DMRG runs.
Finally, in Fig. 12(c) we show results for . Here, accurate long-time evolutions for small systems can be easily performed with the DMRG method. In contrast, the available NEGF approximations are not accurate enough and show poor convergence for long times (results from NEGF simulations not included in the figure).
IV.3 Outlook
After this analysis of the NEGF approach and the illustration of the interesting complementarity with DMRG we briefly discuss questions that will be of interest for future developments. First, it will be very important to extend the arsenal of selfenergies. One important improvement will be achieved by extending the -matrix approximation by including electron-hole contributions and by including dynamical screening effects (FLEX approximation Hermanns (2016)). These choices for the selfenergy will help to extend the interaction range where NEGF properly describes the dynamics. Another way to access larger is to derive novel selfenergies via a perturbation expansion with respect to , i.e., by starting from a Hamiltonian that includes doublons directly. Finally, it would be interesting to further improve the GKBA. While it was found to cure the artificial damping problems of two-time simulations, the dynamics is often too weakly damped. This behavior should improve if one uses instead of HF-propagators correlated propagators Kwong et al. (1998); Latini et al. (2014).
While in this paper, 1D and small quasi-1D systems have been investigated for uncorrelated initial states, it will be interesting to extend the present method comparison to more complex, correlated initial states (including the ground states) as well as to larger 2D and 3D systems. It will also be interesting to analyze the dependence of the dynamics on the sign of the interaction Schneider et al. (2012); Schlünzen et al. (2016), to investigate disordered setups Basko et al. (2006); Gornyi et al. (2005); Nandkishore and Huse (2015); Altman and Vosk (2015) and to compare the fermionic simulations to those for bosonic lattice systems. Finally, the access of long simulation times by NEGF and DMRG for weak and strong coupling, respectively should allow one to study interesting features of the quantum-quench dynamics such as prethermalization Berges et al. (2004); Moeckel and Kehrein (2008); Eckstein et al. (2009).
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Giannetti et al. (2016) C. Giannetti, M. Capone, D. Fausti, M. Fabrizio, F. Parmigiani, and D. Mihailovic, Advances in Physics 65 , 58 (2016).
- 2Gandolfi et al. (unpublished) M. Gandolfi, L. Celardo, F. Borgonovi, G. Ferrini, A. Avella, F. Banfi, and C. Giannetti, p. ar Xiv:1609.07394 (unpublished).
- 3Langen et al. (2015) T. Langen, R. Geiger, and J. Schmiedmayer, Annual Rev. of Condensed Matt. Phys. 6 , 201 (2015).
- 4Eisert et al. (2015) J. Eisert, M. Friesdorf, and C. Gogolin, Nature Phys. 11 , 124 (2015).
- 5Gogolin and Eisert (2016) C. Gogolin and J. Eisert, Rep. Prog. Phys. 79, 056001 (2016) 79 , 056001 (2016).
- 6Schneider et al. (2012) U. Schneider, L. Hackermüller, J. P. Ronzheimer, S. Will, S. Braun, T. Best, I. Bloch, E. Demler, S. Mandt, D. Rasch, et al., Nature Phys. 8 , 213 (2012).
- 7Will et al. (2015) S. Will, D. Iyer, and M. Rigol, Nature Communications 6 , 6009 (2015).
- 8Pertot et al. (2014) D. Pertot, A. Sheikhan, E. Cocchi, L. A. Miller, J. E. Bohn, M. Koschorreck, M. Köhl, and C. Kollath, Phys. Rev. Lett. 113 , 170403 (2014).