P.444 of volume 3 of Gauss's Nachlass includes the following identity for the logarithm of $\vartheta_3(x)$:
$$\mathbb{log}(1+2x+2x^4+2x^9+\cdots) = \frac{2x}{1+x}+\frac{2x^3}{3(1+x^3)}+\frac{2x^5}{5(1+x^5)}+\cdots,$$
and he apparently somehow used this identity to derive the infinite series expression for $\vartheta^4_3(x)$, which is the generating function for the sum of four squares function $r_4(n)$ (this is Jacobi's four squares theorem). As can be seen from the series, it is a kind of Lambert series $\sum_{n=1}^{\infty}a_n\frac{q^n}{1-q^n}$ in which $q=-x$ and $a_{2k+1} = -\frac{2}{2k+1}, a_{2k}=0$.
When I tried to prove Gauss's identity (in order to help myself answer this question), I thought it makes sense to use the Jacobi triple product identity to express $\vartheta_3(x)$ and then to use certain well known relations between the logarithm of infinite products and Lambert series. This led me to the following attempt to prove Gauss's identity:
$$\vartheta_3(q) = \prod_{n=1}^{\infty}(1-q^{2n})\prod_{n=1}^{\infty}(1+q^{2n-1})^2 = \prod_{n=1}^{\infty}(1-q^{2n})(\frac{\prod_{n=1}^{\infty}(1+q^n)}{\prod_{n=1}^{\infty}(1+q^{2n})})^2$$
and by converting the logarithm of infinite product to an infinite sum of logarithms one gets the following linear combination of certain Lambert series:
$$\mathbb{log}(\vartheta_3(q)) = -\sum_{n=1}^{\infty}\frac{1}{n}\frac{q^{2n}}{1-q^{2n}}+2(\sum_{n=1}^{\infty}\frac{(-1)^{n-1}}{n}\frac{q^n}{1-q^n}-\sum_{n=1}^{\infty}\frac{(-1)^{n-1}}{n}\frac{q^{2n}}{1-q^{2n}})$$
and after some simplifications I arrived at the following expression:
$$\mathbb{log}(\vartheta_3(q))= \sum_{n=1}^{\infty}\frac{(-1)^{n-1}}{n}(\frac{2q^n+(-1)^{n}q^{2n}}{1-q^{2n}})$$
I checked the first terms of my expression and Gauss's identity and they coincide, so I guess I did not make mistakes, but I cannot see a way of transforming my linear combination of Lambert series into Gauss form; there is a "trick" here that I have not figured out yet, despite repeating attempts to show equivalence to Gauss's identity.
So my question is, in general, how to prove Gauss's identity for $\mathbb{log}(\vartheta_3(q))$? and in particular, how to move from my final expression to Gauss's form (it will be much appreciated by me if someone will continue my attempt)?
This is not an answer and I've got to run. Take $n$ odd; then $$ \frac{2q^n - q^{2n}}{1-q^{2n}} = \frac{2q^n - 2q^{2n}}{1-q^{2n}} + \frac{q^{2n}}{1-q^{2n}} = 2 \frac{q^n(1-q^n)}{ 1-q^{2n}} + \frac{q^{2n}}{1-q^{2n}} = 2 \frac{q^n}{ 1+q^n} + \frac{q^{2n}}{1-q^{2n}}. $$ Looks like this may work if you rearrange some terms ($n$ and $2n$ etc.)