# Quantum fluctuations approach to the nonequilibrium _GW_ approximation

E. Schroedter, J. -P. Joost  $|$ , M. Bonitz  $|$ Corresponding author:
<EMAIL_ADDRESS>
(Received April 15, 2022, in final form June 03, 2022)

###### Abstract

The quantum dynamics of fermionic or bosonic many-body systems following
external excitation can be successfully studied using two-time nonequilibrium
Green’s functions (NEGF) or single-time reduced density matrix methods.
Approximations are introduced via a proper choice of the many-particle self-
energy or decoupling of the BBGKY hierarchy. These approximations are based on
Feynman’s diagram approaches or on cluster expansions into single-particle and
correlation operators. Here, we develop a different approach where, instead of
equations of motion for the many-particle NEGF (or density operators), single-
time equations for the correlation functions of fluctuations are analyzed. We
present a derivation of the first two equations of the alternative hierarchy
of fluctuations and discuss possible decoupling approximations. In particular,
we derive the polarization approximation (PA) which is shown to be equivalent
to the single-time version [following by applying the generalized
Kadanoff–Baym ansatz (GKBA)] of the nonequilibrium $GW$ approximation with
exchange effects of NEGF theory, for weak coupling. The main advantage of the
quantum fluctuations approach is that the standard ensemble average can be
replaced by a semiclassical average over different initial realizations, as
was demonstrated before by Lacroix and co-workers [see e.g. D. Lacroix et al.,
Phys. Rev. B, 2014, 90, 125112]. Here, we introduce the stochastic $GW$ (SGW)
approximation and the stochastic polarization approximation (SPA) which are
demonstrated to be equivalent to the single-time $GW$ approximation without
and with exchange, respectively, in the weak coupling limit. In addition to
the standard stochastic approach to sample initial configurations, we also
present a deterministic approach. Our numerical tests confirm that our
approach has the same favorable linear scaling with the computation time as
the recently developed G1-G2 scheme [Schluenzen et al., Phys. Rev. Lett.,
2020, 124, 076601]. At the same time, the SPA and SGW approaches scale more
favorably with the system size than the G1-G2 scheme, allowing to extend
nonequilibrium $GW$ calculations to bigger systems.

Key words: quantum dynamics, nonequilbirium Green’s functions, quantum
fluctuations, Hubbard model, GW approximation

###### Abstract

Квантову динаміку ферміонних та бозонних багаточастинкових систем під дією
зовнішнього збудження можна успішно вивчати за допомогою методів нерівноважних
функцій Ґріна (НФҐ) або приведеної матриці густини. Апроксимації вводяться
шляхом правильного вибору багаточастинкової власної енергії або розчепленням
ланцюжка рівнянь ББҐКІ. Ці наближення ґрунтуються на методі діаграм Фейнмана
або на кластерних розвиненнях в одночастинкових та кореляційних операторах. Ми
застосовуємо інший підхід, де замість рівнянь руху для багаточастинкових НФҐ
(або операторів густини) аналізуються рівняння для кореляційних функцій
флуктуацій. Ми отримуємо перші два рівняння альтернативної ієрархії флуктуацій
та обговорюємо можливі наближення для розчеплення цих рівнянь. Зокрема,
отримано поляризаційне наближення (ПН), яке є еквівалентним до нерівноважного
_GW_ наближення з обмінними ефектами теорії НФҐ у межах узагальненого
замикання Каданова-Бейма для випадку слабкого зв’язку. Основна перевага
підходу квантових флуктуацій полягає в тому, що стандартне усереднення за
ансамблем можна замінити напівкласичним середнім за різними початковими
реалізаціями, як це було показано раніше Лакруа та співавторами [див. D.
Lacroix et al., Phys. Rev. B, 2014, 90, 125112]. Ми виконуємо стохастичне _GW_
наближення (SGW) а також наближення стохастичної поляризації (НСП), які в
границі слабкого зв’язку, як показано в цій статті, є еквівалентними до _GW_
апроксимації як з врахуванням обмінних ефектів, так і без. Крім стандартного
стохастичного методу у формалізмі початкових конфігурацій ми також
представляємо точний підхід. Наші числові розрахунки підтверджують, що
запропонований метод має таке ж сприятливе лінійне масштабування за часами
обчислення, як і нещодавно розроблена схема G1–G2 [Schluenzen et al., Phys.
Rev. Lett., 2020, 124, 076601]. У той же час підходи НСП і SGW краще
масштабуються за розміром, ніж схема G1-G2, що дозволяє застосовувати
нерівноважні _GW_ розрахунки і для більших систем.

Ключов слова: квантова динаміка, нерівноважні функції Ґріна, квантові
флуктуації, модель Габбарда, _GW_ наближення

##  1 Introduction

The dynamics of quantum many-body systems following external excitation is of
great interest in many areas such as dense plasmas, nuclear matter, ultracold
atoms or correlated solids. There is a large variety of methods available to
simulate such systems which include real-time quantum Monte Carlo, density
matrix renormalization group approaches, time-dependent density functional
theory and quantum kinetic theory.

At present, only the formalisms of reduced density operators, e.g., [1] and
nonequilibrium Green’s functions (NEGF), e.g., [2, 3, 4, 5] can rigorously
describe the quantum dynamics of correlated systems in more than one
dimension, e.g., [6]. However, NEGF simulations are computationally expensive,
among other things, due to their cubic scaling with the simulation time
$N_{t}$ (the number of time steps). Only recently, linear scaling with $N_{t}$
could be achieved within the G1-G2 scheme [7, 8] which could be demonstrated
even for advanced self-energies including the $GW$ and the $T$-matrix
approximation. Even the nonequilibrium dynamically screened ladder
approximation, which selfconsistently combines dynamical screening and strong
coupling, is now feasible, at least for lattice models [9].

The advantage of time linear scaling of the G1-G2 scheme comes at a price: a
simultaneous propagation of the single-particle and correlated two-particle
Green’s functions, $G_{1}$ and $\mathcal{G}_{2}$, requires a large
computational effort for computing and storing all matrix elements of
$\mathcal{G}_{2}$. For example, the CPU time of $GW$-G1-G2 simulations scales
as $N_{b}^{6}$, where $N_{b}$ is the basis dimension. Even though this
difficulty can be relieved using massively parallel computer hardware, it is
well worth to look for alternative formulations of the problem that are more
suitable for computations, ideally without loss of accuracy.

Here, we consider such an alternative formulation of the quantum many-body
problem that is based on a stochastic approach to the (single-time) dynamics
of quantum fluctuations. The idea is that correlation effects can be
equivalently formulated in terms of correlation functions of fluctuations.
This has been studied in detail, for classical systems by Klimontovich, see
e.g., [10, 11, 12], for an overview on his work see references [13, 14]. For
quantum systems, similar concepts were developed by Ayik, Lacroix and many
others, e.g., [15, 16, 17]. Before outlining our concept, we briefly mention
how our approach differs from other stochastic methods that, of course, have a
long history in physics.

Fluctuations have been investigated in many areas of classics physics,
including plasmas [18], gases, liquids and the theory of liquid or plasma
turbulence, e.g., [19, 20]. Many-particle systems in thermodynamic equilibrium
with interactions so strong that perturbation theories fail are successfully
described using classical and quantum Monte Carlo methods, e.g., [21, 22, 23].
Stochastic approaches have been also broadly applied to treat the time
evolution of interacting classical systems, e.g., via molecular dynamics which
are averaged over initial conditions, kinetic equations including Monte Carlo
collision techniques (PIC-MCC) or kinetic Monte Carlo methods. Similarly
diverse is the spectrum of stochastic approaches to the dynamics of quantum
systems which includes real-time quantum Monte Carlo, various versions of
quantum molecular dynamics, e.g., [24, 25], path integral molecular dynamics,
e.g., [26], generalized Bohmian dynamics, e.g., [27], quantum corrections to
semiclassical dynamics [28], or the truncated Wigner approximation (TWA) [29,
30, 31] and its recently introduced improvements for spin models [32, 33, 34].
Such an extension of the TWA is presented, for example, in reference [35],
where the initial discrete Wigner distribution is sampled in phase space in
combination with a cluster expansion for the phase space points. The
realizations of the expansion coefficients are then propagated by a set of
semiclassical equations based on the BBGKY hierarchy. A different approach for
the description of correlation functions of fluctuations is given by
Wetterich, e.g., [36], where the time evolution of the generating functionals
for the correlation functions is calculated from an exact functional
differential equation instead of directly considering the correlation
functions. For nonequilibrium situations, approximations in the form of an
appropriate truncation of the generating functional are then applied in order
to find solutions.

Here, we use an equation of motion approach that is similar to the single-time
BBGKY hierarchy for reduced density operators, e.g., [1], and to the many-time
Martin–Schwinger hierarchy of NEGF theory [37]. However, instead of many-
particle Green’s functions or density operators, we consider an approach that
is based on fluctuations of the single-particle Green’s function operator on
the time diagonal, $\delta\hat{G}(t)$, and products thereof, giving rise to
single-time two-particle and higher fluctuations. We derive the first
equations for the fluctuation hierarchy and compare them to the BBGKY
hierarchy of the single-time reduced density operators to which they are fully
equivalent. Further, we discuss possible closures to the hierarchy and
concentrate on two approximations that are equivalent to the single-time $GW$
approximation of many-particle theory with and without exchange which plays a
fundamental role for ground state and nonequilibrium properties of quantum
systems.

While the fluctuation hierarchy has a similar complexity as the BBGKY
hierarchy and the G1-G2 scheme, it has a major advantage: the correlated
dynamics can be reformulated and mapped onto an effective single-particle
problem, as was shown by Klimontovich, for classical systems [10], and by Ayik
and Lacroix, for quantum systems, within their stochastic mean field
approximation [15, 17, 16]. The main result of our paper is the identification
of two equations of motion for the single-particle fluctuation $\delta\hat{G}$
— the stochastic polarization approximation (SPA) and stochastic $GW$ (SGW) —
that are equivalent to single-time version of $GW$ with and without exchange,
respectively, in the weak coupling limit. By performing extensive tests for
small Hubbard clusters, we verify this correspondence. Furthermore, we test
the efficiency of this approach by utilizing various stochastic sampling
schemes.

This paper is structured as follows. In section 2 we briefly recall the
formalism of second quantization and the equations of the G1-G2-scheme. In
section 3 we introduce our concept of fluctuations of the nonequilibrium
Green’s functions and the resulting hierarchy of equations for the
fluctuations. In section 4 we discuss the most important approximations to
decouple the hierarchy. Section 5 is devoted to the stochastic mean field
theory and section 6 to the different sampling approaches. Finally, section 7
is devoted to the application to the Hubbard model and contains our numerical
results, and a summary is presented in section 8.

##  2 Theoretical framework

Throughout this work we consider general nonequilibrium situations that are
produced by an external excitation giving rise to complex time dependencies of
all relevant quantities.

###  2.1 Notation and definitions

Although the general theory of NEGF is formulated in terms of functions
depending on multiple time arguments, here we focus on the special case of
equal times ($t=t^{\prime}$) and use a simplified notation. For example, the
single-particle correlation function will be denoted by $G^{<}_{ij}\equiv
G^{<}_{ij}(t)\equiv G^{<}_{ij}(t,t)$, where we will suppress the time
dependence.

In the following, the formalism of the second quantization is used, which is
characterized by the bosonic/fermionic creation ($\hat{c}_{i}^{\dagger}$) and
annihilation ($\hat{c}_{i}$) operators and a single-particle orbital basis
$\\{|1\rangle,\dots,|n\rangle\\}$ of the underlying single-particle Hilbert
space $\mathcal{H}$, which induces the so-called Fock space $\mathcal{F}$ and
allows for the treatment of states with a varying number of particles. Basic
properties of the mentioned operators are the following 111In this work, the
upper/lower sign, i.e., “$\pm$” and “$\mp$”, refers to bosons/fermions,
respectively.:

$\displaystyle\Big{[}\hat{c}_{i},\hat{c}_{j}^{\dagger}\Big{]}_{\mp}=\delta_{ij},\quad\Big{[}\hat{c}_{i},\hat{c}_{j}\Big{]}_{\mp}=\Big{[}\hat{c}^{\dagger}_{i},\hat{c}^{\dagger}_{j}\Big{]}_{\mp}=0.$
(2.1)

Here, a generic Hamiltonian is used to describe a quantum many-particle
system, which is given by

$\hat{H}=\sum_{ij}h_{ij}\hat{c}^{\dagger}_{i}\hat{c}_{j}+\frac{1}{2}\sum_{ijkl}w_{ijkl}\hat{c}_{i}^{\dagger}\hat{c}^{\dagger}_{j}\hat{c}_{l}\hat{c}_{k},$
(2.2)

where the single-particle contribution is defined by $h_{ij}\coloneqq\langle
i|\hat{T}+\hat{V}|j\rangle$. The operator $\hat{T}$ describes the kinetic
energy and $\hat{V}$ an external potential that may be time dependent. For the
second term we defined $w_{ijkl}\coloneqq\langle ij|\hat{W}|kl\rangle$, where
$\hat{W}$ describes two-particle interactions. Note that $\hat{W}$ may also be
time-dependent, e.g., when we compute a correlated initial state from an
uncorrelated state, via the adiabatic switching procedure [38, 9].

The starting point for NEGF theory is the one-body Green’s function that
depends on two time arguments situated on the Keldysh contour $\mathcal{C}$
which is defined as

$G_{ij}(z,z^{\prime})=\frac{1}{\mathrm{i}\hbar}\left\langle\mathcal{T}_{\mathcal{C}}\left\\{\hat{c}_{i}(z)\hat{c}^{\dagger}_{j}(z^{\prime})\right\\}\right\rangle,$
(2.3)

where $\mathcal{T}_{\mathcal{C}}$ is the time-ordering operator on the
contour, and averaging is performed with the correlated unperturbed density
operator of the system. Details of the formalism can be found in textbooks,
e.g., [5, 4] but will not be needed here. For us, it will be sufficient to
consider the correlation functions $G^{\gtrless}$ for real time arguments on
the time-diagonal $(t=t^{\prime})$. The definition of these functions is again
given by an ensemble average,

$\displaystyle G^{<}_{ij}$
$\displaystyle=\pm\frac{1}{\mathrm{i}\hbar}\left\langle\hat{c}^{\dagger}_{j}\hat{c}_{i}\right\rangle,$
$\displaystyle G^{>}_{ij}$
$\displaystyle=\frac{1}{\mathrm{i}\hbar}\left\langle\hat{c}_{i}\hat{c}^{\dagger}_{j}\right\rangle,$
(2.4) $\displaystyle\hat{G}^{<}_{ij}$
$\displaystyle=\pm\frac{1}{\mathrm{i}\hbar}\hat{c}^{\dagger}_{j}\hat{c}_{i},$
$\displaystyle\hat{G}^{>}_{ij}$
$\displaystyle=\frac{1}{\mathrm{i}\hbar}\hat{c}_{i}\hat{c}^{\dagger}_{j}.$
(2.5)

On the real-time-diagonal, the lesser component is, therefore, proportional to
the single-particle density matrix,
$n_{ij}\coloneqq\langle\hat{c}^{\dagger}_{j}\hat{c}_{i}\rangle=\pm\mathrm{i}\hbar
G^{<}_{ij}$. Note that, in the second line, equation (2.5), we omit the
ensemble averaging. The corresponding Green’s function operators will be
denoted ${\hat{G}}^{\gtrless}$ and are the starting point of our theory. Here,
we will not consider bosons in the condensate, but an extension of the
approach is straightforward and would involve also anomalous correlators.

Analogously to the single-particle NEGF, higher order Green’s functions can be
defined, but since only the limit of equal times is of interest in this work,
only one special case of the two-particle Green’s function on the time-
diagonal will be considered here,

$G^{(2),<}_{ijkl}\coloneqq-\frac{1}{\hbar^{2}}\left\langle\hat{c}^{\dagger}_{k}\hat{c}^{\dagger}_{l}\hat{c}_{j}\hat{c}_{i}\right\rangle.$
(2.6)

This quantity deviates from the standard two-particle density matrix by an
additional prefactor of $(-1/\hbar^{2})$. Furthermore, we consider a
decomposition of the two-particle Green’s function on the time-diagonal into a
mean-field (Hartree–Fock) and a correlation contribution,

$G^{(2),<}_{ijkl}=G^{<}_{ik}G^{<}_{jl}\pm
G^{<}_{il}G^{<}_{jk}+\mathcal{G}_{ijkl},$ (2.7)

where the first two terms are the Hartree and Fock contributions, whereas
$\mathcal{G}$ describes two-particle correlations. For reduced density
matrices, this decomposition is called cluster expansion [39]. The correlation
contribution is the central quantity of the so-called G1-G2 scheme, as
introduced in [7, 8], and obeys the following (pair-)exchange symmetries,

$\displaystyle\mathcal{G}_{ijkl}$
$\displaystyle=\mathcal{G}_{jilk}=\Big{[}\mathcal{G}_{klij}\Big{]}^{*}=\pm\mathcal{G}_{jikl}=\pm\mathcal{G}_{ijlk}.$
(2.8)

###  2.2 Time-diagonal Keldysh–Kadanoff–Baym equation

The equations of motion (EOMs) of the NEGF are the Keldysh–Kadanoff–Baym
equations (KBE), which contain the single-particle as well as the two-particle
NEGF, where the latter is usually approximated by various self-energy
approximations [4, 5, 40]. As was shown in references [7, 8], the EOM for the
lesser component of the NEGF on the time-diagonal is given by the time-
diagonal KBE and can be written as 222For two single-particle quantities $A$
and $B$, we define the the commutator as
$[A,B]_{ij}\coloneqq\sum_{k}\\{A_{ik}B_{kj}-B_{ik}A_{kj}\\}$.

$\mathrm{i}\hbar\partial_{t}G^{<}_{ij}=\Big{[}h^{\mathrm{HF}},G^{<}\Big{]}_{ij}+\left[{I}^{(\mathcal{G})}+{I}^{(\mathcal{G})\dagger}\right]_{ij},$
(2.9)

where the Hartree–Fock contribution of the two-particle Green’s function is
included in an effective single-particle Hartree–Fock Hamiltonian, which is
defined as

$h^{\mathrm{HF}}_{ij}\coloneqq
h_{ij}\pm\mathrm{i}\hbar\sum_{kl}w^{\pm}_{ikjl}G_{lk}^{<},$ (2.10)

where we introduced an (anti-)symmetrized interaction tensor

$w^{\pm}_{ijkl}\coloneqq w_{ijkl}\pm w_{ijlk}.$ (2.11)

The last term in equation (2.9) corresponds to the collision integral and
couples to the correlation part of the two-particle Green’s function. It is
given by

${I}^{(\mathcal{G})}_{ij}\coloneqq\pm\mathrm{i}\hbar\sum_{klp}w_{iklp}\mathcal{G}_{lpjk}.$
(2.12)

Equation (2.9) is the first equation of the G1-G2 scheme and requires the
solution of an EOM for $\mathcal{G}$, which can be derived from the known
self-energy approximations and applying the GKBA with Hartree–Fock propagators
[8, 9]. An alternative derivation can be given starting from the single-time
BBGKY hierarchy [39]. A detailed comparison of both apporaches was given in
reference [9]. Here, we reproduce the system of single-time equations that was
derived in reference [9] and explain the terms in the equations as well as the
approximations applied in the present paper,

$\displaystyle\mathrm{i}\hbar\partial_{t}\mathcal{G}_{ijkl}-$
$\displaystyle[h^{(2),\rm
HF},\mathcal{G}]_{ijkl}=\Psi^{\pm}_{ijkl}+L_{ijkl}+P^{\pm}_{ijkl}+C^{(3)}_{ijkl},$
(2.13) $\displaystyle h^{(2),\rm HF}_{ijkl}=$ $\displaystyle h_{ik}^{\rm
HF}\delta_{jl}+h_{il}^{\rm HF}\delta_{ik},$ (2.14)
$\displaystyle\Psi^{\pm}_{ijkl}\coloneqq$
$\displaystyle(\mathrm{i}\hbar)^{2}\sum_{pqrs}w^{\pm}_{pqrs}\Big{\\{}G^{>}_{ip}G^{>}_{jq}G^{<}_{rk}G^{<}_{sl}-(>\leftrightarrow<)\Big{\\}},$
(2.15) $\displaystyle L_{ijkl}=$
$\displaystyle\sum_{rs}\left\\{{V}_{ijrs}\,\mathcal{G}_{rskl}-\mathcal{G}_{ijrs}\,{V}^{*}_{rskl}\right\\},$
(2.16) $\displaystyle P^{\pm}_{ijkl}=$
$\displaystyle\Pi^{\pm}_{ijkl}\pm\Pi^{\pm}_{ijlk},$ (2.17)
$\displaystyle\Pi^{\pm}_{ijkl}\coloneqq$
$\displaystyle\pm\mathrm{i}\hbar\sum_{pqr}\left\\{\mathcal{G}_{rjpl}(w^{\pm}_{ipqr}G^{<}_{qk}-w^{\pm}_{qpkr}G^{<}_{iq})+\mathcal{G}_{iqkp}(w^{\pm}_{pjqr}G^{<}_{rl}-w^{\pm}_{prql}G^{<}_{jr})\right\\},$
(2.18) $\displaystyle C^{(3)}_{ijkl}\coloneqq$
$\displaystyle\pm\mathrm{i}\hbar\sum_{pqr}\left\\{w_{ipqr}\mathcal{G}^{(3)}_{rqjpkl}+w_{pjqr}\mathcal{G}^{(3)}_{iqrkpl}-w_{pqkr}\mathcal{G}^{(3)}_{irjpql}-w_{pqrl}\mathcal{G}^{(3)}_{ijrkqp}\right\\},$
(2.19) $\displaystyle V_{ijkl}=$
$\displaystyle\sum_{rs}(\delta_{ir}\delta_{jl}\pm n_{ir}\delta_{js}\pm
n_{js}\delta_{ir})w_{rskl}.$ (2.20)

Here, $h^{(2),\rm HF}$, equation (2.14), denotes the two-particle Hartree–Fock
Hamiltonian, $\Psi^{\pm}$ is an inhomogeneity (it does not depend on
$\mathcal{G}$) that gives rise to direct and exchange second order Born
scattering terms. The term $L$, equation (2.16), contains the ladder diagrams
that are relevant for strong coupling effects (particle-particle $T$-matrix
approximation, TPP), whereas $\Pi^{\pm}$, equation (2.18), comprises all
polarization diagrams that describe the dynamical screening effects. Note that
the term $P^{\pm}_{ijkl}$, equation (2.17), in addition to $\Pi^{\pm}_{ijkl}$,
contains a term $\Pi^{\pm}_{ijlk}$ which represents particle-hole $T$-matrix
(TPH) diagrams [9]. Finally, $C^{(3)}$ couples to the rest of the system via
three-particle correlations, $\mathcal{G}^{(3)}$.

In what follows, a central role is played by the $GW$ approximation of many-
particle physics. Its single-time version is recovered in the G1-G2-scheme,
equation (2.13), if the three–particle correlations, $C^{(3)}$, the ladder
term $L$ and the particle-hole $T$-matrix contribution, $\Pi^{\pm}_{ijlk}$,
are neglected. The remaining polarization term $\Pi^{\pm}_{ijkl}$ then
corresponds to $GW$ with exchange corrections which will be denoted
$GW^{\pm}$. On the other hand, standard $GW$ (without exchange diagrams)
follows by replacing $\Psi^{\pm}\to\Psi$ and $\Pi^{\pm}_{ijkl}\to\Pi$, which
means that, in these terms, the replacement $w^{\pm}\to w$ is made [39, 8, 9].
In the present paper we present stochastic versions of $GW$ and $GW^{\pm}$.

##  3 Quantum fluctuations of the Green’s functions

Instead of considering the correlations via the correlated part of the two-
particle Green’s function, $\mathcal{G}$, we now take a different approach
that is analogous to the classical fluctuation approach that was developed by
Klimontovich [12]. He considered fluctuations of the microscopic phase space
density [12], $\delta N(\textbf{r},\textbf{p},t)\equiv
N(\textbf{r},\textbf{p},t)-nf(\textbf{r},\textbf{p},t)$, where $n$ is the
density and $f$ is the single-particle probability density in phase space that
obeys a standard kinetic equation. Here, we generalize this approach to
quantum many-particle systems, fully including quantum effects and spin
statistics.

The natural starting point of our quantum fluctuation approach is then to
replace $nf\to G^{<}$ and $N\to\hat{G}^{<}$. Consequently, the defining
quantity is the single-particle fluctuation operator,

$\delta\hat{G}_{ij}\coloneqq\hat{G}^{<}_{ij}-G^{<}_{ij},$ (3.1)

where, obviously, $G^{<}_{ij}=\langle\hat{G}^{<}_{ij}\rangle$ and
$\langle\delta\hat{G}_{ij}\rangle=0$. The attractive feature of this approach
is that two- and three-particle fluctuation operators emerge simply as
products of single-particle fluctuation operators. In what follows, we
consider in detail the two-particle and three-particle fluctuations which are
expectation values of the operator products,

$\displaystyle\gamma_{ijkl}$
$\displaystyle\coloneqq\left\langle\delta\hat{G}_{ik}\delta\hat{G}_{jl}\right\rangle,$
(3.2) $\displaystyle\Gamma_{ijklmn}$
$\displaystyle\coloneqq\left\langle\delta\hat{G}_{il}\delta\hat{G}_{jm}\delta\hat{G}_{kn}\right\rangle.$
(3.3)

Before proceeding with an analysis of the properties and equations of motion
of these quantities, we note that analogous quantities were analyzed in the
theory of reduced density matrices where they are called ‘‘correlation
matrices’’ [41]. Therein, an operator $\hat{\mathcal{Q}}$ is defined as the
projector on the subspace orthogonal to the $N$-particle state $|\Psi\rangle$,

$\hat{\mathcal{Q}}\coloneqq\sum_{\Psi^{\prime}\neq\Psi}|\Psi^{\prime}\rangle\langle\Psi^{\prime}|\coloneqq\hat{\mathrm{I}}-|\Psi\rangle\langle\Psi|,$
(3.4)

where $\hat{\mathrm{I}}$ denotes the identity operator. It is straightforward
to express our two-particle fluctuations in terms of the projector
$\hat{\mathcal{Q}}$:

$\gamma_{ijkl}=\left\langle\hat{G}_{ik}\hat{\mathcal{Q}}\hat{G}_{jl}\right\rangle.$
(3.5)

This alternative expression shows that fluctuations describe correlation
effects of $N$ particles which undergo virtual excitations and de-excitations,
in order to avoid each other. The main drawback of this approach is, however,
that it is incapable of describing single-particle fluctuations, as defined in
equation (3.1). Therefore, we proceed with our approach that is based on the
fluctuation operator (3.1).

###  3.1 Properties of the single-particle and two-particle fluctuations

Since only the limit of equal times is considered here, there exists only a
single one-particle fluctuation operator [cf. equations (2.1), (2.4) and
(2.5)],

$\displaystyle\delta\hat{G}_{ij}=\hat{G}^{>}_{ij}-G^{>}_{ij}=\hat{G}^{<}_{ij}-G^{<}_{ij}.$
(3.6)

Thus, fluctuations of the greater and lesser components of the single-particle
Green’s function coincide on the time-diagonal. This means that, in the case
of fermions, the sum of fluctuations of particles and holes vanishes.

Using the properties of the creation and annihilation operators defined in
equation (2.1), the two-particle Green’s function can be expressed in terms of
two-particle fluctuations [41],

$G^{(2),<}_{ijkl}=G^{<}_{ik}G^{<}_{jl}\mp\frac{1}{\mathrm{i}\hbar}\delta_{il}G^{<}_{jk}+\gamma_{ijkl}.$
(3.7)

Now, using relation (2.7), we establish the relation between two-particle
correlations and fluctuations:

$\displaystyle\gamma_{ijkl}$ $\displaystyle=\mathcal{G}_{ijkl}\pm
G^{>}_{il}G^{<}_{jk}$ $\displaystyle=\mathcal{G}_{ijkl}\pm\gamma^{s}_{ijkl}.$
(3.8)

This shows what the Kronecker delta term in equation (3.7) physically means:
it contains the exchange contribution and additionally a term describing
particle-hole polarization effects [41]. Analogously to the classical case
[12], we define the difference of correlations and fluctuations as ‘‘source
term’’,

$\gamma_{ijkl}^{s}\coloneqq G^{>}_{il}G^{<}_{jk}.$ (3.9)

This source term can be interpreted as simultaneous excitation of a particle-
hole pair: a particle is annihilated in orbital $k$ and created in orbital $j$
whereas a hole is annihilated in orbital $i$ and created in $l$. Such
‘‘particle-hole fluctuations’’ can be understood as the basis for all
fluctuations. They are always present in a quantum system, even if the
particles are non-interacting 333They arise from vacuum fluctuations. In a
relativistic system, the same source term would describe particle-antiparticle
pair excitations.. In the latter case, $\mathcal{G}=0$ and the total
fluctuation $\gamma$ coincides with $\gamma^{s}$. Of course, in the presence
of correlations, the latter will affect the dynamics of the source term. The
standard example for e-h pairs is the electron gas at low temperature.
Therein, spontaneous interaction with a photon of energy $\omega$ and wave
number $q$ may excite an electron from a state $p$ (below the Fermi energy) to
a state $p+q$ above it. The probability of this process is proportional to
$n_{p}(1-n_{p+q})$, ($n_{p}$ is the momentum distribution), in agreement with
the definition (3.9).

By construction, the two-particle correlation function $\mathcal{G}$ obeys the
exchange symmetries (2.8). However, fluctuations in general do not obey these
symmetries. By using the relation between correlations and fluctuations,
equation (3.8), it is straighforward to establish the symmetries of the two-
particle fluctuations under the exchange of particles444Note that, in the
classical case, the fluctuations obey the exchange symmetry.,

$\displaystyle\gamma_{ijkl}\mp\gamma_{ijkl}^{s}$
$\displaystyle=\gamma_{jilk}\mp\gamma^{s}_{jilk}=\pm\gamma_{jikl}-\gamma^{s}_{jikl}=\pm\gamma_{ijlk}-\gamma^{s}_{ijlk}.$
(3.10)

Lastly, we find that two-particle fluctuations obey the following exchange
symmetry,

$\gamma_{ijkl}=\Big{[}\gamma_{lkji}\Big{]}^{*}.$ (3.11)

###  3.2 Single-particle dynamics

Using equation (3.8), it is straightforward to express the time-diagonal KBE
[cf. equation (2.9)] in terms of two-particle fluctuations,

$\mathrm{i}\hbar\partial_{t}G^{<}_{ij}=\Big{[}h^{\mathrm{H}},G^{<}\Big{]}_{ij}+\Big{[}I^{(\gamma)}+I^{(\gamma)\dagger}\Big{]}_{ij},$
(3.12)

where we introduced an effective single-particle Hartree Hamiltonian,

$h^{\mathrm{H}}_{ij}\coloneqq
h_{ij}\pm\mathrm{i}\hbar\sum_{kl}w_{ikjl}G^{<}_{lk},$ (3.13)

and a collision term,

$I^{(\gamma)}_{ij}\coloneqq\pm\mathrm{i}\hbar\sum_{klp}w_{iklp}\gamma_{plkj}.$
(3.14)

Note that equations (2.9) and (3.12) differ with respect to the effective
single-particle Hamiltonians, as defined in equations (2.10) and (3.13). This
is a feature of fluctuations, which, unlike Green’s functions or density
matrices, do not obey symmetries under exchange of particles. Likewise,
effective interactions of fluctuations do not contain any exchange
contributions (no Fock term). As a consequence, quantum fluctuations behave
similar to classical fluctuations [12]. Finally, the collision terms given by
equations (2.12) and (3.14), which lead to the coupling to higher order
contributions, have the same structure but involve different two-particle
quantities, i.e., $\mathcal{G}$ and $\gamma$, respectively.

To close equation (3.12) we now need an equation for $\gamma$. To derive this
equation, we start with deriving an EOM of the single-particle fluctuation
operator. This can be done by subtracting the EOMs for $\hat{G}^{<}$ and
$G^{<}$ from each other, where the EOM for $\hat{G}^{<}$ is given by

$\displaystyle\mathrm{i}\hbar\partial_{t}\hat{G}^{<}_{ij}=\left[\hat{h}^{\mathrm{H}},\hat{G}^{<}\right]_{ij},$
(3.15) $\displaystyle\hat{h}^{\mathrm{H}}_{ij}\coloneqq
h_{ij}\pm\mathrm{i}\hbar\sum_{kl}w_{ikjl}\hat{G}^{<}_{lk}.$ (3.16)

Subtracting equations (3.15) and (3.12) and replacing all operators by their
expectation values and fluctuations, as in equation (3.1), leads to the
following EOM for the single-particle fluctuation operator:

$\displaystyle\mathrm{i}\hbar\partial_{t}\delta\hat{G}_{ij}=$
$\displaystyle\mathrm{i}\hbar\partial_{t}\left\\{\hat{G}^{<}_{ij}-G^{<}_{ij}\right\\}=\left[h^{\mathrm{H}},\delta\hat{G}\right]_{ij}+\left[\delta\hat{U}^{\mathrm{H}},G^{<}\right]_{ij}+\left[\delta\hat{U}^{\mathrm{H}},\delta\hat{G}\right]_{ij}-\Big{[}I^{(\gamma)}+I^{(\gamma)\dagger}\Big{]}_{ij},$
(3.17)

where we defined the operator of an effective single-particle Hartree
potential induced by fluctuations,

$\delta\hat{U}^{\mathrm{H}}_{ij}\coloneqq\pm\mathrm{i}\hbar\sum_{kl}w_{ikjl}\delta\hat{G}^{<}_{lk}.$
(3.18)

Therefore, the corresponding commutator represents the interaction of $G^{<}$
and a ‘‘fluctuation mean-field’’. Furthermore, the last two terms in equation
(3.17) correspond to fluctuations of fluctuations (second-order fluctuations),
which can be seen by identically rewriting the collision term, equation
(3.14), as

$\Big{[}I^{(\gamma)}+I^{(\gamma)\dagger}\Big{]}_{ij}=\left\langle\left[\delta\hat{U}^{\mathrm{H}},\delta\hat{G}\right]\right\rangle_{ij}.$
(3.19)

Thus, the last two terms in equation (3.17) can be expressed as

$\displaystyle\delta\left(\left[\delta\hat{U}^{\mathrm{H}},\delta\hat{G}\right]_{ij}\right)=\pm\mathrm{i}\hbar\sum_{klp}\left\\{w_{iklp}\delta\hat{\gamma}_{plkj}-w_{kljp}\delta\hat{\gamma}_{ipkl}\right\\}.$
(3.20)

While expressing these terms as second-order fluctuations is helpful for
physical interpretation, this form turns out to be mathematically not
advantageous. Instead of an EOM for $\delta\gamma$, we consider equations for
$\gamma$ and $\delta\hat{G}$. Remarkably, by using just the EOM of
$\delta\hat{G}$ [cf. equation (3.17)], it is possible to derive the EOM of any
$N$-particle fluctuation, which makes this equation the cornerstone of our
formalism.

###  3.3 Dynamics of two-particle fluctuations $\gamma$

To close equation (3.12), we need an EOM for $\gamma$. It follows directly
from equation (3.17) which is multiplyied by another factor $\delta\hat{G}$
and by a subsequent ensemble average555For two two-particle quantities $A$ and
$B$, we define the commutator as
$[A,B]_{ijkl}\coloneqq\sum_{pq}\left\\{A_{ijpq}B_{pqkl}-B_{ijpq}A_{pqkl}\right\\}$.:

$\displaystyle\mathrm{i}\hbar\partial_{t}\gamma_{ijkl}$
$\displaystyle=\left[h^{(2),\rm{H}},\gamma\right]_{ijkl}+\pi_{ijkl}+C^{(\Gamma)}_{ijkl},$
(3.21) $\displaystyle h^{(2),\mathrm{H}}_{ijkl}$ $\displaystyle\coloneqq
h^{\mathrm{H}}_{ik}\delta_{jl}+h^{\mathrm{H}}_{jl}\delta_{ik},$ (3.22)

where we introduced an effective two-particle Hartree Hamiltonian. The second
term on the r.h.s. in equation (3.21), is a polarization contribution
involving fluctuations,

$\displaystyle\pi_{ijkl}\coloneqq$
$\displaystyle\pm\mathrm{i}\hbar\sum_{pqr}\gamma_{rjpl}\left\\{w_{ipqr}G^{<}_{qk}-w_{qpkr}G^{<}_{iq}\right\\}$
$\displaystyle\pm\mathrm{i}\hbar\sum_{pqr}\gamma_{iqkp}\left\\{w_{pjqr}G^{<}_{rl}-w_{prql}G^{<}_{jr}\right\\},$
(3.23)

which is similar to the polarization term $\Pi^{\pm}$ of the G1-G2-scheme, cf.
equation (2.17). However, in contrast to the latter, equation (3.23) does not
contain any exchange contributions (it contains the pair potential $w$ instead
of $w^{\pm}$) and also misses the particle-hole ladder term. Using the
property (3.10), it is possible to rewrite $\pi$ such that $G^{<}$ is
eliminated:

$\displaystyle\pi_{ijkl}=$
$\displaystyle(\mathrm{i}\hbar)^{2}\sum_{pqrs}w_{pqrs}\left\\{\gamma_{iskq}\gamma_{jrlp}-\gamma_{rjpl}\gamma_{siqk}\right\\}.$
(3.24)

Finally, the last term on the r.h.s. in equation (3.21), couples to three-
particle fluctuations,

$\displaystyle C^{(\Gamma)}_{ijkl}\coloneqq$
$\displaystyle\pm\mathrm{i}\hbar\sum_{pqr}\left\\{w_{ipqr}\Gamma_{rqjpkl}+w_{pjqr}\Gamma_{iqrkpl}\right\\}$
$\displaystyle\mp\mathrm{i}\hbar\sum_{pqr}\left\\{w_{pqkr}\Gamma_{irjpql}+w_{pqrl}\Gamma_{ijrkqp}\right\\},$
(3.25)

which is again similar to the corresponding term $C^{(3)}$, equation (2.19) in
equation (2.13) that involves the three-particle correlation function. Thus,
similar to the equations for the reduced density operators which form the
BBGKY hierarchy, we again obtain a hierarchy of equations for the
fluctuations. We underline that the system of equations (3.12) and (3.21) is
still exact provided that the three-particle fluctuation $\Gamma$ is known.
Interestingly, equation (3.21) does not explicitly contain exchange, second-
order Born and $T$-matrix contributions giving rise to a simpler formal
structure. Those missing terms are, of course, ‘‘hidden’’ in the three-
particle fluctuations $\Gamma$. The relation between $\gamma$ and
$\mathcal{G}$ is further illustrated in section 4 and requires further
investigation of the three-particle contributions.

##  4 Approximations

Similar to the case of the single-particle BBGKY hierarchy, we now explore
possible approximations to decouple the fluctuation hierarchy. We also analyze
how the approximations known from the BBGKY formalism, e.g., [9], can be
recovered in the present scheme and what new approximations might be offered
by the fluctuations approach. For the latter, it is interesting to compare to
the classical case.

###  4.1 Approximation of first moments

The simplest approximations that were introduced in the classical case are the
so-called ‘‘approximations of moments’’ [12], where either two- or three-
particle fluctuations are neglected. In the case $\gamma\rightarrow 0$, i.e.,
the approximation of first moments, we find a closed equation for $G^{<}$,
which is a mean-field (Hartree or nonlinear Vlasov) equation,

$\mathrm{i}\hbar\partial_{t}G_{ij}^{<}=\Big{[}h^{\mathrm{H}},G^{<}\Big{]}_{ij}.$
(4.1)

Since equation (4.1) contains no Fock term, this approximation is only
reasonable if the exchange effects are of minor importance such as for weakly
degenerate systems.

###  4.2 Time-dependent Hartree–Fock approximation

The missing Fock term in equation (4.1) is restored if the collision integral
is not neglected, as in the approximation of first moments, but treated in the
following approximation: we assume that two-particle correlations,
$\mathcal{G}$, are negligible, compared to the source contributions, i.e.,
$\gamma_{ijkl}\approx\pm\gamma_{ijkl}^{s}$. Inserting this expression into
equation (3.14) leads to the well-known Hartree–Fock approximation (TDHF),

$\mathrm{i}\hbar\partial_{t}G^{<}_{ij}=\Big{[}h^{\mathrm{HF}},G^{<}\Big{]}_{ij}.$
(4.2)

This further shows that fluctuations contain a variety of additional effects
that are not related to correlations and, thus, fluctuations are generally not
negligible, even in the case of weak coupling. Further note that this is
different from the classical case. Therein, the source term $\gamma^{s}$ does
not contribute to the collision integral at all. This difference is not
surprising because exchange (Fock term) is a pure quantum effect.

###  4.3 Approximation of second moments (2M)

The approximation of second moments follows with setting $\Gamma\rightarrow
0$, with the result

$\displaystyle\mathrm{i}\hbar\partial_{t}\gamma^{\mathrm{2M}}_{ijkl}=\left[h^{(2),\rm{H}},\gamma^{\mathrm{2M}}\right]_{ijkl}+\pi_{ijkl}.$
(4.3)

In order to investigate this approximation further, we need to reconsider the
relation between two-particle correlations and fluctuations and find, for the
polarization term [cf. equation (3.23)],

$\pi_{ijkl}=\Pi_{ijkl}+\Psi_{ijkl}.$ (4.4)

Here, $\Pi$ denotes the polarization contribution (2.18), whereas $\Psi$ is
the inhomogeneity, equation (2.15), in the EOM of $\mathcal{G}$. Note that the
superscript ‘‘$\pm$’’ is missing here which indicates that no
(anti-)symmetrized potential appears, i.e., $w^{\pm}\to w$. While still
accounting for Pauli blocking, exchange effects are not fully incorporated,
similar to the case of the approximation of first moments.

The appearance of the $\Psi$ and $\Pi$ terms in equation (4.3) indicates that
the 2M-approximation already contains parts of the standard $GW$
contributions, as given within the G1-G2 scheme (for the correspondence to the
$GW$ approximation, see the end of section 2.2). To identify the missing
contributions, we consider the source term and its dynamics which are derived
from equation (3.12) by using $\partial_{t}G^{<}=\partial_{t}G^{>}$,

$\displaystyle\mathrm{i}\hbar\partial_{t}\gamma_{ijkl}^{{\rm
2M},s}=\left[h^{(2),\rm{H}},\gamma^{{\rm
2M},s}\right]_{ijkl}+R^{(\gamma)}_{ijkl},$ (4.5) $\displaystyle
R^{(\gamma)}_{ijkl}\coloneqq
G^{>}_{il}\Big{[}I^{(\gamma)}+I^{(\gamma)\dagger}\Big{]}_{jk}+\Big{[}I^{(\gamma)}+I^{(\gamma)\dagger}\Big{]}_{il}G^{<}_{jk},\quad$
(4.6)

where $R^{(\gamma)}$ is a residual term. Alternatively, equation (4.5) can be
rewritten with the Hartree–Fock Hamiltonian in the commutator by using two-
particle correlations instead of fluctuations, taking into account equation
(2.9):

$\mathrm{i}\hbar\partial_{t}\gamma_{ijkl}^{{\rm
2M},s}=\left[h^{(2),\rm{HF}},\gamma^{{\rm
2M},s}\right]_{ijkl}+{R}^{(\mathcal{G})}_{ijkl},$ (4.7)

where the residual term $R^{(\mathcal{G})}$ follows from $R^{(\gamma)}$ by
replacing $\gamma\to\mathcal{G}$, in the collision integrals. It is important
to note that this equation accounts for the otherwise missing Fock term.

The final step for the comparison of the 2M approximation, equation (4.3),
with the $GW$ approximation, is to add the equations for $\gamma^{\rm 2M}$ and
$\gamma^{{\rm 2M},s}$, taking into account relation (3.8) and denoting
$\mathcal{G}^{\rm 2M}=\gamma^{\rm 2M}\mp\gamma^{{\rm 2M},s}$,

$\displaystyle\mathrm{i}\hbar\partial_{t}\mathcal{G}^{\mathrm{2M}}_{ijkl}=$
$\displaystyle\left[h^{(2),\mathrm{H}},\mathcal{G}^{\mathrm{2M}}\right]_{ijkl}+\Psi_{ijkl}+\Pi_{ijkl}\mp{R}^{(\mathcal{G})}_{ijkl}\mp\left[U^{(2),\mathrm{F}},\gamma^{{\rm
2M},s}\right]_{ijkl},$ (4.8)

where we defined $U^{(2),\mathrm{F}}\coloneqq
h^{(2),\mathrm{HF}}-h^{(2),\mathrm{H}}$. We conclude that the EOM (4.3) for
$\gamma^{\mathrm{2M}}$ is missing exchange contributions to the effective two-
particle Hartree–Fock Hamiltonian and terms that compensate for
${R}^{(\mathcal{G})}$ in order to be fully equivalent to the $GW$
approximation.

Whether the 2M approximation has an area of application by itself is not known
at the moment and requires further analysis. However, in what follows we
concentrate on the $GW$ approximation. In the next paragraph, we identify the
terms that are missing so far, to reproduce the $GW$ and the $GW^{\pm}$
approximations.

###  4.4 Quantum polarization approximation

Next, we consider the quantum analogue of the classical polarization
approximation (PA) which is known to be equivalent to the Balescu–Lenard
equation [12]. This, in turn, is the classical limit of the nonequilibrium
$GW$ approximation [39].

####  4.4.1 Approximation for the three-particle fluctuations

In the classical case, the PA is derived by assuming that three-particle
correlations are negligible, and two-particle correlations are small compared
to a product of two single-particle distributions. This motivates us to use
the following three weak coupling approximations:

i)

$|\mathcal{G}^{(3)}|\rightarrow 0$,

ii)

$|\mathcal{G}|\ll|\gamma|$,

iii)

$|\mathcal{G}|\ll|\gamma^{s}|$.

The task is now to identify the relevant approximation for the three-particle
term $C^{(\Gamma)}$ in the equation for $\gamma$, since this term is the
source of the terms we are missing so far.

We start by considering the equation of motion of the two-particle source
term, equation (4.7). Applying condition iii) leads to a neglect of the
residual term $R^{(\mathcal{G})}$,

$\mathrm{i}\hbar\partial_{t}\gamma_{ijkl}^{\mathrm{PA},s}=\left[h^{(2),\rm{HF}},\gamma^{\mathrm{PA},s}\right]_{ijkl}.$
(4.9)

Next, we return to the full equation for the two-particle fluctuation,
equation (3.21), and establish a relation between three-particle correlations,
$\mathcal{G}^{(3)}$, and fluctuations, $\Gamma$, which is analogous to
equation (3.8):

$\displaystyle\mathcal{G}^{(3)}_{ijklmn}$
$\displaystyle=\Gamma_{ijklmn}-G^{>}_{im}G^{>}_{jn}G^{<}_{kl}-G^{>}_{in}G^{<}_{jl}G^{<}_{km}$
$\displaystyle\mp G^{<}_{jl}\mathcal{G}_{ikmn}\mp
G^{<}_{km}\mathcal{G}_{ijln}-G^{<}_{kl}\mathcal{G}_{ijmn}$ $\displaystyle\mp
G^{>}_{im}\mathcal{G}_{jkln}\mp
G^{>}_{jn}\mathcal{G}_{iklm}-G^{>}_{in}\mathcal{G}_{jklm},$ (4.10)

where the detailed derivation is given in Appendix A.

Equation (4.10) is now simplified taking into account approximations i), ii)
and iii). Here, it is important to note that the single-particle Green’s
functions in the first line of equation (4.10) allow for different source
terms to be considered, e.g., in the first term,
$G^{>}_{im}G^{>}_{jn}G^{<}_{kl}=G^{>}_{im}\gamma_{jkln}^{s}=G^{>}_{jn}\gamma_{iklm}^{s}$,
and similar for the second term. Therefore, various terms containing two-
particle correlations can be rewritten in terms of fluctuations, i.e.,
$G^{>}_{im}\mathcal{G}_{jkln}$ or $G^{>}_{jn}\mathcal{G}_{iklm}$. This leads
to four possibilities to approximate $\Gamma$, but only two of these will be
considered here:

$\displaystyle\Gamma_{ijklmn}^{(1)}$ $\displaystyle\coloneqq\pm
G^{>}_{im}\gamma_{jkln}\pm G^{<}_{jl}\gamma_{ikmn},$ (4.11)
$\displaystyle\Gamma_{ijklmn}^{(2)}$ $\displaystyle\coloneqq\pm
G^{>}_{jn}\gamma_{iklm}\pm G^{<}_{km}\gamma_{ijln}.$ (4.12)

Inserting equations (4.11) and (4.12) into equation (2.19) leads to the
following approximation for the three-particle term,

$\displaystyle C^{\mathrm{(\Gamma),PA}}_{ijkl}\coloneqq$
$\displaystyle\pm\mathrm{i}\hbar\sum_{pqr}\left\\{w_{ipqr}\Gamma^{(1)}_{rqjpkl}+w_{pjqr}\Gamma^{(2)}_{iqrkpl}\right\\}$
$\displaystyle\mp\mathrm{i}\hbar\sum_{pqr}\left\\{w_{pqkr}\Gamma^{(1)}_{irjpql}+w_{pqrl}\Gamma^{(2)}_{ijrkqp}\right\\}$
(4.13) $\displaystyle=$
$\displaystyle\pi^{\mathrm{F}}_{ijkl}+\left[U^{(2),\mathrm{F}},\gamma\right]_{ijkl}.$
(4.14)

The commutator, $\left[U^{(2),\mathrm{F}},\gamma\right]$, contains the
exchange contribution that is missing in the 2M approximation, equation (4.3).
Together with the Hartree term in this equation, this yields the Hartree–Fock
Hamiltonian,

$\left[h^{(2),\rm{H}},\gamma\right]_{ijkl}+\left[U^{(2),\mathrm{F}},\gamma\right]_{ijkl}=\left[h^{(2),\rm{HF}},\gamma\right]_{ijkl}.$
(4.15)

Further, the term $\pi^{\mathrm{F}}$ in equation (4.14) is given by

$\pi^{\mathrm{F}}_{ijkl}\coloneqq\pi^{\pm}_{ijkl}-\pi_{ijkl},$ (4.16)

where $\pi^{\pm}$ is given by the replacement $w\rightarrow w^{\pm}$. This
term thus denotes exchange contributions to the polarization term which is
also present in the G1-G2 equations, cf. the term $\Pi^{\pm}$, equation
(2.18).

####  4.4.2 Polarization and $GW$ approximations for $\gamma$

Collecting all terms together, we can write down the EOM for $\gamma$:

$\mathrm{i}\hbar\partial_{t}\gamma^{\mathrm{PA}}_{ijkl}=\left[h^{(2),\rm{HF}},\gamma^{\mathrm{PA}}\right]_{ijkl}+\pi^{\pm}_{ijkl}.$
(4.17)

This is already a significant improvement compared to the 2M approximation,
equation (4.3), because we recovered the missing exchange corrections.

Let us now investigate to what approximation in the G1-G2-scheme equation
(4.17) corresponds. Proceeding as in section 4.3, the corresponding equation
for $\mathcal{G}$ is found by adding the equations for $\gamma^{\rm PA}$ and
$\gamma^{\rm PA,s}$, equation (4.9), and introducing $\mathcal{G}^{\rm
PA}\coloneqq\gamma^{\mathrm{PA}}\mp\gamma^{\mathrm{PA},s}$:

$\mathrm{i}\hbar\partial_{t}\mathcal{G}^{\mathrm{PA}}_{ijkl}=\left[h^{(2),\mathrm{HF}},\mathcal{G}^{\mathrm{PA}}\right]_{ijkl}+\Psi^{\pm}_{ijkl}+\Pi^{\pm}_{ijkl},$
(4.18)

where $\Psi^{\pm}$ and $\Pi^{\pm}$ are defined in equations (2.15) and (2.18).
This equation can be considered the fluctuation equivalent of the $GW$
approximation with exchange corrections included, i.e., to $GW^{\pm}$. Thus,
we conclude that the polarization approximation, $\gamma^{\rm PA}$, is
equivalent to $GW^{\pm}$, as long as the weak coupling conditions i)–iii) are
fulfilled.

However, this is not a conserving approximation because complete
(anti-)symmetrization requires to also include the particle-hole ladder (TPH)
terms which leads to the replacement $\Pi^{\pm}\to P^{\pm}$, cf. equation
(2.17) [9]. So far, the TPH term has not been identified in the fluctuation
approach. Thus, the polarization approximation, equation (4.17) is directly
applicable to systems for which the exchange corrections vanish.

In order to derive the $GW$ approximation without exchange, we now look for
modifications to equation (4.17) which do not contain the exchange terms,
i.e., $\Psi^{\pm}\to\Psi$ and $\Pi^{\pm}\to\Pi$. Such an equation for $\gamma$
can be indeed derived using condition iii), $|\mathcal{G}|\ll|\gamma^{(s)}|$,
and additionally neglecting the exchange contributions to the polarization
term,

$\displaystyle C^{(\Gamma)}_{ijkl}$
$\displaystyle\approx\left[U^{(2),\mathrm{F}},\gamma\right]_{ijkl},$ (4.19)
$\displaystyle\pi^{\mathrm{F}}_{ijkl}$ $\displaystyle\rightarrow 0.$ (4.20)

This leads to the following EOM for the two-particle fluctuation which we
denote by $\gamma^{\mathrm{GW}}$:

$\mathrm{i}\hbar\partial_{t}\gamma^{\mathrm{GW}}_{ijkl}=\Big{[}h^{(2),\mathrm{HF}},\gamma^{\mathrm{GW}}\Big{]}_{ijkl}+\pi_{ijkl},$
(4.21)

which, together with equation (4.9) for the source term, corresponds to the
familiar $GW$-G2 equation

$\mathrm{i}\hbar\partial_{t}\mathcal{G}^{\mathrm{GW}}_{ijkl}=\Big{[}h^{(2),\mathrm{HF}},\mathcal{G}^{\mathrm{GW}}\Big{]}_{ijkl}+\Psi_{ijkl}+\Pi_{ijkl}.$
(4.22)

The equivalence of the approximations of the fluctuation approach, equation
(4.21) and the G1-G2 scheme on the $GW$ level holds for weak to moderate
coupling, due to the assumptions i)–iii) and analogous approximations in the
G1-G2 scheme (the neglect of TPP terms). To extend the fluctuation approach to
strong coupling, one has to consider additional contributions from the
residual term $R^{(\gamma)}$ [cf. equation (4.6)] and from three-particle
fluctuations, $\Gamma$.

####  4.4.3 Polarization and $GW$ approximations expressed in terms of the
single-particle fluctuation $\delta\hat{G}$

Equivalently to approximating three-particle fluctuations, it is possible to
derive the quantum polarization approximation solely on the basis of single-
particle fluctuations. This can be done by again considering equation (3.17).
In any approximation, it must be noted that the expectation value of single-
particle fluctuations needs to be zero and thus also their time derivative,
i.e.,
$\langle\partial_{t}\delta\hat{G}\rangle=\partial_{t}\langle\delta\hat{G}\rangle=0$.
This means that any approximation of
$[\delta\hat{U}^{\mathrm{H}},\delta\hat{G}]$ needs to also properly account
for the collision integrals, $I^{(\gamma)}+I^{(\gamma)\dagger}$, in order to
ensure this condition.

As it was seen in the derivation of the quantum polarization approximation,
the missing contributions are given by the commutator of the effective two-
particle Fock-potential with two-particle fluctuations and the exchange part
of the polarization term. Similarly to the relation between correlations and
fluctuations at the two-particle level [cf. equation (3.8)], we now define an
analogous relation for the fluctuation operator $\delta\hat{\gamma}$, where we
take into account the linearity of equation (3.8)

$\displaystyle\delta\hat{\gamma}_{ijkl}$
$\displaystyle\equiv\delta\hat{G}_{ik}\delta\hat{G}_{jl}-\langle\delta\hat{G}_{ik}\delta\hat{G}_{jl}\rangle$
$\displaystyle\equiv\pm\delta\hat{\gamma}^{s}_{ijkl}+\delta\hat{\mathcal{G}}_{ijkl},$
(4.23) $\displaystyle\delta\hat{\gamma}^{s}_{ijkl}$ $\displaystyle\coloneqq
G^{>}_{il}\delta\hat{G}_{jk}+\delta\hat{G}_{il}G^{<}_{jk}.$ (4.24)

Analogously to the derivation of the TDHF approximation, the quantum
polarization approximation can be regarded as an approximation that neglects
‘‘fluctuations’’ of two-particle correlations, i.e.,
$\delta\hat{\mathcal{G}}_{ijkl}\equiv 0$,

$\delta\hat{\gamma}^{\mathrm{PA}}_{ijkl}\coloneqq\pm\delta\hat{\gamma}^{s}_{ijkl}=\pm\left\\{G^{>}_{il}\delta\hat{G}_{jk}+\delta\hat{G}_{il}G^{<}_{jk}\right\\}.$
(4.25)

This leads to the following EOM for $\delta\hat{G}^{\mathrm{PA}}$,

$\mathrm{i}\hbar\partial_{t}\delta\hat{G}^{\mathrm{PA}}_{ij}=\left[h^{\mathrm{HF}},\delta\hat{G}^{\mathrm{PA}}\right]_{ij}+\left[\delta\hat{U}^{\mathrm{HF}},G^{<}\right]_{ij}.$
(4.26)

We now verify that this equation is equivalent to EOM for
$\gamma^{\mathrm{PA}}$, by differentiating the definition

$\mathrm{i}\hbar\partial_{t}\gamma^{\mathrm{PA}}_{ijkl}=\mathrm{i}\hbar\partial_{t}\left\langle\delta\hat{G}^{\mathrm{PA}}_{ik}\delta\hat{G}^{\mathrm{PA}}_{jl}\right\rangle,$
(4.27)

which contains the following terms,

$\displaystyle\underbrace{\bigg{\langle}\left[h^{\mathrm{HF}},\delta\hat{G}^{\mathrm{PA}}\right]_{ik}\delta\hat{G}^{\mathrm{PA}}_{jl}\bigg{\rangle}+\left\langle\delta\hat{G}^{\mathrm{PA}}_{ik}\left[h^{\mathrm{HF}},\delta\hat{G}^{\mathrm{PA}}\right]_{jl}\right\rangle}_{\left[h^{(2),\rm{HF}},\gamma^{\mathrm{PA}}\right]_{ijkl}}$
$\displaystyle+$
$\displaystyle\underbrace{\left\langle\left[\delta\hat{U}^{\mathrm{HF}},G^{<}\right]_{ik}\delta\hat{G}^{\mathrm{PA}}_{jl}\right\rangle+\left\langle\delta\hat{G}^{\mathrm{PA}}_{ik}\left[\delta\hat{U}^{\mathrm{HF}},G^{<}\right]_{jl}\right\rangle}_{\pi_{ijkl}^{\pm}}.$
(4.28)

We have thus proven that approximation (4.26) is equivalent to the quantum
polarization approximation for $\gamma^{\rm PA}$, equation (4.17). Recall
that, to recover the polarization ($GW^{\pm}$) approximation for
$\mathcal{G}$, additionally, in the equation for $\gamma^{s}$ the $R$ term
should be neglected, giving rise to $\gamma^{\rm PA,s}$, equation (4.9).

In the same way, the standard $GW$ approximation can also be recovered. This
is achieved by assuming $\delta\hat{\gamma}_{ijkl}\approx\pm
G^{<}_{jk}\delta\hat{G}_{il}$, leading to the EOM for
$\delta\hat{G}^{\mathrm{GW}}$,

$\mathrm{i}\hbar\partial_{t}\delta\hat{G}^{\mathrm{GW}}_{ij}=\Big{[}h^{\mathrm{HF}},\delta\hat{G}^{\mathrm{GW}}\Big{]}_{ij}+\Big{[}\delta\hat{U}^{\mathrm{H}},G^{<}\Big{]}_{ij}.$
(4.29)

The two equations (4.26) and (4.29) are remarkable results. They contain the
complex physics corresponding to the nonequilibrium $GW$ approximation (with
and without exchange, respectively) in an equation for a single-particle
quantity $\delta\hat{G}$. This equation is closed and, obviously, much simpler
than the corresponding pair of equations in the G1-G2-scheme or the
corresponding KBE of two-particle NEGF theory.

The main problem, however, is that this EOM for the single-particle
fluctuations is an operator equation, for which a direct solution is not
possible. This problem can be solved approximately by using a semi-classical
averaging approach. The result provides an approximation for $GW$ and
$GW^{\pm}$, respectively, at the single-particle level, for weak to moderate
coupling. This will be discussed more in detail in section 5.5. However, it is
not trivially possible to include contributions to compensate for the residual
term [cf. equation (4.6)] at the single-particle level. Thus, these
approximations cannot be directly extended to the strong coupling regime in
this framework at the moment.

##  5 Stochastic mean-field theory

###  5.1 General concept

The basic ideas of the stochastic mean-field theory (SMF) [16] are the
replacement of quantum mechanical operators by stochastic quantities (random
realizations ‘‘$\lambda$’’), and of the quantum mechanical expectation value
by an arithmetic mean,

$\displaystyle\hat{A}$ $\displaystyle\longrightarrow A^{\lambda},$ (5.1)
$\displaystyle\langle\hat{A}\rangle$
$\displaystyle\longrightarrow\overline{A^{\lambda}}.$ (5.2)

In practice, random realizations ‘‘$\lambda$’’ of the initial state are
generated, according to a statistical ensemble, and then propagated in time.
This allows the operator equations to be solved approximately, and thus makes
possible a solution of the many-body problem at the single-particle level, by
sampling single-particle density matrices and propagating them using simple
mean-field dynamics. Each realization of the initial state can now be
propagated independently, thus reducing the solution of the many-body problem
to a) the construction of the probability distribution describing the initial
state, b) mean-field dynamics and c) semiclassical averaging [16].

However, by replacing operators with random variables, important properties of
operators are not properly accounted for, and the procedure requires special
care. For example, the expectation value of products of two non-commuting
operators, $\hat{A},\hat{B}$, should be symmetrized, i.e.,

$\hat{A}\hat{B}=\frac{1}{2}\left[\hat{A},\hat{B}\right]_{+}+\frac{1}{2}\left[\hat{A},\hat{B}\right]_{-}.$
(5.3)

Now, the commutator can be evaluated, and the symmetric expression (anti-
commutator) can be replaced with the random ensemble, i.e.,
$[\hat{A},\hat{B}]_{+}/2\rightarrow A^{\lambda}B^{\lambda}$. This procedure is
already known from other (semiclassical) approaches such as the symmetric
representation of phase space operators [42].

We now apply this procedure to our fluctuation approach and define symmetric
two-particle fluctuations as

$\tilde{\gamma}_{ijkl}\coloneqq\frac{\gamma_{ijkl}+\gamma_{jilk}}{2}=\frac{1}{2}\left\langle\left[\delta\hat{G}_{ik},\delta\hat{G}_{jl}\right]_{+}\right\rangle,$
(5.4)

and find, on the time-diagonal,

$\displaystyle\gamma_{ijkl}$
$\displaystyle=\tilde{\gamma}_{ijkl}+\frac{1}{2}\left\langle\left[\delta\hat{G}_{ik},\delta\hat{G}_{jl}\right]\right\rangle$
$\displaystyle=\tilde{\gamma}_{ijkl}\pm\frac{1}{2\mathrm{i}\hbar}\left\\{\delta_{il}G^{<}_{jk}-\delta_{jk}G^{<}_{il}\right\\}.$
(5.5)

Similarly to fluctuation operators, semiclassical fluctuations are defined as
deviations from the average,

$\Delta G^{\lambda}_{ij}\coloneqq
G^{<,\lambda}_{ij}-\overline{G^{<,\lambda}_{ij}}\equiv
G^{<,\lambda}_{ij}-\tilde{G}^{<}_{ij}.$ (5.6)

We now apply this approach to equation (3.15) and solve an EOM for
$G^{<,\lambda}$ or construct a system of differential equations for
semiclassical fluctuations, $\Delta G^{<,\lambda}$, and $\tilde{G}^{<}$.
Analogously to the quantum fluctuation approach, one can now derive an EOM for
$\overline{\Delta G^{\lambda}\Delta G^{\lambda}}$ that couples to
semiclassical three-particle fluctuations, therefore leading to a hierarchy of
coupled equations [43]. Thus, the solution of an EOM for $G^{<,\lambda}$ is
equivalent to the solution of the entire hierarchy. While this can be
considered as an advantage of the SMF theory, it is necessary to find a
suitable probability distribution which correctly describes all statistical
moments that are known from quantum fluctuations as well as their relation to
correlations. This will be discussed more in detail in section 5.3. Any
probability distribution that fails to correctly reproduce all moments leads
to the propagation of wrong contributions for higher order fluctuations.
However, constructing the correct distribution can be shown to be impossible
due to the differences of quantum mechanical and statistical expectation
values [44].

###  5.2 Observables

The SMF approach also allows one to calculate expectation values of
observables. For an arbitrary one-body observable $A^{(1)}$, the expectation
value can be computed as

$\big{\langle}\hat{A}^{(1)}\big{\rangle}=\pm\mathrm{i}\hbar\sum_{ij}A^{(1)}_{ij}G^{<}_{ji}.$
(5.7)

For a two-body observable $A^{(2)}$, we first have to consider the expectation
value in the framework of $\gamma$,

$\displaystyle\big{\langle}\hat{A}^{(2)}\big{\rangle}=-\hbar^{2}\sum_{ijkl}A^{(2)}_{ijkl}G^{(2),<}_{klij}=-\hbar^{2}\sum_{ijkl}A^{(2)}_{ijkl}\left(G^{<}_{ki}G^{<}_{lj}\mp\frac{1}{\mathrm{i}\hbar}\delta_{jk}G^{<}_{li}+\gamma_{klij}\right).$
(5.8)

The expectation value of $A^{(2)}$ is then computed within the SMF approach by
the symmetrized analog of equation (5.8)

$\displaystyle\big{\langle}\hat{A}^{(2)}\big{\rangle}=-\hbar^{2}\sum_{ijkl}A^{(2)}_{ijkl}\left(G^{<}_{ki}G^{<}_{lj}\mp\frac{1}{2\mathrm{i}\hbar}(\delta_{jk}G^{<}_{li}+\delta_{li}G^{<}_{jk})+\overline{\Delta
G^{\lambda}_{ki}\Delta G^{\lambda}_{lj}}\right).$ (5.9)

This can be generalized to an arbitrary $N$-particle observable.

###  5.3 Probability distribution

To construct the probability distribution describing the initial state at
$t=t_{0}$, we consider an ideal (uncorrelated) state, i.e.,
$\mathcal{G}(t_{0})=0$. We underline that this is not a restriction because,
starting from an ideal initial state a correlated initial state, can be
produced via the adiabatic switching protocol, e.g., [40], which we also apply
to the present fluctuation approach. This has the advantage that we only need
to consider the properties of the initially ideal system and the corresponding
probability distribution.

The expectation values of quantum fluctuation operators and their products are
then given by

$\displaystyle\left\langle\delta\hat{G}_{ij}(t_{0})\right\rangle$
$\displaystyle=0,$ (5.10)
$\displaystyle\left\langle\delta\hat{G}_{ik}(t_{0})\delta\hat{G}_{jl}(t_{0})\right\rangle$
$\displaystyle=-\frac{1}{\hbar^{2}}\delta_{il}\delta_{jk}n_{i}(1\pm n_{j}),$
(5.11)

with $n_{i}\coloneqq\pm\mathrm{i}\hbar G^{<}_{ij}(t_{0})\delta_{ij}$. Applying
the SMF approach to equations (5.10) and (5.11), we obtain

$\displaystyle\overline{\Delta G^{\lambda}_{ij}(t_{0})}$ $\displaystyle=0,$
(5.12) $\displaystyle\overline{\Delta G^{\lambda}_{ik}(t_{0})\Delta
G^{\lambda}_{jl}(t_{0})}$
$\displaystyle=-\frac{1}{2\hbar^{2}}\delta_{il}\delta_{jk}\\{n_{i}(1\pm
n_{j})+n_{j}(1\pm n_{i})\\}.$ (5.13)

Analogously, it is possible to derive expressions for higher moments of
fluctuations, although this becomes increasingly more difficult [44]. The
construction of a suitable probability distribution to best describe the
initial state proves to be one of the main challenges of the SMF approach. The
inability to construct a classical probability distribution can be seen when
considering the second moment at zero temperature for fermions. Here, it
follows

$\overline{\Delta G^{\lambda}_{ik}(t_{0})\Delta
G^{\lambda}_{jl}(t_{0})}=-\frac{1}{2\hbar^{2}}\delta_{il}\delta_{jk}\overline{\delta}_{n_{i}n_{j}},$
(5.14)

where we defined $\overline{\delta}_{ij}\coloneqq 1-\delta_{ij}$. This means,
in particular, that $\overline{|\Delta G^{\lambda}_{ii}|^{2}}=0$ and thus
$\Delta G^{\lambda}_{ii}\equiv 0$, for all realizations. However, the third
and fourth moments are nonzero for expectation values that contain terms with
matching indices, i.e., $\Delta G^{\lambda}_{ii}$. Further problems were
highlighted in reference [44] when considering these higher moments. It was
shown there that no standard probability distribution can simulate a quantum
system exactly.

There exist a variety of approaches to approximately construct the initial
state. Originally, standard probability distributions, e.g., a Gaussian
distribution, are chosen so that the first moments are correctly reproduced
without considering higher moments. This can be generalized by also
considering quasiprobability distributions such as the Husimi distribution
[45]. We also explore Gaussian sampling, in section 6.1. In addition, we
introduce a new possibility to construct the initial state, which is not based
on random sampling from a (quasi-)probability distribution but on a
deterministic procedure. This is shown in section 6.2.

###  5.4 Dynamics

Since two-particle fluctuations appear only in the collision term, we now
concentrate on the latter. Inserting equation (5.5) into equation (3.14) leads
to [we drop the superscript $(\gamma)$]

$\Big{[}I+I^{\dagger}\Big{]}_{ij}=\Big{[}I^{S}+I^{S\dagger}\Big{]}_{ij}+\Big{[}S+S^{\dagger}\Big{]}_{ij},$
(5.15)

where we used the definitions,

$\displaystyle I^{S}_{ij}\coloneqq$
$\displaystyle\pm\mathrm{i}\hbar\sum_{klp}w_{iklp}\tilde{\gamma}_{plkj},$
(5.16) $\displaystyle S_{ij}\coloneqq$
$\displaystyle\frac{1}{2}\sum_{kl}w_{kljk}G^{<}_{il}.$ (5.17)

Therefore, using symmetric fluctuations, yields the following EOM for $G^{<}$:

$\displaystyle\mathrm{i}\hbar\partial_{t}G^{<}_{ij}=$
$\displaystyle\Big{[}h^{\mathrm{H}},G^{<}\Big{]}_{ij}+\Big{[}S+S^{\dagger}\Big{]}_{ij}+\Big{[}I^{S}+I^{S\dagger}\Big{]}_{ij}.$
(5.18)

Replacing $\tilde{\gamma}$ in equation (5.18) with $\overline{\Delta
G^{\lambda}\Delta G^{\lambda}}$ results in the following equation for the
symmetrized expectation value $\tilde{G}^{<}$:

$\displaystyle\mathrm{i}\hbar\partial_{t}\tilde{G}^{<}_{ij}=$
$\displaystyle\left[h^{\mathrm{H}},\tilde{G}^{<}\right]_{ij}+\Big{[}S+S^{\dagger}\Big{]}_{ij}+\Big{[}I^{\mathrm{SMF}}+I^{\mathrm{SMF}\dagger}\Big{]}_{ij},$
(5.19)

with the SMF-collision term defined as

$I^{\mathrm{SMF}}_{ij}\coloneqq\pm\mathrm{i}\hbar\sum_{klp}w_{iklp}\overline{\Delta
G^{\lambda}_{pk}\Delta G^{\lambda}_{lj}},$ (5.20)

where the semi-classical averaging, equation (5.2), is used. Instead of using
an additional equation of motion for the expectation value on the right and
solving the entire fluctuation hierarchy, we can now proceed differently: we
construct an EOM for $\Delta G^{\lambda}$ that is analogous to equation (3.17)
and solve it. The expectation value is computed at the final step.

To accomplish this task, we first construct an EOM for $G^{<,\lambda}$ which
is itself symmetric:

$\mathrm{i}\hbar\partial_{t}G^{<,\lambda}_{ij}=\Big{[}h^{\mathrm{H},\lambda},G^{<,\lambda}\Big{]}_{ij}+\Big{[}S^{\lambda}+S^{\lambda\dagger}\Big{]}_{ij},$
(5.21)

where we introduced the definition

$S^{\lambda}_{ij}\coloneqq\frac{1}{2}\sum_{kl}w_{kljk}G^{<,\lambda}_{il}.$
(5.22)

The EOM for $\Delta G^{<,\lambda}$ becomes,

$\displaystyle\mathrm{i}\hbar\partial_{t}\Delta G^{\lambda}_{ij}=$
$\displaystyle\Big{[}h^{\mathrm{H}},\Delta
G^{<,\lambda}\Big{]}_{ij}+\Big{[}\Delta
U^{\mathrm{H,\lambda}},\tilde{G}^{<}\Big{]}_{ij}+\Big{[}\Delta
S^{\lambda}+\Delta S^{\lambda\dagger}\Big{]}_{ij}+\Big{[}\Delta
U^{\mathrm{H,\lambda}},\Delta G^{<,\lambda}\Big{]}_{ij}$
$\displaystyle-\Big{[}I^{\mathrm{SMF}}+I^{\mathrm{SMF}\dagger}\Big{]}_{ij},$
(5.23)

where we used the definitions

$\displaystyle\Delta S^{\lambda}_{ij}\coloneqq$
$\displaystyle\frac{1}{2}\sum_{kl}w_{kljk}\Delta G^{\lambda}_{il},$ (5.24)
$\displaystyle\Delta U^{\mathrm{H}}_{ij}\coloneqq$
$\displaystyle\pm\mathrm{i}\hbar\sum_{kl}w_{ikjl}\Delta G^{<}_{lk}.$ (5.25)

When three-particle contributions are neglected, as in our case,
$\tilde{\gamma}$ obeys the same EOM as $\gamma$, due to the linearity of the
equation.

For completeness, we mention that, when symmetric three-particle contributions
are taken into account, there arise differences in the equations of motion
that are analogous to the case of the single-particle EOM. This leads to very
complex equations, which are beyond the scope of the present work.

###  5.5 Stochastic GW approximation

The SMF approach can now be used to find new approximations for quantum
fluctuations. Applying the SMF approach to the GW approximation [cf. equation
(4.29)] allows for a treatment of this two-particle approximation at the
single-particle level, therefore significantly reducing computational effort.

Additionally, the neglect of any coupling to higher order fluctuations causes
the neglect of contributions of higher moments, thus potentially minimizing
the impact of choosing an approximate probability distribution. The stochastic
GW approximation (SGW) [stochastic version of equation (4.29)] takes the
following form:

$\displaystyle\mathrm{i}\hbar\partial_{t}\tilde{G}^{<,\mathrm{SGW}}_{ij}=$
$\displaystyle\Big{[}h^{\mathrm{H}},\tilde{G}^{<,\mathrm{SGW}}\Big{]}_{ij}+\Big{[}S+S^{\dagger}\Big{]}_{ij}+\Big{[}I^{\mathrm{SGW}}+I^{\mathrm{SGW}\dagger}\Big{]}_{ij},$
(5.26) $\displaystyle\mathrm{i}\hbar\partial_{t}\Delta
G^{\mathrm{SGW},\lambda}_{ij}=$ $\displaystyle\Big{[}h^{\mathrm{HF}},\Delta
G^{\mathrm{SGW},\lambda}\Big{]}_{ij}+\Big{[}\Delta
U^{\mathrm{H,SGW},\lambda},G^{<,\mathrm{SGW}}\Big{]}_{ij},$ (5.27)

where the following definitions were introduced,

$\displaystyle I^{\mathrm{SGW}}_{ij}$
$\displaystyle\coloneqq\pm\mathrm{i}\hbar\sum_{klp}w_{iklp}\overline{\Delta
G^{\mathrm{SGW},\lambda}_{pk}\Delta G^{\mathrm{SGW},\lambda}_{lj}},$ (5.28)
$\displaystyle\Delta U^{\mathrm{H,SGW},\lambda}_{ij}$
$\displaystyle\coloneqq\pm\mathrm{i}\hbar\sum_{kl}w_{ikjl}\Delta
G^{\mathrm{SGW},\lambda}_{lk}.$ (5.29)

Instead of considering equation (4.29), we can also derive an analogous
expression for equation (4.26) by replacing $\Delta U^{\mathrm{H,PA},\lambda}$
with $\Delta U^{\mathrm{HF,PA},\lambda}$. This gives rise to a set of
equations that is equivalent to the quantum polarization approximation which
will be called stochastic polarization approximation (SPA).

The SGW fluctuation equations, (5.26) and (5.27), as well as the SPA equations
should be solved together with the EOM for the single-particle Green’s
function. The resulting system of equations has got the same complexity as
standard TDHF equations where the numerical scaling (CPU time) of this
approximation is mostly determined by the number of samples $K$. Therefore,
the total numerical scaling of the SPA and SGW models is given by
$\mathcal{O}(KN_{b}^{4}N_{t})$, where $N_{b}$ is the number of basis states
and $N_{t}$ is the number of time steps. It is interesting to compare this to
the scaling of the $GW$ approximation within the G1-G2-scheme. Therein, the
CPU-time scaling is [8] of order $\mathcal{O}(N_{b}^{6}N_{t})$. Thus, we find
that the stochastic approximations are advantageous for (up to coefficients
that may be different) $K<N_{b}^{2}$. The specific choice of sampling methods
provides flexibility and allows for further optimization depending on the
system.

##  6 Sampling

Since it is shown to be impossible to construct a classical probability
distribution in order to correctly describe the initial quantum state, the
standard approaches to SMF only consider the known probability distributions
such as a Gaussian distribution or a uniform distribution. In what follows, we
discuss some of the possibilities to construct the initial state, though, we
only consider a fermionic system at $T=0$ with a finite basis.

###  6.1 Stochastic sampling

The first approach we want to discuss is the standard approach of
approximating the quantum state with a known probability distribution. For the
construction of a suitable probability distribution, we define the variance of
the distribution as

$\sigma_{ij}\coloneqq\frac{1}{2}\overline{\delta}_{n_{i}n_{j}}.$ (6.1)

In the original formulation of the SMF approach, a complex Gaussian
distribution was assumed [44],

$P^{N}_{ij}(x,y)=\frac{1}{2\piup\sigma^{2}_{ij}}\exp\left(-\frac{x^{2}+y^{2}}{2\sigma^{2}_{ij}}\right),$
(6.2)

where $x,y$ denote the real and imaginary part and are independently sampled.
This distribution correctly reproduces the first and second moments, but not
the higher moments. An alternative is given by a complex generalization of the
two-point distribution,

$\displaystyle P^{2p}_{ij}(x,y)=$
$\displaystyle\frac{1}{4}\Big{[}\delta(y)\\{\delta(x+\sqrt{2}\sigma_{ij})+\delta(x-\sqrt{2}\sigma_{ij})\\}+\delta(x)\\{\delta(y+\sqrt{2}\sigma_{ij})+\delta(y-\sqrt{2}\sigma_{ij})\\}\Big{]}.$
(6.3)

While this distribution also fails to correctly reproduce higher moments, it
approximates them better than the Gaussian distribution. In reference [44],
several distributions were compared for the exactly solvable
Lipkin–Meshkov–Glick model, where it was shown that the two-point distribution
had the best agreement with the exact solution. This further illustrates the
importance of an appropriate distribution and motivates the construction of
model distributions beyond the simple Gaussian and two-point assumptions. In
reference [45], the possibility of a quasiprobability distribution was
discussed and compared to previous approaches, where it was shown to further
improve the accuracy of the SMF theory.

In practical application, instead of calculating the expectation value, one
usually generates samples from the (quasi-)probability distribution and
calculates the average over all $K$ samples, i.e.,

$\overline{A^{\lambda}}\approx\frac{1}{K}\sum_{\lambda=1}^{K}A^{\lambda},$
(6.4)

which becomes exact in the limit $K\rightarrow\infty$. Thus, to achieve the
desired accuracy for the expectation value, the number $K$ should be chosen
sufficiently large.

###  6.2 Deterministic sampling

As mentioned in section 6.1, in most practical applications, random
realizations of the initial state are generated, and the expectation value of
the ensemble is approximated by the average over all samples, cf. equation
(5.2). By contrast, the idea of the present sampling method is to construct
realizations of the initial state so that the average is exactly equal to the
expectation value. To this end, we consider the following system of nonlinear
equations

$\displaystyle\sum_{\lambda=1}^{M}\Delta n^{\lambda}_{ij}$ $\displaystyle=0,$
(6.5) $\displaystyle\sum_{\lambda=1}^{M}\Delta n^{\lambda}_{ik}\Delta
n^{\lambda}_{jl}$
$\displaystyle=\frac{M}{2}\delta_{il}\delta_{jk}\overline{\delta}_{n_{i}n_{j}},$
(6.6)

where we defined $\Delta n^{\lambda}_{ij}\coloneqq-\mathrm{i}\hbar\Delta
G^{\lambda}_{ij}(t_{0})$. By finding a solution to these equations, we can
construct an ensemble of the initial state that exactly fulfills the
constraints given by the first two moments. Here, $M$ is a finite parameter
that should be chosen such that the system of nonlinear equations has a
solution and thus depends on the size of the basis.

In what follows, we consider a fermionic system with the spin configurations
$\uparrow,\downarrow$, though this approach can be generalized to arbitrary
spins. Additionally, we assume spin symmetry of the initial state, i.e.,
$n_{i}^{\uparrow}=n_{i}^{\downarrow}$. Let $N_{b}\in\mathbb{N}$ denote the
size of the basis for one spin component and $N_{p},N_{h}\in\mathbb{N}$ denote
the number of occupied and unoccupied orbitals, respectively, for one spin
component, with $N_{p}+N_{h}=N_{b}$. Without loss of generality, we can assume
$n_{i}^{\uparrow}=n_{i}^{\downarrow}=1$ for $i=1,\dots,N_{p}$ and
$n_{i}^{\uparrow}=n_{i}^{\downarrow}=0$ for $i=N_{p}+1,\dots,N_{b}$.

The goal of this approach is not to find all solutions of the system of
equations given by (6.5) and (6.6), for a given $M$, but rather to construct a
single solution while minimizing the number of samples. The simplest way to
fulfill equation (6.5) is to sample also $-\Delta n^{\lambda}$, for every
$\Delta n^{\lambda}$, thus ensuring that the sum over all samples for every
component is equal to zero. Additionally, equation (6.6) can be partially
fulfilled by setting $\Delta n^{\lambda}_{ij}=0$ for
$(i,j)\in\\{1,\dots,N_{p}\\}\times\\{1,\dots,N_{h}\\}\cup\\{N_{p}+1,\dots,N_{b}\\}\times\\{N_{h}+1,\dots,N_{b}\\}$.
The matrix representing a single sample can then be written as

$\Delta n^{\lambda,\sigma}=\begin{pmatrix}0&\mathcal{A}^{\lambda,\sigma}\\\
\mathcal{A}^{\lambda,\sigma\dagger}&0\end{pmatrix}$ (6.7)

with $\mathcal{A}^{\lambda,\sigma}\in\mathbb{C}^{N_{h}\times N_{p}}$. Here, it
was used that $\Delta n^{\lambda}$ is self-adjoint due to the properties of
the single-particle density matrix. We can define a bijection666Any bijection
$\varphi$ can be chosen here for the construction of the solution.
$\varphi:\,\\{1,\dots,N_{p}\\}\times\\{1,\dots,N_{h}\\}\rightarrow\\{1,\dots,M\\}$
with $M\coloneqq N_{p}N_{h}$. If we now set $K=8M$, we can find a solution to
equations (6.5) and (6.6) of the form

$\displaystyle\mathcal{A}_{ij}^{\alpha\beta,\uparrow}$
$\displaystyle=\begin{cases}\mathrm{i}^{\beta}M^{1/2},&\mbox{for
}\varphi(i,j)=\alpha,\,\beta=1,\dots,4,\\\ 0,&\mbox{else},\end{cases}$ (6.8)
$\displaystyle\mathcal{A}_{ij}^{\alpha\beta,\downarrow}$
$\displaystyle=\begin{cases}\mathrm{i}^{\beta}M^{1/2},&\mbox{for
}\varphi(i,j)=\alpha,\,\beta=5,\dots,8,\\\ 0,&\mbox{else},\end{cases}$ (6.9)

for $(\alpha,\beta)\in\\{1,\dots,M\\}\times\\{1,\dots,8\\}$. Using this
approach for the construction of a solution, we require a total of
$K=8N_{p}N_{h}$ samples, for each spin component, which means that this
approach becomes numerically favorable for systems that are either mostly
occupied or unoccupied.

For example, for a two-state system with $N_{p}=N_{h}=1$, this leads to 16
samples of the form

$\displaystyle\Delta n^{1,\uparrow}=\Delta n^{5,\downarrow}=-\Delta
n^{2,\uparrow}=-\Delta n^{6,\downarrow}$
$\displaystyle=\begin{pmatrix}0&\mathrm{i}\\\ -\mathrm{i}&0\end{pmatrix},$
(6.10) $\displaystyle\Delta n^{3,\uparrow}=\Delta n^{7,\downarrow}=-\Delta
n^{4,\uparrow}=-\Delta n^{8,\downarrow}$ $\displaystyle=\begin{pmatrix}0&1\\\
1&0\end{pmatrix},$ (6.11) $\displaystyle\Delta n^{5,\uparrow}=\Delta
n^{6,\uparrow}=\Delta n^{7,\uparrow}=\Delta n^{8,\uparrow}$
$\displaystyle=\begin{pmatrix}0&0\\\ 0&0\end{pmatrix},$ (6.12)
$\displaystyle\Delta n^{1,\downarrow}=\Delta n^{2,\downarrow}=\Delta
n^{3,\downarrow}=\Delta n^{4,\downarrow}$ $\displaystyle=\begin{pmatrix}0&0\\\
0&0\end{pmatrix}.$ (6.13)

Although the first two moments are exactly reproduced, higher moments are not
well approximated. In fact, the $n$-th moment can be proportional to
$M^{({n-2})/{2}}$, for $n>2$, thus making this approach barely applicable to
standard SMF theory. However, this approach turns out to be well suited for
the SPA and SGW, because only the first two moments directly contribute here.
Using this approach, we find that the total numerical scaling for SGW is given
by $\mathcal{O}(N_{p}N_{h}N_{b}^{4}N_{t})$. This sampling method leads to the
same numerical scaling as for $GW$ if we consider a system with
$N_{p}=N_{h}=N_{b}/2$. However, we find that for systems away from half-
filling, this approach provides a significant improvement.

##  7 Application to the Fermi–Hubbard model

In what follows, we consider the Hubbard model since it allows for simple
numerical tests of our theoretical results. It is not only a key model of
correlated electrons in condensed matter but also directly applicable to cold
atoms in optical lattices. Furthermore, here extensive NEGF and G1-G2 results
are available, e.g., references [6, 46, 7] that can be used for comparison.

###  7.1 Hubbard Hamiltonian

For the Fermi–Hubbard model, the general pair-interaction of equation (2.2)
transforms into

$w_{ijkl}^{\alpha\beta\gamma\delta}=U\delta_{ij}\delta_{ik}\delta_{il}\delta_{\alpha\gamma}\delta_{\beta\delta}\overline{\delta}_{\alpha\beta},$
(7.1)

with the on-site interaction $U$ and the spin components denoted by greek
indices. Additionally, the kinetic energy is replaced by a hopping Hamiltonian

$h_{ij}=-\delta_{\langle i,j\rangle}J,$ (7.2)

which describes nearest-neighbor hopping. The total Hamiltonian is then given
by

$\hat{H}=-J\sum_{\langle
i,j\rangle}\sum_{\sigma\in\\{\uparrow,\downarrow\\}}\hat{c}^{\dagger}_{i\sigma}\hat{c}_{j\sigma}+U\sum_{i}\hat{n}_{i}^{\uparrow}\hat{n}_{i}^{\downarrow}.$
(7.3)

###  7.2 Implementation

The EOMs for the single-particle Green’s function in SGW, equations (5.26) and
(5.27), take the following form, where we drop all superscripts
‘‘$\mathrm{SGW}$’’ and ‘‘$<$’’, for simplicity, e.g.,
$G^{<,\mathrm{SGW}}_{ij}\equiv G_{ij}$,

$\displaystyle\mathrm{i}\hbar\partial_{t}G^{\sigma}_{ij}=$
$\displaystyle\Big{[}h^{\sigma},G^{\sigma}\Big{]}_{ij}+\Big{[}I+I^{\dagger}\Big{]}^{\sigma}_{ij},$
(7.4) $\displaystyle\mathrm{i}\hbar\partial_{t}\Delta
G^{\lambda,\sigma}_{ij}=$ $\displaystyle\Big{[}h^{\sigma},\Delta
G^{\lambda,\sigma}\Big{]}_{ij}+\Big{[}\Delta
U^{\lambda,\sigma},G^{\sigma}\Big{]}_{ij},$ (7.5)

where the Hartree(–Fock) Hamiltonian and fluctuation Hartree-potential in
equations (7.4) and (7.5) in the Hubbard basis become [cf. equations (2.10)
and (5.29)]

$\displaystyle h_{ij}^{\sigma}$ $\displaystyle\equiv
h_{ij}^{\mathrm{HF},\sigma}\equiv h_{ij}^{\mathrm{H},\sigma}=-\delta_{\langle
i,j\rangle}J-\mathrm{i}\hbar\delta_{ij}UG^{{\overline{\sigma}}}_{ii},$ (7.6)
$\displaystyle\Delta U^{\lambda,\sigma}_{ij}$ $\displaystyle\equiv\Delta
U_{ij}^{\mathrm{H},\lambda,\sigma}\equiv\Delta
U_{ij}^{\mathrm{HF},\lambda,\sigma}=-\mathrm{i}\hbar\delta_{ij}U\Delta
G^{\lambda,{\overline{\sigma}}}_{ii}.$ (7.7)

Here, $\sigma=\uparrow(\downarrow)$ implies
$\overline{\sigma}=\downarrow(\uparrow)$. This shows that all exchange
contributions vanish, due to the specific choice of the pair-interaction so
that the SGW is equivalent to SPA.

The collision term in equation (7.4) takes the form [cf. equation (5.28)]

$I^{\sigma}_{ij}=-\mathrm{i}\hbar U\overline{\Delta
G^{\lambda,{\overline{\sigma}}}_{ii}\Delta G^{\lambda,\sigma}_{ij}},$ (7.8)

and the contributions due to symmetrization take the form [cf. equation
(5.17)]

$S_{ij}^{\sigma}=\frac{1}{2}UG^{\sigma}_{ij}=-S_{ij}^{\dagger,\sigma}.$ (7.9)

Thus, symmetrization also does not lead to any additional contributions in the
Hubbard basis. The initial state of the system in the natural orbital basis of
$n(t_{0})$ is chosen such that

$\displaystyle G_{ij}^{\sigma}(t_{0})$
$\displaystyle=-\frac{1}{\mathrm{i}\hbar}\delta_{ij}n^{\sigma}_{i},$ (7.10)
$\displaystyle\overline{\Delta G_{ij}^{\lambda,\sigma}(t_{0})}$
$\displaystyle=0,$ (7.11) $\displaystyle\overline{\Delta
G^{\lambda,\sigma}_{ij}(t_{0})\Delta G^{\lambda,\sigma^{\prime}}_{kl}(t_{0})}$
$\displaystyle=-\frac{1}{2\hbar^{2}}\delta_{il}\delta_{jk}\delta_{\sigma\sigma^{\prime}}\overline{\delta}_{n^{\sigma}_{i}n_{j}^{\sigma^{\prime}}}.$
(7.12)

Depending on the configuration of the system, it becomes necessary to perform
a transformation from the natural orbital basis to the Hubbard basis. This can
be achieved by diagonalization of the Hamiltonian and the transformation of
$G$ and $\Delta G^{\lambda}$ using its eigenvectors.

In general, it is necessary to compute a nontrivial interacting ground state
from which the externally driven dynamics start. Here, this is done using the
so-called ‘‘adiabatic switching method’’ [6] by replacing the on-site
interaction $U$ with a time dependent interaction $U(t)$. Calculations start
at $t_{s}$ with an uncorrelated ground state, with $U(t_{s})=0$. The on-site
interaction is increased monotonously and rather slowly such that at $t_{0}$
the system is in a fully correlated ground state with $U(t_{0})=U$.

The observables we consider are the densities on all sites, $n_{i}(t)$, and,
in addition, the kinetic energy ($E_{\mathrm{kin}}$), Hartree–Fock
($E_{\mathrm{HF}}$) and correlation energies ($E_{\mathrm{cor}}$) which follow
from equations (5.7) and (5.9),

$\displaystyle E_{\mathrm{kin}}$ $\displaystyle=\mathrm{i}\hbar J\sum_{\langle
i,j\rangle}\sum_{\sigma}G^{\sigma}_{ji},\qquad
E_{\mathrm{HF}}=-\hbar^{2}U\sum_{i}G^{\uparrow}_{ii}G^{\downarrow}_{ii},\qquad
E_{\mathrm{cor}}=-\hbar^{2}U\sum_{i}\overline{\Delta
G^{\lambda,\uparrow}_{ii}\Delta G^{\lambda,\downarrow}_{ii}}.$ (7.13)

In Appendix B it is shown that the approximations presented in section 4.4.2
conserve the total energy for the Fermi–Hubbard model. This directly implies
that total energy conservation is also given for SGW/SPA due to the structure
of the underlying set of equations.

###  7.3 Numerical results

Figure 1:  (Colour online) Panel (a) corresponds to the density evolution of a
half-filled $8$-site chain at $U=0.1J$ from SGW for initial states generated
from deterministic sampling (section 6.2) and random sampling (section 6.1)
using a complex four-point and Gaussian distribution with $10^{4}$ samples for
each spin configuration. Starting point of the simulations was a
noninteracting initial state. The deviation displayed in (b) is computed as
$\Delta n=n^{\mathrm{det}}-n^{\mathrm{random}}$.

Figure 2: (Colour online) Same as figure 2, but for $U/J=1.0$. Figure 3:
(Colour online) Density evolution of a half-filled $20$-site chain at $U=0.1J$
(a) and $U=1.0J$ (c) from SGW, for initial states generated from deterministic
and random sampling using a complex four-point distribution with $10^{4}$
samples for each spin configuration. Starting point of the simulations was an
uncorrelated initial state. The deviation for $U=0.1J$ (b) and $U=1.0J$ (d) is
computed as $\Delta n=n^{\mathrm{det}}-n^{\mathrm{random}}$.

####  7.3.1 Comparison of sampling methods

First, we compare three different sampling methods. These are the complex
four-point and Gaussian distributions, as well as the deterministic sampling
method. In particular, we compare numerical results for 1D chains without
periodic boundary conditions for different system sizes and interaction
strengths. The test system is a chain with eight sites, of which the leftmost
four sites are fully occupied and the right-hand half is empty, with a weak
on-site interaction ($U=0.1J$). For the two stochastic sampling methods
$10^{4}$ random realizations are used, for each spin component. The
deterministic sampling method has $128$ samples for each spin component.

Figure 4: (Colour online) Mean absolute deviation $\varepsilon$ [cf. equation
(7.14)] for half-filled chains with $N_{b}=8,20,30$ using the complex four-
point distribution and a fit, which is proportional to $1/\sqrt{K}$.

Figure 2 shows that all distributions display good agreement with each other.
In particular, the distributions have absolute deviations of less than $0.02$
for all times, where the deviation is defined as the density difference
between the deterministic sampling method and the stochastic methods. Thus,
the results of the deterministic sampling method can be interpreted as the
limit of the dynamics of stochastic approaches and that all distributions with
the same first two moments lead to the same dynamics. This is a significant
difference compared to the standard SMF approach, since there the choice of
the distribution for generating the initial state has a considerable influence
on the dynamics of the system.

Next, we extend the simulations for this system to an interaction strength of
$U=1.0J$. Despite the same number of samples, we see in figure 2 significant
deviations in the density dynamics for times $t\gtrsim 50\hbar/J$, showing
that a stronger coupling leads to greater sensitivity with respect to
deviations from the ideal initial state. Yet, for the cases considered, there
does not appear to be any significant difference between the two stochastic
sampling methods.

Therefore, it is of interest to investigate this dependence on the interaction
strength also for other systems. To this end, we consider a chain with twenty
sites, where again the ten leftmost sites are fully occupied and the remaining
ones are empty. Again, we consider the cases of $U=0.1J$ and $U=1.0J$ and use
$10^{4}$ samples for each spin component for the stochastic sampling method.
Here, we only use the complex four-point distribution, since all considered
probability distributions lead to equivalent results within this
approximation. For the deterministic sampling method, we now use $800$ samples
for each spin component. In figure 3 we see the behavior for the twenty-site
chain similar to the behaviour of the eight-site chain for a weak on-site
interaction: the absolute value of the deviation is smaller than $0.02$, for
all times. Differences between the two system sizes are more apparent for an
interaction strength of $U=1.0J$. Here, the density deviations are slightly
larger than for $U=0.1J$, but they do not exceed 0.04 in magnitude. The
dynamics of the two considered systems differ mainly in their speed.
Deviations therefore seem to occur earlier the faster the dynamics proceed.
This suggests that for a strong coupling, large systems are less sensitive to
deviations from the exact initial fluctuations than smaller systems. Hence,
the deterministic sampling method is advantageous over the stochastic one for
small systems, especially for a stronger coupling. For the systems we
considered, the required number of samples for a converged calculation was
found to depend only weakly on the system size. This is further highlighted
when considering the mean absolute error of the initial samples, which we
define by

$\varepsilon\coloneqq\frac{1}{4N_{b}^{4}}\sum_{ijkl}\sum_{\sigma\sigma^{\prime}}\left|\overline{\Delta
G^{\lambda,\sigma}_{ik}(t_{0})\Delta
G^{\lambda,\sigma^{\prime}}_{jl}(t_{0})}-\tilde{\gamma}^{\sigma\sigma^{\prime}\sigma\sigma^{\prime}}_{ijkl}(t_{0})\right|.$
(7.14)

Using this expression, we can quantify the deviations from the ideal initial
state given by equation (7.12). Figure 4 shows that the mean absolute error is
nearly independent of the system size $N_{b}$. Additionally, our fit with a
function proportional to $1/\sqrt{K}$ illustrates the expected trend due to
the central limit theorem. Assuming $10^{4}$ samples for each spin component
and half-filling so that $N_{p}=N_{h}=N_{b}/2$, stochastic sampling appears to
be advantageous, especially for system sizes larger than $100$ sites.

####  7.3.2 Comparison of SGW to _GW_ -G1-G2 and exact results

We now compare the results obtained from SGW with the exact solution using
exact diagonalization and the $GW$ approximation of the G1-G2 scheme. All
further calculations for SGW are done using the deterministic sampling method.
We again consider 1D chains without periodic boundary conditions for different
on-site interactions and systems sizes.

Figure 5: (Colour online) Evolution of the density from SGW, the $GW$-G1-G2
scheme and exact diagonalization for a half-filled $8$-site chain at $U=0.1J$.
The initial state of the system was uncorrelated.

Figure 6: (Colour online) Evolution of the kinetic, HF and correlation energy
from SGW, the $GW$-G1-G2 scheme and exact diagonalization for a half-filled
$8$-site chain at $U=0.1J$. The initial state of the system was uncorrelated.

Figure 7: (Colour online) Same as figure 6, but for $U=0.5J$.

Figure 8: (Colour online) Evolution of the kinetic, HF and correlation energy
from SGW, the $GW$-G1-G2 scheme and exact diagonalization for a half-filled
$8$-site chain at $U=0.5J$. The initial state of the system was uncorrelated.
Figure 9: (Colour online) Same as figure 6, but for $U=1.0J$.

First, we consider a half-filled eight-site chain with weak on-site
interaction ($U=0.1J$). The initial state is chosen so that the leftmost half
of the sites is fully occupied. It can be seen in figure 6 that the density
dynamics from SGW and the $GW$-G1-G2 scheme show good agreement with each
other. Both reproduce the exact result reasonably well for times $t\lesssim
40\hbar/J$, after which SGW qualitatively displays a better agreement for
times $t\lesssim 50\hbar/J$. After this point, deviations between the two
approximations and the exact result start to increase, with the $GW$
approximation showing slightly better density dynamics than SGW with respect
to the exact solution. However, this trend is only observable at the beginning
and, already for times $t\lesssim 70\hbar/J$, there are only minimal
deviations between the two approximations, although both fail to correctly
reproduce the exact dynamics of the system. While the exact dynamics displays
some sort of damping behavior with increasing time, we observe a revival of
the dynamics for both approximations which is due to the missing damping in
the HF propagators of the GKBA approach.

Next, consider the energy dynamics of this system. As seen in figure 6 the
same overall trend as for the density dynamics can be observed. However, the
different energies display a different behavior, i.e., the HF energy shows
good agreement for all calculations for times $t\lesssim 50\hbar/J$, after
which a revival of the dynamics can be observed for SGW and $GW$. Both
approximations show a qualitative good agreement and only small deviations
arise. Similar results can be seen when considering the kinetic energy.
However, already for times $t\gtrsim 20\hbar/J$, there are significant
differences between the two approximations and the exact result. However, SGW
and $GW$ show a very good agreement between each other, for the entire time
propagation. For times $t\gtrsim 70\hbar/J$, we can see a good agreement for
all three calculations. Finally, consider the correlation energy. Here, we
observe the same overall behavior as for the kinetic energy, although SGW and
$GW$ show deviations for later times that compensate the difference in HF
energy. In conclusion, the energy dynamics shows an analogous behavior to the
density dynamics. Both SGW and $GW$ fail to correctly reproduce the kinetic
and correlation energy by significantly underestimating correlations.
Additionally, we see that SGW conserves the total energy within the accuracy
of the numerical calculation, thus confirming the results presented in
Appendix B.

Figure 10: (Colour online) Density evolution of a half-filled 20-site chain at
$U=0.1J$ (a) and $U=1.0J$ (b) from SGW and the $GW$-G1-G2 scheme. (c) and (d)
show the density evolution for a 30-site chain for $U=0.1J$ and $U=1.0J$,
respectively. Starting point of the simulations was a uncorrelated initial
state.

Figure 11: (Colour online) Evolution of the kinetic, HF and correlation energy
from SGW and the $GW$-G1-G2 scheme for a half-filled 30-site chain at $U=0.1J$
(a) and $U=1.0J$ (b). The initial state of the system was uncorrelated.

Figure 12: (Colour online) Density dynamics on the first site of a 30-site
chain with the leftmost five sites fully occupied for $U=1.0J$ from SGW and
the $GW$-G1-G2 scheme. Starting point of the calculations was an uncorrelated
state. The $GW$ results diverge for $t\geqslant 55\hbar/J$.

Figure 13: (Colour online) Dynamics of the kinetic, HF and correlation
energies of the setup described in figure 13.  Figure 14: (Colour online)
Ground state calculations of the kinetic, HF, correlation and total energies
for $U=0.0J,\dots,2.0J$ for a half-filled $8$-site chain. SGW and $GW$
energies are calculated using the adiabatic switching method and the exact
energies by direct calculation of the correlated ground state using exact
diagonalization.

Next, we consider the same system but with an increased on-site interaction of
$U=0.5J$. In figure 8 we see that the SGW and $GW$ show a very good agreement
between each other for times $t\lesssim 40\hbar/J$, whereas after this time
small differences show up. For times $t\gtrsim 60\hbar/J$, the oscillations of
the density for both approximations have a similar amplitude. At the same
time, both approximations reproduce the exact dynamics only for $t\lesssim
10\hbar$. For larger times, the amplitude of the oscillations has the correct
order of magnitude, but the dominant frequencies are not captured. Again, the
deviations of SGW and $GW$ from the exact benchmark are comparable.

Figure 8 shows that the Hartree–Fock energy of the SGW and $GW$ agree well
with each other. For the kinetic and correlation energies, we observe
deviations similar to the density dynamics. Again, both approximations
significantly underestimate the correlations and thus overestimate the kinetic
energy of the system.

One of the main assumptions in the derivation of SGW (as well as the $GW$
approximation) was that correlations are small. It is therefore interesting to
consider the density dynamics of a system for stronger coupling, i.e.,
$U\geqslant 1.0J$, in order to estimate the range of validity. In figure 9 we
again consider a half-filled eight-site chain where the leftmost sites are
fully occupied. For $U=1.0J$, we observe significant deviations of both
approximations compared to the exact results, and the deviations start earlier
when the coupling is increased, in agreement with earlier stochastic mean
field calculations [16]. Interestingly, the overall trend of the density —
relaxation towards $1.0$ at long times — is correctly reproduced by SGW and
$GW$. Additionally, we observe increasing differences between SGW and $GW$
with increasing coupling. For $U=1.0J$, both approximations show a good
agreement for times $t\lesssim 40\hbar/J$.

Finally, let us compare SGW to $GW$ for larger systems with $20$ and $30$
sites, even though we do not have exact benchmarks available. We again
consider an initial state where the leftmost half of the sites is completely
filled and the remaining sites are empty. In figure 11, we find that for
$U=0.1J$, both approximations show a good agreement between each other for
$N=20$ as well as $N=30$, showing that for a slower dynamics in larger systems
and weak coupling the equivalence of the SGW and $GW$ holds for longer times.
The agreement of the density dynamics between both approximations, for
$U=1.0J$, now extends to longer times, i.e., $t\sim 40\hbar/J$ ($t\sim
20\hbar/J$), for $30$ ($20$) sites. At the same time, figure 11 shows that the
energy dynamics are in very good agreement, between both approximations, not
only for $U=0.1J$ but even for $U=1.0J$. Likewise, in figures 13 and 13 the
density and energy dynamics of SGW and $GW$ are compared for a thirty-site
chain with $U=1.0J$ where in the initial state only the first five sites are
fully occupied. Up to a time of $t\sim 50\hbar/J$, both methods are in good
agreement except for a short period around $t\sim 35\hbar/J$ where $GW$ shows
an unphysical negative density. Further, at $t\sim 55\hbar/J$, the $GW$
dynamics become unstable which is not reproduced by SGW. Instability is a well
known issue of the time-linear approach. In order to fix it, in general one
has to resort to the numerically expensive procedure of purification [9].
Therefore, the increased numerical stability of SGW compared to $GW$ is a huge
advantage of the stoachastic approach.

After comparing SGW and $GW$ for systems far from equilibrium, we now inquire
how well ground state properties agree with each other. To this end, we again
consider a half-filled eight-site chain. We consider half filling, i.e., each
site is singly occupied. Using the adiabatic switching method, we calculate
the ground state energies for varying on-site interaction strength for SGW and
$GW$. As can be seen in figure 14, the HF energies for SGW, $GW$ agree with
exact diagonalization results over the entire $U$ range. Slightly worse is the
agreement of the total energies. On the other hand, kinetic and correlation
energies are found to be in good agreement with the benchmark up to $U/J\sim
1.25$. As observed before, cf. figure 6, $GW$ displays slightly smaller
deviations from the exact result than SGW.

##  8 Discussion and outlook

In this paper we have developed a fluctuation approach to quantum many-
particle systems in and out of equilibrium. Our approach was motivated by the
classical fluctuation approach of Klimontovich [12]. Indeed, our theory
contains Klimontovich’s results as a limiting case. In particular, the
classical limit of the two-particle fluctuations, cf. equation (3.8), is
directly transformed into the correlation function of the fluctuations of the
microscopic phase space density $N(x)$[12],

$\displaystyle\gamma_{ijkl}$ $\displaystyle\equiv\mathcal{G}_{ijkl}\pm
G^{>}_{il}G^{<}_{jk}$ $\displaystyle\to\mathcal{G}^{\rm
cl}_{ijkl}+\delta_{il}\,n_{jk}$ $\displaystyle\to
n^{2}g_{2}(x,x^{\prime})+\delta(x-x^{\prime})nf(x)\equiv\langle\delta
N(x)\delta N(x^{\prime})\rangle,$

where, in the second line we used the (anti-)commutation rules of the field
operators, $G^{>}_{ij}=\frac{1}{\mathrm{i}\hbar}(\delta_{ij}\pm
n_{ij})\approx\frac{1}{\mathrm{i}\hbar}\delta_{ij}$, and, in the third line,
transformed from a general basis to the Wigner representation with
$x=\textbf{r},\textbf{p}$. The second terms in the expression for the
fluctuations are the respective source terms, $\gamma^{s}$, and their
classical limit. A similar correspondence holds between the different
approximations. In particular, our result for $\gamma^{\rm GW}$, which is
equivalent to the $GW$ approximation in the classical limit, goes over to the
non-Markovian Balescu–Lenard equation. While in the classical case this
equation has never been solved numerically, in this paper we included
extensive numerical results for the more general quantum case.

On the other hand, our work has been motivated by the stochastic mean field
concept of Ayik and Lacroix. Our theory was formulated at the most fundamental
level — fluctuations of the single-particle Green’s functions operators
$\delta\hat{G}$. This allowed us to establish the connections between
fluctuations and correlations and to derive the first equations of the
fluctuations hierarchy. The equations of this hierarchy are formally simpler
than the BBGKY hierarchy, containing fewer terms. However, a number of
physical effects, such as exchange diagrams and strong coupling contributions,
become entangled with other terms in a non-trivial manner. The reason is that
the two-particle fluctuations, $\gamma$, contain contributions from both two-
particle correlations, $\mathcal{G}$, and particle-hole fluctuations,
$\gamma^{s}$, cf. equation (3.8).

In order to decouple the fluctuation hierarchy, we considered various
approximations, in particular those that were introduced before, for classical
systems, such as the approximations of first and second moments [12]. While
the approximation of first moments — corresponding to the time-dependent
Hartree (quantum Vlasov) approximation — has a clear range of applicability,
the approximation of second moments does not. We, therefore, investigated
additional approximations and worked out which of them correspond to the known
many-body approximations of the BBGKY or NEGF formalisms. Of particular
interest was to identify the $GW$ approximation with and without exchange
corrections because it is a cornerstone for many-body theory for weakly and
moderately coupled systems, both in the ground state and out of equilibrium.
Here, we were able to identify the polarization approximation,
$\gamma\to\gamma^{\rm PA}$, and the $GW$ approximation, $\gamma\to\gamma^{\rm
GW}$, as the proper equivalent to the $GW$-G1-G2 equations with and without
exchange, respectively. Both of them are valid for weak to moderate coupling
such that the residual term $R^{(\mathcal{G})}$ could be neglected.

While the effort to solve the fluctuation equations (propagate $G_{1}$ and
$\gamma$) is comparable to that of the G1-G2-scheme, the fluctuation approach
has a large potential advantage: the computationally costly propagation of
two-particle quantity $\gamma$ can be avoided entirely if the corresponding
equation of motion for $\delta\hat{G}$ is solved instead. Here, we were able
to derive the equation of motion for $\delta\hat{G}$ that is equivalent
[again, up to the term $R^{(\mathcal{G})}$] to the polarization approximation
and, thus, to the $GW$ approximation with and without exchange corrections for
weak and moderate coupling. Applying the concept of the semiclassical SMF
approach, we mapped this problem on the stochastic PA (SPA) and stochastic
$GW$ (SGW) models that can be solved by averaging over a large number $K$ of
trajectories that start from different initial states. The numerical scaling
of SGW vs $GW$-G1-G2 was found to be $\mathcal{O}(KN_{b}^{4}N_{t})$ vs
$\mathcal{O}(N_{b}^{6}N_{t})$ which will be advantageous if the number of
samples that is required is below $N_{b}^{2}$, i.e., for large systems. This
is exactly the case where the G1-G2-scheme becomes computationally
challenging, so the fluctuation approach gives access to situations that are
presently out of reach with other methods.

We tested the SGW model for small Hubbard clusters (here SPA and SGW coincide)
in the ground state and far from equilibrium. Moreover, we explored a variety
of stochastic sampling methods including Gaussian sampling, the four point
approach and a deterministic method. Our simulations confirm that all sampling
methods that reproduce the first two moments correctly are equivalent, in the
framework of the SGW approximation. Therefore, there is a particularly large
freedom to minimize the numerical effort by an appropriate choice of the
sampling method. The deterministic sampling method provides advantages
especially for systems which are away from half filling. On the other hand,
stochastic sampling is favorable for large systems (100 sites or larger), due
to the weak dependence of the number of samples $K$ on the system size
$N_{b}$. In that case, the advantage of the stochastic approach, that scales
in the Hubbard case as $\mathcal{O}(KN_{b}^{2}N_{t})$, over the $GW$-G1-G2
scheme [$\mathcal{O}(N_{b}^{4}N_{t})$] is maximal.

Our numerical SGW results showed excellent agreement with $GW$-G1-G2 results,
and both, overall, agree well with the exact results, where available, within
the applicability range of the approximation, i.e., for moderate coupling,
$U/J\lesssim 1$. This confirms that the fluctuation approach constitutes an
interesting and efficient alternative to NEGF or density operator methods,
provided the relevant physical approximation has been identified. To extend
the applicability range of the fluctuation approach in the future, it would be
particularly interesting to identify the strong coupling ($T$-matrix)
approximation, e.g., [47], as well as the dynamically screened ladder
approximation [9].

Finally, three fundamental issues concerning the limitations of our
fluctuation approach should be addressed. The first issue is the validity
range of replacing the quantum-mechanical average by a semiclassical mean, cf.
equation (5.2). This procedure is only exact if there are no quantum coherence
effects between different samples. This is the case if the individual samples
are selected for an ideal system, as was done in our simulations, or in a
basis of natural orbitals. In that case, the density matrix is diagonal and no
quantum coherence effects occur. The second issue is the use of moments of the
probability distribution that correspond to a non-interacting state. However,
this is not an approximation if correlations are subsequently selfconsistently
created via the adiabatic switching procedure, as in our approach. The third
issue is the use of an ensemble of fluctuations of the initial state alone
without any update of the ensemble at later times. Even if the density matrix
is diagonal in the initial state, non-diagonal elements (quantum coherences)
may develop in the course of the dynamics. This is an issue in the standard
stochastic mean field approach, e.g., [16] where an ensemble of realizations
of the single-particle Green’s function operator, $\hat{G}_{1}$, is propagated
and added semiclassically. However, our SGW approach is different: instead of
many realizations of $G_{1}$ we propagate its expectation value quantum-
mechanically. The main difference to the G1-G2 scheme lies in the treatment of
the collision term $I$. While in the G1-G2 scheme, $I$ is computed from a
trace over the quantum two-particle correlation function $\mathcal{G}$, cf.
equation (2.12), in our fluctuations approach we compute the trace over
$\overline{\Delta G^{\mathrm{SGW},\lambda}_{pk}\Delta
G^{\mathrm{SGW},\lambda}_{lj}}$, cf. equation (5.28). This latter expression,
in general, does not capture all quantum coherence effects that build up
during the evolution. However, the computation of the trace gives rise to
dephasing effects that tend to reduce the differences to the exact result.

Fortunately, we are able to perform the ultimate test of the relevance of
these three issues. In fact, our extensive numerical results that compare SGW
and $GW$-G1-G2 simulations are an independent strong confirmation that no
systematic errors are introduced by our stochastic approach for the systems
and parameters under consideration. This gives strong confirmation for the
power of the present fluctuation approach which should be very valuable, in
particular for studying the dynamics of large quantum systems.

## Acknowledgements

We acknowledge helpful comments from Stefan Kehrein and Niclas Schlünzen. This
work is supported by the Deutsche Forschungsgemeinschaft via project
BO1366/16.

##  Appendix A Three-particle correlations and fluctuations

The relation between three-particle correlations, $\mathcal{G}^{(3)}$ and
fluctuations, $\Gamma$, can be found in similar fashion as for two-particle
correlations and fluctuations in equation (3.8). First, the full three-
particle Green’s function is expressed in terms of fluctuations (all
quantities depend on a single time),

$\displaystyle G^{(3),<}_{ijklmn}=$
$\displaystyle\pm\frac{1}{(\mathrm{i}\hbar)^{3}}\langle\hat{c}^{\dagger}_{l}\hat{c}^{\dagger}_{m}\hat{c}^{\dagger}_{n}\hat{c}_{k}\hat{c}_{j}\hat{c}_{i}\rangle=\pm\frac{1}{(\mathrm{i}\hbar)^{3}}\left\\{\pm\langle\hat{c}^{\dagger}_{l}\hat{c}^{\dagger}_{m}\hat{c}^{\dagger}_{n}\hat{c}_{i}\hat{c}_{j}\hat{c}_{k}\rangle\right\\}$
$\displaystyle=$
$\displaystyle\pm\frac{1}{(\mathrm{i}\hbar)^{3}}\left\\{\pm\langle\hat{c}^{\dagger}_{l}\hat{c}^{\dagger}_{m}(\mp\delta_{in}\pm\hat{c}_{i}\hat{c}^{\dagger}_{n})\hat{c}_{j}\hat{c}_{k}\rangle\right\\}$
$\displaystyle=$
$\displaystyle\pm\frac{1}{(\mathrm{i}\hbar)^{3}}\Big{\\{}-\delta_{in}\langle\hat{c}^{\dagger}_{l}(\mp\delta_{jm}\pm\hat{c}_{j}\hat{c}^{\dagger}_{m})\hat{c}_{k}\rangle+\langle\hat{c}^{\dagger}_{l}(\mp\delta_{im}\pm\hat{c}_{i}\hat{c}^{\dagger}_{m})(\mp\delta_{jn}\pm\hat{c}_{j}\hat{c}^{\dagger}_{n})\hat{c}_{k}\rangle\Big{\\}}$
$\displaystyle=$
$\displaystyle\pm\frac{1}{(\mathrm{i}\hbar)^{2}}\delta_{jm}\delta_{in}G^{<}_{kl}-\frac{1}{\mathrm{i}\hbar}\delta_{in}(G^{<}_{jl}G^{<}_{km}+\gamma_{jklm})+\frac{1}{(\mathrm{i}\hbar)^{2}}\delta_{im}\delta_{jn}G^{<}_{kl}\mp\frac{1}{\mathrm{i}\hbar}\delta_{im}(G^{<}_{jl}G^{<}_{kn}+\gamma_{jkln})$
$\displaystyle\mp\frac{1}{\mathrm{i}\hbar}\delta_{jn}(G^{<}_{il}G^{<}_{km}+\gamma_{iklm})+G^{<}_{il}G^{<}_{jm}G^{<}_{kn}+G^{<}_{il}\gamma_{jkmn}+G^{<}_{jm}\gamma_{ikln}+G^{<}_{kn}\gamma_{ijlm}+\Gamma_{ijklmn}$
$\displaystyle=$ $\displaystyle\pm G^{>}_{in}G^{>}_{jm}G^{<}_{kl}\mp
G^{>}_{in}G^{<}_{kl}G^{<}_{jm}\mp G^{>}_{jm}G^{<}_{in}G^{<}_{kl}$
$\displaystyle\pm
G^{<}_{in}G^{<}_{kl}G^{<}_{jm}-G^{>}_{in}G^{<}_{jl}G^{<}_{km}+G^{<}_{in}G^{<}_{jl}G^{<}_{km}$
$\displaystyle+G^{>}_{jn}G^{>}_{im}G^{<}_{kl}-G^{>}_{jn}G^{<}_{im}G^{<}_{kl}-G^{>}_{im}G^{<}_{jn}G^{<}_{kl}$
$\displaystyle+G^{<}_{kl}G^{<}_{im}G^{<}_{jn}\mp
G^{>}_{im}G^{<}_{jl}G^{<}_{kn}\pm G^{<}_{im}G^{<}_{jl}G^{<}_{kn}$
$\displaystyle\mp G^{>}_{jn}G^{<}_{il}G^{<}_{km}\pm
G^{<}_{jn}G^{<}_{il}G^{<}_{km}+G^{<}_{il}G^{<}_{jm}G^{<}_{kn}$
$\displaystyle-G^{>}_{in}\gamma_{jklm}+G^{<}_{in}\gamma_{jklm}\mp
G^{>}_{im}\gamma_{jkln}$ $\displaystyle\pm G^{<}_{im}\gamma_{jkln}\mp
G^{>}_{jn}\gamma_{iklm}\pm G^{<}_{jn}\gamma_{iklm}$
$\displaystyle+G^{<}_{im}\gamma_{jkln}+G^{<}_{jm}\gamma_{ikln}+G^{<}_{kn}\gamma_{ijlm}$
$\displaystyle+\Gamma_{ijklmn}.$ (A.1)

On the other hand, the three-particle Green’s function is expressed via the
cluster expansion by

$\displaystyle G^{(3),<}_{ijklmn}=$ $\displaystyle
G^{<}_{il}G^{<}_{jm}G^{<}_{kn}+G^{<}_{im}G^{<}_{jn}G^{<}_{kl}+G^{<}_{in}G^{<}_{jl}G^{<}_{km}$
$\displaystyle\pm G^{<}_{im}G^{<}_{jl}G^{<}_{kn}\pm
G^{<}_{in}G^{<}_{jm}G^{<}_{kl}\pm G^{<}_{il}G^{<}_{jn}G^{<}_{km}$
$\displaystyle+G^{<}_{il}\mathcal{G}_{jkmn}\pm G^{<}_{im}\mathcal{G}_{jkln}\pm
G^{<}_{in}\mathcal{G}_{jkml}$ $\displaystyle+G^{<}_{jm}\mathcal{G}_{ikln}\pm
G^{<}_{jl}\mathcal{G}_{ikmn}\pm G^{<}_{jn}\mathcal{G}_{iklm}$
$\displaystyle+G^{<}_{kn}\mathcal{G}_{ijlm}\pm G^{<}_{kl}\mathcal{G}_{ijnm}\pm
G^{<}_{km}\mathcal{G}_{ijln}$ $\displaystyle+\mathcal{G}^{(3)}_{ijklmn}.$
(A.2)

Comparing equations (A.1) and (A.2), we identify the general connection
between the correlated three-particle Green’s function and the three-particle
and two-particle fluctuations:

$\displaystyle\mathcal{G}^{(3)}_{ijklmn}=$
$\displaystyle\Gamma_{ijklmn}-G^{>}_{im}G^{>}_{jn}G^{<}_{kl}-G^{>}_{in}G^{<}_{jl}G^{<}_{km}$
$\displaystyle\mp G^{<}_{jl}\mathcal{G}_{ikmn}\mp
G^{<}_{km}\mathcal{G}_{ijln}-G^{<}_{kl}\mathcal{G}_{ijmn}$ $\displaystyle\mp
G^{>}_{im}\mathcal{G}_{jkln}\mp
G^{>}_{jn}\mathcal{G}_{iklm}-G^{>}_{in}\mathcal{G}_{jklm}$ (A.3)
$\displaystyle=$
$\displaystyle\Gamma_{ijklmn}-G^{>}_{in}G^{<}_{km}G^{<}_{jl}-G^{>}_{jn}G^{>}_{im}G^{<}_{kl}$
$\displaystyle\mp G^{>}_{in}G^{<}_{jm}G^{<}_{kl}\mp
G^{>}_{in}G^{>}_{jm}G^{<}_{kl}$ $\displaystyle\pm G^{<}_{jl}\gamma_{ikmn}\pm
G^{<}_{km}\gamma_{ijln}+G^{<}_{kl}\gamma_{ijmn}$ $\displaystyle\pm
G^{>}_{im}\gamma_{jkln}\pm G^{>}_{jn}\gamma_{iklm}+G^{>}_{in}\gamma_{jklm}.$
(A.4)

##  Appendix B Total energy conservation of the SGW approximation

The total energy of a system is given by the sum of the single-particle and
two-particle contributions [cf. equations (5.7) and (5.8)]. We start with a
representation in terms of two-particle correlations, rather than
fluctuations, so that it follows

$\displaystyle
E_{\mathrm{tot}}=\pm\mathrm{i}\hbar\sum_{ij}h_{ij}G^{<}_{ji}-\frac{\hbar^{2}}{2}\sum_{ijkl}w_{ijkl}\left(G^{\mathrm{HF},(2),<}_{klij}+\mathcal{G}_{klij}\right),$
(B.1)

where $G^{\mathrm{HF},(2),<}$ denotes the Hartree–Fock contribution to the
two-particle Green’s function [cf. equation (2.7)]. Differentiating equation
(B.1) with respect to time leads to

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}E_{\mathrm{tot}}=$
$\displaystyle\pm\sum_{ij}h_{ij}\left(\left[h^{\mathrm{HF}},G^{<}\right]_{ji}+\left[I^{(\mathcal{G})}+I^{(\mathcal{G})\dagger}\right]_{ji}\right)+\mathrm{i}\hbar\sum_{ijkl}w^{\pm}_{ijkl}G^{<}_{ki}\left(\left[h^{\mathrm{HF}},G^{<}\right]_{lj}+\left[I^{(\mathcal{G})}+I^{(\mathcal{G})\dagger}\right]_{lj}\right)$
$\displaystyle-\frac{\hbar^{2}}{2}\sum_{ijkl}w_{ijkl}\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{G}_{klij},$
(B.2)

where we used that
$w_{ijkl}G^{\mathrm{HF},(2),<}_{klij}=w^{\pm}_{ijkl}G^{<}_{ki}G^{<}_{lj}$ and
that, using the exchange symmetries of the interaction tensor, allows us to
rewrite the following derivative as

$\sum_{ijkl}w^{\pm}_{ijkl}\frac{\mathrm{d}}{\mathrm{d}t}\left(G^{<}_{ki}G^{<}_{lj}\right)=2\sum_{ijkl}w^{\pm}_{ijkl}G^{<}_{ki}\frac{\mathrm{d}}{\mathrm{d}t}G^{<}_{lj}.$
(B.3)

Obviously, total energy conservation depends on the applied approximation of
two-particle correlations (or fluctuations) and on the underlying system. For
the GW approximation, as given in section 4.4.2 [cf. equation (4.22), with
included residual term $R$, cf. equation (4.6)] and for the case of the
Fermi–Hubbard system it follows, due to spin symmetry,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}E_{\mathrm{tot}}=$
$\displaystyle\,2J\sum_{\langle
i,j\rangle}\left(\left[h^{\uparrow},G^{<,\uparrow}\right]_{ji}+\left[I^{(\mathcal{G})}+I^{(\mathcal{G})\dagger}\right]^{\uparrow}_{ji}\right)+2\mathrm{i}\hbar
U\sum_{i}G^{<,\uparrow}_{ii}\left(\left[h^{\downarrow},G^{<,\downarrow}\right]_{ii}+\left[I^{(\mathcal{G})}+I^{(\mathcal{G})\dagger}\right]^{\downarrow}_{ii}\right)$
$\displaystyle-\mathrm{i}\hbar
U\sum_{i}\left(\left[h^{(2),\uparrow\downarrow},\mathcal{G}^{\uparrow\downarrow\uparrow\downarrow}\right]_{iiii}+\Psi^{\uparrow\downarrow\uparrow\downarrow}_{iiii}+\Pi^{\uparrow\downarrow\uparrow\downarrow}_{iiii}-R^{(\mathcal{G}),\uparrow\downarrow\uparrow\downarrow}_{iii}\right).$
(B.4)

In reference [9] it was shown that the GW-G1-G2 approximation is energy
conserving so that

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}E_{\mathrm{tot}}=\mathrm{i}\hbar
U\sum_{i}R^{(\mathcal{G}),\uparrow\downarrow\uparrow\downarrow}_{iiii}.$ (B.5)

By inserting the definition of the residual term [cf. equation (4.6)] and
using that $G^{\lessgtr,\uparrow\downarrow}\equiv 0$, we find

$R^{(\mathcal{G}),\uparrow\downarrow\uparrow\downarrow}_{iiii}=G^{>,\uparrow\downarrow}_{ii}\left[I^{(\mathcal{G})}+I^{(\mathcal{G})\dagger}\right]^{\downarrow\uparrow}_{ii}+G^{<,\downarrow\uparrow}_{ii}\left[I^{(\mathcal{G})}+I^{(\mathcal{G})\dagger}\right]^{\uparrow\downarrow}_{ii}\equiv
0.$ (B.6)

Therefore, $\frac{\mathrm{d}}{\mathrm{d}t}E_{\mathrm{tot}}\equiv 0$, so the
total energy is conserved for the SGW approximation.

## References

  * [1] Bogolyubov N. N., Problems of the Dynamcial Theory in Statistical Physics, Gostekhizdat, Moscow–Leningrad, 1946, (in Russian).
  * [2] Keldysh L., Zh. Eksp. Teor. Fiz., 1964, 47, 1515–1527.
  * [3] Bonitz M., Kremp D., Scott D. C., Binder R., Kraeft W. D., Köhler H. S., J. Phys.: Condens. Matter, 1996, 8, No. 33, 6057, doi:10.1088/0953-8984/8/33/012.
  * [4] Stefanucci G., van Leeuwen R., Nonequilibrium Many-Body Theory of Quantum Systems, Cambridge University Press, New York, 2013.
  * [5] Balzer K., Bonitz M., Nonequilibrium Green’s Functions Approach to Inhomogeneous Systems, Springer, Berlin, 2013.
  * [6] Schlünzen N., Hermanns S., Bonitz M., Verdozzi C., Phys. Rev. B, 2016, 93, 035107,
doi:10.1103/PhysRevB.93.035107.

  * [7] Schlünzen N., Joost J. P., Bonitz M., Phys. Rev. Lett., 2020, 124, No. 7, 076601,
doi:10.1103/PhysRevLett.124.076601.

  * [8] Joost J. P., Schlünzen N., Bonitz M., Phys. Rev. B, 2020, 101, 245101, doi:10.1103/PhysRevB.101.245101.
  * [9] Joost J. P., Schlünzen N., Ohldag H., Bonitz M., Lackner F., Březinová I., Phys. Rev. B, 2022, 105, 165155,
doi:10.1103/PhysRevB.105.165155.

  * [10] Klimontovich Y., Zh. Eksp. Teor. Fiz., 1957, 33, 982.
  * [11] Klimontovich Y., Zh. Eksp. Teor. Fiz., 1972, 62, 1770.
  * [12] Klimontovich Y. L., Kinetic Theory of Nonideal Gases and Nonideal Plasmas, Pergamon Press, 1982.
  * [13] Bonitz M., Ebeling W., Romanovsky Yu. M., Contrib. Plasma Phys., 2003, 43, No. 5–6, 247–251,
doi:10.1002/ctpp.200310019.

  * [14] Dufty J. W., Bonitz M., J. Phys. Conf. Ser., 2005, 11, No. 1, 47, doi:10.1088/1742-6596/11/1/005.
  * [15] Ayik S., Phys. Lett. B, 2008, 658, 174–179, doi:10.1016/j.physletb.2007.09.072.
  * [16] Lacroix D., Hermanns S., Hinz C. M., Bonitz M., Phys. Rev. B, 2014, 90, 125112,
doi:10.1103/PhysRevB.90.125112.

  * [17] Lacroix D., Ayik S., Eur. Phys. J. A, 2014, 50, 95, doi:10.1140/epja/i2014-14095-8.
  * [18] Rostoker N., Nucl. Fusion, 1961, 1, 101, doi:10.1088/0029-5515/1/2/004.
  * [19] Kadomtsev B., Plasma Turbulence, Academic Press, New York, 1968.
  * [20] Sitenko A. G., Electromagnetic Fluctuations in Plasma, Academic Press, Cambridge, 1967.
  * [21] Foulkes W. M. C., Mitas L., Needs R. J., Rajagopal G., Rev. Mod. Phys., 2001, 73, 33–83,
doi:10.1103/RevModPhys.73.33.

  * [22] Dornheim T., Groth S., Bonitz M., Phys. Rep., 2018, 744, 1–86, doi:10.1016/j.physrep.2018.04.001.
  * [23] Filinov V. S., Bonitz M., Ebeling W., Fortov V. E., Plasma Phys. Controlled Fusion, 2001, 43, No. 6, 743,
doi:10.1088/0741-3335/43/6/301.

  * [24] Filinov V., Mol. Phys., 1986, 88, 1517, doi:10.1080/00268979609484533.
  * [25] Filinov V., Thomas P., Varga I., Meier T., Bonitz M., Fortov V., Koch S. W., Phys. Rev. B, 2002, 65, 165124,
doi:10.1103/PhysRevB.65.165124.

  * [26] Braams B. J., Manolopoulos D. E., J. Chem. Phys., 2006, 125, No. 12, 124105, doi:10.1063/1.2357599.
  * [27] Larder B., Gericke D. O., Richardson S., Mabey P., White T. G., Gregori G., Sci. Adv., 2019, 5, No. 11,
doi:10.1126/sciadv.aaw1634.

  * [28] Polkovnikov A., Ann. Phys., 2010, 325, No. 8, 1790–1852, doi:10.1016/j.aop.2010.02.006.
  * [29] Polkovnikov A., Phys. Rev. A, 2003, 68, 053604, doi:10.1103/PhysRevA.68.053604.
  * [30] Blakie P. B., Bradley A. S., Davis M. J., Ballagh R. J., Gardiner C. W., Adv. Phys., 2008, 57, No. 5, 363–455,
doi:10.1080/00018730802564254.

  * [31] Martin A. D., Ruostekoski J., New J. Phys., 2010, 12, No. 5, 055018, doi:10.1088/1367-2630/12/5/055018.
  * [32] Pinna D., Ryan C. A., Ohki T., Kent A. D., Phys. Rev. B, 2016, 93, 184412, doi:10.1103/PhysRevB.93.184412.
  * [33] Sundar B., Wang K. C., Hazzard K. R. A., Phys. Rev. A, 2019, 99, 043627, doi:10.1103/PhysRevA.99.043627.
  * [34] Kunimi M., Nagao K., Goto S., Danshita I., Phys. Rev. Res., 2021, 3, 013060,
doi:10.1103/PhysRevResearch.3.013060.

  * [35] Pucci L., Roy A., Kastner M., Phys. Rev. B, 2016, 93, 174302, doi:10.1103/PhysRevB.93.174302.
  * [36] Wetterich C., Phys. Rev. Lett., 1997, 78, 3598–3601, doi:10.1103/PhysRevLett.78.3598.
  * [37] Martin P. C., Schwinger J., Phys. Rev., 1959, 115, No. 6, 1342–1373, doi:10.1103/PhysRev.115.1342.
  * [38] Hermanns S., Schlünzen N., Bonitz M., Phys. Rev. B, 2014, 90, 125111, doi:10.1103/PhysRevB.90.125111.
  * [39] Bonitz M., Quantum Kinetic Theory, Teubner-Texte zur Physik, 2 edn., Springer, Cham, 2016.
  * [40] Schlünzen N., Hermanns S., Scharnke M., Bonitz M., J. Phys.: Condens. Matter, 2020, 32, 103001,
doi:10.1088/1361-648X/ab2d32.

  * [41] Valdemoro C., Alcoba D. R., Tel L. M., Pérez-Romero E., Int. J. Quantum Chem., 2009, 109, No. 11, 2622–2638,
doi:10.1002/qua.22042.

  * [42] Kubo R., J. Phys. Soc. Jpn., 1964, 19, No. 11, 2127–2139, doi:10.1143/JPSJ.19.2127.
  * [43] Lacroix D., Tanimura Y., Ayik S., Yilmaz B., Eur. Phys. J. A, 2016, 52, 94, doi:10.1140/epja/i2016-16094-1.
  * [44] Ulgen I., Yilmaz B., Lacroix D., Phys. Rev. C, 2019, 100, 054603, doi:10.1103/PhysRevC.100.054603.
  * [45] Yilmaz B., Lacroix D., Curebal R., Phys. Rev. C, 2014, 90, 054617, doi:10.1103/PhysRevC.90.054617.
  * [46] Schlünzen N., Joost J. P., Heidrich-Meisner F., Bonitz M., Phys. Rev. B, 2017, 95, 165139,
doi:10.1103/PhysRevB.95.165139.

  * [47] Schlünzen N., Bonitz M., Contrib. Plasma Phys., 2016, 56, No. 1, 5–91, doi:10.1002/ctpp.201610003.

Ukrainian 0

Квантово-флуктуаційний підхід до нерівноважного _GW_ наближення

Е. Шрьоедтер, Я. -Ф. Йоост, М. Боніц

Інститут теоретичної фізики та астрофізики, Кільський університет ім.
Крістіана Альбрехта, D-24098 Кіль, Німеччина