# Modelling excitonic Mott transitions in two-dimensional semiconductors

A. Kudlis  I. Iorsh  Department of Physics and Technology, ITMO University,
St. Petersburg 197101, Russia

###### Abstract

We analyze the many-particle correlations that affect the optical properties
of two-dimensional semiconductors. These correlations manifest themselves
through the specific optical resonances such as excitons, trions, etc.
Starting from the generic electron-hole Hamiltonian and employing the
microscopic Heisenberg equation of motion the infinite hierarchy of
differential equations can be obtained. In order to decouple the system we
address the cluster expansion technique which provides a regular procedure of
consistent accounting of many-particle correlation contributions into the
interband polarization dynamics. In particular, the partially taken into
account three-particle correlations modify the behavior of absorption spectra
with the emergence of a trion-like peak additional to excitonic ones. In
contrast to many other approaches, the proposed one allows us to model the
optical response of 2d semiconductors in the regime when the Fermi energies
are of the order of the exciton and trion binding energies, thus allowing us
to rigorously model the onset of the excitonic Mott transition, the regime
being recently studied in various 2d semiconductors, such as transition metal
dichalcogenides.

##  I Introduction

Two-dimensional semiconductors, such as monolayers of transition metal
dichalcogenides (TMDs), appear to be an ideal platform for the exploration of
the excitonic complexes Wang  _et al._ (2018a). Peculiar and appealing
properties of TMDs are largely dictated by their two-dimensional nature:
suppressed screening leads to the emergence of the tightly bound
Cheiwchanchamnangij and Lambrecht (2012); Komsa and Krasheninnikov (2012);
Chernikov  _et al._ (2014) and at the same time strongly interacting excitons
in these structures Ramasubramaniam (2012); He  _et al._ (2014). The former
property allows for the efficient optical probing of the exciton structure Qiu
_et al._ (2013); Berkelbach  _et al._ (2013); Ugeda  _et al._ (2014);
Berghäuser and Malic (2014); Wang  _et al._ (2018b), and latter leads to the
emergence of the pronounced many-particle correlations. Besides that, the 2d
nature of TMDs allows for the efficient doping of these structures by means of
an external gating Xu  _et al._ (2015). The most pronounced effect arising in
doped TMDs is the formation of the additional peak in photoluminescence
redshifted with respect to the excitonic one Wang  _et al._ (2015); Mak  _et
al._ (2012); Ross  _et al._ (2013); Shang  _et al._ (2015); You  _et al._
(2015).

From the theoretical side, there still exists an ambiguity in the
interpretation of the origin of this peak. While some of the researchers
describe the peak as trions, tightly bound complexes of two electrons and a
hole (or vice versa) Zhang  _et al._ (2014); Kidd  _et al._ (2016); Ganchev
_et al._ (2015); Velizhanin and Saxena (2015); Mayers  _et al._ (2015); Jones
_et al._ (2015); Vaclavkova  _et al._ (2018); Durnev and Glazov (2018), the
other part of the community exploits the approach based on the Fermi polarons:
excitons, dressed by the sea of the residual electrons Sidler  _et al._
(2016); Efimkin and MacDonald (2017); Back  _et al._ (2017); Roch  _et al._
(2019); Tan  _et al._ (2020). Although it has been shown, that the two
approaches are equivalent in the low doping limit Glazov (2020), they produce
qualitatively different results in the large doping limit Klawunn and Recati
(2011). Moreover, at elevated doping and thus electron concentration the
composite nature of the exciton should be taken into account: ultimately, as
the average separation between the electrons is of the order of the exciton
Bohr radius, the excitons can be no longer treated as bound pairs which
results in the onset of the exciton Mott transition, where the exciton gas
transforms to the electron-hole plasma. The Mott transition in TMD monolayers
(MLs) has been recently explored experimentally Yu  _et al._ (2020). The
quantitative model of high density behavior of the trions in TMDs is still yet
to appear. Moreover, the experimentally relevant regime of the excitonic Mott
transition Lin  _et al._ (2019) still lacks adequate quantitative theoretical
description.

Figure 1: Schematic representation of the procedure employed in the paper. The
key quantity within the work is the interband polarization $P(\omega)$ which
allows us to calculate the optical absorption. The former in turn can be
expressed via two-operator expectation values ($\braket{2}_{t}$) which enter
the hierarchy of differential equations. The Heisenberg equation together with
the cluster expansion technique allow to regularly truncate a such infinite
system within the chosen order of correlations
($\braket{N}_{t}^{c}=\braket{\hat{a}_{1}\dots\hat{a}_{N}\hat{a}_{N}^{\dagger}\dots\hat{a}_{1}^{\dagger}}_{t}^{c}$).
In this paper we partly take into account the dynamics of three-particle
correlations ($\braket{3}_{t}^{c}$) thereby expanding the solution given by
the semiconductor Bloch equation.

One of the methods proved powerful for the description of the excitons in the
vicinity of the Mott transition is based on the so-called cluster expansion
technique Hall (1975); Schoeller (1994); Fricke (1996). Within this approach,
we start with the many-particle Hamiltonian of electrons interacting with the
classical time-dependent electromagnetic field. This quantum problem allows
for the reformulation in terms of an infinite series of coupled differential
equations for many-particle correlations. The cluster expansion defines the
scheme of the truncation of this system to obtain the many-body dynamics with
controlled accuracy. The cluster expansion technique has been used to
quantitatively model the excitonic Mott transition in GaAs quantum wells Kira
and Koch (2006); Haug and Koch (2009).

It has however never been used to model the high-density trion dynamics due to
the computational complexity of the resulting equations and the weakness of
the trion response in GaAs quantum wells.

In this paper, we apply the cluster expansion technique for the TMD MLs to
explore the effect of the high-order many-particle correlations on the optical
response of these structures (see Fig. 1). By extending the system of
equations to capture three-particle contributions we were able to extract the
emergence of the trion peak in the optical polarizability of TMDs starting
barely from the electron-hole Hamiltonian. We believe that the results of the
paper prove that the cluster expansion methods could be extremely useful to
study the dynamical many-body correlations in TMDs and their effect on the
transient optical properties.

The paper is organized as follows. In Sec. II we describe the model and a set
of approximations adopted within this work. In Sec. III we show the approach
used and also the equation of motion (EOM) for different orders. The final
analytical expressions and the corresponding numerical results are presented
and discussed in Sec. IV. Finally, in Sec. V we will draw a conclusion.

Throughout the work we try to use indices $\bm{i},\bm{j}$ as summation ones,
while primed indices – $\bm{i}^{\prime},\bm{j}^{\prime}$ – will be used as
external ones. Also, if some quantity is time-dependent, we can demonstrate it
as $Q(t)$, $Q_{t}$ or $Q^{t}$. In this case, the corresponding Fourier
transform is denoted by $Q(\omega)$, $Q_{\omega}$ or $Q^{\omega}$
respectively.

##  II Description of the model

In what follows we limit ourselves to the case of just two bands. The
generalization to the case of the multiband structure is straightforward. We
also neglect the spin degree of freedom and consider the single valley
dynamics, thus neglecting the intervalley scattering. Finally, we do not
consider the electron-phonon interaction in this work. Thus, the model
Hamiltonian reads:

$\displaystyle\hat{H}$
$\displaystyle=\hat{H}_{0}+\hat{V},\quad\hat{V}\equiv\hat{H}_{el,p},\quad\hat{H}_{0}=\hat{H}_{el,k}+\hat{H}_{I},$
(1)

where

$\displaystyle\hat{H}_{el,k}=$
$\displaystyle\sum\limits_{\lambda,\bm{k}}E_{\lambda,\bm{k}}\hat{a}^{\dagger}_{\lambda,\bm{k}}\hat{a}_{\lambda,\bm{k}}=\hbar\sum\limits_{\bm{i}}\varepsilon_{\bm{i}}\hat{a}^{\dagger}_{\bm{i}}\hat{a}_{\bm{i}},\quad\hat{H}_{I}\approx-\sum\limits_{\bm{k}}\mathcal{E}(t)d_{cv}\left(\hat{a}^{\dagger}_{c,\bm{k}}\hat{a}_{v,\bm{k}}+\hat{a}^{\dagger}_{v,\bm{k}}\hat{a}_{c,\bm{k}}\right)=\hbar\sum\limits_{\bm{i},\bm{j}}h^{ext}_{\bm{i}\bm{j}}(t)\hat{a}^{\dagger}_{\bm{i}}\hat{a}_{\bm{j}},$
(2) $\displaystyle\hat{V}=$
$\displaystyle\frac{1}{2}\sum\limits_{\begin{subarray}{c}\,\bm{k}_{1},\bm{k}_{2}\\\
\bm{q}\neq
0\end{subarray}}V_{q}\left[\hat{a}^{\dagger}_{c,\bm{k}_{1}+\bm{q}}\hat{a}^{\dagger}_{c,\bm{k}_{2}-\bm{q}}\hat{a}_{c,\bm{k}_{2}}\hat{a}_{c,\bm{k}_{1}}+\hat{a}^{\dagger}_{v,\bm{k}_{1}+\bm{q}}\hat{a}^{\dagger}_{v,\bm{k}_{2}-\bm{q}}\hat{a}_{v,\bm{k}_{2}}\hat{a}_{v,\bm{k}_{1}}+2\hat{a}^{\dagger}_{c,\bm{k}_{1}+\bm{q}}\hat{a}^{\dagger}_{v,\bm{k}_{2}-\bm{q}}\hat{a}_{v,\bm{k}_{2}}\hat{a}_{c,\bm{k}_{1}}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{4}\hbar\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}\,\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}\hat{a}_{\bm{j}_{1}},\
\ \text{where}\quad\bm{i}=(\lambda,\bm{k}),\quad\lambda\in\\{c,v\\}.$ (3)

Indices $c$ and $v$ correspond to conduction and valence bands respectively.
For single-particle kinetic energies we turn to the parabolic isotropic
dispersion relation. $\mathcal{E}(t)$ is the electric field strength, while
$d_{cv}$ is effective dipole matrix element of interband absorption. As usual,
the creation and annihilation operators obey the following commutation
relations:

$\begin{split}&\big{[}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}},\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\big{]}_{+}=0,\quad\big{[}\hat{a}_{\bm{i}^{\prime}_{1}},\hat{a}_{\bm{i}^{\prime}_{2}}\big{]}_{+}=0,\\\
&\big{[}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}},\hat{a}_{\bm{i}^{\prime}_{1}}\big{]}_{+}=\delta_{\bm{i}^{\prime}_{1},\bm{j}^{\prime}_{1}}\equiv\delta_{\lambda_{1},\bm{k}_{1};\lambda_{2},\bm{k}_{2}}\equiv\delta_{\bm{k}_{1},\bm{k}_{2}}\delta_{\lambda_{1},\lambda_{2}}.\end{split}$
(4)

For a greater versatility for each term of the Hamiltonian one can introduce
the tensor functions $h^{ext}_{\bm{i}\bm{j}}$ and
$v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}$. The connection formulas
between them and standard notation can be found in Appendix B. Note only here
that small potential – $v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}$ – has
to satisfy the following symmetry conditions:

$\displaystyle
v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}=-v_{\bm{i}_{2},\bm{i}_{1},\bm{j}_{1},\bm{j}_{2}},\
\
v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}=-v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}}.\
\ $ (5)

It should be marked also that during the derivation of all basic expressions
we do not specify the form of these functions. Hence, any model reducible to
this structure can be treated by the approach used in this work. Moreover,
this formalism can be applied to systems with bosons.

Speaking about the TMD MLs, when it comes to specific numerical calculations,
the conventional Coulomb potential in two dimensions ($\sim 1/q$) is no longer
applicable Rösner  _et al._ (2016); Stier  _et al._ (2016); Raja  _et al._
(2017). There are two reasons of its deviation and both of them relate to
screening phenomena. The first modification is associated with the dynamical
screening caused by the presence of free carriers. The frequency and momentum
dependencies of this effect often lead to difficulties of theoretical
description. In a such situation it is necessary to exploit some
approximation. In particular, within the random-phase approximation (RPA) the
standard Coulomb potential is replaced with some effective one given by the
Lindhard formula Stern (1967); Haug and Koch (2009). Schematically, it can be
expressed as:

$\displaystyle V_{\bm{q}}\rightarrow
V_{\bm{q}}^{eff}=\frac{V_{\bm{q}}}{1-V_{\bm{q}}\Pi(\bm{q},\omega)},$ (6)

where $\Pi(\bm{q},\omega)$ is free-particle polarization, described by:

$\displaystyle\Pi(\bm{q},\omega)=\sum\limits_{\bm{k}}\frac{f_{\bm{k}-\bm{q}}-f_{\bm{k}}}{\hbar(\omega+i\delta+\varepsilon_{\bm{k}-\bm{q}}-\varepsilon_{\bm{k}})},$
(7)

where $f_{\bm{k}}$ is the Fermi-Dirac distribution function. Within the
present paper we expect that effects connected with screening caused by free
carriers are incorporated automatically by taking into account of many-
particle correlations and its artificial injection into the theory via
different effective potentials is excessive. The only thing which we assume is
that the dynamical screening is fully developed and we deal with stationary
systems.

The second type of screening effects is associated with dielectric properties
of both substrate, superstrate and environment which are always present in
realistic physical experiments. Due to the dimensional confinement in MLs any
inhomogeneity of surroundings lead to a significant distortion of the 2d
Coulomb law in the layer. By solving the Poisson equation in such compound
structure one can obtain some potential for MLs of a finite thickness Van Tuan
_et al._ (2018). Within the strict two-dimensional limit, however, a such
potential tends to the well-known Rytova-Keldysh form Rytova (1967); Keldysh
(1979); Schmitt-Rink and Ell (1985); Cudazzo  _et al._ (2011). The Fourier
transform of it is as follows:

$V_{\bm{q}}=\frac{2\pi e^{2}}{L^{2}}\frac{1}{\varepsilon(q)\,q},$ (8)

where $\varepsilon(q)=\varepsilon_{0}(1+q\,r_{0})$. The characteristics
$\varepsilon_{0}$ and $r_{0}$ can be considered as phenomenological parameters
of the theory. They correspond to the average dielectric constant of
surroundings and effective screening length. As was just mentioned, the
expression (8) was obtained taking the strict 2d confinement into account.
Such approximation is applicable if the typical radius of the bonded electron-
hole pairs exceeds the lattice constant, which is realized in the case of TMD
MLs. Moreover, the Rytova-Keldysh potential turns out to be well applicable in
order to analyze the formation of dipolar excitons and accompanied phenomena
in double layer heterostructures with using the different TMD MLs Berman and
Kezerashvili (2017).

Thus, as a final potential used in our numerical calculations we take the
expression (8).

As mentioned above, we aim at calculating the interband polarization, which
can be written as

$\bm{P}(t)=\sum\limits_{\begin{subarray}{c}\lambda,\lambda^{\prime}\\\
\bm{k}\end{subarray}}\braket{\hat{a}^{\dagger}_{\lambda,\bm{k}}\hat{a}_{\lambda^{\prime},\bm{k}}}_{t}\bm{d}_{\lambda,\lambda^{\prime}}=\sum\limits_{\begin{subarray}{c}\lambda,\lambda^{\prime}\\\
\bm{k}\end{subarray}}P_{\lambda\lambda^{\prime},\bm{k}}(t)\bm{d}_{\lambda\lambda^{\prime}}.$
(9)

Taking into account the two-band approximation adopted in this work, it is
necessary to consider only the following quantity:

$\tilde{\mathfrak{P}}_{\bm{k}}^{t}\equiv\tilde{\mathfrak{P}}(\bm{k},t)\equiv
P_{vc,\bm{k}}(t)=\braket{\hat{a}^{\dagger}_{v,\bm{k}}\hat{a}_{c,\bm{k}}}_{t}.$
(10)

However, further only due to the computational reasons we prefer to work with
the following variable:

$\mathfrak{P}_{\bm{k}}^{t}\equiv\mathfrak{P}(\bm{k},t)=\braket{\hat{a}_{c,\bm{k}}\hat{a}^{\dagger}_{v,\bm{k}}}_{t}\equiv\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}=-\tilde{\mathfrak{P}}_{\bm{k}}^{t},$
(11)

where we used compound indices: $\bm{j}^{\prime}=(v,\bm{k}),\
\bm{i}^{\prime}=(c,\bm{k})$. The transition from one variable to another is
trivial.

##  III Cluster expansion

As was previously mentioned, during the derivation of dynamical equations for
quantities of interest we are faced with the infinite system of equations. In
this section, we describe the procedure of the proper truncation which makes
the system closed. Also, we present the series of approximations for the part
which we treat numerically.

###  III.1 Equation of motion

In order to describe the time evolution of the arbitrary operator $\hat{A}$
one can address the microscopic Heisenberg EOM. In terms of the model (1) it
reads

$\hbar\frac{d}{dt}\hat{A}+i\big{[}\hat{A},\hat{H}_{el,k}\big{]}+i\big{[}\hat{A},\hat{H}_{I}\big{]}=-i\big{[}\hat{A},\hat{V}\big{]}.$
(12)

The same is true for the expectation value of operator $\hat{A}$ with an
initial statistical operator $\rho_{0}$ at an initial time:
$\braket{A}_{t}=Tr[\rho_{0}A(t)]$. Speaking about $\rho_{0}$ there are no
special restrictions on its structure.

Further, we focus on the particular type of the operator $\hat{A}$ for which
all the expressions presented in this paper are valid. The corresponding
general form can be written as

$\hat{A}=\hat{a}_{\bm{i}^{\prime}_{1}}\dots\hat{a}_{\bm{i}^{\prime}_{n}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{n}}\dots\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}.$
(13)

It is clear, that with such operator in hand in the non-interacting case
($V=0$) from (12) we obtain a closed system of dynamical equations which
contains only one type of expectation values. A completely different
situation, however, is observed if we are dealing with nonzero potential. In
particular, the dynamics of the two-operator expectation value –
$\braket{\hat{a}^{\dagger}_{\bm{i}}\hat{a}_{\bm{i}^{\prime}}}_{t}$ – is
coupled with four-operator ones. The latter, in their turn, depend already on
six-operator expectation value dynamics and so on. These steps lead us to
infinite hierarchy of differential equations. In order to solve this system,
one should find a proper way of truncation of this system. One of the options
is factorization of many-operator expectation values into the product of
dominant two-operator terms. This procedure results in the RPA that,
unfortunately, does not allow one to tackle multi-particle effects. To
overcome this problem in Fricke (1996) some approach on the basis of the
cluster expansion technique (CET) presented in Schoeller (1994) was suggested.
The main idea is to rewrite all differential equations on expectation values
in terms of correlations. By means of the CET the expectation value of the
product of an arbitrary combination of creation and annihilation operators
$\hat{b}_{i}$ can be expanded as follows:

$\begin{split}\braket{b_{1}}_{t}&=\braket{b_{1}}_{t}^{c},\\\
\braket{b_{1}b_{2}}_{t}&=\braket{b_{1}b_{2}}_{t}^{c}+\braket{b_{1}}_{t}^{c}\braket{b_{2}}_{t}^{c},\\\
\braket{b_{1}b_{2}b_{3}}_{t}&=\braket{b_{1}b_{2}b_{3}}_{t}^{c}+\braket{b_{1}b_{2}}_{t}^{c}\braket{b_{3}}_{t}^{c}+\braket{b_{2}}_{t}^{c}\braket{b_{1}b_{3}}_{t}^{c}\\\
&\quad+\braket{b_{1}}_{t}^{c}\braket{b_{2}b_{3}}_{t}^{c}+\braket{b_{1}}_{t}^{c}\braket{b_{2}}_{t}^{c}\braket{b_{3}}_{t}^{c},\\\
&\dots\ .\end{split}$ (14)

For the $n$-th order correlations the sum extends over all disjoint partitions
of the set $\\{b_{1},\dots,b_{n}\\}$. These expressions can be considered as
the definition of correlations. This means that in order to obtain the
$n$-operator correlation we have to subtract from the $n$-operator expectation
value all the lower-order correlations ($n-1$, $n-2$, $\dots$). The operators
in each correlation retain their order. The sign of each term is defined by
the number of permutations of fermionic operators in order to coincide with
the initial one.

Within this work we require the conservation of the fermionic occupation
numbers. This automatically causes that correlations for odd number of
operators vanish as well as for combinations where numbers of annihilation and
creation operators do not coincide. In particular, for two-operator
expectation values we have

$\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}=\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c},\quad\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}_{\bm{i}^{\prime}}}_{t}=\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}_{\bm{i}^{\prime}}}_{t}^{c},$
(15)

while for four-operator ones the corresponding expansion reads:

$\displaystyle\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}=\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}+\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}$
$\displaystyle\qquad\qquad\qquad\qquad-\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}.$
(16)

Through the canonical commutation relations (4) it is not difficult to see
that the correlations must meet the following symmetry conditions:

$\begin{split}\braket{\dots\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}\dots}^{c}_{t}&=-\braket{\dots\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}_{\bm{i}^{\prime}}\dots}^{c}_{t}\\\
\braket{\dots\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\dots}^{c}_{t}&=-\braket{\dots\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}_{\bm{i}^{\prime}_{1}}\dots}^{c}_{t}\\\
\braket{\dots\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\dots}^{c}_{t}&=-\braket{\dots\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\dots}^{c}_{t}.\end{split}$
(17)

Thus, by means of (15) the polarization components (11) can be expressed via
correlations in a trivial way:

$\displaystyle\mathfrak{P}_{\bm{k}}^{t}=\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c},\
\bm{i}^{\prime}=(c,\bm{k}),\ \bm{j}^{\prime}=(v,\bm{k}).$ (18)

In addition, there is another important quantity worth rewriting in terms of
correlations – particle density operators:

$\begin{split}\mathfrak{n}_{\lambda,\bm{k}}^{t}=\braket{\hat{n}_{\lambda,\bm{k}}}_{t}^{c}=&\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}_{\bm{i}^{\prime}}}_{t}^{c},\
\ \bm{i}^{\prime}=\bm{j}^{\prime}=(\lambda,\bm{k}),\\\
\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}&=1-\mathfrak{n}_{\lambda,\bm{k}}^{t}.\end{split}$
(19)

Further, we will see that taking into account some other physical assumptions
only these two functions enter the EOM. Also, within the present study we do
not consider the correlations higher than three-particle ones. In such
situation it looks very reasonable to introduce a specific variables. For two-
and three-particle correlations they are as follows:

$\displaystyle\mathfrak{D}^{t}_{\bm{i}^{\prime}_{1},\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2},\bm{j}^{\prime}_{1}}$
$\displaystyle=\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c},$
(20)
$\displaystyle\mathfrak{T}^{t}_{\bm{i}^{\prime}_{1},\bm{i}^{\prime}_{2},\bm{i}^{\prime}_{3},\bm{j}^{\prime}_{3},\bm{j}^{\prime}_{2},\bm{j}^{\prime}_{1}}$
$\displaystyle=\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}_{\bm{i}^{\prime}_{3}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{3}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c},$
(21)

where all $\bm{i}^{\prime}$ and $\bm{j}^{\prime}$ are compound indices. Here
one should make a remark regarding the features of notation. The expressions
$\braket{\hat{a}_{1}\dots\hat{a}_{N}\hat{a}^{\dagger}_{N}\dots\hat{a}^{\dagger}_{1}}_{t}$
and $\braket{2N}_{t}$ are used for the average of the $2N$-operator product,
while
$\braket{\hat{a}_{1}\dots\hat{a}_{N}\hat{a}^{\dagger}_{N}\dots\hat{a}^{\dagger}_{1}}_{t}^{c}$
and $\braket{N}_{t}^{c}$ for correlations of $N$-th order.

Having obtained an idea about correlations, let us finally figure out what is
the benefit to work with them instead of expectation values. The understanding
could be best achieved by comparing the structures of the EOM in both cases.
Based on the interaction form (3), following the notation in Kira and Koch
(2006) schematically for $2N$-operator expectation values the differential
equation can be written as follows:

$\hbar\frac{d}{dt}\braket{2N}_{t}=\tilde{T}_{N}\big{[}\braket{2N}_{t}\big{]}+\tilde{V}_{2,N}\big{[}\braket{2N+2}_{t}\big{]},$
(22)

while for $N$-particle correlations the corresponding EOM have the following
form:

$\displaystyle\hbar\frac{d}{dt}\braket{N}^{c}_{t}$
$\displaystyle=T_{N}\big{[}\braket{N}^{c}_{t}\big{]}+V_{2,N}\big{[}\braket{N+1}^{c}_{t}\big{]}$
(23)
$\displaystyle\quad+V_{1,N}\big{[}\braket{N+1}_{\braket{N}^{c}_{t},\braket{N-1}^{c}_{t},\dots,\braket{1}^{c}_{t}}\big{]}.$

From (22) one can see that in order to obtain a closed system one has to omit
the term $\tilde{V}_{2,N}$, but in this case we are totally losing the
information about contributions into dynamics from the interaction. In the
case of correlations, however, the presence of $V_{1,N}$ which contains
correlation of order no higher than $N$ allows us to construct a closed system
of differential equations without loosing the interaction information within
the given order of approximation. Thus, in order to obtain a closed system of
the EOM within an $N$-particle correlation approximation, we shall neglect all
correlations of order $N+1$ and higher. In the next section in terms of this
schematic equation we describe all the approximations for which analytical and
partly numerical results with some simplification are presented in this paper.

###  III.2 Series of approximations

####  III.2.1 Free system dynamics

First, for the completeness of the study we reproduce the expression for free-
particle polarization. The corresponding schematic equation by means of (23)
reads

$\hbar\frac{d}{dt}\braket{1}^{c}_{t}=T_{1}\big{[}\braket{1}^{c}_{t}\big{]},$
(24)

where we omit all the contributions connected with Coulomb interaction.

####  III.2.2 One-particle dynamics

Following the adopted strategy, within the one-particle correlation
approximation the EOM looks like

$\hbar\frac{d}{dt}\braket{1}^{c}_{t}=T_{1}\big{[}\braket{1}^{c}_{t}\big{]}+V_{1,1}\big{[}\braket{2}_{\braket{1}^{c}_{t}}\big{]}.$
(25)

The notation $\braket{2}_{\braket{1}^{c}_{t}}$ means that two-particle
correlations enter the corresponding equation only via the product of one-
particle correlations. Further, based on this expression, the well-known
semiconductor Bloch equation is restored.

####  III.2.3 Two-particle dynamics via one-particle correlations

As the next step we include into consideration two-particle correlation
dynamics, however, within this approximation only one-particle contributions
enter the corresponding EOM. These corrections can be associated with
scattering processes. Let us note that the EOM for one-particle correlations
is already exact. Recapitulating what was said above we come to the following
system:

$\begin{split}&\hbar\frac{d}{dt}\braket{1}^{c}_{t}=T_{1}\big{[}\braket{1}^{c}_{t}\big{]}+V_{1,1}\big{[}\braket{2}_{\braket{1}^{c}_{t}}\big{]}+V_{2,1}\big{[}\braket{2}^{c}_{t}\big{]},\\\
&\hbar\frac{d}{dt}\braket{2}^{c}_{t}=K\left[T_{2}\big{[}\braket{2}^{c}_{t}\big{]}\right]+V_{1,2}\big{[}\braket{3}_{\braket{1}^{c}_{t}}\big{]}.\end{split}$
(26)

where operator $K$ discards from $T_{2}$ all the explicitly field-dependent
parts. In terms of two-particle correlations it formally can be expressed as:

$K\left[\braket{\big{[}A,H_{el,k}\big{]}}_{t}+\braket{\big{[}A,H_{I}\big{]}}_{t}\right]=\braket{\big{[}A,H_{el,k}\big{]}}_{t}.$
(27)

This operation allows us to select from the right-hand side of the second
equation in (26) only the one pure two-particle term which totally coincides
with those from the left-hand side.

####  III.2.4 Two-particle dynamics

Here, in addition to the previous case, we include into the second line of
(26) two-particle correlation terms themselves. Thus, the corresponding system
reads:

$\begin{split}&\hbar\frac{d}{dt}\braket{1}^{c}_{t}=T_{1}\big{[}\braket{1}^{c}_{t}\big{]}+V_{1,1}\big{[}\braket{2}_{\braket{1}^{c}_{t}}\big{]}+V_{2,1}\big{[}\braket{2}^{c}_{t}\big{]},\\\
&\hbar\frac{d}{dt}\braket{2}^{c}_{t}=T_{2}\big{[}\braket{2}^{c}_{t}\big{]}+V_{1,2}\big{[}\braket{3}_{\braket{2}^{c}_{t},\braket{1}^{c}_{t}}\big{]}.\end{split}$
(28)

The second line of system (LABEL:app_4_system_body_of_article) as will be
further demonstrated has a very different form for different types of two-
particle correlations. This feature stems from the desire to omit all the
quadratically and higher field-dependent contributions into polarization
dynamics. One of the effects appearing with including into consideration the
pure two-particle correlations is connected with screening of the Coulomb
interaction which affects the one-particle correlations dynamics by means of a
coupled system of differential equations.

####  III.2.5 Three-particle dynamics via one- and two-particle correlations

In this approximation the three-particle terms are taken into account. First,
however, we omit pure three-particle contributions. Here, within the linear in
field approximation the EOMs for one- and two-particle correlations are exact.
Thus, the corresponding system reads

$\displaystyle\hbar\frac{d}{dt}\braket{1}^{c}_{t}=T_{1}\big{[}\braket{1}^{c}_{t}\big{]}+V_{1,1}\big{[}\braket{2}_{\braket{1}^{c}_{t}}\big{]}+V_{2,1}\big{[}\braket{2}^{c}_{t}\big{]},$
$\displaystyle\hbar\frac{d}{dt}\braket{2}^{c}_{t}=T_{2}\big{[}\braket{2}^{c}_{t}\big{]}+V_{1,2}\big{[}\braket{3}_{\braket{2}^{c}_{t},\braket{1}^{c}_{t}}\big{]}+V_{2,2}\big{[}\braket{3}^{c}_{t}\big{]},$
$\displaystyle\hbar\frac{d}{dt}\braket{3}^{c}_{t}=K\big{[}T_{3}\big{[}\braket{3}^{c}_{t}\big{]}\big{]}+V^{\prime}_{1,3}\big{[}\braket{4}_{\braket{2}^{c}_{t},\braket{1}^{c}_{t}}\big{]}.$
(29)

The primed function $V^{\prime}_{1,3}$ means that only terms without momentum
summation are considered. This will be discussed later.

Within this paper we limit ourselves only by analyzing three-particle
contributions. This is motivated by the will to obtain trion-like behavior of
the absorption spectrum. Moreover, due to the noticeable numerical complexity
we deviate from the presented set of approximations at the last step by
introducing some simplifications which will be discussed further.

##  IV Results

###  IV.1 Analytics

Before we proceed with analytics and numerics, let us make some assumptions,
which allows to dramatically decrease the complexity of further computations.
From now, following the problem statement in Kira and Koch (2006) we suppose
that the analyzed systems are excited only by homogeneous electric field with
polarization lying in the sample plane. This results in coincidence of total
momenta of annihilation and creation operators in expectation values. In terms
of one-particle correlations this requirement reads as:

$\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}=\delta_{\bm{k},\bm{k}^{\prime}}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c},$
(30)

with compound indices $\bm{i}^{\prime}=(\lambda,\bm{k})$ and
$\bm{j}^{\prime}=(\lambda^{\prime},\bm{k}^{\prime})$. This simplification
leads to the fact that among all one-particle correlations only
$\mathfrak{P}(\bm{k},t)$ and $\mathfrak{n}(\lambda,\bm{k},t)$ survive in the
EOM. In the general case, the homogeneity condition is expressed as follows:

$\displaystyle\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\dots\hat{a}_{\bm{i}^{\prime}_{n}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{n}}\dots\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}=\delta_{\bm{k}_{1}+\dots+\bm{k}_{n},\bm{k}^{\prime}_{1}+\dots+\bm{k}^{\prime}_{n}}$
$\displaystyle\times\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\dots\hat{a}_{\bm{i}^{\prime}_{n}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{n}}\dots\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c},$
(31)

with similar index structure to (30). In fact, this assumption allows one to
reduce a dimension of all integrals appearing within calculations. Within the
body of the paper we present only the final expressions for low-order
approximations. Due to the declared interest in description of optical
spectra, the results will be presented for the $\bm{k}$-component of
susceptibility $\chi(\omega)$ entering the following relation:

$\tilde{\mathfrak{P}}_{\bm{k}}^{\omega}=\chi(\bm{k},\omega)\mathcal{E}(\omega).$
(32)

It should be noted that electric susceptibility $\chi(\omega)$ is one of the
most calculated quantities due to the fact that it contains a lot of
information about the optical properties of materials including oscillator
strength, absorption, refractive index, etc. The TMDs are not an exception;
there are plenty of works where the susceptibility was analyzed by means of
different theoretical approaches (see, e.g. Refs. Brunetti  _et al._ (2018);
Scholz  _et al._ (2013) and references therein).

In this work, due to the equilibrium system requirements, we also replace all
$\mathfrak{n}_{\lambda,\bm{k}}^{t}$ by their equilibrium values, i.e., Fermi-
Dirac distribution functions. Moreover, from the computational point of view
it is quite useful to work with holes instead of valence electrons. The
corresponding relations are as follows:

$\displaystyle
f_{(c,\bm{k})}=\frac{1}{e^{\beta\left(E_{g}+\hbar^{2}\bm{k}^{2}/2m_{c}-\mu_{c}\right)}+1},$
(33) $\displaystyle
f_{(v,\bm{k})}=\frac{1}{e^{\beta\left(\hbar^{2}\bm{k}^{2}/2m_{v}-\mu_{v}\right)}+1},$
(34) $\displaystyle
f_{(h,\bm{k})}=\frac{1}{e^{\beta\left(\hbar^{2}\bm{k}^{2}/2m_{h}-\mu_{h}\right)}+1},$
(35) $\displaystyle f_{(v,\bm{k})}=1-f_{(h,\bm{k})},$ (36)

where $E_{g}$ is the band gap energy, $\beta=1/kT$, $m_{h}$ and $m_{v}$ are
the hole and valence electron masses ($m_{h}=-m_{v}$, with $m_{h}>0$), while
$\mu_{h}$ and $\mu_{v}$ are the chemical potentials of holes and valence
electrons, respectively ($\mu_{h}=-\mu_{v}$) and $m_{c}$ are the conduction
band electron mass. Also, we use further $f_{(e,\bm{k})}\equiv f_{(c,\bm{k})}$
and $m_{e}\equiv m_{c}$. All the details of calculation can be found in
Appendix C.

####  IV.1.1 Free system dynamics

Within this approximation the cluster expansion technique is unnecessary. We
set the potential equal to zero and derive from the Eq. (12) the following
expression:

$\chi_{I}(\bm{k},\omega)=-\frac{d_{cv}\big{[}1-f_{(e,\bm{k})}-f_{(h,\bm{k})}\big{]}}{\hbar\big{[}\omega+i\delta-(\varepsilon_{(e,\bm{k})}+\varepsilon_{(h,\bm{k})})\big{]}},$
(37)

where $\varepsilon_{(e,\bm{k})}$ and $\varepsilon_{(h,\bm{k})}$ are one-
particle energies, defined in (72). This result coincides with the well-known
formula for susceptibility of non-interacting systems.

####  IV.1.2 One-particle dynamics and semiconductor Bloch equation

From (25) by means of the cluster expansion (14) we find:

$\chi_{II}(\bm{k},\omega)=\Gamma_{2}(\bm{k})\chi_{I}^{R}(\bm{k},\omega),$ (38)

where functions $\Gamma_{II}$ and $\chi_{I}^{R}$ are as follows:

$\displaystyle\chi_{I}^{R}(\bm{k},\omega)$ $\displaystyle=$
$\displaystyle-\frac{d_{cv}\big{[}1-f_{(e,\bm{k})}-f_{(h,\bm{k})}\big{]}}{\hbar\big{[}\omega+i\delta-(\epsilon_{(e,\bm{k})}+\epsilon_{(h,\bm{k})})\big{]}},$
(39) $\displaystyle\Gamma_{II}(\bm{k})$ $\displaystyle=$ $\displaystyle
1+\frac{1}{d_{cv}}\sum\limits_{\bm{q}\neq\bm{k}}\,\chi_{I}^{R}(\bm{q},\omega)V_{\bm{k}-\bm{q}}\Gamma_{II}(\bm{q}).$
(40)

The superscript $R$ denotes that the energies in the denominator of (39)
contains renormalized energies in contrast to (37). It can also be noted that
$\Gamma_{II}(\bm{k})$ coincides with the generalized Rabi frequency up to a
factor $d_{cv}\mathcal{E}/\hbar$. The details of calculation are presented in
Appendix C.2. This result coincides with answer obtained by means of the well-
known semiconductor Bloch equation Haug and Koch (2009).

####  IV.1.3 Multiparticle analysis

Starting from the approximation (26) the equation for dynamics of one-particle
correlations is already exact. However, it contains the terms which we do not
know explicitly. Fortunately, for these contributions some equations can be
derived.

For the susceptibility component from the first line of (26) one can obtain by
means of (92) and (21) the following expression:

$\displaystyle\chi(\bm{k},\omega)$ $\displaystyle=$
$\displaystyle\chi_{I}^{R}(\bm{k},\omega)\Bigg{[}1+\frac{1}{d_{cv}}\sum\limits_{\bm{q}\neq\bm{k}}\,V_{\bm{k}-\bm{q}}\chi(\bm{q},\omega)\Bigg{]}+\frac{1}{\mathcal{E}(\omega)\hbar\big{[}\omega+i\delta-(\epsilon_{(c,\bm{k})}-\epsilon_{(v,\bm{k})})\big{]}}$
(41) $\displaystyle\times$
$\displaystyle\sum\limits_{\bm{k}^{\prime}_{2},\bm{q}^{\prime}\neq
0}\,V_{\bm{q}^{\prime}}\Big{[}\mathfrak{D}^{\omega}_{(c,\bm{k}-\bm{q}^{\prime}),(c,\bm{k}^{\prime}_{2}),(c,\bm{k}^{\prime}_{2}-\bm{q}^{\prime}),(v,\bm{k})}+\mathfrak{D}^{\omega}_{(c,\bm{k}-\bm{q}^{\prime}),(v,\bm{k}^{\prime}_{2}),(v,\bm{k}^{\prime}_{2}-\bm{q}^{\prime}),(v,\bm{k})}$
$\displaystyle\qquad\qquad\,\,-\mathfrak{D}^{\omega}_{(c,\bm{k}),(v,\bm{k}^{\prime}_{2}),(v,\bm{k}^{\prime}_{2}-\bm{q}^{\prime}),(v,\bm{k}+\bm{q}^{\prime})}-\mathfrak{D}^{\omega}_{(c,\bm{k}),(c,\bm{k}^{\prime}_{2}),(c,\bm{k}^{\prime}_{2}-\bm{q}^{\prime}),(v,\bm{k}+\bm{q}^{\prime})}\Big{]}.$

Within all the further calculations the structure of this equation remains the
same. Only the functions $\mathfrak{D}$ will be calculated within the
different approximation orders. In order to get an idea about the structure of
functions $\mathfrak{D}$, from the second line of the system (26) within the
assumptions adopted in this paper we find the following equation for one of
the terms in (41):

$\displaystyle\mathfrak{D}^{\omega,III}_{(c,\bm{k}-\bm{q}^{\prime}),(c,\bm{k}^{\prime}_{2}),(c,\bm{k}^{\prime}_{2}-\bm{q}^{\prime}),(v,\bm{k})}=F_{\mathfrak{D},1}(\chi_{II},\bm{k}-\bm{q}^{\prime},\bm{k}^{\prime}_{2},\bm{k}^{\prime}_{2}-\bm{q}^{\prime},\bm{k},\omega,\\{cccv\\}),\
\text{where function $F_{\mathfrak{D},1}$ defines as:}$ (42) $\displaystyle
F_{\mathfrak{D},1}(\chi,\bm{k}-\bm{q}^{\prime},\bm{k}^{\prime}_{2},\bm{k}^{\prime}_{2}-\bm{q}^{\prime},\bm{k},\omega,\\{cccv\\})=\frac{\mathcal{E}(\omega)}{\hbar\Big{[}\omega+i\delta-\big{[}\varepsilon_{(c,\bm{k}-\bm{q}^{\prime})}+\varepsilon_{(c,\bm{k}^{\prime}_{2})}-\varepsilon_{(c,\bm{k}^{\prime}_{2}-\bm{q}^{\prime})}-\varepsilon_{(v,\bm{k})}\big{]}\Big{]}}$
$\displaystyle\quad\times\Bigg{[}\chi(\bm{k},\omega)V_{\bm{k}-\bm{k}^{\prime}_{2}}\Big{[}-f_{(c,\bm{k}^{\prime}_{2})}f_{(c,\bm{k}-\bm{q}^{\prime})}+f_{(c,\bm{k}^{\prime}_{2}-\bm{q}^{\prime})}\big{[}-1+f_{(c,\bm{k}^{\prime}_{2})}+f_{(c,\bm{k}-\bm{q}^{\prime})}\big{]}\Big{]}$
$\displaystyle\qquad+\chi(\bm{k}^{\prime}_{2},\omega)V_{\bm{k}-\bm{k}^{\prime}_{2}}\Big{[}f_{(c,\bm{k}-\bm{q}^{\prime})}\big{[}1-f_{(c,\bm{k}^{\prime}_{2}-\bm{q}^{\prime})}-f_{(v,\bm{k})}\big{]}+f_{(c,\bm{k}^{\prime}_{2}-\bm{q}^{\prime})}f_{(v,\bm{k})}\Big{]}$
$\displaystyle\qquad+\chi(\bm{k},\omega)V_{\bm{q}^{\prime}}\Big{[}f_{(c,\bm{k}^{\prime}_{2})}f_{(c,\bm{k}-\bm{q}^{\prime})}+f_{(c,\bm{k}^{\prime}_{2}-\bm{q}^{\prime})}\big{[}1-f_{(c,\bm{k}^{\prime}_{2})}-f_{(c,\bm{k}-\bm{q}^{\prime})}\big{]}\Big{]}$
$\displaystyle\qquad+\chi(\bm{k}-\bm{q}^{\prime},\omega)V_{\bm{q}^{\prime}}\Big{[}f_{(c,\bm{k}^{\prime}_{2})}\big{[}-1+f_{(c,\bm{k}^{\prime}_{2}-\bm{q}^{\prime})}+f_{(v,\bm{k})}\big{]}-f_{(c,\bm{k}^{\prime}_{2}-\bm{q^{\prime}})}f_{(v,\bm{k})}\Big{]}\Bigg{]}.$
(43)

As can be seen, we introduce the superscripts for two-particle correlations:
$\mathfrak{D}^{\omega,3}$. This is dictated by the needs of numerical
calculations. For the approximate solution of the emergent system of algebraic
equations we develop an iterative procedure. Let us briefly describe it. All
the functions which we introduce here can be found in Appendix C.3. As the
first basic step we take the solution obtained by means of Eq. (38):

$\displaystyle\chi_{II}(\bm{k},\omega)=F_{\chi,1}(\chi_{II},\bm{k},\omega).\quad$
(44)

As was mentioned above, this equation can be solved by the matrix inversion
approach. On the basis of this initial point, one can iteratively find further
approximations presented in the previous section by means of the following
expressions. For $N\in\\{III,IV\\}$ we have:

$\displaystyle\mathfrak{D}^{\omega,N}_{(\lambda_{1},\bm{p}_{1}),(\lambda_{2},\bm{p}_{2}),(\lambda_{3},\bm{p}_{3}),(\lambda_{4},\bm{p}_{4})}$
$\displaystyle\quad=F_{\mathfrak{D},1}(\chi_{II},\bm{p}_{1},\bm{p}_{2},\bm{p}_{3},\bm{p}_{4},\omega,\\{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}\\})$
$\displaystyle\quad+F_{\mathfrak{D},2}(\mathfrak{D}^{\omega,N-I},\bm{p}_{1},\bm{p}_{2},\bm{p}_{3},\bm{p}_{4},\omega,\\{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}\\}),$
$\displaystyle\chi_{N}(\bm{k},\omega)$
$\displaystyle\quad=F_{\chi,1}(\chi_{N},\bm{k},\omega)+F_{\chi,2}(\mathfrak{D}^{\omega,N},\bm{k},\omega),$
(45)

where by definition we consider
$\mathfrak{D}^{\omega,I}\equiv\mathfrak{D}^{\omega,II}\equiv 0$. Having
obtained the function $\mathfrak{D}^{\omega,N}$ one can calculate the modified
solution for $\chi_{N}(\bm{k},\omega)$ again by means of the matrix inversion
approach. The analysis of $N=V$ is isolated, due to the needs of applying some
simplification. Owing to the computational requirements, the corresponding
system of equations for $N=V$ reads

$\displaystyle\mathfrak{T}^{\omega,V}_{(\lambda_{1},\bm{p}_{1}),(\lambda_{2},\bm{p}_{2}),(\lambda_{3},\bm{p}_{3}),(\lambda_{4},\bm{p}_{4})}$
$\displaystyle\quad=F_{\mathfrak{T},I}(\chi_{II},\mathfrak{D}^{\omega,III},\bm{p}_{1},\bm{p}_{2},\bm{p}_{3},\bm{p}_{4},\omega,\\{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}\\}),$
$\displaystyle\mathfrak{D}^{\omega,V}_{(\lambda_{1},\bm{p}_{1}),(\lambda_{2},\bm{p}_{2}),(\lambda_{3},\bm{p}_{3}),(\lambda_{4},\bm{p}_{4})}$
$\displaystyle\quad=F_{\mathfrak{D},III}(\mathfrak{T}^{\omega,V},\bm{p}_{1},\bm{p}_{2},\bm{p}_{3},\bm{p}_{4},\omega,\\{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}\\}),$
$\displaystyle\chi_{V}(\bm{k},\omega)$
$\displaystyle\quad=F_{\chi,1}(\chi_{V},\bm{k},\omega)+F_{\chi,2}(\mathfrak{D}^{\omega,V},\bm{k},\omega).$
(46)

As was said previously, all the presented functionals $F$ can be found in
Appendix C.3 and C.4.

Figure 2: The behavior of absorption spectrum in the low-density regime
($n_{e}=n_{h}=10^{8}\,$cm-2). Inset (a) demonstrates the deviation of our
results ($E_{n}$) for exciton binding energies from the hydrogenic series of
an ideal 2d system ($E_{n}^{R}$). The Rydberg constant $R$ is calculated under
the following condition: $E_{5}=E_{5}^{R}$. Inset (b) demonstrates the
deviation of our results ($f_{n}/f_{1}$) for relative oscillator strength from
the 2d hydrogenic model of excitons ($f_{n}^{R}/f_{1}^{R}$).

Also, they are put into separate Mathematica-files, which can be found in the
Supplemental Material.

Figure 3: Change of absorption spectra behavior with increasing of exciton
concentration in the case of three different temperature values: $273\,$K,
$100\,$K, $5\,$K. The circled letter M indicates the vicinity of the Mott
transition. The star sign as a superscript for different $n_{r}$ values means
that the corresponding spectrum has rescaled regions. The shaded areas
demonstrate that the corresponding dependence has a rescaled region. Figure 4:
The dependence of $1s$ ($n=1$) excitonic peak position on the electron-hole
concentration within $\chi_{II}$ approximation. The linear regression gives
the following result for the slope ratio
$|dE_{1}/dn_{e}|=0.0026\,$meV$\cdot$$\text{\,}\mathrm{\SIUnitSymbolMicro m}$2.

###  IV.2 Numerics and discussion

In this section we obtain and compare the absorption spectra for different
approximations and physical parameter values. Due to the relatively low
computational cost of obtaining $\chi_{II}(\bm{k},\omega)$, we present it for
a wide range of carrier densities and temperatures. We note however that since
we did not account for the electron-phonon interaction, the temperature here
only affects the smearing of the Fermi surface. For some relevant values of
temperature and concentration combinations we obtain higher-order corrections
up to $\chi_{V}(\bm{k},\omega)$, where as expected the trion-like resonance
manifest itself. The detailed description of the numerical computations can be
found in Appendix D. Let us present here only the values of some physical
parameters. For the average dielectric constant and effective screening length
we choose the following reasonable numbers: $\varepsilon_{0}=2$ and
$r_{0}=5\,$nm.

Figure 5: The dependence of difference between $1s$ excitonic peak position
obtained within $\chi_{II}$ and $\chi_{III}$ approximations on the electron-
hole concentration ($\Delta=E_{1}^{\chi_{II}}-E_{1}^{\chi_{III}}$).

As for the effective masses of electrons and holes we stopped on the typical
values: $m_{e}=0.4\,m_{0}$ and $m_{h}=0.6\,m_{0}$, where $m_{0}$ is the
electron rest mass. The band gap energy $E_{g}$ is set to $1.8\,$eV. All the
other parameters which are meaningless from the physical point of view and
affect only the efficiency of the calculation procedure are defined and
discussed in Appendix D.

Figure 6: Change of absorption spectra behavior with increasing of electron
concentration when hole density is zero. The dependencies are presented for
three different temperature values: $237\,$K, $100\,$K, $5\,$K. Figure 7: The
dependence of $1s$ excitonic peak position on the electron concentration when
hole density is equated to zero. The linear regression gives the following
result for the slope ratio:
$|dE_{1}/dn_{e}|=\,0.000996\,$meV$\cdot$$\text{\,}\mathrm{\SIUnitSymbolMicro
m}$2. Figure 8: The behavior of peak height for $1s$ excitonic state. The
absorption is normalized to unity, which, however, allows us to understand the
relative change. Figure 9: The dependence of difference between $1s$ excitonic
peak position obtained within $\chi_{II}$, $\chi_{III}$, and $\chi_{IV}$
approximations on the electron-hole concentration
($\Delta_{23}=E_{1}^{\chi_{II}}-E_{1}^{\chi_{III}}$,
$\Delta{24}=E_{1}^{\chi_{II}}-E_{1}^{\chi_{IV}}$).

In order to assess the correctness of the analytical calculations presented in
the previous section, we analyze a number of dependencies. First, the
dependence of resonances position on carrier concentrations is analyzed at
fixed temperatures. Next, we study the excitonic resonance behavior relative
to the doped carrier concentration. Also, some additional dependencies which
are usually obtained via experiments will be presented.

####  IV.2.1 Carriers concentration variation

In this case we fix the temperature and vary the densities. We limit ourselves
to temperatures $5\,$K, $100\,$K, and $273\,$K. The corresponding dependencies
are presented in Fig. 3. It should be noted, that high-temperature behavior
can not be described without taking into account the phonons, therefore in the
cases of $273\,$K and $100\,$K the pictures can not be perceived as genuine
ones. Also, we considered separately the low density regime when the electron-
hole pair concentration is equal to $10^{8}\,$cm-2 (see Fig. 2). We compare
the calculated binding energy series with the analytical formula for exciton
energies of an ideal 2D system: $-R/(n-1/2)^{2}$. Assuming that for large
principal quantum numbers (in our case $n=5$) the both dependencies have to
coincide, we extract the Rydberg constant: $R=256\,$meV. It should be noted
that the estimate obtained by means of analytical expression is much greater:
$e^{4}m_{r}/2\varepsilon_{0}\hbar^{2}\approx 816\,$meV. It is seen that our
results deviate from this simple model dependence which is generally accepted
for TMDs Chernikov  _et al._ (2014); He  _et al._ (2014); Wang  _et al._
(2018b). For the case of $5\,$K, we also obtain the drift of $1s$ exciton peak
position with growth of electron-hole pairs concentration within the
$\chi_{II}$ and $\chi_{III}$ approximations. The corresponding trend is
depicted in Fig. 4. The extracted estimate of slope ration ($\approx
0.0026\,$meV$\cdot$$\text{\,}\mathrm{\SIUnitSymbolMicro m}$2) is of the same
order of magnitude as the number recently obtained within the experiment
Kravtsov  _et al._ (2020) on TMD MLs. This number allows one to extract the
exciton Bohr radius by means of the following relation:
$dE_{1}/dn_{r}=2.07E_{1}a^{2}_{B}$ Shahnazaryan  _et al._ (2017). It leads to
the following estimate: $a_{B}=2.6\,$nm. Also, we observe the onset of the
Mott transition at the concentrations with $n_{r}\approx 10,\,170,\,400$ for
$5\,$K, $100\,$K, and $273\,$K respectively. At these concentrations the
absorption peaks associated with the $1s$ exciton peak disappear and the
system is characterized by negative absorption coefficients. Moreover, at low
temperature $T=5\,$K we observe the negative absorption region in the vicinity
of the second exciton state. This is due to the fact that we disregard the
electron-phonon interaction, the excited exciton states do not relax to the
ground ones. We note also, that since we only consider a single valley in our
model, the results can not be directly compared to the experimentally observed
exciton shifts in TMDs since no contribution of the intervalley scattering can
be accounted for.

In Fig. 5 we plot the dependence of the difference between exciton shifts
obtained within the $\chi_{II}$ and $\chi_{III}$ approximations.

Figure 10: The behavior of absorption spectra calculated by means of
$\chi_{II}$ and $\chi_{V}$ approximations. The ”LF” superscript denotes that
the corresponding dependence is constructed by means of Lorentzian function(s)
which fit the calculated points. The last in turn are denoted with the
”$\star$” symbol. The quantities $X^{\chi_{II}}$, $X^{\chi_{V}}$, and
$T^{\chi_{V}}$ are equal to $1.6199(7)\,$eV, $1.6174(7)\,$eV, and
$1.6116(6)\,$eV respectively.

####  IV.2.2 Doping level variation

In this section as previously we choose the same temperature values set but
fix the number of hole concentration at value $10^{10}\,$cm-2. It is worth
noting here that the picture would change only slightly if we set the hole
concentration to zero. In particular, the peak position would remain almost
unchanged. The corresponding dependencies are depicted in Fig. 6. For zero-
valued hole concentration we calculate the behavior of $1s$ excitonic peak
position with growth of electron density which is presented in Fig. 7. The
value of the slope ratio extracted from this dependence is equal to
$\approx\,$1.0 $\text{\,}\mathrm{\SIUnitSymbolMicro
eV}$$\cdot$$\text{\,}\mathrm{\SIUnitSymbolMicro m}$2. The corresponding
dependence of difference between 1s excitonic peak position obtained within
different approximations on the electron-hole concentration is also presented
in Fig. 9. We also obtain the dependence of the normalized peak height on
electron density which is presented in Fig. 8.

In the case when the electron concentration exceeds the hole one it is
reasonable to perform the search of trionic states. For such analysis we
choose the following set of calculation parameters: $T=5\,$K
$n_{h}=10^{10}\,$cm-2, $n_{e}=40n_{h}$, $\delta=6.72\,$meV. The resuls of the
corresponding computations are presented in Fig.10. The moderate scattering of
calculated points from the expected Lorenz-like behavior of peaks is caused by
the integration accuracy of corrections. The detailed description of the
numerical procedure is presented in Appendix D.2. The fitting of obtained
results by means of Lorentzian functions allows us to extract the numerical
estimate for trion binding energy. It turns out to be about $6\,$meV which is
at least twice smaller than the experimentally observed values Bellus  _et
al._ (2015); Lundt  _et al._ (2018). This discrepancy is mainly due to the
abundance of the multivalley structure of TMDs in our model, since it is known
that the fundamental trion state corresponds to the two electrons filling
different valleys.

##  V Conclusion

To conclude, we have extended the cluster expansion technique to account for
the higher-order correlations and applied it to model the optical absorption
of the TMD monolayer. Our account of the three-particle correlations allowed
us to model the absorption peaks associated with the trion quasiparticles. The
developed technique allows us to directly model experimentally accessible
absorption spectra. The main advantage of the proposed formalism is that
despite the computational complexity, it allows us to address the regime of
large electron-hole densities and large doping, when the Fermi energy becomes
comparable to the exciton and trion binding energies and, ultimately, the
exciton Mott transition.

The natural development of the presented formalism would be to include the
electron-phonon interaction as well as spin and valley degrees of freedom. The
cluster expansion technique proved to be a powerful tool to model the
absorption spectra mediated by the strong many body correlations present in
TMD monolayers and heterostructures.

###### Acknowledgements.

We highly appreciate support of RFBR (Project No. 19-52-51010), at the expense
of which all the numerical computations were fulfilled. The analytical
calculations were performed through support of RFBR and DFG (Project No.
21-52-12038). A.K. also gratefully acknowledges Mathew Colgrove from NVidia
for the help in debugging the code.

##  Appendix A The information about the Supplemental Material

All the functionals presented in this paper in terms of correlations can be
found in Mathematica-file (convoluted_expressions_in_terms_of_arb_mom.nb).

##  Appendix B Connection between notations

Here we present the connection formulas between two types of notation in (2)
and (3). For the model considered in this work these expressions are as
follows:

$\begin{split}&\hbar
h^{ext}_{\bm{i}\bm{j}}(t)=-\mathcal{E}(t)d_{cv}\sum\limits_{\bm{k}}\big{[}\delta_{\bm{j},(v,\bm{k})}\delta_{\bm{i},(c,\bm{k})}+\delta_{\bm{j},(c,\bm{k})}\delta_{\bm{i},(v,\bm{k})}\big{]},\\\
&v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}=\frac{1}{\hbar}\sum\limits_{\begin{subarray}{c}\,\bm{k}_{1},\bm{k}_{2}\\\
\bm{q}\neq
0\end{subarray}}2\,V_{q}\big{[}\mathfrak{Kron}(\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1})-\mathfrak{Kron}(\bm{i}_{2},\bm{i}_{1},\bm{j}_{2},\bm{j}_{1})+\mathfrak{Kron}(\bm{i}_{2},\bm{i}_{1},\bm{j}_{1},\bm{j}_{2})-\mathfrak{Kron}(\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2})\big{]},\\\
&\mathfrak{Kron}(\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1})=\frac{1}{4}\big{[}\delta_{(c,\bm{k_{1}+q}),\bm{i}_{1}}\delta_{(c,\bm{k_{2}-q}),\bm{i}_{2}}\delta_{(c,\bm{k}_{2}),\bm{j}_{2}}\delta_{(c,\bm{k}_{1}),\bm{j}_{1}}+\delta_{(v,\bm{k_{1}+q}),\bm{i}_{1}}\delta_{(v,\bm{k_{2}-q}),\bm{i}_{2}}\delta_{(v,\bm{k}_{2}),\bm{j}_{2}}\delta_{(v,\bm{k}_{1}),\bm{j}_{1}}\\\
&\qquad\qquad\qquad\qquad\,\
+\delta_{(c,\bm{k_{1}+q}),\bm{i}_{1}}\delta_{(v,\bm{k_{2}-q}),\bm{i}_{2}}\delta_{(v,\bm{k}_{2}),\bm{j}_{2}}\delta_{(c,\bm{k}_{1}),\bm{j}_{1}}+\delta_{(v,\bm{k_{1}+q}),\bm{i}_{1}}\delta_{(c,\bm{k_{2}-q}),\bm{i}_{2}}\delta_{(c,\bm{k}_{2}),\bm{j}_{2}}\delta_{(v,\bm{k}_{1}),\bm{j}_{1}}\big{]}.\end{split}$
(47)

These artificially long expressions are caused by the symmetrization
conditions (5).

##  Appendix C EsOM. Series of approximations

Because we are interested in the frequency spectrum, we introduce the Fourier
transform and work further with algebraic equations instead of differential
ones. Within this paper for the Fourier transform we accept the following
definition:

$\tilde{\mathfrak{P}}(\bm{k},\omega)=\int\limits_{-\infty}^{\infty}\,dt\,\tilde{\mathfrak{P}}(\bm{k},t)\,\exp{(i\omega
t)},\quad\tilde{\mathfrak{P}}(\bm{k},t)=\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}\,d\omega\,\tilde{\mathfrak{P}}(\bm{k},\omega)\,\exp{(-i\omega
t)}.$ (48)

Let us derive the general expressions which are necessary for all the
approximation orders. On the basis of (12), we derive the dynamical equations
on correlations. First, the commutator of
$\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}$ and $H_{el,k}$
gives

$i\braket{\big{[}\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}},H_{el,k}\big{]}}_{t}=i\hbar\big{[}\varepsilon_{\bm{i}^{\prime}}-\varepsilon_{\bm{j}^{\prime}}\big{]}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}=i\hbar\big{[}\varepsilon_{\bm{i}^{\prime}}-\varepsilon_{\bm{j}^{\prime}}\big{]}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}.$
(49)

The similar relation is true for the operator (13):

$i\braket{\big{[}\hat{a}_{\bm{i}^{\prime}_{1}}\dots\hat{a}_{\bm{i}^{\prime}_{n}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{n}}\dots\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}},H_{el,k}\big{]}}_{t}=i\hbar\big{[}\varepsilon_{\bm{i}^{\prime}_{1}}+\dots+\varepsilon_{\bm{i}^{\prime}_{n}}-\varepsilon_{\bm{j}^{\prime}_{n}}-\dots-\varepsilon_{\bm{j}_{1}^{\prime}}\big{]}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\dots\hat{a}_{\bm{i}^{\prime}_{n}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{n}}\dots\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}.$
(50)

The next important contribution is caused by the presence of the external
field:

$\displaystyle
i\braket{\big{[}\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}},H_{I}\big{]}}_{t}$
$\displaystyle=$ $\displaystyle
i\hbar\sum\limits_{\bm{i},\bm{j}}h^{ext}_{\bm{i}\bm{j}}\big{[}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}^{\dagger}_{\bm{i}}\hat{a}_{\bm{j}}}_{t}-\braket{\hat{a}^{\dagger}_{\bm{i}}\hat{a}_{\bm{j}}\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}\big{]}=i\hbar\bigg{[}\sum\limits_{\bm{j}}h^{ext}_{\bm{i}^{\prime}\bm{j}}\braket{\hat{a}_{\bm{j}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}-\sum\limits_{\bm{i}}h^{ext}_{\bm{i}\bm{j}^{\prime}}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}}}_{t}\bigg{]}$
(51) $\displaystyle=$ $\displaystyle
i\hbar\bigg{[}\sum\limits_{\bm{j}}h^{ext}_{\bm{i}^{\prime}\bm{j}}\braket{\hat{a}_{\bm{j}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}-\sum\limits_{\bm{i}}h^{ext}_{\bm{i}\bm{j}^{\prime}}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}}}_{t}^{c}\bigg{]}=i\hbar\bigg{[}\sum\limits_{\bm{j}}T^{2}_{\\{\bm{i}^{\prime}\bm{j}\\},\\{\bm{j},\bm{j}^{\prime}\\}}-\sum\limits_{\bm{i}}T^{2}_{\\{\bm{i}\bm{j}^{\prime}\\},\\{\bm{i}^{\prime},\bm{i}\\}}\bigg{]},$

where $T^{2}_{\\{\bm{i}^{\prime}\bm{j}\\},\\{\bm{j},\bm{j}^{\prime}\\}}$ is
defined as follows:

$T^{2}_{\\{\bm{i}^{\prime}\bm{j}\\},\\{\bm{j},\bm{j}^{\prime}\\}}=h^{ext}_{\bm{i}^{\prime}\bm{j}}\braket{\hat{a}_{\bm{j}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}=h^{ext}_{\bm{i}^{\prime}\bm{j}}\braket{\hat{a}_{\bm{j}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}^{c}_{t}=T^{2}_{c,\\{\bm{i}^{\prime}\bm{j}\\},\\{\bm{j},\bm{j}^{\prime}\\}}.$
(52)

The commutator of $H_{I}$ with four- and six-operator combinations can be
written by analogy:

$\displaystyle
i\braket{\big{[}\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}},H_{I}\big{]}}_{t}$
$\displaystyle=$ $\displaystyle
i\hbar\bigg{[}\sum\limits_{\bm{j}}\left(T^{4}_{\\{\bm{i}^{\prime}_{1}\bm{j}\\},\\{\bm{j}\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}+T^{4}_{\\{\bm{i}^{\prime}_{2}\bm{j}\\},\\{\bm{i}^{\prime}_{1}\bm{j},\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}\right)-\sum\limits_{\bm{i}}\left(T^{4}_{\\{\bm{i}\bm{j}^{\prime}_{2}\\},\\{\bm{i}^{\prime}_{1}\bm{i}^{\prime}_{2},\bm{i}\bm{j}^{\prime}_{1}\\}}+T^{4}_{\\{\bm{i}\bm{j}^{\prime}_{1}\\},\\{\bm{i}^{\prime}_{1}\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2}\bm{i}\\}}\right)\bigg{]},$
$\displaystyle
i\braket{\big{[}\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}_{\bm{i}^{\prime}_{3}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{3}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}},H_{I}\big{]}}_{t}$
$\displaystyle=$ $\displaystyle
i\hbar\bigg{[}\sum\limits_{\bm{j}}\left(T^{6}_{\\{\bm{i}^{\prime}_{1}\bm{j}\\},\\{\bm{j}\bm{i}^{\prime}_{2}\bm{i}^{\prime}_{3},\bm{j}^{\prime}_{3}\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}+T^{6}_{\\{\bm{i}^{\prime}_{2}\bm{j}\\},\\{\bm{i}^{\prime}_{1}\bm{j}\bm{i}^{\prime}_{3},\bm{j}^{\prime}_{3}\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}+T^{6}_{\\{\bm{i}^{\prime}_{3}\bm{j}\\},\\{\bm{i}^{\prime}_{1}\bm{i}^{\prime}_{2}\bm{j},\bm{j}^{\prime}_{3}\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}\right)$
(53)
$\displaystyle\quad-\sum\limits_{\bm{i}}\left(T^{6}_{\\{\bm{i}\bm{j}^{\prime}_{3}\\},\\{\bm{i}^{\prime}_{1}\bm{i}^{\prime}_{2}\bm{i}^{\prime}_{3},\bm{i}\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}+T^{6}_{\\{\bm{i}\bm{j}^{\prime}_{2}\\},\\{\bm{i}^{\prime}_{1}\bm{i}^{\prime}_{2}\bm{i}^{\prime}_{3},\bm{j}^{\prime}_{3}\bm{i}\bm{j}^{\prime}_{1}\\}}+T^{6}_{\\{\bm{i}\bm{j}^{\prime}_{1}\\},\\{\bm{i}^{\prime}_{1}\bm{i}^{\prime}_{2}\bm{i}^{\prime}_{3},\bm{j}^{\prime}_{3}\bm{j}^{\prime}_{2}\bm{i}\\}}\right)\bigg{]},$

where
$T^{4}_{\\{\bm{i}^{\prime}_{1}\bm{j}\\},\\{\bm{j}\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}$
and
$T^{6}_{\\{\bm{i}^{\prime}_{1}\bm{j}\\},\\{\bm{j}\bm{i}^{\prime}_{2}\bm{i}^{\prime}_{3},\bm{j}^{\prime}_{3}\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}$
read as:

$T^{4}_{\\{\bm{i}^{\prime}_{1}\bm{j}\\},\\{\bm{j}\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}=h^{ext}_{\bm{i}^{\prime}_{1}\bm{j}}\braket{\hat{a}_{\bm{j}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t},\quad
T^{6}_{\\{\bm{i}^{\prime}_{1}\bm{j}\\},\\{\bm{j}\bm{i}^{\prime}_{2}\bm{i}^{\prime}_{3},\bm{j}^{\prime}_{3}\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}=h^{ext}_{\bm{i}^{\prime}_{1}\bm{j}}\braket{\hat{a}_{\bm{j}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}_{\bm{i}^{\prime}_{3}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{3}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}.$
(54)

Let us note, however, that these quantities, in contrast to $T^{2}$, do not
coincide with their correlation counterparts. For this purpose one should
address the cluster expansion (14).

Based on the derived expressions the logic of the construction of commutators
with many-operator combinations is clear. All the aforementioned relations
allow one to rewrite the left-hand side of (12) for expectation values in
terms of correlations,

$\displaystyle\frac{d}{dt}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}+\frac{1}{\hbar}i\braket{\big{[}\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}},H_{el,k}\big{]}}_{t}+\frac{1}{\hbar}i\braket{\big{[}\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}},H_{I}\big{]}}_{t}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad=\frac{d}{dt}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}+i\big{[}\varepsilon_{\bm{i}^{\prime}}-\varepsilon_{\bm{j}^{\prime}}\big{]}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}^{c}_{t}+i\bigg{[}\sum\limits_{\bm{j}}h^{ext}_{\bm{i}^{\prime}\bm{j}}\braket{\hat{a}_{\bm{j}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}-\sum\limits_{\bm{i}}h^{ext}_{\bm{i}\bm{j}^{\prime}}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}}}_{t}^{c}\bigg{]},$
(55)

or in terms of schematic notation (24),

$\displaystyle\hbar\frac{d}{dt}\braket{1}^{c}_{t}=\hbar\frac{d}{dt}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}^{c}_{t},\quad
T_{1}\big{[}\braket{1}^{c}_{t}\big{]}=-\hbar\Bigg{[}i\big{[}\varepsilon_{\bm{i}^{\prime}}-\varepsilon_{\bm{j}^{\prime}}\big{]}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}^{c}_{t}+i\bigg{[}\sum\limits_{\bm{j}}h^{ext}_{\bm{i}^{\prime}\bm{j}}\braket{\hat{a}_{\bm{j}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}-\sum\limits_{\bm{i}}h^{ext}_{\bm{i}\bm{j}^{\prime}}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}}}_{t}^{c}\bigg{]}\Bigg{]}.$
(56)

In the next sections we include the consideration of Coulomb interaction.

###  C.1 Free system dynamics

First, let us reproduce free-system polarization. We set $V=0$ and obtain from
(55) the following equation:

$\frac{d}{dt}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}+i\big{[}\varepsilon_{\bm{i}^{\prime}}-\varepsilon_{\bm{j}^{\prime}}\big{]}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}^{c}_{t}+i\bigg{[}\sum\limits_{\bm{j}}h^{ext}_{\bm{i}^{\prime}\bm{j}}\braket{\hat{a}_{\bm{j}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}-\sum\limits_{\bm{i}}h^{ext}_{\bm{i}\bm{j}^{\prime}}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}}}_{t}^{c}\bigg{]}=0.$
(57)

For the purposes of further computations it is quite instructive to introduce
the diagrammatic rules and to accompany all the analytic expressions by their
graphical counterparts. Moreover, the three-particle correlation dynamics will
be derived only by means of the powerful method of Feynman diagrams. The
detailed description of the technique with explanation of prefactor and sign
choice of a diagram was presented in Fricke (1996). In this work we add only
some customization. Also, we limit ourselves only by consideration of
connected graphs which appear in equations on correlations in contrast to
differential equations on expectation values where unconnected diagrams are
also presented. The diagram elements are as follows:

$\displaystyle=$
$\displaystyle\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}},\quad\begin{gathered}\includegraphics[height=48.36958pt]{diag/rules/field_diag_dest.pdf}\end{gathered}=\sum\limits_{\bm{j}}h^{ext}_{\bm{i}^{\prime}\bm{j}},\quad\begin{gathered}\includegraphics[height=48.36958pt]{diag/rules/field_diag_cret.pdf}\end{gathered}=-\sum\limits_{\bm{i}}h^{ext}_{\bm{i}\bm{j}^{\prime}},$
(61) $\displaystyle=$
$\displaystyle\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\dots\hat{a}_{\bm{i}^{\prime}_{n}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{n}}\dots\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c},\quad\begin{gathered}\includegraphics[height=48.36958pt]{diag/rules/contraction.pdf}\end{gathered}=\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c},\quad\begin{gathered}\includegraphics[height=48.36958pt]{diag/rules/external_line_cr.pdf}\end{gathered},\begin{gathered}\includegraphics[height=48.36958pt]{diag/rules/external_line_ds.pdf}\end{gathered},$
(66)

where the last two graphs – external vertices – are involved in different
contractions. In this notation Eq. (57) is presented in Fig. 11. The result of
applying operator $[\hat{D}+\hat{E}]$ is obvious. Here and further in the
brackets one can see the number of terms which correspond to a presented
unlabeled diagram. For $\mathfrak{P}_{\bm{k}}^{t}$ from (18) one derives the
following:

$\displaystyle\frac{d}{dt}\braket{\hat{a}_{(c,\bm{k})}\hat{a}^{\dagger}_{(v,\bm{k})}}_{t}^{c}+i\big{[}\varepsilon_{(c,\bm{k})}$
$\displaystyle-$
$\displaystyle\varepsilon_{(v,\bm{k})}\big{]}\braket{\hat{a}_{(c,\bm{k})}\hat{a}^{\dagger}_{(v,\bm{k})}}_{t}^{c}+i\bigg{[}\sum\limits_{\bm{j}}h^{ext}_{(c,\bm{k})\bm{j}}\braket{\hat{a}_{\bm{j}}\hat{a}^{\dagger}_{(v,\bm{k})}}_{t}^{c}-\sum\limits_{\bm{i}}h^{ext}_{\bm{i}(v,\bm{k})}\braket{\hat{a}_{(c,\bm{k})}\hat{a}^{\dagger}_{\bm{i}}}_{t}^{c}\bigg{]}=0,$
(67) $\displaystyle
i\sum\limits_{\bm{j}}h^{ext}_{(c,\bm{k})\bm{j}}\braket{\hat{a}_{\bm{j}}\hat{a}^{\dagger}_{(v,\bm{k})}}_{t}^{c}$
$\displaystyle=$ $\displaystyle
i\left(-\frac{\mathcal{E}(t)d_{cv}}{\hbar}\right)\sum\limits_{\bm{k}^{\prime}}\sum\limits_{\bm{j}}\big{[}\delta_{(c,\bm{k}),(c,\bm{k}^{\prime})}\delta_{\bm{j},(v,\bm{k^{\prime}})}+\delta_{\bm{j},(c,\bm{k}^{\prime})}\delta_{(c,\bm{k}),(v,\bm{k^{\prime}})}\big{]}\braket{\hat{a}_{\bm{j}}\hat{a}^{\dagger}_{(v,\bm{k})}}_{t}^{c}=$
(68) $\displaystyle=$ $\displaystyle
i\left(-\frac{\mathcal{E}(t)d_{cv}}{\hbar}\right)\big{[}1-\mathfrak{n}_{v,\bm{k}}^{t}\big{]},$
$\displaystyle
i\sum\limits_{\bm{i}}h^{ext}_{\bm{i}(v,\bm{k})}\braket{\hat{a}_{(c,\bm{k})}\hat{a}^{\dagger}_{\bm{i}}}_{t}^{c}$
$\displaystyle=$ $\displaystyle
i\left(-\frac{\mathcal{E}(t)d_{cv}}{\hbar}\right)\big{[}1-\mathfrak{n}_{c,\bm{k}}^{t}\big{]}.$
(69)

Taking into account the relation (11), in terms of polarization component we
obtain:

$\hbar\frac{d}{dt}\tilde{\mathfrak{P}}_{\bm{k}}^{t}+i\hbar\big{[}\varepsilon_{(c,\bm{k})}-\varepsilon_{(v,\bm{k})}\big{]}\tilde{\mathfrak{P}}_{\bm{k}}^{t}+i\mathcal{E}(t)d_{cv}\big{[}\mathfrak{n}_{c,\bm{k}}^{t}-\mathfrak{n}_{v,\bm{k}}^{t}\big{]}=0.$
(70)

Figure 11: The differential equation on one-particle correlations within free-
system approximation ($N=I$).

Thus, we derive the dynamical equation on the polarization component in the
noninteracting case. The corresponding algebraic equation for the
susceptibility component (32) is as follows:

$\chi_{I}(\bm{k},\omega)=\frac{d_{cv}\big{[}f_{(c,\bm{k})}-f_{(v,\bm{k})}\big{]}}{\hbar\big{[}\omega+i\delta-(\varepsilon_{(c,\bm{k})}-\varepsilon_{(v,\bm{k})})\big{]}}=-\frac{d_{cv}\big{[}1-f_{(e,\bm{k})}-f_{(h,\bm{k})}\big{]}}{\hbar\big{[}\omega+i\delta-(\varepsilon_{(e,\bm{k})}+\varepsilon_{(h,\bm{k})})\big{]}}.$
(71)

Also, instead of $\mathfrak{n}$ we consider their equilibrium values – Fermi-
Dirac distributions of carriers in valence and conductivity bands. The one-
particle energies are defined as follows:

$\hbar\varepsilon_{(e,\bm{k})}\equiv\hbar\varepsilon_{(c,\bm{k})}=E_{g}+\frac{\hbar^{2}\bm{k}^{2}}{2m_{c}},\quad\hbar\varepsilon_{(v,\bm{k})}=\frac{\hbar^{2}\bm{k}^{2}}{2m_{v}},\quad\hbar\varepsilon_{(h,\bm{k})}=\frac{\hbar^{2}\bm{k}^{2}}{2m_{h}}.$
(72)

###  C.2 One-particle dynamics

In this part we aim to derive the function
$V_{1,1}\big{[}\braket{2}_{\braket{1}^{c}_{t}}\big{]}$ from (25). For this
purpose, however, we have to rewrite commutator of
$\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}$ and $\hat{V}$
in terms of correlations. For the right-hand side of (12) we obtain

$\displaystyle-i\braket{\big{[}\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}},V\big{]}}_{t}=-i\hbar\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}CM_{1}[\bm{i}^{\prime},\bm{j}^{\prime},\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}],$
(73)

where for the convenience of the further calculations the two-operator-
potential commutator was introduced:

$CM_{1}[\bm{i}^{\prime},\bm{j}^{\prime},\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}]=\frac{1}{4}\Bigg{[}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}-\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}\hat{a}_{\bm{j}_{1}}\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}\Bigg{]}.$
(74)

By means of the cluster expansion technique (14) the appeared difference of
expectation values can be expanded as follows:

$\displaystyle\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}$
$\displaystyle-$
$\displaystyle\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}\hat{a}_{\bm{j}_{1}}\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}$
(75) $\displaystyle=$
$\displaystyle\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}-\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}-\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}$
$\displaystyle+$
$\displaystyle\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}-\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}}}_{t}^{c}+\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}$
$\displaystyle+$
$\displaystyle\braket{\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}}}_{t}^{c}-\braket{\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}+\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}$
$\displaystyle+$
$\displaystyle\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\braket{\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}-\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}-\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}$
$\displaystyle+$
$\displaystyle\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}+\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}-\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}$
$\displaystyle-$
$\displaystyle\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}.$

Due to the symmetry properties of potential (5) and correlations (17), the
expression (75) can be reduced to the combination of 3 unique terms (grouped
by means of square brackets):

$\displaystyle\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}$
$\displaystyle-$
$\displaystyle\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}\hat{a}_{\bm{j}_{1}}\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}$
(76) $\displaystyle=$
$\displaystyle\,4\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}\bigg{[}-\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}+\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}}}_{t}^{c}\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}\bigg{]}$
$\displaystyle+$
$\displaystyle\,2\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\bigg{[}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}+\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}\bigg{]}$
$\displaystyle-$
$\displaystyle\,2\braket{\hat{a}_{\bm{j}_{1}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}\bigg{[}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}+\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}}}_{t}^{c}\bigg{]}.$

This grouping will be clear later. Thus, taking into account (4) and (15) for
commutator (73) we find:

$\displaystyle-i\braket{\big{[}\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}},V\big{]}}_{t}=-i\hbar\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}$
$\displaystyle\Bigg{[}\frac{1}{2}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\delta_{\bm{j}^{\prime}\bm{j}_{1}}-\frac{1}{2}\braket{\hat{a}_{\bm{j}_{1}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}\delta_{\bm{i}^{\prime}\bm{i}_{1}}$
$\displaystyle+\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}$
$\displaystyle\bigg{[}-\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}+\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}}}_{t}^{c}\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}\bigg{]}\Bigg{]}.$
(77)

In terms of the schematic equation (23) we obtain the following parts:

$\displaystyle V_{1,1}\big{[}\braket{2}_{\braket{1}^{c}_{t}}\big{]}$
$\displaystyle=$
$\displaystyle-i\hbar\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}\Bigg{[}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}\bigg{[}-\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}+\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}}}_{t}^{c}\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}\bigg{]}\Bigg{]},$
(78) $\displaystyle V_{2,1}\big{[}\braket{2}^{c}_{t}\big{]}$ $\displaystyle=$
$\displaystyle-i\hbar\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}\Bigg{[}\frac{1}{2}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\delta_{\bm{j}^{\prime}\bm{j}_{1}}-\frac{1}{2}\braket{\hat{a}_{\bm{j}_{1}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}\delta_{\bm{i}^{\prime}\bm{i}_{1}}\Bigg{]}.$
(79)

The expression (78) allows us to write the equation on polarization dynamics
within a one-particle approximation. Hence, for (25) we find

$\displaystyle\frac{d}{dt}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}+i\big{[}\varepsilon_{\bm{i}^{\prime}}-\varepsilon_{\bm{j}^{\prime}}\big{]}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}^{c}_{t}+i\bigg{[}\sum\limits_{\bm{j}}h^{ext}_{\bm{i}^{\prime}\bm{j}}\braket{\hat{a}_{\bm{j}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}-\sum\limits_{\bm{i}}h^{ext}_{\bm{i}\bm{j}^{\prime}}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}}}_{t}^{c}\bigg{]}$
$\displaystyle\qquad\qquad=-i\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}\Bigg{[}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}\bigg{[}-\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}+\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}}}_{t}^{c}\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}\bigg{]}\Bigg{]},$
(80)

where the corresponding diagrammatic counterpart is depicted in Fig. 12.

Figure 12: The differential equation on one-particle correlations within one-
particle approximation ($N=II$).

Let us obtain now the algebraic equation for susceptibility. For
$\mathfrak{P}(\bm{k},t)$ from (18) one has to treat
$V_{1,1}\big{[}\braket{2}_{\braket{1}^{c}_{t}}\big{]}$:

$\displaystyle V_{1,1}\big{[}\braket{2}_{\braket{1}^{c}_{t}}\big{]}$
$\displaystyle=$ $\displaystyle-i\sum\limits_{\bm{q}\neq
0}V_{\bm{q}}\bigg{[}(\mathfrak{n}_{v,\bm{k}+\bm{q}}^{t}-\mathfrak{n}_{c,\bm{k}+\bm{q}}^{t})\mathfrak{P}_{\bm{k}}^{t}+(\mathfrak{n}_{c,\bm{k}}^{t}-\mathfrak{n}_{v,\bm{k}}^{t})\mathfrak{P}_{\bm{k}+\bm{q}}^{t}\bigg{]}$
(81) $\displaystyle=$ $\displaystyle
i\Big{[}\mathfrak{P}_{\bm{k}}^{t}\sum\limits_{\bm{q}\neq\bm{k}}\,V_{\bm{k}-\bm{q}}(\mathfrak{n}_{c,\bm{q}}^{t}-\mathfrak{n}_{v,\bm{q}}^{t})+(\mathfrak{n}_{v,\bm{k}}^{t}-\mathfrak{n}_{c,\bm{k}}^{t})\sum\limits_{\bm{q}\neq\bm{k}}\,V_{\bm{k}-\bm{q}}\mathfrak{P}_{\bm{q}}^{t}\Big{]}.$

Combining the previously obtained results and keeping in mind that
$\tilde{\mathfrak{P}}_{\bm{k}}^{t}=-\mathfrak{P}_{\bm{k}}^{t}$, we obtain the
well-known semiconductor Bloch equation Haug and Koch (2009):

$i\hbar\frac{d}{dt}\tilde{\mathfrak{P}}_{\bm{k}}^{t}-\Bigg{[}\hbar\epsilon_{(c,\bm{k})}-\hbar\epsilon_{(v,\bm{k})}\Bigg{]}\tilde{\mathfrak{P}}_{\bm{k}}^{t}=\big{[}\mathfrak{n}_{c,\bm{k}}^{t}-\mathfrak{n}_{v,\bm{k}}^{t}\big{]}\Bigg{[}\mathcal{E}(t)d_{cv}+\sum\limits_{\bm{q}\neq\bm{k}}\,V_{\bm{k}-\bm{q}}\tilde{\mathfrak{P}}_{\bm{q}}^{t}\Bigg{]}.$
(82)

where the renormalized energies were introduced:

$\displaystyle\hbar\epsilon_{(c,\bm{k})}$
$\displaystyle=\hbar\varepsilon_{(c,\bm{k})}-\sum\limits_{\bm{q}\neq\bm{k}}\,V_{\bm{k}-\bm{q}}\mathfrak{n}_{c,\bm{q}}^{t},$
(83) $\displaystyle\hbar\epsilon_{(v,\bm{k})}$
$\displaystyle=\hbar\varepsilon_{(v,\bm{k})}-\sum\limits_{\bm{q}\neq\bm{k}}\,V_{\bm{k}-\bm{q}}\mathfrak{n}_{v,\bm{q}}^{t},$
(84) $\displaystyle\hbar\epsilon_{(h,\bm{k})}$
$\displaystyle=\hbar\varepsilon_{(h,\bm{k})}-\sum\limits_{\bm{q}\neq\bm{k}}\,V_{\bm{k}-\bm{q}}\mathfrak{n}_{h,\bm{q}}^{t}.$
(85)

The corresponding algebraic equation for susceptibility component
$\chi(\bm{k},\omega)$ reads:

$\hbar\Bigg{[}\omega+i\delta-\bigg{[}\epsilon_{(c,\bm{k})}-\epsilon_{(v,\bm{k})}\bigg{]}\Bigg{]}\chi_{II}(\bm{k},\omega)=\big{[}f_{(c,\bm{k})}-f_{(v,\bm{k})}\big{]}\Bigg{[}d_{cv}+\sum\limits_{\bm{q}\neq\bm{k}}\,V_{\bm{k}-\bm{q}}\chi_{II}(\bm{q},\omega)\Bigg{]}.$
(86)

This equation can be rewritten:

$\chi_{II}(\bm{k},\omega)=\Gamma_{II}(\bm{k})\chi_{I}^{R}(\bm{k},\omega),$
(87)

where in terms of electrons and holes we have:

$\chi_{I}^{R}(\bm{k},\omega)=-\frac{d_{cv}\big{[}1-f_{(e,\bm{k})}-f_{(h,\bm{k})}\big{]}}{\hbar\big{[}\omega+i\delta-(\epsilon_{(e,\bm{k})}+\epsilon_{(h,\bm{k})})\big{]}},\quad\Gamma_{II}(\bm{k})=1+\frac{1}{d_{cv}}\sum\limits_{\bm{q}\neq\bm{k}}\,\chi_{I}^{R}(\bm{q},\omega)V_{\bm{k}-\bm{q}}\Gamma_{II}(\bm{q}).$
(88)

The introducing of the function $\Gamma$ allows one to solve the integral
equation (87) by means of the matrix inversion approach.

###  C.3 Two-particle dynamics

The first line of the system (26) can be simply derived by means of the
expression (79). Adding this term to (C.2), we obtain exact equation for one-
particle correlations dynamics:

$\displaystyle\frac{d}{dt}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}+i\big{[}\varepsilon_{\bm{i}^{\prime}}-\varepsilon_{\bm{j}^{\prime}}\big{]}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}^{c}_{t}+i\bigg{[}\sum\limits_{\bm{j}}h^{ext}_{\bm{i}^{\prime}\bm{j}}\braket{\hat{a}_{\bm{j}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}-\sum\limits_{\bm{i}}h^{ext}_{\bm{i}\bm{j}^{\prime}}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}}}_{t}^{c}\bigg{]}=-i\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}$
$\displaystyle\times\Bigg{[}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}\bigg{[}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}}}_{t}^{c}\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}-\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\bigg{]}+\frac{1}{2}\bigg{[}\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\delta_{\bm{j}^{\prime}\bm{j}_{1}}-\braket{\hat{a}_{\bm{j}_{1}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}_{t}^{c}\delta_{\bm{i}^{\prime}\bm{i}_{1}}\bigg{]}\Bigg{]},$
(89)

or in terms of the previously introduced quantities:

$\displaystyle\Big{[}\frac{d}{dt}+i\big{[}\varepsilon_{\bm{i}^{\prime}}-\varepsilon_{\bm{j}^{\prime}}\big{]}\Big{]}$
$\displaystyle\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}^{\dagger}_{\bm{j}^{\prime}}}^{c}_{t}$
(90)
$\displaystyle=-i\bigg{[}\sum\limits_{\bm{j}}T^{2}_{c,\\{\bm{i}^{\prime}\bm{j}\\},\\{\bm{j},\bm{j}^{\prime}\\}}-\sum\limits_{\bm{i}}T^{2}_{c,\\{\bm{i}\bm{j}^{\prime}\\},\\{\bm{i}^{\prime},\bm{i}\\}}\bigg{]}-i\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}\,CM_{1}[\bm{i}^{\prime},\bm{j}^{\prime},\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}],$

where the corresponding graphical equation is presented in Fig. 13.

Figure 13: The exact differential equation on one-particle correlations
dynamics.

Taking into account the potential form (47), for polarization component from
(C.3) the dynamical equation looks as follows:

$\displaystyle\hbar\frac{d}{dt}\mathfrak{P}_{\bm{k}}^{t}+i\hbar\Bigg{[}\epsilon_{(c,\bm{k})}-\epsilon_{(v,\bm{k})}\Bigg{]}\mathfrak{P}_{\bm{k}}^{t}-i\big{[}\mathfrak{n}_{c,\bm{k}}^{t}-\mathfrak{n}_{v,\bm{k}}^{t}\big{]}\Bigg{[}\mathcal{E}(t)d_{cv}-\sum\limits_{\bm{q}\neq\bm{k}}\,V_{\bm{k}-\bm{q}}\mathfrak{P}_{\bm{q}}^{t}\Bigg{]}$
$\displaystyle=i\sum\limits_{\bm{k}_{2},\bm{q}\neq
0}\,V_{\bm{q}}\Big{[}\braket{\hat{a}_{(c,\bm{k}-\bm{q})}\hat{a}_{(c,\bm{k}_{2})}\hat{a}^{\dagger}_{(c,\bm{k}_{2}-\bm{q})}\hat{a}^{\dagger}_{(v,\bm{k})}}_{t}^{c}+\braket{\hat{a}_{(c,\bm{k}-\bm{q})}\hat{a}_{(v,\bm{k}_{2})}\hat{a}^{\dagger}_{(v,\bm{k}_{2}-\bm{q})}\hat{a}^{\dagger}_{(v,\bm{k})}}_{t}^{c}$
$\displaystyle\qquad\qquad\qquad-\braket{\hat{a}_{(c,\bm{k})}\hat{a}_{(v,\bm{k}_{2})}\hat{a}^{\dagger}_{(v,\bm{k}_{2}-\bm{q})}\hat{a}^{\dagger}_{(v,\bm{k}+\bm{q})}}_{t}^{c}-\braket{\hat{a}_{(c,\bm{k})}\hat{a}_{(c,\bm{k}_{2})}\hat{a}^{\dagger}_{(c,\bm{k}_{2}-\bm{q})}\hat{a}^{\dagger}_{(v,\bm{k}+\bm{q})}}_{t}^{c}\Big{]},$
(91)

while the corresponding algebraic equation for susceptibility component reads
as:

$\displaystyle\chi(\bm{k},\omega)=\chi_{I}^{R}(\bm{k},\omega)\Bigg{[}1+\frac{1}{d_{cv}}\sum\limits_{\bm{q}\neq\bm{k}}\,V_{\bm{k}-\bm{q}}\chi(\bm{q},\omega)\Bigg{]}+\frac{1}{\mathcal{E}(\omega)\hbar\big{[}\omega+i\delta-(\epsilon_{(e,\bm{k})}+\epsilon_{(h,\bm{k})})\big{]}}$
$\displaystyle\times\sum\limits_{\bm{k}_{2},\bm{q}\neq
0}\,V_{\bm{q}}\Big{[}\braket{\hat{a}_{(c,\bm{k}-\bm{q})}\hat{a}_{(c,\bm{k}_{2})}\hat{a}^{\dagger}_{(c,\bm{k}_{2}-\bm{q})}\hat{a}^{\dagger}_{(v,\bm{k})}}_{\omega}^{c}+\braket{\hat{a}_{(c,\bm{k}-\bm{q})}\hat{a}_{(v,\bm{k}_{2})}\hat{a}^{\dagger}_{(v,\bm{k}_{2}-\bm{q})}\hat{a}^{\dagger}_{(v,\bm{k})}}_{\omega}^{c}$
$\displaystyle\qquad\qquad\qquad-\braket{\hat{a}_{(c,\bm{k})}\hat{a}_{(v,\bm{k}_{2})}\hat{a}^{\dagger}_{(v,\bm{k}_{2}-\bm{q})}\hat{a}^{\dagger}_{(v,\bm{k}+\bm{q})}}_{\omega}^{c}-\braket{\hat{a}_{(c,\bm{k})}\hat{a}_{(c,\bm{k}_{2})}\hat{a}^{\dagger}_{(c,\bm{k}_{2}-\bm{q})}\hat{a}^{\dagger}_{(v,\bm{k}+\bm{q})}}_{\omega}^{c}\Big{]}.$
(92)

All the multiparticle effects are encoded in the second and third lines of
(92). Unfortunately, the exact form of two-particle correlations
$\braket{\hat{a}_{\bm{i}^{\prime}}\hat{a}_{\bm{j}^{\prime}}\hat{a}^{\dagger}_{\bm{l}^{\prime}}\hat{a}^{\dagger}_{\bm{s}^{\prime}}}_{\omega}^{c}$
is unknown. However, by means of the Heisenberg equation of motion (12) and
cluster expansion technique (14) we can derive the exact dynamical equation
for two-particle correlations, which schematically looks like

$\hbar\frac{d}{dt}\braket{2}^{c}_{t}=T_{2}\big{[}\braket{2}^{c}_{t}\big{]}+V_{1,2}\big{[}\braket{3}_{\braket{2}^{c}_{t},\braket{1}^{c}_{t}}\big{]}+V_{2,2}\big{[}\braket{3}^{c}_{t}\big{]}.$
(93)

Having obtained Eq. (93), we can simply obtain the result for the currently
considered (26) and further approximations by neglecting the corresponding
contributions from (93).

There is no the equation of motion for the correlations themselves. However,
one can apply the Heisenberg EOM (12) to four-operator expectation value and
after that by means of the cluster expansion (14) derives the similar equation
for two-particle correlations. From (16) we have:

$\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}=\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}-\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}+\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}.$
(94)

It is reasonable to consider the following expression:

$\displaystyle\Big{[}\frac{d}{dt}+i\big{[}\varepsilon_{\bm{i}^{\prime}_{1}}+\varepsilon_{\bm{i}^{\prime}_{2}}-\varepsilon_{\bm{j}^{\prime}_{2}}-\varepsilon_{\bm{j}^{\prime}_{1}}\big{]}\Big{]}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}=\Big{[}\frac{d}{dt}+i\big{[}\varepsilon_{\bm{i}^{\prime}_{1}}+\varepsilon_{\bm{i}^{\prime}_{2}}-\varepsilon_{\bm{j}^{\prime}_{2}}-\varepsilon_{\bm{j}^{\prime}_{1}}\big{]}\Big{]}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}$
$\displaystyle\qquad-\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}\Big{[}\frac{d}{dt}+i\big{[}\varepsilon_{\bm{i}^{\prime}_{2}}-\varepsilon_{\bm{j}^{\prime}_{2}}\big{]}\Big{]}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}-\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}\Big{[}\frac{d}{dt}+i\big{[}\varepsilon_{\bm{i}^{\prime}_{1}}-\varepsilon_{\bm{j}^{\prime}_{1}}\big{]}\Big{]}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}$
$\displaystyle\qquad+\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}\Big{[}\frac{d}{dt}+i\big{[}\varepsilon_{\bm{i}^{\prime}_{2}}-\varepsilon_{\bm{j}^{\prime}_{1}}\big{]}\Big{]}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}+\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}\Big{[}\frac{d}{dt}+i\big{[}\varepsilon_{\bm{i}^{\prime}_{1}}-\varepsilon_{\bm{j}^{\prime}_{2}}\big{]}\Big{]}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}.$
(95)

In (95) we know everything about the second and third lines [instead of
derivatives for one-particle correlations we can use (90)]. Thus, we have to
handle the following part:

$\Big{[}\frac{d}{dt}+i\big{[}\varepsilon_{\bm{i}^{\prime}_{1}}+\varepsilon_{\bm{i}^{\prime}_{2}}-\varepsilon_{\bm{j}^{\prime}_{2}}-\varepsilon_{\bm{j}^{\prime}_{1}}\big{]}\Big{]}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}.$
(96)

To understand what we have to obtain as output let us look back at Heisenberg
equation:

$\hbar\frac{d}{dt}\braket{A}_{t}+i\braket{\big{[}A,H_{el,k}\big{]}}_{t}=-i\braket{\big{[}A,H_{I}\big{]}}_{t}-i\braket{\big{[}A,V\big{]}}_{t}.$
(97)

In fact, the left-hand side of(97) coincides with (96) up to $\hbar$.
Therefore, we have to calculate now two expectation values of the commutators
$\big{[}A,H_{I}\big{]}$ and $\big{[}A,V\big{]}$. The first one is trivial(we
already found it for the considered type of quantity $A$):

$-i\braket{\big{[}\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}},H_{I}\big{]}}_{t}=-i\hbar\bigg{[}\sum\limits_{\bm{j}}\left(T^{4}_{\\{\bm{i}^{\prime}_{1}\bm{j}\\},\\{\bm{j}\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}+T^{4}_{\\{\bm{i}^{\prime}_{2}\bm{j}\\},\\{\bm{i}^{\prime}_{1}\bm{j},\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}\right)-\sum\limits_{\bm{i}}\left(T^{4}_{\\{\bm{i}\bm{j}^{\prime}_{2}\\},\\{\bm{i}^{\prime}_{1}\bm{i}^{\prime}_{2},\bm{i}\bm{j}^{\prime}_{1}\\}}+T^{4}_{\\{\bm{i}\bm{j}^{\prime}_{1}\\},\\{\bm{i}^{\prime}_{1}\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2}\bm{i}\\}}\right)\bigg{]},$
(98)

where $T^{4}$ can be expressed via correlations as follows:

$T^{4}_{\\{\bm{i}^{\prime}_{1}\bm{j}\\},\\{\bm{j}\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}=h^{ext}_{\bm{i}^{\prime}_{1}\bm{j}}\braket{\hat{a}_{\bm{j}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}=h^{ext}_{\bm{i}^{\prime}_{1}\bm{j}}\bigg{[}\braket{\hat{a}_{\bm{j}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}+\braket{\hat{a}_{\bm{j}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}-\braket{\hat{a}_{\bm{j}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\bigg{]},$
(99)

while to calculate the second average within the used formalism we have to
address to cluster expansion. For the two-particle case this gives:

$\displaystyle-i\braket{\big{[}A,V\big{]}}_{t}=-i\braket{\big{[}\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}},V\big{]}}_{t}=-i\frac{1}{4}\hbar\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}\Bigg{[}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}-\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}\hat{a}_{\bm{j}_{1}}\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}\Bigg{]}$
$\displaystyle=-i\frac{1}{4}\hbar\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}\Bigg{[}4\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{i}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}+4\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+4\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}-4\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-4\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}+4\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-4\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{2}}}_{t}^{c}-4\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{i}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-4\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}+4\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{i}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-4\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}}_{t}^{c}+4\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-4\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}}_{t}^{c}+4\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+2\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}}_{t}^{c}-2\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-4\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}+4\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-4\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}+4\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+2\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{i}_{2}}}_{t}^{c}\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}+2\braket{\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{2}}}_{t}^{c}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+2\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}\delta_{\bm{i}_{1},\bm{i}^{\prime}_{1}}-4\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\delta_{\bm{i}_{1},\bm{i}^{\prime}_{2}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+2\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{2}}\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\delta_{\bm{i}^{\prime}_{1},\bm{i}_{2}}-4\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\delta_{\bm{i}^{\prime}_{1},\bm{i}_{2}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-2\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{2}}\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\delta_{\bm{i}_{2},\bm{i}^{\prime}_{2}}-2\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\delta_{\bm{i}_{2},\bm{i}^{\prime}_{2}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-2\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}\delta_{\bm{j}_{1},\bm{j}^{\prime}_{1}}+2\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}\delta_{\bm{j}_{1},\bm{j}^{\prime}_{2}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-4\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\delta_{\bm{j}^{\prime}_{1},\bm{j}_{2}}+2\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\delta_{\bm{j}^{\prime}_{1},\bm{j}_{2}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+4\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\delta_{\bm{j}_{2},\bm{j}^{\prime}_{2}}-2\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\delta_{\bm{j}_{2},\bm{j}^{\prime}_{2}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-2\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\delta_{\bm{i}^{\prime}_{1},\bm{i}_{2}}-2\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\delta_{\bm{i}_{2},\bm{i}^{\prime}_{2}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad-2\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\delta_{\bm{j}^{\prime}_{1},\bm{j}_{2}}-2\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\delta_{\bm{j}_{2},\bm{j}^{\prime}_{2}}\Bigg{]}$
$\displaystyle=-i\hbar\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}CM_{2}[\bm{i}^{\prime}_{1},\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2},\bm{j}^{\prime}_{1},\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}].$
(100)

Thus, we know all the necessary expressions in terms of correlations for
equation (95). Let us present for completeness of the study, despite the
cumbersomeness of all expressions, all terms entering (95) via correlations:

$\displaystyle-\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}\Big{[}\frac{d}{dt}+i\big{[}\varepsilon_{\bm{i}^{\prime}_{2}}-\varepsilon_{\bm{j}^{\prime}_{2}}\big{]}\Big{]}$
$\displaystyle\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}=+i\bigg{[}\sum\limits_{\bm{j}}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}T^{2}_{\\{\bm{i}^{\prime}_{2}\bm{j}\\},\\{\bm{j},\bm{j}^{\prime}_{2}\\}}-\sum\limits_{\bm{i}}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}T^{2}_{\\{\bm{i}\bm{j}^{\prime}_{2}\\},\\{\bm{i}^{\prime}_{2},\bm{i}\\}}\bigg{]}$
(101)
$\displaystyle+i\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}\,\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}\,CM_{1}[\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2},\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}],$
$\displaystyle-\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}\Big{[}\frac{d}{dt}+i\big{[}\varepsilon_{\bm{i}^{\prime}_{1}}-\varepsilon_{\bm{j}^{\prime}_{1}}\big{]}\Big{]}$
$\displaystyle\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}=+i\bigg{[}\sum\limits_{\bm{j}}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}T^{2}_{\\{\bm{i}^{\prime}_{1}\bm{j}\\},\\{\bm{j},\bm{j}^{\prime}_{1}\\}}-\sum\limits_{\bm{i}}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}T^{2}_{\\{\bm{i}\bm{j}^{\prime}_{1}\\},\\{\bm{i}^{\prime}_{1},\bm{i}\\}}\bigg{]}$
(102)
$\displaystyle+i\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}\,\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}\,CM_{1}[\bm{i}^{\prime}_{1},\bm{j}^{\prime}_{1},\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}],$
$\displaystyle\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}\Big{[}\frac{d}{dt}+i\big{[}\varepsilon_{\bm{i}^{\prime}_{2}}-\varepsilon_{\bm{j}^{\prime}_{1}}\big{]}\Big{]}$
$\displaystyle\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}=-i\bigg{[}\sum\limits_{\bm{j}}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}T^{2}_{\\{\bm{i}^{\prime}_{2}\bm{j}\\},\\{\bm{j},\bm{j}^{\prime}_{1}\\}}-\sum\limits_{\bm{i}}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}T^{2}_{\\{\bm{i}\bm{j}^{\prime}_{1}\\},\\{\bm{i}^{\prime}_{2},\bm{i}\\}}\bigg{]}$
(103)
$\displaystyle-i\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}\,\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}\,CM_{1}[\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{1},\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}],$
$\displaystyle\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}\Big{[}\frac{d}{dt}+i\big{[}\varepsilon_{\bm{i}^{\prime}_{1}}-\varepsilon_{\bm{j}^{\prime}_{2}}\big{]}\Big{]}$
$\displaystyle\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}=-i\bigg{[}\sum\limits_{\bm{j}}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}T^{2}_{\\{\bm{i}^{\prime}_{1}\bm{j}\\},\\{\bm{j},\bm{j}^{\prime}_{2}\\}}-\sum\limits_{\bm{i}}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}T^{2}_{\\{\bm{i}\bm{j}^{\prime}_{2}\\},\\{\bm{i}^{\prime}_{1},\bm{i}\\}}\bigg{]}$
(104)
$\displaystyle-i\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}\,\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}\,CM_{1}[\bm{i}^{\prime}_{1},\bm{j}^{\prime}_{2},\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}],$

and the last one is

$\displaystyle\Big{[}\frac{d}{dt}+i\big{[}\varepsilon_{\bm{i}^{\prime}_{1}}+\varepsilon_{\bm{i}^{\prime}_{2}}-\varepsilon_{\bm{j}^{\prime}_{2}}-\varepsilon_{\bm{j}^{\prime}_{1}}\big{]}\Big{]}$
$\displaystyle\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}=\frac{1}{\hbar}\Bigg{[}-i\braket{\big{[}A,H_{I}\big{]}}_{t}-i\braket{\big{[}A,V\big{]}}_{t}\Bigg{]}$
$\displaystyle=$
$\displaystyle-i\bigg{[}\sum\limits_{\bm{j}}\left(T^{4}_{\\{\bm{i}^{\prime}_{1}\bm{j}\\},\\{\bm{j}\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}+T^{4}_{\\{\bm{i}^{\prime}_{2}\bm{j}\\},\\{\bm{i}^{\prime}_{1}\bm{j},\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}\right)-\sum\limits_{\bm{i}}\left(T^{4}_{\\{\bm{i}\bm{j}^{\prime}_{2}\\},\\{\bm{i}^{\prime}_{1}\bm{i}^{\prime}_{2},\bm{i}\bm{j}^{\prime}_{1}\\}}+T^{4}_{\\{\bm{i}\bm{j}^{\prime}_{1}\\},\\{\bm{i}^{\prime}_{1}\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2}\bm{i}\\}}\right)\bigg{]}$
(105)
$\displaystyle\qquad\qquad-i\bigg{[}\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}CM_{2}[\bm{i}^{\prime}_{1},\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2},\bm{j}^{\prime}_{1},\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}]\bigg{]}.$

Thus, for two-particle correlations one obtains:

$\displaystyle\Big{[}\frac{d}{dt}+i\big{[}\varepsilon_{\bm{i}^{\prime}_{1}}+\varepsilon_{\bm{i}^{\prime}_{2}}-\varepsilon_{\bm{j}^{\prime}_{2}}-\varepsilon_{\bm{j}^{\prime}_{1}}\big{]}\Big{]}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}$
$\displaystyle=$
$\displaystyle\,-i\Bigg{[}\sum\limits_{\bm{j}}\bigg{[}-\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}T^{2}_{\\{\bm{i}^{\prime}_{2}\bm{j}\\},\\{\bm{j},\bm{j}^{\prime}_{2}\\}}-\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}T^{2}_{\\{\bm{i}^{\prime}_{1}\bm{j}\\},\\{\bm{j},\bm{j}^{\prime}_{1}\\}}+\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}T^{2}_{\\{\bm{i}^{\prime}_{2}\bm{j}\\},\\{\bm{j},\bm{j}^{\prime}_{1}\\}}+\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}T^{2}_{\\{\bm{i}^{\prime}_{1}\bm{j}\\},\\{\bm{j},\bm{j}^{\prime}_{2}\\}}$
(106)
$\displaystyle\qquad\,\,+T^{4}_{\\{\bm{i}^{\prime}_{1}\bm{j}\\},\\{\bm{j}\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}+T^{4}_{\\{\bm{i}^{\prime}_{2}\bm{j}\\},\\{\bm{i}^{\prime}_{1}\bm{j},\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}\bigg{]}-\sum\limits_{\bm{i}}\bigg{[}-\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}T^{2}_{\\{\bm{i}\bm{j}^{\prime}_{2}\\},\\{\bm{i}^{\prime}_{2},\bm{i}\\}}-\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}T^{2}_{\\{\bm{i}\bm{j}^{\prime}_{1}\\},\\{\bm{i}^{\prime}_{1},\bm{i}\\}}$
$\displaystyle\qquad\,\,+\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}T^{2}_{\\{\bm{i}\bm{j}^{\prime}_{1}\\},\\{\bm{i}^{\prime}_{2},\bm{i}\\}}+\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}T^{2}_{\\{\bm{i}\bm{j}^{\prime}_{2}\\},\\{\bm{i}^{\prime}_{1},\bm{i}\\}}+T^{4}_{\\{\bm{i}\bm{j}^{\prime}_{2}\\},\\{\bm{i}^{\prime}_{1}\bm{i}^{\prime}_{2},\bm{i}\bm{j}^{\prime}_{1}\\}}+T^{4}_{\\{\bm{i}\bm{j}^{\prime}_{1}\\},\\{\bm{i}^{\prime}_{1}\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2}\bm{i}\\}}\bigg{]}\Bigg{]}$
$\displaystyle-i\Bigg{[}\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}\,\bigg{[}-\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}\,CM_{1}[\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2},\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}]-\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}\,CM_{1}[\bm{i}^{\prime}_{1},\bm{j}^{\prime}_{1},\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}]$
$\displaystyle+\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}\,CM_{1}[\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{1},\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}]+\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}\,CM_{1}[\bm{i}^{\prime}_{1},\bm{j}^{\prime}_{2},\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}]$
$\displaystyle\qquad\qquad\qquad\qquad\,+CM_{2}[\bm{i}^{\prime}_{1},\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2},\bm{j}^{\prime}_{1},\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}]\bigg{]}\Bigg{]}.$

One can simplify each term in (106):

$\displaystyle\sum\limits_{\bm{j}}\bigg{[}-\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}T^{2}_{\\{\bm{i}^{\prime}_{2}\bm{j}\\},\\{\bm{j},\bm{j}^{\prime}_{2}\\}}-\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}T^{2}_{\\{\bm{i}^{\prime}_{1}\bm{j}\\},\\{\bm{j},\bm{j}^{\prime}_{1}\\}}+\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}T^{2}_{\\{\bm{i}^{\prime}_{2}\bm{j}\\},\\{\bm{j},\bm{j}^{\prime}_{1}\\}}+\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}T^{2}_{\\{\bm{i}^{\prime}_{1}\bm{j}\\},\\{\bm{j},\bm{j}^{\prime}_{2}\\}}$
$\displaystyle\qquad\,\,+T^{4}_{\\{\bm{i}^{\prime}_{1}\bm{j}\\},\\{\bm{j}\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}+T^{4}_{\\{\bm{i}^{\prime}_{2}\bm{j}\\},\\{\bm{i}^{\prime}_{1}\bm{j},\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}\bigg{]}=\sum\limits_{\bm{j}}\bigg{[}h^{ext}_{\bm{i}^{\prime}_{1}\bm{j}}\braket{\hat{a}_{\bm{j}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}+h^{ext}_{\bm{i}^{\prime}_{2}\bm{j}}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{j}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}\bigg{]}$
$\displaystyle=\sum\limits_{\bm{j}}\bigg{[}T^{4}_{c,\\{\bm{i}^{\prime}_{1}\bm{j}\\},\\{\bm{j}\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}+T^{4}_{c,\\{\bm{i}^{\prime}_{2}\bm{j}\\},\\{\bm{i}^{\prime}_{1}\bm{j},\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}\bigg{]},$
(107)

where we introduce:

$\displaystyle
T^{4}_{c,\\{\bm{i}^{\prime}_{1}\bm{j}\\},\\{\bm{j}\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2}\bm{j}^{\prime}_{1}\\}}$
$\displaystyle=h^{ext}_{\bm{i}^{\prime}_{1}\bm{j}}\braket{\hat{a}_{\bm{j}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}.$
(108)

This definition is being expanded to arbitrary correlations. Similarly for
$i$-summation one can find:

$\displaystyle\sum\limits_{\bm{i}}\bigg{[}-\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}T^{2}_{\\{\bm{i}\bm{j}^{\prime}_{2}\\},\\{\bm{i}^{\prime}_{2},\bm{i}\\}}-\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}T^{2}_{\\{\bm{i}\bm{j}^{\prime}_{1}\\},\\{\bm{i}^{\prime}_{1},\bm{i}\\}}+\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}T^{2}_{\\{\bm{i}\bm{j}^{\prime}_{1}\\},\\{\bm{i}^{\prime}_{2},\bm{i}\\}}+\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}T^{2}_{\\{\bm{i}\bm{j}^{\prime}_{2}\\},\\{\bm{i}^{\prime}_{1},\bm{i}\\}}$
$\displaystyle\qquad\,\,+T^{4}_{\\{\bm{i}\bm{j}^{\prime}_{2}\\},\\{\bm{i}^{\prime}_{1}\bm{i}^{\prime}_{2},\bm{i}\bm{j}^{\prime}_{1}\\}}+T^{4}_{\\{\bm{i}\bm{j}^{\prime}_{1}\\},\\{\bm{i}^{\prime}_{1}\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2}\bm{i}\\}}\bigg{]}=\sum\limits_{\bm{i}}\bigg{[}T^{4}_{c,\\{\bm{i}\bm{j}^{\prime}_{1}\\},\\{\bm{i}^{\prime}_{1}\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2}\bm{i}\\}}+T^{4}_{c,\\{\bm{i}\bm{j}^{\prime}_{2}\\},\\{\bm{i}^{\prime}_{1}\bm{i}^{\prime}_{2},\bm{i}\bm{j}^{\prime}_{1}\\}}\bigg{]}.$
(109)

Finally, for a sum with the potential we obtain 15 unique contributions (as
previously, they are grouped by means of square brackets):

$\displaystyle\Bigg{[}\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}\,\bigg{[}-\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}\,CM_{1}[\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2},\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}]-\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}\,CM_{1}[\bm{i}^{\prime}_{1},\bm{j}^{\prime}_{1},\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}]$
$\displaystyle\qquad\,\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}^{c}_{t}\,CM_{1}[\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{1},\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}]+\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}\,CM_{1}[\bm{i}^{\prime}_{1},\bm{j}^{\prime}_{2},\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}]+CM_{2}[\bm{i}^{\prime}_{1},\bm{i}^{\prime}_{2},\bm{j}^{\prime}_{2},\bm{j}^{\prime}_{1},\bm{i}_{1},\bm{i}_{2},\bm{j}_{2},\bm{j}_{1}]\bigg{]}\Bigg{]}$
$\displaystyle=\sum\limits_{\begin{subarray}{c}\,\bm{i}_{1},\bm{i}_{2}\\\
\bm{j}_{1},\bm{j}_{2}\end{subarray}}v_{\bm{i}_{1},\bm{i}_{2},\bm{j}_{1},\bm{j}_{2}}\,\bigg{[}\big{[}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{i}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}-\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{2}}}_{t}^{c}\big{]}$
$\displaystyle\qquad\qquad\qquad\qquad-\big{[}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}\left(\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}}_{t}^{c}-\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\right)\big{]}$
$\displaystyle\qquad\qquad\qquad\qquad-\big{[}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\left(\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}}_{t}^{c}-\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\right)\big{]}$
$\displaystyle\qquad\qquad\qquad\qquad-\big{[}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}\left(\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{1}}}_{t}^{c}-\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\right)\big{]}$
$\displaystyle\qquad\qquad\qquad\qquad-\big{[}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\left(\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{1}}}_{t}^{c}-\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\right)\big{]}$
$\displaystyle\qquad\qquad\qquad\qquad-\frac{1}{2}\big{[}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}}_{t}^{c}\left(\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}_{\bm{j}_{2}}}_{t}^{c}-\braket{\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\right)\big{]}$
$\displaystyle\qquad\qquad\qquad\qquad+\frac{1}{2}\big{[}\braket{\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\left(\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{2}}}_{t}^{c}-\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}}_{t}^{c}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{i}_{2}}}_{t}^{c}\right)\big{]}$
$\displaystyle\qquad\qquad\qquad\qquad-\big{[}\braket{\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\delta_{\bm{i}^{\prime}_{1},\bm{i}_{2}}\big{]}+\big{[}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\delta_{\bm{i}_{2},\bm{i}^{\prime}_{2}}\big{]}$
$\displaystyle\qquad\qquad\qquad\qquad-\big{[}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\delta_{\bm{j}^{\prime}_{1},\bm{j}_{2}}\big{]}+\big{[}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{j}_{1}}}_{t}^{c}\delta_{\bm{j}_{2},\bm{j}^{\prime}_{2}}\big{]}$
$\displaystyle\qquad\qquad\qquad\qquad-\frac{1}{2}\big{[}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\delta_{\bm{i}^{\prime}_{1},\bm{i}_{2}}\big{]}+\frac{1}{2}\big{[}\braket{\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{j}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\delta_{\bm{i}_{2},\bm{i}^{\prime}_{2}}\big{]}$
$\displaystyle\qquad\qquad\qquad\qquad-\frac{1}{2}\big{[}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}}_{t}^{c}\delta_{\bm{j}^{\prime}_{1},\bm{j}_{2}}\big{]}+\frac{1}{2}\big{[}\braket{\hat{a}^{\dagger}_{\bm{i}_{2}}\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}^{\dagger}_{\bm{i}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}_{\bm{j}_{1}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}_{t}^{c}\delta_{\bm{j}_{2},\bm{j}^{\prime}_{2}}\big{]}\bigg{]}.$
(110)

Combining all the above obtained contributions, for each term from (93) we
find:

$\displaystyle\hbar\frac{d}{dt}\braket{2}^{c}_{t}$ $\displaystyle=$
$\displaystyle\hbar\frac{d}{dt}\braket{\hat{a}_{\bm{i}^{\prime}_{1}}\hat{a}_{\bm{i}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{2}}\hat{a}^{\dagger}_{\bm{j}^{\prime}_{1}}}^{c}_{t}$
(111) $\displaystyle T_{2}\big{[}\braket{2}^{c}_{t}\big{]}$ $\displaystyle=$