Use the fact that the approximation in the paper is equal to the inner product of e^{-x^2} with a Dirac comb [1]. The amazing thing about the Gaussian function and the Dirac comb is that they're both preserved by the Fourier transform. So apply the Fourier transform to both of them, observe that the Fourier transform is unitary (essentially a rotation), and therefore doesn't affect the inner product; then expand the inner product of the Fourier transforms.
Essentially, it's the same proof, but not in the language of Theta functions. In fact, I'd argue it's a better proof, because it generalises.
Generalisation: This technique applies to all Riemann sums, as long as you can compute the Fourier transform of the function. The thinner the tails of the function's Fourier transform, the faster its Riemann sums will converge.