Interpretation of a certain general theorem used by Gauss in his work on theta functions. 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.
 A: 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.
A: Let's use the variables $q, z$ with $q=x, y=e^{2iz}$ so that $$P(x, y) =\vartheta_3(z,q),Q(x,y)=\vartheta_4(z,q),R(x,y)=\vartheta_2(z,q)$$ and we can now transcribe the general theorem of Gauss as $$\vartheta_3(z+w,q)\vartheta_3(z-w,q)=\vartheta_3(2z,q^2)\vartheta_3(2w,q^2)+\vartheta_2(2z,q^2)\vartheta_2(2w,q^2)$$ (with $t=e^{2iw}$) as an identity between Jacobi theta functions.
This is one of the most fundamental identities between theta functions and almost all algebraic relations between theta functions can be derived using this. You may have a look at this paper at arXiv for some identities derived via this general theorem of Gauss
The proof of the same can be given by considering the ratio of left and right sides and showing that it is a doubly periodic functions with no poles. And thus is a constant. It requires some effort to show that the constant is $1$ but can be shown with some algebraic manipulation on the series corresponding to these functions with $z=0,w=0$.

At the moment I don't have a direct algebraic proof of the above identity and will need to check Jacobi Fundamenta Nova to see if Jacobi provided any such proof. Also as you have remarked in your question, Jacobi Theta functions are defined for all complex $z, q$ with $|q|<1$.
