matrix $f_{\alpha\beta}({\bf p})= {\rm Tr}\,\varrho a^{\dagger}_{{\bf p}\beta}a_{{\bf p}\alpha}$, where $\varrho$ can be either an equilibrium or non-equilibrium statistical operator, $a^{\dagger}_{{\bf p}\alpha}$, $a_{{\bf p}\alpha}$ are the creation and annihilation operators, with index $\alpha$ running $2S+1$ values. Depending on the spin value, these operators meet the following bosonic commutation (integer spin) or fermionic anticommutation (half-integer spin) relations: $$\begin{gathered}
[a_{{\bf p}\alpha},a^{\dagger}_{{\bf p}'\alpha'}]_{B}=\delta_{{\bf p}{\bf p}'}\delta_{\alpha\alpha'}, \quad [a_{{\bf p}\alpha},a_{{\bf p}'\alpha'}]_{B}=0, \label{eq:2.1} \\
\{a_{{\bf p}\alpha},a^{\dagger}_{{\bf p}'\alpha'}\}_{F}=\delta_{{\bf p}{\bf p}'}\delta_{\alpha\alpha'}, \quad \{a_{{\bf p }\alpha},a_{{\bf p}'\alpha'}\}_{F}=0. \label{eq:2.2}\end{gathered}$$ Since below we study a homogeneous interacting Bose gas, we use the momentum ${\bf p}$ to specify the individual state of a particle. For the lattice models, one should consider a lattice site index instead of ${\bf p}$.

In the case of spin-1/2 system, the density matrix $f_{\alpha\beta}({\bf p})$, being a square matrix of the second order, can be written as a linear combination of the Pauli matrices $\sigma^{i}_{\alpha\beta}$ and unit matrix $I$ which form a basis for the vector space of $2\times 2$ matrices. The scalar part of the single-particle density matrix in such a decomposition defines the density of the system, whereas its vectorial part specifies three components of the magnetization vector. The latter quantity is induced by the spin 