I'm trying to understand the meaning of a general proposition stated by Gauss in a posthomous paper (this paper is in pp. 470-481 of volume 3 of Gauss's werke) on theta functions, a proposition which seems to serve as a guiding and organizing principle of the vast amount of relations among theta functions that he found.
Gauss's notation and definitions
Denote by $P(x,y),Q(x,y),R(x,y)$ the following functions:
$$P(x,y)=1+x(y+\frac{1}{y})+x^4(y^2+\frac{1}{y^2})+x^9(y^3+\frac{1}{y^3})+...$$ $$Q(x,y)= 1-x(y+\frac{1}{y})+x^4(y^2+\frac{1}{y^2})-x^9(y^3+\frac{1}{y^3})+...$$ $$R(x,y)=x^{\frac{1}{4}}(y^{\frac{1}{2}}+y^{-\frac{1}{2}})+x^{\frac{9}{4}}(y^{\frac{3}{2}}+y^{-\frac{3}{2}})+x^{\frac{25}{4}}(y^{\frac{5}{2}}+y^{-\frac{5}{2}})+...$$
These functions include Jacobi theta functions in their usual meaning as special cases; if $y$ is a complex number whose absolute value is $1$, and $z$ is defined to be a real number such that $y = e^{2iz}$, then we have:
$$P(x,y)=1+2cos(2z)x+2cos(4z)x^4+2cos(6z)x^9+...=\vartheta_3(z,x)$$
which follows from the identity $cos(2nz)= \frac{e^{2inz}+e^{-2inz}}{2}$. In paticular, we have:
$$P(x,1)=1+2x+2x^4+2x^9+...=\vartheta_3(0,x)$$, So one can understand $P(x,y),Q(x,y),R(x,y)$ as a generalization of Jacobi theta function $\vartheta(z,x)$ from purely real $z$ to a complex $z$ (non-zero imaginary part of z), so that $|y| \ne 1$.
Remark: I'm not very familiar with Jacobi's publications, so it's quite possible that Jacobi's original definition of his theta functions includes also the case when $z$ is complex, so Gauss's functions $P(x,y),Q(x,y),R(x,y)$ are nothing else than simply Jacobi's theta functions with different notation.
Gauss's theorem
On August 6, 1827, Gauss stated the following "general theorem":
$$P(x,ty)\cdot P(x,\frac{y}{t}) = P(x^2,t^2)P(x^2,y^2) + R(x^2,t^2)R(x^2,y^2) $$
and then goes on to derive a multitude of relations from it.
For more comprehensive background on this question, please look at the answer to HSM stackexchange post https://hsm.stackexchange.com/questions/6256/did-gauss-know-jacobis-four-squares-theorem.
Therefore, i'd like to know how to interpret the general theorem stated by Gauss.
The definition of the Gauss theta functions can be written as
$$ P(x,y) = \sum_{n\in\mathbb{Z}} x^{n^2}y^n,\;\; R(x,y) = \sum_{n\in\mathbb{Z}+\frac12} x^{n^2}y^n. \tag{1} $$
Now consider the product of two theta functions
$$ S := P(x,ty)\cdot P(x,y/t) = \left(\sum_{n\in\mathbb{Z}} x^{n^2}(ty)^n\right) \! \left(\sum_{m\in\mathbb{Z}} x^{m^2}(y/t)^m\right). \tag{2} $$
This can be rewritten as a double sum
$$ S = \sum_{n,m\in\mathbb{Z}} x^{n^2+m^2} y^{n+m}t^{n-m}. \tag{3} $$
Rewrite this using new variables
$$ j = \frac{n+m}2,\;\; k = \frac{n-m}2 \;\; \text{ where } \;\; n = j+k,\;\; m = j-k \tag{4} $$
to get
$$ S = \sum_{n,m\in\mathbb{Z}} x^{2(j^2+k^2)} y^{2j}t^{2k}. \tag{5} $$
The double sum $\,S\,$ splits into two cases. One is $\,S_0\,$ where $\,n,m\,$ have the same parity with $\,j,k\in\mathbb{Z}.\,$ The other is $\,S_1\,$ where $\,n,m\,$ have different parity with $\,j,k\in\mathbb{Z}+\frac12.\,$ Rewrite the sums as products
$$ S_0 = \sum_{j,k\in\mathbb{Z}} (x^2)^{k^2}(t^2)^k \cdot (x^2)^{j^2}(y^2)^j = P(x^2,t^2)P(x^2,y^2) \tag{6} $$
and
$$ S_1 = \sum_{j,k\in\mathbb{Z}+\frac12} (x^2)^{k^2}(t^2)^k \cdot (x^2)^{j^2}(y^2)^j = R(x^2,t^2)R(x^2,y^2). \tag{7} $$
The end result is
$$ S = S_0+S_1 = P(x^2,t^2)P(x^2,y^2) + R(x^2,t^2)R(x^2,y^2). \tag{8} $$
I think that this is similar to what Gauss' original proof was but I have no way to know that. This approach must be very old.