More precisely, $\omega^{2D}_p \approx \sqrt{2\pi ne^2/m} ~ \sqrt{q}$. There is a limit on the validity of these formula only in large $q$.

Thanks for your comments. Here is a reference: rmf.smf.mx/pdf/rmf/39/4/39_4_640.pdf

In quantum mechanics, in order to represent an entity by a matrix, the prerequisite is that this entity must act in the Hilbert space. A Grassmannn variable is not any operator and it does not operate in the Hilbert space. So, it is impossible to put it as a matrix.

This does not answer the question, because still Grassmann variables remain Grassmann variables, which are not any kind of operators and can't be reduced to a set of real numbers!

Definitely you can represent these Fermionic operators by matrices, which is standard quantum mechamics! But you can't do this for Grassmann variables. This is what I claimed in my answer.

In a scalar field like what you wrote, the particles created by $a^{\dagger}_n$ are doomed to be annihilated by the quartic term $\phi^4$. This is no more than the usual phonon case: the number of phonon can not be conserved in the presence of anharmonic terms.