The Euler/Jacobi theta function (using the notation of this question) is $\vartheta_3(\tau) := \sum_{n\in \mathbb Z} q^{n^2}$ where $q = e^{2\pi i\tau}$ is the nome. The square $(\vartheta_3(\tau))^2$ then has the series representation $\sum_{k=0}^\infty r_2(k) q^k$, where $r_2(k)$ is the number of pairs $(m,n) \in \mathbb Z^2$ s.t. the number $k$ can be written as the sum of two squares $m^2+n^2$. This function has a decently compact formula, explained beautifully in this 3b1b video. Regardless, from the very fundamental definition of the area of a circle of radius $R$ having area $\pi R^2$, we see that $\frac 1{N^2} \sum_{k=0}^N r_2(k) \to \pi$ as $N\to\infty$.
Because $r_2(k)$ are the Fourier coefficients of $(\vartheta_3(\tau))^2$, and we're doing a sum of the $r_2(k)$ from $0$ to $N$ (really $-N$ to $N$ since the negative power Fourier coefficients are all $0$ in this case), this screams "Poisson summation formula!" to me. The sum $\frac 1{N^2} \sum_{k=0}^N r_2(k)$ is also highly reminiscent of a Riemann sum expression, but I can't get things to quite work out. Of course, the well known Gaussian integral is $\int_{\mathbb R} e^{-x^2} d x = \sqrt \pi$, and this quadratic exponential is really looking similar to the quadratic exponential in $\vartheta_3$. Moreover, the square root in the Gaussian integral sort of looks like it comes from the fact that we're studying $(\vartheta_3(\tau))^2$.
For these reasons, it is looking like there might be a proof of the value of the Gaussian integral using these building blocks. However, I can't put the pieces together, with some of the difficulties being the fact that $\frac 1{N^2} \sum_{k=0}^N r_2(k)$ has that $\frac 1{N^2}$ factor out front, and the fact that although the Poisson summation formula looks very promising on the Fourier/frequency domain, it doesn't seem to be promising on the time domain.
There's a similar question here Computing the Gaussian integral with Fourier methods?, but does not mention about the "number of sum of squares" function at al.
