NESSi
v1.0.2
The NonEquilibrium Systems Simulation Library

Simulating the time evolution of strongly driven manybody 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, Photoinduced gap closure in an excitonic insulator. Phys. Rev. B 94, 035121 (2016). Fieldtheoretical approaches based on Green's functions provide a versatile framework to derive systematic approximations for quantum manyparticle 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 trARPES. 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 manybody 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 fluctuationdissipation relations. At the same time, it sets the basis for the use of diagrammatic techniques in numerical simulations of the timeevolution or nonequilibrium steady states of quantum systems. This is where NESSi
comes into play:
In very superficial terms, a nonequilibrium propagator \(G_{ij}(t,t')\) describes a twotime correlation between the creation and annihilation of excitations, typically labelled by singleparticle states \(i\) and \(j\) referring to orbital, spin, ... . Formulated in real time, the equation of motion for a twotime 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 onebody Hamiltonian (which includes mean fields), and the selfenergy \(\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 nonMarkovian 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 longtime 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 nonequilibrium 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 manybody 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 Lshaped contour.Nonequilibrium steady states (NESS) are states in which the memory of the initial state is lost, and twotime correlations \(C(t,t')\) do not depend on the average time \(t_{av}=(t+t')/2\) but only on relative time \(tt'\). 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 realtime 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 frequencydomain 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 contourordered expectation values, where the contourordering simply orders operators according to the timedirection 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 twotime 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 timedependent 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 contourordered expectation value into regular timedependent 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 welldefined 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 pumpexcited Mott insulator. Driving a Hubbard model in the Mott insulating regime with an electric field pulse whose frequency is comparable to the onsite interaction \(U\) creates longlived doubleoccupations (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 metalinsulator crossover regime, is qualitatively different.
Doublon relaxation in photoexcited 1TTaS \(_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 photoexcited Mott systems, such a test has been performed on 1T \(TaS_2\), a Mott insulating material with a triangular lattice structure. In this study, timeresolved photoemission was used to track the population of doublons after an ultrashort laser pulse. Photoemission measures the lesser component of the timedependent Greens function, a quantity that can be directly computed within nonequilibrium DMFT. The theoretical modeling was based on a triangular lattice singleband 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 holedoped sample.
Benchmarking quantum simulators. In an other recent theoryexperiment 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 singleband Hubbard model, so that a direct comparison with nonequilibrium DMFT results is meaningful. Within the accuracy of the experiments, the DMFT calculations based on strongcoupling perturbative solvers, as implemented in the PPSC library, produced consistent results for the dynamics, both in the resonant and offresonant driving regime.
previous page Getting started with NESSi – next page Manual