One-particle Greens function methods can model molecular and solid spectra at zero or non-zero temperatures. One-particle Greens functions directly provide electronic energies and one-particle properties, such as dipole moment. However, the evaluation of two-particle properties, such as $langle{S^2}rangle$ and $langle{N^2}rangle$ can be challenging, because they require a solution of the computationally expensive Bethe--Salpeter equation to find two-particle Greens functions. We demonstrate that the solution of the Bethe--Salpeter equation can be complitely avoided. Applying the thermodynamic Hellmann--Feynman theorem to self-consistent one-particle Greens function methods, we derive expressions for two-particle density matrices in a general case and provide explicit expressions for GF2 and GW methods. Such density matrices can be decomposed into an antisymmetrized product of correlated one-electron density matrices and the two-particle electronic cumulant of the density matrix. Cumulant expressions reveal a deviation from ensemble representability for GW, explaining its known deficiencies. We analyze the temperature dependence of $langle{S^2}rangle$ and $langle{N^2}rangle$ for a set of small closed-shell systems. Interestingly, both GF2 and GW show a non-zero spin contamination and a non-zero fluctuation of the number of particles for closed-shell systems at the zero-temperature limit.