  One-particle Green's function methods can model molecular and solid spectra
at zero or non-zero temperatures. One-particle Green's 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 Green's 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
Green's 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.
