Quantum fluctuations approach to the nonequilibrium $GW$ approximation

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


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 (the number of time steps). Only recently, linear scaling with could be achieved within the G1-G2 scheme [7,8] which could be demonstrated even for advanced self-energies including the and the -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]. 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.

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.

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 ( = ) and use a simplified notation. For example, the single-particle correlation function will be denoted by < ≡ < ( ) ≡ < ( , ), 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 (ˆ †) and annihilation (ˆ) operators and a single-particle orbital basis {|1 , . . . , | } of the underlying single-particle Hilbert space H , which induces the so-called Fock space F and allows for the treatment of states with a varying number of particles. Basic properties of the mentioned operators are the following : Here, a generic Hamiltonian is used to describe a quantum many-particle system, which is given bŷ where the single-particle contribution is defined by ℎ |ˆ+ˆ| . The operatorˆdescribes the kinetic energy andˆan external potential that may be time dependent. For the second term we defined |ˆ| , whereˆdescribes two-particle interactions. Note thatˆmay also be time-dependent, e.g., when we compute a correlated initial state from an uncorrelated state, via the adiabatic switching procedure [9,38].
The starting point for NEGF theory is the one-body Green's function that depends on two time arguments situated on the Keldysh contour C which is defined as where T 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., [4,5] but will not be needed here. For us, it will be sufficient to consider the correlation functions ≷ for real time arguments on the time-diagonal ( = ). The definition of these functions is again given by an ensemble average,

23401-3
averaging. The corresponding Green's function operators will be denotedˆ≷ 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, (2),< − 1 ℏ 2 ˆ †ˆ †ˆˆ . (2.6) This quantity deviates from the standard two-particle density matrix by an additional prefactor of (−1/ℏ 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, (2) ,< = < < ± < < + G , where the first two terms are the Hartree and Fock contributions, whereas 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, (2.8)

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 iℏ < = ℎ HF , < + ( G) + ( G) † , (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 where we introduced an (anti-)symmetrized interaction tensor ± ± . (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 ( G) ±iℏ ∑︁ G .
(2.12) Equation (2.9) is the first equation of the G1-G2 scheme and requires the solution of an EOM for 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 For two single-particle quantities and , we define the the commutator as [ , ] { − }.
In what follows, a central role is played by the approximation of many-particle physics. Its single-time version is recovered in the G1-G2-scheme, equation (2.13), if the three-particle correlations, (3) , the ladder term and the particle-hole -matrix contribution, Π ± , are neglected. The remaining polarization term Π ± then corresponds to with exchange corrections which will be denoted ± .
On the other hand, standard (without exchange diagrams) follows by replacing Ψ ± → Ψ and Π ± → Π, which means that, in these terms, the replacement ± → is made [8,9,39]. In the present paper we present stochastic versions of and ± .

Quantum fluctuations of the Green's functions
Instead of considering the correlations via the correlated part of the two-particle Green's function, 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], (r, p, ) ≡ (r, p, ) − (r, p, ), where is the density and 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 → < and →ˆ<. Consequently, the defining quantity is the single-particle fluctuation operator, where, obviously, < = ˆ< and ˆ = 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, 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 operatorQ is defined as the projector on the subspace orthogonal to the -particle state |Ψ , whereÎ denotes the identity operator. It is straightforward to express our two-particle fluctuations in terms of the projectorQ: This alternative expression shows that fluctuations describe correlation effects of 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).

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)], 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 twoparticle Green's function can be expressed in terms of two-particle fluctuations [41], Now, using relation (2.7), we establish the relation between two-particle correlations and fluctuations: 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", This source term can be interpreted as simultaneous excitation of a particle-hole pair: a particle is annihilated in orbital and created in orbital whereas a hole is annihilated in orbital and created in . 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 . In the latter case, G = 0 and the They arise from vacuum fluctuations. In a relativistic system, the same source term would describe particle-antiparticle pair excitations.

23401-6
total fluctuation coincides with . 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 and wave number may excite an electron from a state (below the Fermi energy) to a state + above it. The probability of this process is proportional to (1 − + ), ( is the momentum distribution), in agreement with the definition (3.9).
By construction, the two-particle correlation function 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 particles , Lastly, we find that two-particle fluctuations obey the following exchange symmetry, (3.11)

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, where we introduced an effective single-particle Hartree Hamiltonian, and a collision term, (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., G and , respectively. To close equation (3.12) we now need an equation for . 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ˆ< and < from each other, where the EOM forˆ< is given by 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: (3.17) Note that, in the classical case, the fluctuations obey the exchange symmetry.

23401-7
where we defined the operator of an effective single-particle Hartree potential induced by fluctuations, Therefore, the corresponding commutator represents the interaction of < 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 Thus, the last two terms in equation (3.17) can be expressed as 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 , we consider equations for andˆ. Remarkably, by using just the EOM ofˆ[cf. equation (3.17)], it is possible to derive the EOM of any -particle fluctuation, which makes this equation the cornerstone of our formalism.

Dynamics of two-particle fluctuations
To close equation (3.12), we need an EOM for . It follows directly from equation (3.17) which is multiplyied by another factorˆand by a subsequent ensemble average : 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, which is similar to the polarization term Π ± 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 instead of ± ) and also misses the particle-hole ladder term. Using the property (3.10), it is possible to rewrite such that < is eliminated: Finally, the last term on the r.h.s. in equation (3.21), couples to three-particle fluctuations, For two two-particle quantities and , we define the commutator as [ , ] − .

23401-8
Quantum fluctuations approach to the nonequilibrium GW approximation which is again similar to the corresponding term (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 Γ is known. Interestingly, equation (3.21) does not explicitly contain exchange, second-order Born and -matrix contributions giving rise to a simpler formal structure. Those missing terms are, of course, "hidden" in the three-particle fluctuations Γ. The relation between and G is further illustrated in section 4 and requires further investigation of the three-particle contributions.

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.

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 → 0, i.e., the approximation of first moments, we find a closed equation for < , which is a mean-field (Hartree or nonlinear Vlasov) equation, 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.

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, G, are negligible, compared to the source contributions, i.e., ≈ ± . Inserting this expression into equation (3.14) leads to the well-known Hartree-Fock approximation (TDHF), 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 does not contribute to the collision integral at all. This difference is not surprising because exchange (Fock term) is a pure quantum effect.

Approximation of second moments (2M)
The approximation of second moments follows with setting Γ → 0, with the result 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)],

23401-9
Here, Π denotes the polarization contribution (2.18), whereas Ψ is the inhomogeneity, equation (2.15), in the EOM of G. Note that the superscript "±" is missing here which indicates that no (anti-)symmetrized potential appears, i.e., ± → . 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 Ψ and Π terms in equation (4.3) indicates that the 2M-approximation already contains parts of the standard contributions, as given within the G1-G2 scheme (for the correspondence to the 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 where ( ) 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): where the residual term ( G) follows from ( ) by replacing → 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 approximation, is to add the equations for 2M and 2M, , taking into account relation (3.8) and denoting where we defined (2),F ℎ (2) ,HF − ℎ (2),H . We conclude that the EOM (4.3) for 2M is missing exchange contributions to the effective two-particle Hartree-Fock Hamiltonian and terms that compensate for ( G) in order to be fully equivalent to the 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 approximation. In the next paragraph, we identify the terms that are missing so far, to reproduce the and the ± approximations.

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 approximation [39].

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: The task is now to identify the relevant approximation for the three-particle term (Γ) in the equation for , 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 ( G) , Next, we return to the full equation for the two-particle fluctuation, equation (3.21), and establish a relation between three-particle correlations, G (3) , and fluctuations, Γ, which is analogous to equation (3.8): 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, > > < = > = > , and similar for the second term. Therefore, various terms containing two-particle correlations can be rewritten in terms of fluctuations, i.e., > G or > G . This leads to four possibilities to approximate Γ, but only two of these will be considered here: Inserting equations (4.11) and (4.12) into equation (2.19) leads to the following approximation for the three-particle term, (4.14) The commutator, (2) ,F , , 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, Further, the term F in equation (4.14) is given by where ± is given by the replacement → ± . This term thus denotes exchange contributions to the polarization term which is also present in the G1-G2 equations, cf. the term Π ± , equation (2.18).

Polarization and approximations for
Collecting all terms together, we can write down the EOM for : 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 G is found by adding the equations for PA and PA,s , equation (4.9), and introducing G PA PA ∓ PA, : where Ψ ± and Π ± are defined in equations (2.15) and (2.18). This equation can be considered the fluctuation equivalent of the approximation with exchange corrections included, i.e., to ± . Thus, we conclude that the polarization approximation, PA , is equivalent to ± , 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 Π ± → ± , 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 approximation without exchange, we now look for modifications to equation (4.17) which do not contain the exchange terms, i.e., Ψ ± → Ψ and Π ± → Π. Such an equation for can be indeed derived using condition iii), |G| | ( ) |, and additionally neglecting the exchange contributions to the polarization term, This leads to the following EOM for the two-particle fluctuation which we denote by GW : which, together with equation (4.9) for the source term, corresponds to the familiar -G2 equation The equivalence of the approximations of the fluctuation approach, equation (4.21) and the G1-G2 scheme on the 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 ( ) [cf. equation (4.6)] and from three-particle fluctuations, Γ.

Polarization and
approximations expressed in terms of the single-particle fluctuationÊ quivalently 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 singleparticle fluctuations needs to be zero and thus also their time derivative, i.e., ˆ = ˆ = 0. This means that any approximation of [ˆH,ˆ] needs to also properly account for the collision integrals, ( ) + ( ) † , 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ˆ, where we take into account the linearity of equation (3.8) 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.,Ĝ ≡ 0, This leads to the following EOM forˆP A , We now verify that this equation is equivalent to EOM for PA , by differentiating the definition which contains the following terms, (4.28) We have thus proven that approximation (4.26) is equivalent to the quantum polarization approximation for PA , equation (4.17). Recall that, to recover the polarization ( ± ) approximation for G, additionally, in the equation for the term should be neglected, giving rise to PA,s , equation (4.9).
In the same way, the standard approximation can also be recovered. This is achieved by assuminĝ ≈ ± <ˆ, leading to the EOM forˆG W , (4.29) The two equations (4.26) and (4.29) are remarkable results. They contain the complex physics corresponding to the nonequilibrium approximation (with and without exchange, respectively) in an equation for a single-particle quantityˆ. 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 and ± , 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.

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 " "), and of the quantum mechanical expectation value by an arithmetic mean,ˆ− → , In practice, random realizations " " 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 singleparticle 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,ˆ,ˆ, should be symmetrized, i.e.,ˆ= Now, the commutator can be evaluated, and the symmetric expression (anti-commutator) can be replaced with the random ensemble, i.e., [ˆ,ˆ] + /2 → . 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˜ and find, on the time-diagonal, Similarly to fluctuation operators, semiclassical fluctuations are defined as deviations from the average, We now apply this approach to equation (3.15) and solve an EOM for <, or construct a system of differential equations for semiclassical fluctuations, Δ <, , and˜<. Analogously to the quantum fluctuation approach, one can now derive an EOM for Δ Δ that couples to semiclassical threeparticle fluctuations, therefore leading to a hierarchy of coupled equations [43]. Thus, the solution of an EOM for <, 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].

Observables
The SMF approach also allows one to calculate expectation values of observables. For an arbitrary one-body observable (1) , the expectation value can be computed as For a two-body observable (2) , we first have to consider the expectation value in the framework of , (2) ,< = −ℏ 2 ∑︁ (2) The expectation value of (2) is then computed within the SMF approach by the symmetrized analog of equation (5.8) This can be generalized to an arbitrary -particle observable.

Probability distribution
To construct the probability distribution describing the initial state at = 0 , we consider an ideal (uncorrelated) state, i.e., G( 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 with ±iℏ < ( 0 ) . Applying the SMF approach to equations (5.10) and (5.11), we obtain 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 where we defined 1 − . This means, in particular, that |Δ | 2 = 0 and thus Δ ≡ 0, for all realizations. However, the third and fourth moments are nonzero for expectation values that contain terms with matching indices, i.e., Δ . 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.

23401-15
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.

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 ( )] where we used the definitions, 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 Δ 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 <, which is itself symmetric: where we introduced the definition The EOM for Δ <, becomes, where we used the definitions When three-particle contributions are neglected, as in our case,˜obeys the same EOM as , 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 singleparticle EOM. This leads to very complex equations, which are beyond the scope of the present work.

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: where the following definitions were introduced, 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 . Therefore, the total numerical scaling of the SPA and SGW models is given by O ( 4 ), where is the number of basis states and is the number of time steps. It is interesting to compare this to the scaling of the approximation within the G1-G2-scheme. Therein, the CPU-time scaling is [8] of order O ( 6 ). Thus, we find that the stochastic approximations are advantageous for (up to coefficients that may be different) < 2 . The specific choice of sampling methods provides flexibility and allows for further optimization depending on the system.

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 = 0 with a finite basis.

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 1 2 . (6.1) In the original formulation of the SMF approach, a complex Gaussian distribution was assumed [44], where , 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, 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 samples, i.e., which becomes exact in the limit → ∞. Thus, to achieve the desired accuracy for the expectation value, the number should be chosen sufficiently large.

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 where we defined Δ −iℏΔ ( 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, 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.

23401-18
In what follows, we consider a fermionic system with the spin configurations ↑, ↓, though this approach can be generalized to arbitrary spins. Additionally, we assume spin symmetry of the initial state, i.e., ↑ = ↓ . Let ∈ N denote the size of the basis for one spin component and , ℎ ∈ N denote the number of occupied and unoccupied orbitals, respectively, for one spin component, with + ℎ = . Without loss of generality, we can assume ↑ = ↓ = 1 for = 1, . . . , and ↑ = ↓ = 0 for = + 1, . . . , . 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 , 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 −Δ , for every Δ , 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 Δ = 0 for ( , ) ∈ {1, . . . , } × {1, . . . , ℎ } ∪ { + 1, . . . , } × { ℎ + 1, . . . , }. The matrix representing a single sample can then be written as with A , ∈ C ℎ × . Here, it was used that Δ is self-adjoint due to the properties of the singleparticle density matrix. We can define a bĳection : Although the first two moments are exactly reproduced, higher moments are not well approximated. In fact, the -th moment can be proportional to ( −2)/2 , for > 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 O ( ). This sampling method leads to the same numerical scaling as for if we consider a system with = ℎ = /2. However, we find that for systems away from half-filling, this approach provides a significant improvement.

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 Any bĳection can be chosen here for the construction of the solution.

Hubbard Hamiltonian
For the Fermi-Hubbard model, the general pair-interaction of equation (

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 "SGW" and "<", for simplicity, e.g., <,SGW ≡ , 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 ( 0 ) is chosen such that

23401-20
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 and Δ 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 with a time dependent interaction ( ). Calculations start at with an uncorrelated ground state, with ( ) = 0. The on-site interaction is increased monotonously and rather slowly such that at 0 the system is in a fully correlated ground state with ( 0 ) = .
The observables we consider are the densities on all sites, ( ), and, in addition, the kinetic energy ( kin ), Hartree-Fock ( HF ) and correlation energies ( cor ) which follow from equations (5.7) and (5.9),

∑︁
, 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.

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 ( = 0.1 ). 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 1 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 = 1.0 . Despite the same number of samples, we see in figure 2 significant deviations in the density dynamics for times 50ℏ/ , 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 = 0.1 and = 1.0 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 = 1.0 . Here, the density deviations are slightly larger than for = 0.1 , 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 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 . Additionally, our fit with a function proportional to 1/ √ illustrates the expected trend due to the central limit theorem. Assuming 10 4 samples for each spin component and half-filling so that = ℎ = /2, stochastic sampling appears to be advantageous, especially for system sizes larger than 100 sites.

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 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.
First, we consider a half-filled eight-site chain with weak on-site interaction ( = 0.1 ). The initial state is chosen so that the leftmost half of the sites is fully occupied. It can be seen in figure 5 that the density dynamics from SGW and the -G1-G2 scheme show good agreement with each other. Both  reproduce the exact result reasonably well for times 40ℏ/ , after which SGW qualitatively displays a better agreement for times 50ℏ/ . After this point, deviations between the two approximations and the exact result start to increase, with the 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 70ℏ/ , 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 50ℏ/ , after which a revival of the dynamics can be observed for SGW and . 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 20ℏ/ , there are significant differences between the two approximations and the exact result. However, SGW and show a very good agreement between each other, for the entire time propagation. For times 70ℏ/ , 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 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 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.
Next, we consider the same system but with an increased on-site interaction of = 0.5 . In figure 7 we see that the SGW and show a very good agreement between each other for times 40ℏ/ , whereas after this time small differences show up. For times 60ℏ/ , the oscillations of the density for both approximations have a similar amplitude. At the same time, both approximations reproduce the exact dynamics only for 10ℏ. 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 from the exact benchmark are comparable. Figure 8 shows that the Hartree-Fock energy of the SGW and 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 approximation) was that correlations are small. It is therefore interesting to consider the density dynamics of a system for stronger coupling, i.e., 1.0 , in order to estimate the range of validity. In figure 9 we again consider a halffilled eight-site chain where the leftmost sites are fully occupied. For = 1.0 , 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 . Additionally, we observe increasing differences between SGW and with increasing coupling. For = 1.0 , both approximations show a good agreement for times 40ℏ/ . Finally, let us compare SGW to 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 10, we find that for = 0.1 , both approximations show a good agreement between each other for = 20 as well as = 30, showing that for a slower dynamics in larger systems and weak coupling the equivalence of the SGW and holds for longer times. The agreement of the density dynamics between both approximations, for = 1.0 , now extends to longer times, i.e., ∼ 40ℏ/ ( ∼ 20ℏ/ ), 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 = 0.1 but even for = 1.0 . Likewise, in figures 12 and 13 the density and energy dynamics of SGW and are compared for a thirty-site chain with = 1.0 where in the initial state only the first five sites are fully occupied. Up to a time of ∼ 50ℏ/ , both methods are in good agreement except for a short period around ∼ 35ℏ/ where shows an unphysical negative density. Further, at ∼ 55ℏ/ , the 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 is a huge advantage of the stoachastic approach.
After comparing SGW and 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 . As can be seen in figure 14, the HF energies for SGW, agree with exact diagonalization results over the entire 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 / ∼ 1.25. As observed before, cf. figure 6, displays slightly smaller deviations from the exact result than SGW.

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 ( ) [12], where, in the second line we used the (anti-)commutation rules of the field operators, > = 1 iℏ ( ± ) ≈ 1 iℏ , and, in the third line, transformed from a general basis to the Wigner representation with = r, p. The second terms in the expression for the fluctuations are the respective source terms, , and their classical limit. A similar correspondence holds between the different approximations. In particular, our result for GW , which is equivalent to the 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ˆ. 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, , contain contributions from both two-particle correlations, G, and particle-hole fluctuations, , 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 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, → PA , and the approximation, → GW , as the proper equivalent to the -G1-G2 equations with and without exchange, respectively. Both of them are valid for weak to moderate coupling such that the residual term ( G) could be neglected.
While the effort to solve the fluctuation equations (propagate 1 and ) 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 can be avoided entirely if the corresponding equation of motion forˆis solved instead. Here, we were able to derive the equation of motion forˆthat is equivalent [again, up to the term ( G) ] to the polarization approximation and, thus, to the 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 (SGW) models that can be solved by averaging over a large number of trajectories that start from different initial states. The numerical scaling of SGW vs -G1-G2 was found to be O ( 4 ) vs O ( 6 ) which will be advantageous if the number of samples that is required is below 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 on the system size . In that case, the advantage of the stochastic approach, that scales in the Hubbard case as O ( 2 ), over the -G1-G2 scheme [O ( 4 )] is maximal. Our numerical SGW results showed excellent agreement with -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, / 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 ( -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,ˆ1, is propagated and added semiclassically. However, our SGW approach is different: instead of many realizations of 1 we propagate its expectation value quantum-mechanically. The main difference to the G1-G2 scheme lies in the treatment of the collision term . While in the G1-G2 scheme, is computed from a trace over the quantum two-particle correlation function G, cf. equation (2.12), in our fluctuations approach we compute the trace over Δ SGW, Δ SGW, , 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 -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.
Comparing equations (A.1) and (A.2), we identify the general connection between the correlated threeparticle Green's function and the three-particle and two-particle fluctuations:

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 tot = ±iℏ where we used that HF, (2),< = ± < < and that, using the exchange symmetries of the interaction tensor, allows us to rewrite the following derivative as Obviously, total energy conservation depends on the applied approximation of two-particle correlations (or fluctuations) and on the underlying system. In reference [9] it was shown that the GW-G1-G2 approximation is energy conserving so that Therefore, d d tot ≡ 0, so the total energy is conserved for the SGW approximation.