It is well known that if $p$ is a prime, then $p$ can be written in the form $x^2 + n y^2$ under certain congruences conditions. For example,
\begin{align} p = x^2+y^2 &\Leftrightarrow p = 2 \text{ or } p = 1 \text{ }(\mathrm{mod}\text{ }4) \\ p = x^2+2y^2 &\Leftrightarrow p = 2 \text{ or } p = 1,3 \text{ }(\mathrm{mod}\text{ }8) \\ p = x^2+3y^2 &\Leftrightarrow p = 3 \text{ or } p = 1 \text{ }(\mathrm{mod}\text{ }3) \end{align}
At the same time, one can constructs a generating function for the number of solutions of the equation $x^2 + n y^2 = m$, using the theta function. If one define,
\begin{equation} \phi(q) = \sum_{m=-\infty}^\infty q^{m^2} \end{equation}
Then, we have the following generating function,
\begin{equation} \phi(q)\phi(q^n) = \sum_{m=0}^\infty q^m \text{ nb_solutions}(x^2 + n y^2 = m) \end{equation}
Question :
I do wonder whether it is possible to demonstrate the congruences above using an analytical method through the generating function written above, and the identities of the theta function. This would for example imply verifying that the coefficient is zero for primes satisfying a certain congruence, and verifying it is not for others.
Note :
Being a physics student, I would be really pleased by a solution which only involves symbolic manipulations, and identites of the theta function. I do not understand high mathematical theories concerning modular forms, elliptic functions or whatever. I have come to understand such things only by playing with it alone, and finding some information on the internet.