  The machinery of computing vacuum expectation values of a time-ordered
sequence of position operators of the simple harmonic oscillator is already
well established. It rests on a Wick theorem, which enables one to decompose
such a quantity in terms of products of \emph{pairwise contractions}, which are
vacuum expectation values of a time-ordered sequence of position operators
taken two at a time. This result naturally leads to a diagrammatic approach of
computing such correlators, and is already well known in the form of Feynman
diagrams. We generalise this setup to encompass expectation values of a general
ordered sequence of position operators (Wightman sequences) in general density
matrices of the simple harmonic oscillator. A Wick theorem is first developed
for this situation and consequently a diagrammatics is laid down. Wightman
correlators in general density matrices of the \emph{anharmonic} oscillator are
also analysed and a diagrammatic formalism is developed for them too.
