NESSi
v1.0.2
The NonEquilibrium Systems Simulation Library
|
Simulating the time evolution of strongly driven many-body quantum systems is challenging because behavior distinct from the equilibrium properties can emerge on vastly different timescales.
Such simulations are however of interest in a broad variety of contexts:
NESSi
library. For details, see: Denis Golež, Philipp Werner, and Martin Eckstein, Photo-induced gap closure in an excitonic insulator. Phys. Rev. B 94, 035121 (2016). Field-theoretical approaches based on Green's functions provide a versatile framework to derive systematic approximations for quantum many-particle systems. Green's functions directly measure spectra and occupations of electronic quasiparticles and collective excitations in the solid, and therefore give access to spectroscopic quantities like tr-ARPES. This is complementary to exact methods which give access to the full wave function, but suffer from the exponential scaling of the Hilbert space with system size. The NEGF approach, as pioneered by Keldysh, Kadanoff and Baym, extends concepts from equilibrium many-body physics (Feynman diagrams, path integrals) to the study of nonequilibrium phenomena. For a basic introduction to the Keldysh formalism, one may consider, e.g., the following books:
The Keldysh formalism provides the analytical foundation for important concepts in nonequilibrium theory, such as the quantum Boltzmann equation or fluctuation-dissipation relations. At the same time, it sets the basis for the use of diagrammatic techniques in numerical simulations of the time-evolution or nonequilibrium steady states of quantum systems. This is where NESSi
comes into play:
In very superficial terms, a non-equilibrium propagator \(G_{ij}(t,t')\) describes a two-time correlation between the creation and annihilation of excitations, typically labelled by single-particle states \(i\) and \(j\) referring to orbital, spin, ... . Formulated in real time, the equation of motion for a two-time NEGF has the typical form
\begin{equation} \label{dyson_schematic} i\partial_t G(t,t') - H_{mf}(t) G(t,t') -\int_{\text{previous time}} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!d\bar t\,\,\, \Sigma(t,\bar t)G(\bar t,t') = \delta(t,t'). \end{equation}
Here \(H_{mf}\) is an effective one-body Hamiltonian (which includes mean fields), and the self-energy \(\Sigma(t,t')\), which typically is some functional of the Green's function itself, incorporates effects due to interactions and the coupling to the environment. All quantities are matrices in orbital space. Equation \eqref{dyson_schematic} is a non-Markovian integral equation: The propagation forward in time involves a memory integration over all previous times, as illustrated in the figure below. [Note that Eq. \eqref{dyson_schematic} and the figure below is rather schematic; the precise parametrization of times used to propagate NEGFs depends on the implementation and is described in Mathematical formulation as well as Ref. [6].]
The long-time memory can make NEGF simulations computationally costly if many Green's functions with large orbital dimensions are involved. A number of different ways have been reported to precisely parametrize and solve the KB equations on a discrete time grid [7] [5] [3]. The NESSi
software package follows the implementation used successfully in the context of non-equilibrium DMFT [2], and was in part described in [4]. It provides classes representing the various types of Green’s functions, implements the basic operations on these functions and allows to solve integral equations of motion such as \eqref{dyson_schematic}. A detailed description of the implementation is provided in [6].
It was noted in the early 60'es by Keldysh, Schwinger, Kadanoff, Baym, and others that many-body theory can be formulated on a generalized time contour.
libcntr
library of NESSi
is tailored to deal with numerical problems for Green's functions on the L-shaped contour.Non-equilibrium steady states (NESS) are states in which the memory of the initial state is lost, and two-time correlations \(C(t,t')\) do not depend on the average time \(t_{av}=(t+t')/2\) but only on relative time \(t-t'\). To describe NESSs, a time contour in which the initial state is shifted to time \(t=-\infty\) is more appropriate, and this contour retains only the real-time branches \(\mathcal{C}_{1,2}\). libcntr
can describe such states only if one would simulate the relaxation from a given initial state to the NESS. While this is possible in some cases, a (much simpler) code based on a frequency-domain representation of NESS Green's functions is more appropriate. A corresponding library may be added in a later release of NESSi
.
This section summarizes a few relevant definitions, to fix the notation. For a more detailed description, see [6], Chapter 3.
All relevant expectation values and correlations functions can be naturally represented as contour-ordered expectation values, where the contour-ordering simply orders operators according to the time-direction on the contour (as indicated by the arrows in the figure above),
\begin{align} T_\mathcal{C} \big\{\hat{A}(t_1) \hat{B}(t_2) \big\} = \begin{cases} \hat{A}(t_1) \hat{B}(t_2) & : t_1 \succ t_2 \\ \xi\hat{B}(t_2) \hat{A}(t_1) & : t_2 \succ t_1 \ . \end{cases} \end{align}
Here the sign \(\xi\) is +1 (-1) if the operators are bosonic (fermionic). The basic two-time correlator used in NESSi
is
\begin{equation} \label{twotimegreens} C_{A,B}(t,t') = -i\langle T_\mathcal{C} A(t) B(t') \rangle_{\mathcal{S}} \equiv -i\frac{\text{tr} \Big[ T_\mathcal{C} e^{\mathcal{S}}\,A(t) B(t')\Big]}{\text{tr} \Big[ T_\mathcal{C}\,e^{\mathcal{S}}\Big]}. \end{equation}
Here \(\mathcal{S}\) is a time-dependent action. In the simplest case of a closed quantum system it would be the evolution with a given Hamiltonian, \(\mathcal{S}=-i\int_\mathcal{C} dt H(t)\). In this case, it is straightforward to unfold the contour-ordered expectation value into regular time-dependent expectation values, see figure below for an illustration. A more general action, such as \(\mathcal{S}=-i\int_\mathcal{C} dt dt' c(t)^\dagger \Delta(t,t')c(t') \) describes coupling to an environment, such as a particle reservoir or a bath of oscillators.
Contour correlators are typically parametrized in terms of the Keldysh components, which correspond to choosing the time arguments on the various branches of \(\mathcal{C}\), and have a well-defined meaning in terms of spectral functions and occupation functions of elementary excitations, or response functions to external perturbations. Important components for the documentation of libcntr
are:
In libcntr
, contour correlators \(C(t,t^\prime)\) are represented by the minimal set of Keldysh components \(\{C^<,C^{\mathrm{R}},C^\rceil,C^{\mathrm{M}}\}\) on \(\mathcal{C}[h,N_t,h_\tau,N_\tau]\). The main data structure is the class herm_matrix
, which is described in detail in Green functions.
Dyson integral equations are reduced to Volterra Integral (differential) equations. Our implementation is based on a numerical solution of Volterra equations which uses Gregory quadrature rules, as described in Brunner and v Houven, The numerical solution of Volterra equations (North Holland, Amsterdam 1986). A Gregory rule of degree \(k\) for an integral \(I=\int_0^{Nh} dt\ f(t)\) is an integration rule
\begin{equation} I \approx \sum_{j=0}^{m(N,k)} w_j\ f(jh) \end{equation}
on an equidistant grid which has the following property:
The detailed explanation of the numerical solution is given in Ref. [6]
cntr::convolution
, cntr::vie2
and cntr::dyson
. The numerical error of the \(k\)-th order Gregory rule is \({\cal O}(N^{-(k+2)})\). For example, \(k=0\) is the simple trapezoidal quadrature rule, with has an error \({\cal O}(N^{-2})\). The error scaling of the implemented integral equation is demonstrated in Test of accuracy of Dyson solver: equilibrium.Prior to its public release in 2019, NESSi
had been used extensively for more than a decade in mostly nonequilibrium DMFT related studies of correlated lattice systems. More than 60 publications report numerical simulations based on NESSi
or its core constituent libcntr
. In the following, we just provide a small glimpse of the type of problems which have been addressed.
Thermalization of a pump-excited Mott insulator. Driving a Hubbard model in the Mott insulating regime with an electric field pulse whose frequency is comparable to the on-site interaction \(U\) creates long-lived double-occupations (doublons) and holes (holons). As the system thermalizes the density of doublons changes and eventually reaches the value which corresponds to the temperature \(T_{eff}\), which can be independently calculated by measuring the energy injected into the system by the pulse. An interesting question is how fast the doublons will thermalize. Using nonequilibrium DMFT simulations based on NESSi
it was shown that the thermalization time increases exponentially with \(U\) (or the gap size), so that even in a purely electronic model, thermalization times can be orders of magnitude longer than the intrinsic timescale (inverse hopping). Note that the relaxation pathway in at \(U=1.5\), in the metal-insulator crossover regime, is qualitatively different.
Doublon relaxation in photo-excited 1T-TaS \(_2\). Benchmarking nonequilibrium DMFT results is not a simple task, since for many problems, there is no better method available. In such a situation, it is useful to resort to the ultimate test - comparison to experiments! For the doublon dynamics in photo-excited Mott systems, such a test has been performed on 1T- \(TaS_2\), a Mott insulating material with a triangular lattice structure. In this study, time-resolved photoemission was used to track the population of doublons after an ultrashort laser pulse. Photoemission measures the lesser component of the time-dependent Greens function, a quantity that can be directly computed within nonequilibrium DMFT. The theoretical modeling was based on a triangular lattice single-band Hubbard model and the simulations were performed using NESSi
in combination with the PPSC library. Theory and experiment agree almost quantitatively under the assumption of a slightly hole-doped sample.
Benchmarking quantum simulators. In an other recent theory-experiment comparison, NESSi
was used to benchmark the production or suppression of doublons in periodically driven cold atom systems. In these “quantum simulators’’ the Hamiltonian is well defined and essentially given by a single-band Hubbard model, so that a direct comparison with nonequilibrium DMFT results is meaningful. Within the accuracy of the experiments, the DMFT calculations based on strong-coupling perturbative solvers, as implemented in the PPSC library, produced consistent results for the dynamics, both in the resonant and off-resonant driving regime.
previous page Getting started with NESSi – next page Manual