While reviewing the integral test, I remembered a formula for bounding sums the test applies to (non-negative, monotone decreasing $f(x)$): $$\sum_{N+1}^M f(n) \leq \int_N^M f(x)dx$$
Specifically, for $S_N = \sum_{k=1}^N e^{-k^2}$, $$S_\infty = \sum_{k=1}^\infty e^{-k^2} \leq \int_0^\infty e^{-x^2}dx = \frac{\sqrt{\pi}}{2}\operatorname{erf}(\infty) = \frac{\sqrt{\pi}}{2}$$
I computed $\frac{\sqrt{\pi}}{2}\operatorname{erf}(N) - S_N$ for large $N$ ($10^6$), and it seemed to converge to ~$0.49990832303943183$. This implies that $$\sum_{k=1}^\infty e^{-k^2} = \frac{\sqrt{\pi}}{2} - C,$$ where $C$ is ~$0.49990832303943183$. ries gives a few ideas for $C$, and also for $S_{10^6}$.
Can these estimates get any better? Or, fingers crossed, is there a closed form for this sum? Specifically without Jacobi Theta functions.