Manual for steady-state Kadanoff-Baym equations

Steady state Keldysh formalism 

Another framework for solving the KB equations is the Keldysh formalism for non-equilibrium steady states (NESSs). In comparison to two-time functions defined on the three-legged KB contour \(\mathcal C = \mathcal C_1 \cup \mathcal C_2 \cup \mathcal C_3\), the main assumption in the steady-state formalism is the absence of correlations with the initial noninteracting state on the imaginary time branch \(\mathcal C_3\), such that the mixed self-energies \(\Sigma^{\rceil}\) and \(\Sigma^{\lceil}\) vanish. The vertical branch can then be shifted to \(t = - \infty\) and is eliminated from the equations, such that time arguments are restricted to the two-branch contour \(\mathcal C_K = \mathcal C_1 \cup \mathcal C_2\). Second, any contour-ordered two-time Green’s function \(G(t,t') = -i \langle \mathcal{T}_{\mathcal C} c(t) c^{\dagger}(t') \rangle\) in a NESS exhibits time-translational invariance \(G(t,t') = G(t-t')\), so that a numerical treatment is possible in Fourier representation.

Green’s functions

To represent NESS Green’s functions, we keep the retarded component \(G^R(t)\) and lesser component \(G^<(t)\) on an equidistant real-time grid as the two non-redundant components. Frequency-dependent Green’s functions are given by the spectral representation

(27)\[G^R(t) = -i\theta(t) \int d\omega \, A(\omega) e^{-i\omega t},\]

(28)\[G^R(\omega+i0) = \int_0^\infty dt \, G^R(t) e^{i(\omega+i0) t},\]

(29)\[G^<(t) = \int \frac{d\omega}{2\pi} \, G^<(\omega) e^{-i\omega t},\]

(30)\[G^<(\omega) = \int dt \, G^<(t) e^{i\omega t},\]

with the spectral function

\[A(\omega) =-\frac{1}{2\pi} \left(G^R(\omega+i0)-G^R(\omega+i0)^\dagger\right).\]

In equilibrium at temperature \(1/\beta\), the fluctuation dissipation function implies

(31)\[G^<(\omega)= -\xi 2\pi i F_\mu(\omega) A(\omega),\]

where the sign \(\xi\) is \(\xi=\pm1\) for bosonic (fermionic) Green’s functions, and the distribution function is \(F_\mu(\omega) =1/(e^{\beta(\omega-\mu)}-\xi)\).

Dyson equation

In order to solve the Dyson equation in the steady state, we rewrite the Kadanoff-Baym equations in the time-translationally invariant case. The equation for the retarded Green’s function becomes

(32)\[G^R(\omega+i0^{+}) = \big[\omega+i0^{+} - \epsilon-\Sigma^R(\omega+i0^{+})\big]^{-1}\]

in Fourier representation and the lesser component is given by

(33)\[G^<(\omega) = G^R(\omega+i0) \Sigma^<(\omega) G^R(\omega+i0)^\dagger.\]

Overview 

The current implementation of the steady state extension to NESSi has two versions, which are structured into two different namespaces ness and ness2. We suggest to the user to start with the newer version ness2, it covers and extends all functionalities of the older version. For the sake of completeness, the old version ness is also documented on this page.

`ness2`	Namespace for current implementation of the NESS code, recommended version.
`ness`	Namespace for old implementation of the NESS code.

ness2 

This is a quick overview over the most important classes and functions in the ness2 namespace:

`ness2::fft_array`	Basic class containing the data structure for holding arrays of matrices for frequency and time domain, supports element access and algebraic manipulations.	Section ness2::fft_array
`ness2::herm_matrix_ness`	Class to represent steady state Green’s function (lesser and retarded component), supports element access, Fourier transforms, equilibrium Green’s function construction and reading and writing to file and two-time functions and other utilities.	Section ness2::herm_matrix_ness
`ness2::dyson`	Solve the Dyson equation in the steady state.	Section Dyson equation
`ness2::green_equilibrium_ness`	Construct equilibrium Green’s function from a given density of states.	Section Equilibrium Green’s functions
`ness2::Bubble1_ness`	Calculate particle-hole bubble diagram for two given Green’s functions.	Section Diagram utilities
`ness2::Bubble2_ness`	Calculate particle-particle bubble diagram for two given Green’s functions.	Section Diagram utilities
`ness2::cntr2ness`	Write from two-time function into a steady state Green’s function, also vice versa is supported.	Section Data exchange with two-time functions

ness 

This is a quick overview over the most important classes and functions in the ness namespace:

`ness::GF`	Class to represent steady state Green’s function (lesser and retarded component), supports element access, Fourier transforms, equilibrium Green’s function construction and reading and writing to file and two-time functions and other utilities (either frequency or time domain).	Section ness::GF
`ness::GF_pair`	Class to represent steady state Green’s function (lesser and retarded component), both in frequency and time domain, helper class `fft_solver` provides Fourier transform.	Section ness::GF_pair
`ness::fft_solver`	Helper class to provide Fourier forward and back transform.	Section ness::fft_solver
`ness::dyson_ness`	Solve the Dyson equation in the steady state.	Section Dyson equation
`ness::green_equilibrium_ness`	Construct equilibrium Green’s function from a given density of states.	Section Equilibrium Green’s functions
`ness::Bubble1_ness`	Calculate particle-hole bubble diagram for two given Green’s functions.	Section Diagram utilities
`ness::Bubble2_ness`	Calculate particle-particle bubble diagram for two given Green’s functions.	Section Diagram utilities

Green’s functions 

`ness2::fft_array`

Overview

class

ness2::fft_array

This class is a basic data structure for an object C containing two arrays C.time_ and C.freq_, which hold Nft square matrices stored as consecutive complex numbers. The class provides plans to perform a (non-normalized) discrete Fourier transform

(34)\[\mathcal{F}\left\{ C \right\}[l] = \sum_{j=0}^{N_{\rm ft}-1} e^{i 2 \pi j l / N_{\rm ft}} \, C[j],\]

(35)\[\bar{\mathcal{F}}\left\{ C \right\}[l] = \sum_{j=0}^{N_{\rm ft}-1} e^{-i 2 \pi j l / N_{\rm ft}} \, C[j].\]

using the FFT algorithm. In addition, the class offers routines for standard element access and basic algebraic operations. The ness2::fft_array models an array that is periodic in both time and frequency dimensions, with a period of Nft. This periodic wrapping is handled automatically when accessing array elements. For most efficient performance of the FFT, the size Nft should be a power of two. An object of the class is characterized by the following two arguments:

Nft (integer): FFT array length.
size (integer): matrix size (size x size).

Constructors

`fft_array()`	Default constructor, with `Nft = 4`, `size = 1` and `FFTW_FLAG = FFTW_ESTIMATE`.
`fft_array(int Nft, int size, unsigned FFTW_FLAG = FFTW_ESTIMATE)`	Standard constructor.
`fft_array(const fft_array& other)`	Copy constructor given a `ness2::fft_array` object `other`.

The constructor has an optional third argument FFTW_FLAG, which chooses the plan for the construction of the plan in the FFTW library used to do the Fourier transforms (default is FFTW_ESTIMATE).

Accessing individual elements

General element access:

The following member functions read and set elements of the time (ness2::fft_array::time_) or frequency (ness2::fft_array::freq_) part of the ness2::fft_array, specified by the domain argument, which can take the values ness2::fft_domain::freq or ness2::fft_domain::time. The parameter M is a Matrix type (complex eigen matrices are supported). For [set|get]_element, time or frequency arguments outside the interval [0, Nft-1] are mapped back into this interval by adding an integer multiple of Nft.

`C.set_element(i,M,domain)`	\((C_{\rm time\|freq})[i]\) set to matrix `M`.
`C.get_element(i,M,domain)`	Matrix `M` set to \(C_{\rm time\|freq}[i]\), `M` is automatically resized as needed.

To set a matrix element of \((C_{\rm time|freq})\) one can use the following utility, which operates on all Nft entries:

C.set_matrixelement(i1,i2, fft array &B, j1,j2, domain)

\((C_{\rm time|freq})_{i_1,i_2}\) set to \((B_{\rm time|freq})_{j_1,j_2}\).

Simple Algebra

The following methods cover basic algebraic operations on the ness2::fft_array object. The domain argument is set to ness2::fft_domain::time or ness2::fft_domain::freq to specify whether the operation applies to C.time_ or C.freq_. All routines operate on all Nft entries. The parameter M is a Matrix type (complex eigen matrices are supported).

`C.incr(fft_array &B, cplx a, domain)`	\(C_{\rm time\|freq}\) set to \(C_{\rm time\|freq}+aB_{\rm time\|freq}\).
`C.smul(cplx a, domain)`	\(C_{\rm time\|freq}\) set to \(aC_{\rm time\|freq}\).
`C.set_zero(domain)`	\(C_{\rm time\|freq}\) set to 0 for all entries.
`C.left_multiply(M, domain)`	\(C_{\rm time\|freq}\) set to \(MC_{\rm time\|freq}\).
`C.right_multiply(M, domain)`	\(C_{\rm time\|freq}\) set to \(C_{\rm time\|freq}M\).
`C.left_multiply_hermconj(M, domain)`	\(C_{\rm time\|freq}\) set to \(M^\dagger C_{\rm time\|freq}\).
`C.right_multiply_hermconj(M, domain)`	\(C_{\rm time\|freq}\) set to \(C_{\rm time\|freq}M^\dagger\).

Fourier transform

The following member functions perform the discrete Fourier transforms (34) and (35) using FFT:

`C.fft_to_time()`	Set \(C_{\rm time}\) to \(\bar{\mathcal{F}}\{C_{\rm freq}\}\) using (35) and FFT.
`C.fft_to_freq()`	Set \(C_{\rm freq}\) to \({\mathcal{F}}\{C_{\rm time}\}\) using (34) and FFT.

Example:

int N_grid = 100;
int size = 2;
double h = 0.1;
int ex_id = 25;

// create arrays
fft_array new_array(N_grid, size);
fft_array new_array2(N_grid, size);

// matrices for operations
cdmatrix A(size,size);
cdmatrix B(size,size);
cdmatrix C(size,size);
A(1,0) = 1;
A(0,1) = II; // fill A with values

// element access
new_array.set_element(ex_id,A,ness2::fft_domain::time); // read into new_array from A
new_array.get_element(ex_id,B,ness2::fft_domain::time); // read into B from new_array
new_array2.set_matrixelement(1,1, new_array, 0, 1, ness2::fft_domain::time); // write one matrix element in new_array2 from new_array

// example operations
new_array.incr(new_array2, 0.5, ness2::fft_domain::time);
new_array.smul(2.0, ness2::fft_domain::time);
new_array.left_multiply(A, ness2::fft_domain::time);
new_array.right_multiply(B, ness2::fft_domain::time);
new_array2.set_zero(ness2::fft_domain::time);
new_array.fft_to_freq(); // FFT to frequency
new_array.get_element(ex_id,C,ness2::fft_domain::freq); // e.g. read into C

Grid info

class

ness2::grid_info

This helper class describes the time and frequency grids with steps h and dw and length Nft used in ness2::fft_array and ness2::herm_matrix_ness.

`grid_info(Nft, h)`	Constructor to initialize grid with time step `h` and length `Nft`.
`grid.set_h(new_h)`	Update time step `h` to `new_h` and adjust frequency step `dw` accordingly.
`grid.set_dw(new_dw)`	Update frequency step `dw` to `new_dw` and adjust time step `h` accordingly.
`grid.time_at(n)`	Get time value \(t\) at grid index `n`.
`grid.freq_at(k)`	Get frequency value \(\omega\) at grid index `k`.
`grid.time_grid()`	Generate the full wrapped time grid (`std::vector<double>`).
`grid.freq_grid()`	Generate the full wrapped frequency grid (`std::vector<double>`).

Example:

int N_grid = 100;
double h = 0.1;
double new_h = 0.2;
double new_dw = 0.1;
int ex_id = 25;

// manipulate and show grid info
grid_info new_grid(N_grid, h); // create a grid
cout << "Number of grid points: " << new_grid.Nft_ << endl;
cout << "Time step: " << new_grid.h_ << endl;
cout << "Frequency step: " << new_grid.dw_ << endl;

new_grid.set_h(new_h); // update time step
cout << "New time step: " << new_grid.h_ << endl;
cout << "New frequency step: " << new_grid.dw_ << endl;

new_grid.set_dw(new_dw); // update frequency step
cout << "New time step: " << new_grid.h_ << endl;
cout << "New frequency step: " << new_grid.dw_ << endl;
cout << "Time at example index: " << new_grid.time_at(ex_id) << endl;
cout << "Frequency at example index: " << new_grid.freq_at(ex_id) << endl;

std::vector<double> full_time_grid = new_grid.time_grid(); // store full time grid in vector
std::vector<double> full_freq_grid = new_grid.freq_grid(); // store full frequency grid in vector

`ness2::herm_matrix_ness`

Overview

class

ness2::herm_matrix_ness

The data type ness2::herm_matrix_ness is used to represent a steady-state function G with hermitian symmetry (non-hermitian Green’s functions are not supported in the current extension of libcntr). An object of this class is initialized by:

Nft (integer): FFT array length.
size (integer): matrix size (size x size).

An object G of type ness2::herm_matrix_ness contains a pair of two members ness2::herm_matrix_ness::ret_ and ness2::herm_matrix_ness::les_ of type ness2::fft_array, representing \(G^R\) and \(G^<\) in time and frequency. Time arguments thereby correspond to an equidistant grid

(36)\[t\in \{jh: j=-N_{\rm ft}/2,...,N_{\rm ft}/2-1\},\]

and the frequencies represent the dual FFT grid

(37)\[\omega \in \{ \Delta_\omega j: j= -N_{\rm ft}/2,...,N_{\rm ft}/2-1,\Delta_\omega= \frac{2\pi}{hN_{\rm ft}}\}.\]

We require \(N_{\rm ft}\) to be even, and use the conventional layout of the arrays where \(G^{<|R}(jh)\) is represented by the \(j\)-th element of G.[les|ret]_.time_, if \(0\le j< N_{\rm ft}/2\), and by the \((j+N_{\rm ft})\)-th element if \(-N_{\rm ft}/2\le j< 0\) (analogous for frequency). For convenience we provide a small helper class ness2::fft_grid, see Grid info.

Constructors

`herm_matrix_ness()`	Default constructor, with `Nft = 4`, `size = 1`.
`herm_matrix_ness(Nft, size)`	Standard constructor.
`herm_matrix_ness(const herm_matrix_ness& other)`	Copy constructor given a `herm_matrix_ness` object `other`.

The constructor again has an optional third argument FFTW_FLAG, for the plan of the FFT construction on FFTW.

Accessing individual elements

General element access:

The following member functions read and set elements of the time or frequency part (domain argument set to ness2::fft_domain::freq or ness2::fft_domain::time) of the retarded or lesser component of the ness2::herm_matrix_ness object. The parameter M is a Matrix type (complex eigen matrices are supported). As for ness2::fft_array [set|get]_element, time or frequency (i) arguments outside the interval [0, Nft-1] are mapped back into this interval by adding an integer multiple of Nft.

`G.set_[les\|ret](i,M,domain)`	\((G_{\rm time\|freq}^{<\|R})[i]\) set to matrix `M`.
`G.get_[les\|ret](i,M,domain)`	Matrix `M` set to \(G_{\rm time\|freq}^{<\|R}[i]\).
`G.retarded()`	Return reference to `G.ret_`.
`G.lesser()`	Return reference to `G.les_`.

Simple Algebra

The same functions as for ness2::fft_array are also defined for ness2::herm_matrix_ness, the operation is performed for both retarded and lesser component.

`G.incr(herm_matrix_ness &B, cplx a, domain)`	\(G^{R,<}_{\rm time\|freq}\) set to \(G^{R,<}_{\rm time\|freq}+aB^{R,<}_{\rm time\|freq}\).
`G.smul(cplx a, domain)`	\(G^{R,<}_{\rm time\|freq}\) set to \(a G^{R,<}_{\rm time\|freq}\).
`G.set_zero(domain)`	\(G^{R,<}_{\rm time\|freq}\) set to 0 for all entries.
`G.left_multiply(M, domain)`	\(G^{R,<}_{\rm time\|freq}\) set to \(M G^{R,<}_{\rm time\|freq}\).
`G.right_multiply(M, domain)`	\(G^{R,<}_{\rm time\|freq}\) set to \(G^{R,<}_{\rm time\|freq}M\).
`G.left_multiply_hermconj(M, domain)`	\(G^{R,<}_{\rm time\|freq}\) set to \(M^\dagger G^{R,<}_{\rm time\|freq}\).
`G.right_multiply_hermconj(M, domain)`	\(G^{R,<}_{\rm time\|freq}\) set to \(G^{R,<}_{\rm time\|freq}M^\dagger\).

Integral transforms

For the Fourier transform of Green’s functions, we must distinguish between plain DFT, which can be accessed by the call [ret|les]_.to_[time|freq]() to the members of ness2::herm_matrix_ness, and approximations to the GFs integral transforms. The transformation G.transform_to_freq(h, METHOD) computes an approximation to the integral transforms (28), (30) on the frequency grid (37), assuming a time grid with timestep h. The integrals are computed with boundaries [0,T] (for \(G^R\)) and [-T, T] (for \(G^<\)), with \(T = (Nft/2 - 1)h\). If the argument METHOD is omitted or given by the keyword FFT_TRAPEZ (default), we compute the Fourier integral in the trapezoidal rule. The reverse transformation G.transform_to_time(h, METHOD) computes an approximation to the Fourier integrals (27), (29) on the time grid (36). For both directions, we also provide an implementation METHOD=FFT_CUBIC, where the integrals are approximated by the exact Fourier transform of a piecewise cubic interpolating function.

`G.integral_transform_to_freq(h, METHOD)`	Compute the frequency-domain Fourier integrals (28), (30) of the Green’s function, assuming a timestep `h`.
`G.integral_transform_to_time(h, METHOD)`	Compute the time-domain Fourier integrals (27), (29) of the Green’s function, assuming a timestep `h`.

Density matrix

The density matrix in the steady state is directly evaluated by calculating the lesser component of the Green’s function:

\[\rho_{A} = \xi i A^<(0),\]

with \(\xi=\pm1\) for bosonic (fermionic) Green’s functions. The equal-time convolution

\[\rho_{AB} = \xi_A \xi_B i (A\ast B)^<(0),\]

can be obtained from the Fourier transform of \(C=A\ast B\). The corresponding function calls are ness2::density_matrix and ness2::convolution_density_matrix:

`density_matrix(result, bosefermi, A)`	`result` (complex matrix) is set to \(\rho_{A}\), for a function `A` of type `ness2::herm_matrix_ness`.
`convolution_density_matrix( result, bosefermi, A, B, double h)`	`result` (complex matrix) is set to the steady-state \(\rho_{AB}\), for the convolution of two functions `A` and `B` of type `ness2::herm_matrix_ness`.

File I/O

File input and output from text files and also HDF5 files is supported.

Reading and writing to text files:

Text files start with an initial header line # Nft size precision containing apart from the size properties also the optional accuracy argument. The data is structured into blocks for the retarded and lesser component. Reading expects the same format as writing to text files.

`G.print_to_file(filename, precision = 12)`	Print the `ness2::herm_matrix_ness` object `G` (`ret_`, `les_` and attributes `Nft`, `size`) to a text file. Precision argument is optional.
`G.read_from_file(filename, FFTW_FLAG = FFTW_ESTIMATE)`	Read `ret_`, `les_` and attributes `Nft`, `size` from a text file and store in `ness2::herm_matrix_ness` object `G`, expects same format as `print_to_file`. Second argument `FFTW_FLAG` optional to set the FFT plan of object `G` manually.

Reading and writing to HDF5 files:

To store a ness2::herm_matrix_ness or read it from a HDF5 file, one can use the member functions ness2::herm_matrix_ness::write_to_hdf5 and ness2::herm_matrix_ness::read_from_hdf5. This stores/reads the attributes Nft, size1, and the Green’s function component’s data sorted in groups ret, les.

`G.write_to_hdf5(hid_t group_id)`	Stores the `ness2::herm_matrix_ness` (attributes and data) to a given HDF5 group.
`G.write_to_hdf5(hid_t group_id, const char *groupname)`	Stores the `ness2::herm_matrix_ness` (attributes and data) to a given HDF5 group with given groupname.
`G.write_to_hdf5(const char filename, const char groupname)`	Write data and attributes from the `ness2::herm_matrix_ness` to an HDF5 file under the given group name.
`G.read_from_hdf5(hid_t group_id)`	Reads the `ness2::herm_matrix_ness` (attributes and data) from a given HDF5 group handle.
`G.read_from_hdf5(hid_t group_id, const char *groupname)`	Reads the `ness2::herm_matrix_ness` (attributes and data) from a given HDF5 group handle with given group name.
`G.read_from_hdf5(const char filename, const char groupname)`	Read all data and attributes from an HDF5 file from a given group into the `ness2::herm_matrix_ness`.

Note

The HDF5 reading functions have a third optional argument FFTW_FLAG (default is FFTW_ESTIMATE) to set the FFT plan of object G manually.

Example:

int N_grid = 100;
int size = 2;
double h = 0.1;
int ex_id = 25;

// Green's function object
herm_matrix_ness G(N_grid, size);

// complex matrices
cdmatrix M1(size,size);
cdmatrix M2(size,size);

// Do something with G ...

G.set_zero(ness2::fft_domain::freq);
G.get_les(ex_id,M1,ness2::fft_domain::time); // read in from G into M1
density_matrix(M2, -1, G); // get (fermionic) density matrix from G, save in M2
G.integral_transform_to_freq(h, FFT_TRAPEZ); // evaluate frequency-domain Fourier integrals

// save G in hdf5 file
hid_t file_id = open_hdf5_file(flout); // char *flout points to the filename
store_int_attribute_to_hid(file_id, "N_grid", N_grid); // write N_grid to group "/N_grid"
// same for size, h ...
G.write_to_hdf5(file_id, "G"); // new group "/G" in file
close_hdf5_file(file_id);

`ness::GF`

Overview

class

ness::GF

This class contains the data structures for representing one-time contour Green’s functions \(G(t)\) on the equidistantly discretized real-time Keldysh contour \(\mathcal{C}_K\) with points \(\{i h: i=0,1,2,...,{\tt nt-1}\}\) or equivalently on an equidistant real-frequency grid \(\{ i d \omega: i=0,1,2,...,{\tt nfreq}-1\}\) with spacing \(h\) or \(d \omega\) respectively. Note that the time argument of a Green’s function object \(G(t)\) is understood as a time difference in Keldysh steady state theory and the class ness::GF stores the positive time part. The Green’s functions are initialized by the following parameters:

dgrid (double): spacing of the time or frequency discretization.
ngrid (integer): number of discretization points on the time or frequency axis.
size1 (integer): first orbital or any quantum number dimension.
size2 (integer): second orbital or any quantum number dimension. Each Green’s function \(G(t)\) is a matrix of dimension size1 \(\times\) size2.
sig (integer): takes the values -1 (fermion) or +1 (boson).
gf_type (integer): documents the type of the Green’s function and takes 0 (time type: time_gf) or 1 (frequency type: freq_gf).

Constructors

`GF()`	Default constructor, does not allocate memory, sets all parameters to `0`.
`GF(double dgrid, int ngrid, int size1, int size2, int sign, int gf_type)`	Allocates memory and sets all entries to `0`.
`GF(double dgrid, int ngrid, int size1, int sign, int gf_type)`	Allocates memory and sets all entries to `0` with `size1=size2`.
`GF(const GF& g)`	Constructs a `GF` object as a copy of a given `GF` object `g`.

Accessing individual elements

`G.Retarded`	Return reference to retarded component.
`G.Lesser`	Return reference to lesser component.

The following member functions read and set components or submatrices of a Green’s function G from another Green’s function gf:

`G.set_element(int i1, int i2, const GF& gf, int s1, int s2)`	Set `G.Retarded[n](i1,i2) = gf.Retarded[n](s1,s2)` and `G.Lesser[n](i1,i2) = gf.Lesser[n](s1,s2)` for all `n`.
`G.get_element(int i1, int i2, const GF& gf, int s1, int s2)`	Set `gf.Retarded[n](s1,s2) = G.Retarded[n](i1,i2)` and `gf.Lesser[n](s1,s2) = G.Lesser[n](i1,i2)` for all `n`.
`G.set_submatrix(const std::vector<int> & subid1, const std::vector<int> & subid2, const GF& gf)`	Set `G.Retarded[n](subid1[i1],subid2[i2]) = gf.Retarded[n](i1,i2)` and corresponding lesser for all `i1,i2` and `n`.
`G.get_submatrix(const std::vector<int> & subid1, const std::vector<int> & subid2, const GF& gf)`	Set `gf.Retarded[n](i1,i2) = G.Retarded[n](subid1[i1],subid2[i2])` and corresponding lesser for all `i1,i2` and `n`.

Scalar multiplication and incrementation

`G.smul(cplx weight)`	`G` set to `weight` \(\cdot\) `G` for all times / frequencies.
`G.incr(GF &gf, cplx weight)`	`G` set to `G` + `weight` \(\cdot\) `gf` for all times / frequencies.

Size consistency G.ngrid_==gf.ngrid_ and G.size1_ \(\cdot\) G.size2_==gf.size1_ \(\cdot\) gf.size2_ as well as G.gf_type_==gf.gf_type_ required.

Density matrix

G.Lesser[0]

Density matrix \(\rho_{ab} = i \xi G^<_{a,b}(0)\) (\(\xi\) is the bosonic/fermionic sign).

File I/O

`void ness::outputGF(GF &G, std::string filename)`	Writes the specified Green’s function `G` into a text file `filename`.
`void ness::readGF(GF &G, std::string filename)`	Reads the contents of a text file `filename` into a Green’s function object `G`.

The output or input text file is ordered: Frequency / time, \(\Re \{G^{\mathrm{R}}\}\), \(\Im \{G^{\mathrm{R}}\}\), \(\Im \{G^{<}\}\), \(\Re \{G^{<}\}\), … (row-major order for orbital indices size1_, size2_).

`ness::fft_solver`

Overview

class

ness::fft_solver

This class contains the data structures for transforming a Green’s function \(G(t)\) on the equidistantly discretized real-time Keldysh contour \(\mathcal{C}_K\) with points \(\{i h: i=0,1,2,...,{\tt nt-1}\}\) to a Green’s function \(G(\omega)\) on an equidistant real-frequency grid \(\{ i d \omega: i=0,1,2,...,{\tt nfreq}-1\}\) with spacing \(h\) or \(d \omega\) respectively. We use the fast Fourier transform (FFT) library FFTW3. The class is defined by two parameters:

nt (integer): number of time points of the GF object.
FFTW_FLAG (size_t): The FFTW3 flag, this argument is usually either FFTW_MEASURE or FFTW_ESTIMATE.

Constructor

fft_solver(int nt, size_t FFTW_FLAG)

Default constructor, creates necessary FFTW3 structures (input and output array, FFT plans).

Member functions

`fft.to_time(GF &outG, const GF &inG, int set_dt = 0)`	Transforms the contents of the Green’s function `inG` of type `freq_gf` to time and writes them into the Green’s function `outG` of type `time_gf`. If `set_dt` is set to 1, `outG.dgrid_` is automatically set from `inG.dgrid_`. It is important that `inG.ngrid_ = fft.nfreq_` and `outG.ngrid_ = 3 fft.nfreq_ / 2 + 1`.
`fft.to_freq(GF &outG, const GF &inG, int set_dw = 0)`	Transforms the contents of the Green’s function `inG` of type `time_gf` to frequency and writes them into the Green’s function `outG` of type `freq_gf`. If `set_dw` is set to 1, `outG.dgrid_` is automatically set from `inG.dgrid_`. It is important that `inG.ngrid_ = 3 fft.nfreq_ / 2 + 1` and `outG.ngrid_ = fft.nfreq_`.

`ness::GF_pair`

Overview

class

ness::GF_pair

This class contains the data structures for simultaneously handling a time Green’s function object and frequency Green’s function object. The advantage of this class is that both objects (freq and time) are constructed with the necessary amount of data points for a Fourier transform. The Fourier transform can be called as a member function (to_time|to_freq) on the members of the class. The class is thus defined by all the parameters of a time GF object.

Constructors

`GF_pair()`	Default constructor, does not allocate memory.
`GF_pair(double h, int nt, int size, int sign, fft_solver & fft, int long_freq = 0)`	Allocates memory and sets all entries to `0` with `size1=size2=size`.
`GF_pair(const GF_pair & gp)`	Constructs a `GF_pair` object as a copy of a given `GF_pair` object `gp`.

Accessing members

`gp.freq`	Frequency `GF` object with correct spacing and number of points.
`gp.time`	Time `GF` object with correct spacing and number of points for the FFT.

Member functions

`gp.to_time()`	Fourier transforms `gp.freq` to time and writes the result into `gp.time`.
`gp.to_freq()`	Fourier transforms `gp.time` to frequency and writes the result into `gp.freq`.

Equilibrium Green’s functions 

Also in the steady state we provide tools for calculating the Green’s function from a given spectral function solving the integrals (27), (29) or to force a given Green’s function into equilibrium using relation (31). The provided functions also work in the namespace ness.

Constructing a Green’s function from a given DOS 

As defined in the manual, we can construct a Green’s function from a given density of states \(A(\omega)\) in the steady state:

\[G(t) = \int d \omega \,A(\omega) g_{\omega}(t) \ ,\]

where \(g_\omega(t)\) is the single-orbital noninteracting Green’s function with energy \(\omega\). The latter is defined as (\(F_\xi(x)=1/(e^{\beta x}-\xi )\) is the Fermi/Bose distribution):

\(g_\omega^R (t) = -i e^{-i t \omega}\)
\(g_\omega^< (t) = -i\xi e^{-i t \omega} F_\xi(\omega-\mu)\)

For the class ness2::herm_matrix_ness this is accomplished by the following call of ness2::green_equilibrium_ness:

ness2::green_equilibrium_ness(sign, G, DOS &dos, beta, mu,h, METHOD, limit=100, nn=20)

Set \(G^{<|R}(t)\) to equilibrium Green’s function for Bosons (sign=+1) or Fermions (sign=-1).

The user provides the inverse temperature beta, the time step h and the chemical potential mu (set to 0.0 by default). DOS is a class representing \(A(\omega)\), which provides an operation dos(double omega) to return \(A(\omega)\), and the numbers dos.lo_ and dos.hi_ representing the lower and upper bound of the support of \(A(\omega)\). One may use keywords BOSON and FERMION for the sign \(\pm 1\). METHOD can be FFT_TRAPEZ or FFT_ADAPTIVE, for FFT_TRAPEZ arguments nn and limit are ignored.

If the argument METHOD is FFT_TRAPEZ the integral is computed as in transform_to_time, after initializing the imaginary part of the frequency-dependent components ret_.freq_ and les_.freq_ using the dos function on the frequency grid (37). Note that the routine is implemented only for scalar \(A(\omega)\); if \(G\) is matrix-valued, it is set to a diagonal matrix.
Alternatively, if METHOD = FFT_ADAPTIVE, integrals are computed using the adaptive Fourier integral of libcntr, with a subdivision of the domain of \(A(\omega)\) in at most limit intervals of nn points (default is limit=100 and nn=20). The adaptive Fourier integration is more accurate but slower, because it does not exploit the FFT algorithm.

Example:

int N_grid = 100;
int size = 2;
double h = 0.1;
double beta= 0.1;
double mu= 0.0;
int sign= -1;
herm_matrix_ness G(N_grid, size);

cntr::bethedos dos;
dos.V_ = 0.25;
dos.lo_ = -0.5;
dos.hi_ = 0.5;

ness2::green_equilibrium_ness(-1, G, dos, beta, mu, h, FFT_TRAPEZ);

Setting a Green’s function to thermal equilibrium 

Also for this class, we provide a method ness2::herm_matrix_ness::force_equilibrium, which does a Fourier transform of \(G^R(t)\) to frequency (analogous to ness2::herm_matrix_ness::integral_transform_to_freq with METHOD = FFT_TRAPEZ) to initialize \(G^<(\omega)\) on the frequency grid (37) according to the equilibrium dissipation fluctuation relation (31), and transforms \(G^<(\omega)\) to \(G^<(t)\) analogous to ness2::herm_matrix_ness::integral_transform_to_time.

G.force_equilibrium(sign, beta, mu, h)

Set \(G^<(t)\) according to (29) and the equilibrium dissipation fluctuation relation (31) for Bosons (sign=+1) or Fermions (sign=-1), with \(A(\omega)\) from \(G^R(t)\).

Diagram utilities 

Bubble diagrams 

Apart from algebraic manipulations there is also a steady state implementation of the bubble products \(C(t)=iA_{a,a'}(t) B_{b',b}(-t)\) and \(C(t)=iA_{a,a'}(t) B_{b,b'}(t)\), see the Diagrams section in the main code. Applying the Langreth rules we can provide the two basic functions, ness2::Bubble1_ness and ness2::Bubble2_ness, as in the main code, which work also in the namespace ness.

`ness2::Bubble1_ness(herm_matrix_ness &C,int c1,int c2,herm_matrix_ness &A,int a1,int a2,herm_matrix_ness &B,int b1,int b2)`	\(C_{c1,c2}(t)=i A_{a1,a2}(t) B_{b2,b1}(-t)\) is calculated for given `herm_matrix_ness` objects `A,B` with matrix indices `a1`, `a2` and `b1`, `b2`, respectively.
`ness2::Bubble1_ness(herm_matrix_ness &C,herm_matrix_ness &A,herm_matrix_ness &B)`	\(C(t)=i A(t) B(-t)\) is calculated for given scalar `herm_matrix_ness` objects `A,B` with hermitian symmetry (indices default to zero).
`ness2::Bubble2_ness(herm_matrix_ness &C,int c1,int c2,herm_matrix_ness &A,int a1,int a2,herm_matrix_ness &B,int b1,int b2)`	\(C_{c1,c2}(t)=i A_{a1,a2}(t) B_{b1,b2}(t)\) is calculated for given `herm_matrix_ness` objects `A,B` with matrix indices `a1`, `a2` and `b1`, `b2`, respectively.
`ness2::Bubble2_ness(herm_matrix_ness &C,herm_matrix_ness &A,herm_matrix_ness &B)`	\(C(t)=i A(t) B(t)\) is calculated for given scalar `herm_matrix_ness` objects `A,B` with hermitian symmetry (indices default to zero).

Example:

The following example calculates the second-order self-energy for a general time-dependent interaction \(U(t)\) out of a scalar local Green’s function, \(\Sigma(t,t') = U(t) G(t,t') G(t',t) U(t') G(t,t')\). The diagram is factorized in a particle-hole bubble \(\chi\) and a particle-particle bubble:

\[\Sigma(t,t') = -i W(t,t') G(t,t'), \quad W(t,t')=U(t) \chi(t,t') U(t'), \quad \chi(t,t') = iG(t,t') G(t',t).\]

int N_grid = 100;
int size = 1;
double h = 0.1;
int sign= -1;
cdmatrix U(size,size);
herm_matrix_ness G(N_grid, size);
herm_matrix_ness Sigma(N_grid, size);
herm_matrix_ness W(N_grid, size);

// set G and U ...

// compute Sigma
ness2::Bubble1_ness(W,G,G); // W(t) is set to chi=ii*G(t)G(-t);
W.left_multiply(U, ness2::fft_domain::time); // set W  to W(t)=U chi(t) U:
W.right_multiply(U, ness2::fft_domain::time);
ness2::Bubble2_ness(Sigma,G,W); // Sigma(t) is set to Sigma=ii*G(t)W(t);
Sigma.smul(-1.0, ness2::fft_domain::time); // the final -1 sign

Dyson equation 

At present, for the solution of the Dyson equation we provide only one method METHOD=FFT_TRAPEZ, which is based on a straightforward discrete Fourier transform: The self energy is transformed to the frequency grid (37) using the integral transform ness2::herm_matrix_ness::integral_transform_to_freq(), corresponding to a trapezoidal evaluation of (28) and (30). Then the frequency-dependent functions \(G^{<|R}(\omega)\) are computed on the grid (37) using (32) and (33). For the back-transform to (27) and (29), we again use the implementation ness2::herm_matrix_ness::integral_transform_to_time().

Note

In order to reduce Fourier artifacts from the large frequency region we only use \(G^{<|R}(\omega)\) at frequency grid points \(-(N_{\rm freq}/2-1),...,N_{\rm freq}/2-1\), with \(N_{\rm freq}=N_{\rm ft}/3\), and set \(G^{<|R}(\omega)\) to zero outside this interval. With the stepsize \(\Delta_\omega\), the maximum frequency is therefore \(\omega_{\rm max} = \frac{2\pi}{6h}\). An accurate solution of the Dyson equation requires \(h\) to be small enough such that spectral functions are sufficiently decayed outside \([-\omega_{\rm max},\omega_{\rm max}]\). Moreover, \(N_{\rm ft}\) must be sufficiently large such that the functions \(G^{R|<}(t)\) decay at the largest time \(t_{\rm max} = \frac{hN_{\rm ft}}{2}\), and \(\Delta_\omega\) can resolve the most narrow structures in frequency.

ness2::dyson(herm_matrix_ness &G, double mu, cdmatrix &epsilon, herm_matrix_ness &Sigma, double h, fft_integral_method method=FFT_TRAPEZ, [ETA])

Solve the frequency Dyson equation with timestep h for G, for a given self-energy Sigma, local Hamiltonian epsilon, and further parameters [ETA]. For METHOD, the only currently implemented option is FFT_TRAPEZ.

G and Sigma are the Green’s functions of type ness2::herm_matrix_ness which must be defined on the same gridsize Nft and must have the same matrix dimension.
epsilon is a complex Eigen matrix.
[ETA] summarizes further optional parameters which introduce a long-time regularization of the Dyson equation. Possible options are REG_CONST, REG_GAUSS, REG_OHMIC
This routine is also implemented in the namespace ness using the function call ness::dyson_ness(GF & G, GF & Sigma, const cdmatrix & H, double eta).

Example:

int N_grid = 100;
int size = 1;
double h = 0.1;
double mu=0;
cdmatrix Ham(size,size);
herm_matrix_ness G(N_grid, size);
herm_matrix_ness Sigma(N_grid, size);

// Initialize G, Sigma, Ham, ...

ness2::dyson(G, mu, Ham, Sigma, h, FFT_TRAPEZ);

VIE2

In the old namespace ness, there is also an explicit implementation for solving a VIE2, the involved Green’s functions need to match in their size1_, size2_ and ngrid_ properties and need to be hermitian.

`ness::vie2_ness(GF & G, const GF & F, const GF & Q)`	\((1 + F) \ast G = Q\), `F`, `Q` given `GF` objects in frequency space, `G` is output `GF` in frequency space.
`ness::vie2_ness(GF_pair & G, const GF_pair & F, const GF_pair & Q)`	\((1 + F) \ast G = Q\), `F`, `Q` given `GF_pair` objects in time, `G` is output `GF_pair` in time.

Parallelization with OMP 

class

ness2::FFT_OMP_Manager

The Fourier transforms can be parallelized using a shared-memory parallelization, using the built in functionalities of the FFTW library. In order to make this functionality available, we provide a small helper class ness2::FFT_OMP_Manager. ness2::FFT_OMP_Manager::initialize must be called for initialization at the beginning of the program. Before creating an ness2::fft_array or ness2::herm_matrix_ness object (which contains an fft_plan of the FFTW library), one needs to call ness2::FFT_OMP_Manager::set_threads_for_new_plans. All plans which are created after this call will be generated such that their execution spans over a specified number of nomp threads. Note that new plans are also created once a new ness2::fft_array or ness2::herm_matrix_ness is created via an assignment and or copy assignment. At the end of the program, call ness2::FFT_OMP_Manager::finalize.

`ness2::FFT_OMP_Manager::initialize()`	Initialize parallelization helper class, call at beginning of program.
`ness2::FFT_OMP_Manager::set_threads_for_new_plans(int nomp)`	Fix amount of threads `nomp` for all plans to be created (in `ness2::fft_array` or `ness2::herm_matrix_ness`).
`ness2::FFT_OMP_Manager::finalize()`	Finalize, call at end of program.

Inside the code, you can access the number of threads currently used in the construction of plans from the variable FFT_OMP_Manager::current_threads. To test the parallel setup, we provide a notebook utils/test_ness2_omp.ipynb.

Note

This functionality is added automatically when the library is built with omp=ON and ness=ON. In addition it requires the FFTW library to be compiled with --enable-threads. If only the non-threaded FFTW library is available (or if omp=OFF), the calls to the FFT OMP Manager have simply no effect, and all Fourier transforms are single threaded.
The generation of plans (i.e., the creation of any fft_herm_matrix_ness or fft_array) is not thread-safe, and should be called from a serial region of the code.
In general, one should avoid nesting outer parallelization and parallelization of the Fourier transforms. Using outer parallelization in combination with FFT plans with more than one thread may lead to oversubscription. If outer parallelization is used, plans should therefore simply be created with nomp=1 threads.

Utilities 

Comparing Green’s functions 

To compare the data of two Green’s functions, e.g. for convergence checks, one can analyze the difference using the \(L_2\)-norm difference \(|| M ||\) for \(M=A-B\) of the individual elements. The difference is defined as the \(L_2\)-norm difference, summed over all time or frequency points:

(38)\[\Delta[A,B] = \sum_{i=0}^{\tt ngrid-1} \big( ||A^\mathrm{R}(i \cdot \Delta x)-B^\mathrm{R}(i \cdot \Delta x)|| + ||A^<(i\cdot \delta x)-B^<(i\cdot \delta x)|| \big).\]

\(|| M ||\) is the standard matrix \(L_2\) norm, \(\Delta x = \Delta \omega\) for domain=ness2::fft_domain::freq and \(\Delta x = h\) for domain=ness2::fft_domain::time. The corresponding function call is ness2::distance_norm2:

ness2::distance_norm2(herm_matrix_ness &A,herm_matrix_ness &B,domain)

Returns difference norm \(\|A-B\|_2\) (double) for ness2::herm_matrix_ness or ness2::fft_array objects, on the time or frequency grid (domain=ness2::fft_domain::[time|freq]).

In the old namespace the function call is double ness::GF2norm(GF &g1, GF &g2), the type of g1 and g2 has to match and they must have equal size1_, size2_ and ngrid_.

Example:

int N_grid = 100;
int size = 1;
double h = 0.1;
herm_matrix_ness G(N_grid, size);

// some iterative procedure to determine G
{
    herm_matrix_ness tG(N_grid, size);   //temporary variable
    double convergence_error;
    int iter_max=100;
    for(int iter=0;iter<=iter_max;iter++){
        tG = G;        // store values of G before iteration
        // ... some code to update G ...
        convergence_error = ness2::distance_norm2(tG,G,ness2::fft_domain::time);
        if( convergence_error < some_sufficiently_small_number) break;
    }
    if(iter>iter_max){
        cout << "no convergence!" << endl;
    }
}

Data exchange with two-time functions 

Data exchange between two-time and steady state functions is enabled by the following functions: The function ness2::cntr2ness allows to read the ness2::fft_domain::time data of a steady state function \(G_{\rm ness}\) from a given timeslice \(t_0\) of a two-time object \(G\), such that \(G_{\rm ness}^{R|<}(t)\leftarrow G_{\rm cntr}^{R|<}(t_0,t_0-t)\) for all \(t\). At negative times \(t\), \(\bar G\) is set assuming hermitian symmetry; if tstp is smaller than the maximum time \(N_{\rm ft}/2-1\) in \(G_{\rm ness}\), the remaining entries in \(G_{\rm ness}\) are left zero. The reverse function ness2::ness2cntr sets \(G_{\rm cntr}^{R|<}(t,t')\leftarrow G_{\rm ness}^{R|<}(t-t')\) for all arguments \((t,t')\) where \(t-t'\) is in the domain of \(G_{\rm ness}\), and zero otherwise. The two function calls are:

`ness2::cntr2ness(herm_matrix_ness& Gness, cntr::herm_matrix<T>& Gcntr, int tstp = -1)`	Set values \(G^{<\|R}\) of a `ness2::herm_matrix_ness` object `Gness` from timeslice `tstp` of a two-time function `Gcntr`.
`ness2::ness2cntr(cntr::herm_matrix<T>& Gcntr, const herm_matrix_ness& Gness)`	Set \(G^{<\|R}\) of a two-time function `Gcntr` from a steady state function `Gness`, assuming time-translational invariance.

Gcntr can be cntr::herm_matrix, cntr::herm_matrix_timestep, cntr::herm_matrix_moving or cntr::herm_matrix_timestep_moving.
For the latter three types, the argument tstp in cntr2ness is ignored and can be omitted, because always the leading physical timestep is addressed.
If the argument tstp is omitted for Gcntr of type cntr::herm_matrix, it defaults to the largest physical time (Gcntr.nt).

Resampling 

In the steady state, calculations on much finer time grids are possible than in two-time calculations. Therefore downsampling and upsampling routines are necessary to read and write Green’s functions defined on different grids effectively between the two interfaces. The following routines resample steady state ness2::herm_matrix_ness objects, so that they can be brought into the same shape as their cntr counterpart and be processed with the ness2::cntr2ness and ness2::ness2cntr functions.

Upsampling:

Upsampling a ness2::herm_matrix_ness object to one with finer time grid is equivalent to increasing the maximum frequency on frequency domain. One can thus transform the ness2::herm_matrix_ness object of length Nft to frequency using integral_transform_to_freq and enlarge the frequency interval by adding zeros until the desired length Nft*factor is achieved. Then a back-transform to time using integral_transform_to_time on the new grid Nft*factor returns the upsampled object in time.

ness2::upsample(double h_in,herm_matrix_ness &in, int factor)

Upsample a ness2::herm_matrix_ness object in with initial timestep h_in to a finer grid with a factor of factor more points and spacing h_in/factor.

Downsampling:

Downsampling a ness2::herm_matrix_ness object to one with coarser time grid is equivalent to keeping only every factor-th element in a Green’s function of length Nft in time domain to reduce it to the length Nft/factor. In frequency space this corresponds to dropping all high frequency elements outside the length Nft/factor of the resampled object. Returns ness2::herm_matrix_ness or ness2::fft_array.

ness2::downsample(herm_matrix_ness &in, int factor)

Downsample a ness2::herm_matrix_ness or ness2::fft_array object in to a coarser grid with a factor of factor fewer points. Returns ness2::herm_matrix_ness or ness2::fft_array.

HDF5 python utilities 

Also for the steady state code we provide python tools, which include scripts for reading and post-processing Green’s functions from HDF5 format via the h5py python package. The python tools can be found in libcntr/python3 in the file ReadNESS.py. An example usage is demonstrated in the example programs Steady-state Dyson equation: Anderson Impurity model and DMFT in the steady state.

Note

To make the python module available to a python script, the module needs to be in the python include path. This is achieved by setting export PYTHONPATH=<path>/nessi/libcntr/python3.