NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library
Physics background

Non-equilibrium quantum many-particle physics

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.

Time scales for different light-induced phenomena in lattice systems:
Timescales001v01.png

Such simulations are however of interest in a broad variety of contexts:

  • By driving condensed matter with tailored light one can engineer novel quantum phases. Light-induced superconductivity or Floquet states are among the tantalizing examples. See, e.g., D. N. Basov, R. D. Averitt, and D. Hsieh. Towards properties on demand in quantum materials. Nature Materials, 16 1077 (2017).
  • Analog quantum simulation platforms allow to explore genuine non-equilibrium phenomena, such as the nontrivial dynamics of quantum systems at the boundary between integrable and ergodic behavior.
  • Dissipative driven quantum systems are relevant for quantum transport and nanotechnology, including quantum computing architectures.
  • Time-resolved pump-probe spectroscopy can be used to analyze the interplay of electronic quasiparticles and collective excitations on their microscopic timescale. See. e.g., Claudio Giannetti, Massimo Capone, Daniele Fausti, Michele Fabrizio, Fulvio Parmigiani, and Dragan Mihailovic. Ultrafast optical spectroscopy of strongly correlated materials and high-temperature superconductors: a non-equilibrium approach. Advances in Physics, 65 58, (2016).
Example: Simulated time- and angular-resolved photoemission spectrum (tr-ARPES) of an excitonic insulator:
exciton_pes.gif
The movie illustrates the possibility to obtain the electronic structure out of equilibrium using tr-ARPES: The system is excited with an electric field pulse E(t), and the ARPES spectrum is recorded with a probe pulse at different delays with respect to the pump (upper panel). tr-ARPES reveals multiple dynamical processes: A filling of the conduction bands is followed by a broadening of the bands due to a change of the quasiparticle lifetime, internal relaxation, a closing of the gap (insulator-to-metal transition) due to a non-thermal melting of the exciton condensate, and thermalization. — Time-dependent GW simulations performed using the 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).

↑ to the top of the page ↑

Keldysh Formalism and nonequilibrium Green's functions

Non-equilibrium Green's functions (NEGF)

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:

Kadanoff-Baym (KB) equations: Equations of motion with memory

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].]

Causal time-propagation of a two-time NEGF on an equally spaced time-grid:
In order to solve Eq. \eqref{dyson_schematic} at a given time-slice in the two-time plane (red time-points), the Green's function and the self energy at all previous time-points (shaded yellow) is needed.

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].

↑ to the top of the page ↑

Mathematical formulation

The Keldysh contour

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.

  • The L-shaped time-contour \(\mathcal{C}\) shown in the figure below is best adapted to describe situations in which the time-evolution starts from a mixed or pure state which is representable by a thermal density matrix \(\rho=e^{-\beta H(0_-)}/Z\). By default, this refers to a system which is initially in thermal equilibrium and then subject to perturbations. Because the operator \(H(0_-)\) can be different from the Hamiltonian \(H(t)\) for \(t\ge 0\) which determines the time evolution, this representation still allows for a large variety of nonequilibrium set-ups. The 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.

    L-shaped Kadanoff-Baym contour:
    contour_discrete.png
    L-shaped Kadanoff-Baym contour \(\mathcal{C}\) in the complex time plane, running from over the forward branch \(\mathcal{C}_1\), backward branch \(\mathcal{C}_2\) to the imaginary (Matsubara) branch \(\mathcal{C}_3\). Grey dots sketch the discretization on the real-time branches, while orange dots represent the discretization of the Matsubara branch. For numerical calculations, the contour arguments are discretized according to the sketch in the figure. The real-time branches are divided into \((N_t+1)\) equidistant points \(t_n=n h\), \(n=0\), \(\dots\), \(N_t\), while \(\tau_m= m h_\tau\), \(m=0\), \(\dots\), \(N_\tau\) with \(\tau_0=0^+\), \(\tau_{N_\tau}=\beta^-\) samples the Matsubara branch. The corresponding discretized contour is denoted by \(\mathcal{C}[h\), \(N_t\), \(h_\tau\), \(N_\tau]\)
Contour-ordered Green's functions

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.

Unfolding the contour-ordered expectation value into regular time-dependent expectation values:
For an action \(\mathcal{S}=-i\int_\mathcal{C} dt H(t)\), Eq. \eqref{twotimegreens} can be interpreted as a regular time-dependent correlation function by discretizing the contour-ordered expectation value \(e^{\mathcal{S}}\) into the time-evolution over infinitesimal intervals.

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:

Important Keldysh components (For a detailed list and explanation, see [6], Chapter 3).
  • Lesser component: \(C_{A,B}^<(t,t')= C_{A,B}(t_+,t'_-) = -i\xi\langle B(t')A(t)\rangle\)
  • Greater component: \(C_{A,B}^>(t,t')= C_{A,B}(t_-,t'_+) = -i\langle A(t)B(t')\rangle\)
  • Retarded component: \(C_{A,B}^R(t,t')= \theta(t-t') [C_{A,B}^>(t,t')-C_{A,B}^<(t,t')(t_-,t'_+)] = -i\langle [A(t),B(t')]_\xi\rangle\). The retarded component is related to response functions (Kubo formula) or spectral functions. Lesser and greater functions correspond to the occupied and unoccupied density of states, respectively. (For example, for electrons, with \(A=c_j\) and \(B=c_j^\dagger\), \(C^<\) encodes the result of an idealized photoemission experiment.
  • The imaginary-time component \(C(-i\tau,-t\tau')\) depends on the time difference \(\tau-\tau'\) only, and is related to the standard imaginary time-ordered Green's function \(C(\tau)=-\langle T_\tau A(\tau)B(0)\rangle\) (Matsubara function) like \(C_{A,B}^M(\tau)= -iC_{A,B}(-i\tau,0)\)

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.

↑ to the top of the page ↑

Numerical solution and accuracy

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:

  • For \(N<k, m(N,k)=k\); the integration is based on polynomial interpolation of the function in the interval \([0,kh]\).
  • For \(N>k\), \(m(N,k)=N\); only the weights \(w_j\) for \(j=0,\dots k\) and \(w_j=N,\dots N-k\) differ from 1.

The detailed explanation of the numerical solution is given in Ref. [6]

Accuracy:
We define a routine to be of order \(k\) ( \(k=0,1,2,3,\dots\)) if it uses a Gregory rule of degree \(k\). In practice, this implies that the first \(k\) time steps must be treated differently, and explained in the documentation of 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.

Previous use of NESSi

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.

    Dynamics and relaxation times of doublons:
    thermalization_small.png
    The left panel shows the difference between the doublon density of the pulse-excited system and the value expected after thermalization, for indicated values of \(U\). The right hand panel plots the relaxation times extracted from these curves as a function of \(U\), and demonstrates the exponential dependence. For details, see M. Eckstein and P. Werner, Thermalization of a pump-excited Mott insulator, Phys. Rev. B 84, 035122 (2011).


  • 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.

    Time-resolved photoemission signal:
    1TTaS2small.png
    The left panel shows the measured photoemission spectra of pumped \(1T-TaS_2\) (offset for clarity), with the positions of the lower Hubbard band (LHB) and upper Hubbard band (UHB) indicated. The right panels show the simulated results for a half-filled Mott insulator ( \(n=1\)) and a slightly hole-doped Mott insulator ( \(n=0.95\)). For details, see M. Ligges et al., Ultrafast Doublon Dynamics in Photoexcited \(1T-TaS_2\), Phys. Rev. Lett. 120, 166401 (2018).


  • 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.

    Doublon production for near-resonant driving:
    coldatom_small.png
    The two panels show the experimental (left) and DMFT (right) results for the doublon production in a shaken brickwall lattice, with \(K_0\) indicating the dimensionless driving amplitude. Dynamics beyond 10 ms is influenced by atom losses and is not accessible in the simulations. Error bars on the theoretical curves reflect the uncertainty in the experimental entropy estimation. For details, see K. Sandholzer et al., Quantum simulation meets nonequilibrium DMFT: Analysis of a periodically driven strongly correlated Fermi-Hubbard model, arxiv:18011.12826

↑ to the top of the page ↑

previous page Getting started with NESSinext page Manual