This is a sequel to my comments to the question which was too long to fit in another comment.
We have the formulas for $\vartheta_{3}^{10}(q), \vartheta_{3}^{12}(q)$ from Topics in Analytic Number Theory by Rademacher (famous for proving an infinite series formula to calculate the number of partitions of a positive integer) on page 198:
\begin{align}
\vartheta_{3}^{10}(q) &= 1 + \frac{4}{5}\left\{\sum_{n = 1}^{\infty}\frac{2n^{4}q^{n}}{1 + q^{2n}} + \sum_{n = 1}^{\infty}(-1)^{n - 1}\frac{(2n - 1)^{4}q^{2n - 1}}{1 - q^{2n - 1}}\right\} + \frac{2}{5}\vartheta_{3}^{2}(q)\vartheta_{2}^{4}(q)\vartheta_{4}^{4}(q)\tag{1}\\
\vartheta_{3}^{12}(q) &= 1 + 8\sum_{n = 1}^{\infty}\frac{n^{5}q^{n}}{1 - q^{2n}} - 8\sum_{n = 1}^{\infty}(-1)^{n}\frac{n^{5}q^{2n}}{1 - q^{2n}} + \vartheta_{2}^{4}(q)\vartheta_{3}^{4}(q)\vartheta_{4}^{4}(q)\tag{2}
\end{align}
Finding a general formula for $\vartheta_{3}^{k}(q)$ for even positive integer $k$ is a difficult problem but using the methods given in Rademacher's book it looks like it is possible to obtain such formulas at the cost of heavy symbolic manipulation for a specific $k$.
Update: I found one pattern in your formulas by using the substitution $x = K'(k)/K(k)$ so that when $x = 0$ then $k = 1$ and when $x = \infty$ then $k = 0$ and moreover $$\frac{dx}{dk} = -\frac{\pi}{2kk'^{2}K^{2}}$$ so that the integral of $\vartheta_{4}^{2n}(e^{-\pi x})/(1 + x^{2})$ is transformed into $$\int_{0}^{1}\left(\frac{2k'K}{\pi}\right)^{n}\frac{1}{K^{2} + K'^{2}}\frac{\pi}{2kk'^{2}}\,dk = \left(\frac{2}{\pi}\right)^{n - 1}\int_{0}^{1}\frac{k^{-1}k^{'(n - 2)}K^{n}}{K^{2} + K'^{2}}\,dk$$ and that explains (at least to some extent) the occurrence of $\dfrac{1}{\pi^{n - 1}}$ in your formulas.
Next it is easy to prove one of the formulas in $(14)$. We have $$\vartheta_{2}^{2}\vartheta_{4}^{2} = kk'(2K/\pi)^{2}$$ and hence $$\int_{0}^{\infty}\vartheta_{2}^{2}\vartheta_{4}^{2}\,dx = \int_{0}^{1}kk'\cdot\frac{4K^{2}}{\pi^{2}}\cdot\frac{\pi}{2kk'^{2}K^{2}}\,dk = \frac{2}{\pi}\int_{0}^{1}\frac{dk}{\sqrt{1 - k^{2}}} = 1$$ I wonder if similar technique can be applied to prove other formulas.
If $q = e^{-\pi x}$ then $dx = -\dfrac{dq}{\pi q}$ and interval $(0, \infty)$ changes to $(0, 1)$ and hence we can express the first integral of $(14)$ as $$\frac{1}{\pi}\int_{0}^{1}\vartheta_{2}^{4}(q)\vartheta_{4}^{2}(q)\,\frac{dq}{q} = \frac{16}{\pi}\int_{0}^{1}\psi^{4}(q^{2})\phi^{2}(-q)\,dq$$ Next $$\psi^{4}(q^{2}) = \sum_{n = 0}^{\infty}\frac{(2n + 1)q^{2n}}{1 - q^{4n + 2}}, \phi^{2}(-q) = 1 + 4\sum_{n = 1}^{\infty}\frac{(-1)^{n}q^{n}}{1 + q^{2n}}$$ I wonder if you can utilize the above Lambert series to prove that the desired integral is equal to $1$. It appears that if we express the integrand as a Lambert series then it can also be expressed as the logarithmic derivative of some product of theta functions and the integral can be evaluated. See this paper regarding some integrals related to theta functions (all of it was given by Ramanujan in his lost notebook).
The key is the link between the Dedekind eta function and elliptic integrals.
Let $\tau$ be purely imaginary and in upper half of complex plane and let $$q=\exp(2\pi i\tau) \in(0,1)$$ be the corresponding nome. Consider the elliptic modulus $k\in(0,1)$ corresponding to nome $q$ given in terms of $q$ via Jacobi theta functions $$k=\frac{\vartheta_{2}^{2}(q)}{\vartheta _{3}^{2}(q)},\,\vartheta_{2}(q)=\sum_{n\in\mathbb {Z}} q^{(n+(1/2))^{2}},\,\vartheta _{3}(q)=\sum_{n\in\mathbb {Z}} q^{n^2}\tag{1}$$ Let $k'=\sqrt {1-k^2}$ and we further define elliptic integrals $$K=K(k) =\int_{0}^{\pi/2}\frac{dx} {\sqrt{1-k^2\sin^2 x}}, \, K'=K(k') \tag{2}$$ The circle of these definitions is finally completed by the formula $$\frac{K'} {K} =-2i\tau\tag{3}$$ Let $\tau'$ be another purely imaginary number in upper half of the complex plane such that $$\frac{\tau'} {\tau} =r\in\mathbb {Q} ^{+} \tag{4}$$ Let the corresponding nome be $q'=\exp(2\pi i\tau') $ and the elliptic moduli be $l, l'=\sqrt{1-l^2}$ and the elliptic integrals based on these moduli be denoted by $L, L'$. Then from the relation $\tau'=r\tau$ we get via $(3)$ the modular equation $$\frac{L'} {L} =r\frac{K'} {K}, r\in\mathbb {Q} ^{+} \tag{5}$$ Under these circumstances Jacobi proved using the transformation of elliptic integrals that the relation between moduli $k, l$ is algebraic and the ratio $K/L$ is an algebraic function of $k, l $.
The Dedekind's eta function is related to elliptic integrals via the relation $$\eta(\tau) =q^{1/24}\prod_{n=1}^{\infty} (1-q^n)=2^{-1/6}\sqrt{\frac{2K}{\pi}}k^{1/12}k'^{1/3}\tag{6}$$ Now let $\tau=i/2$ so that $q=e^{-\pi} $ and then from $(3)$ we have $K=K'$ so that $k=k'=1/\sqrt{2}$ and it is well known that for this value of $k$ we have $$K(k) =\frac{\Gamma^{2}(1/4)} {4\sqrt{\pi}} \tag{7}$$ From $(6)$ it now follows that $\eta(\tau) =\eta(i/2)$ is an algebraic multiple of $\Gamma (1/4)\pi^{-3/4}$.
Let $\tau'=ri, r\in \mathbb {Q} ^{+} $ so that $\tau'/\tau=2r$ is a positive rational number. As noted above if $l, L$ correspond to $\tau'$ then the relation between $l$ and $k=1 /\sqrt{2}$ is algebraic so that $l$ is an algebraic number and the ratio $K/L$ is an algebraic function of $k, l $ and thus $K/L$ is also an algebraic number. Thus from equation $(6)$ it follows that $\eta(ri) $ is an algebraic multiple of $\Gamma (1/4)\pi^{-3/4}$.
More generally it can be proved that if $r$ is a positive rational number then the value of $\eta(i\sqrt{r}) $ can be expressed in terms of values of Gamma function at rational points and $\pi$ and certain algebraic numbers.
Also let me complete the link between $\Gamma (1/4)$ and elliptic integrals starting with your approach. We have $$\frac{\Gamma ^2(1/4)}{2\sqrt{\pi}}=\int_{0}^{\pi}\frac{dx}{\sqrt{\sin^4 x+\cos^4 x}}=\int_{0}^{\pi}\frac{dx}{\sqrt{1-2\sin^2 x\cos^2 x}}$$ and the integral can further be written as $$\int_{0}^{\pi}\frac{dx}{\sqrt{1-(1/2)\sin^2 2x}}$$ Putting $2x =t$ we can see that it reduces to $$\frac{1}{2}\int_{0}^{2\pi}\frac{dt}{\sqrt{1-(1/2)\sin^2 t}}=2\int_{0}^{\pi/2}\frac{dx}{\sqrt{1-(1/2)\sin^2 x}}=2K(1/\sqrt{2})$$ and we are done.
Best Answer
You can reduce your problem to $\sum_{n\in \mathbb{Z}} e^{-n^2\pi}$ using the Fourier transform and Poisson summation formula. First a few observations: setting $f(x) = e^{-\pi x^2}$, we have $\hat{f}(k)=f(k)$, and also $x^2 f(x)= \frac{1}{4\pi^2}(fââ(x)+2\pi f(x))$. By the Poisson summation formula,
$$\sum_{n\in\mathbb{Z}} n^2 f(n) = \sum_{n\in\mathbb{Z}} \mathcal{F}(x^2f(x))(n)= \sum_{n\in\mathbb{Z}} \frac{1}{4\pi^2}[(2\pi i n)^2\hat{f}(n)+2\pi \hat{f}(n)] $$
Moving the first term to the other side shows that
$$ \sum_{n\in\mathbb{Z}}n^2 f(n) = \frac{1}{4\pi}\sum _{n\in\mathbb{Z}} f(n) $$