Euler-Maclaurin summation for $e^{-x^2}$

I want to approximate the sum $$\sum_{k=0}^\infty e^{-k^2}$$ using the Euler-Maclaurin formula $$\sum_{k=0}^\infty f(k) = \int_0^\infty f(x) \, dx + \frac{1}{2}(f(0) + f(\infty)) + \frac{1}{12}(f'(\infty) - f'(0)) - \frac{1}{720}(f'''(\infty) - f'''(0)) + \ldots$$

where $f(\infty)$ means $\lim_{x \to \infty} f(x)$. But putting $f(x) = e^{-x^2}$, it's not too hard to see that all derivatives of odd order are an odd polynomial times $f$, implying they vanish at $x = 0$. But this means all terms of the Euler-Maclaurin formula vanish except for the first two! But I know this is wrong because I evaluated the sum numerically in Mathematica.

What's wrong with my calculations?

Nothing is wrong with the details of the computation, but you are treating the Euler-Maclaurian formula as an equality, when it isn't.

Before continuing, I will switch the sum a little bit: You are working with $$S:=\sum_{k=0}^{\infty} e^{-k^2}.$$ The formulas are a little cleaner in terms of $$T:=\sum_{k=-\infty}^{\infty} e^{-k^2}.$$ These are related by $T=2S-1$, so it is easy to switch between them.

The Euler-Macluarin formula with remainder term, combined with your correct computation that $f^{(k)}(x)$ goes to $0$ as $x \to \pm \infty$, tells us that we have $$\sum_{k=-\infty}^{\infty} e^{-k^2} = \int_{x=-\infty}^{\infty} e^{-x^2} dx + \int_{-\infty}^{\infty} \frac{B_N(x-\lfloor x \rfloor)}{N!} \frac{d^N e^{-x^2}}{(dx)^N} dx$$ for any positive integer $N$. Here $B_N$ is the $N$-th Bernouli polynomial and $\lfloor x \rfloor$ is $x$ rounded down to an integer. For some functions, the remainder $\int \frac{B_N(x-\lfloor x \rfloor)}{N!} \frac{d^N f}{(dx)^N} dx$ goes to $0$ as $N \to \infty$, but it often doesn't, and it doesn't in your case.

To give some other examples, if you compute the asymptotics of $\sum_{k \leq M} \frac{1}{k}$ using Euler-Maclaurin, the Euler-Mascheroni constant occurs as the limit of the remainder term. If you derive Stirling's formula by Euler-Maclaurin summation of $\sum \log k$, then the $\log (2 \pi)$ occurs as the limit of the remainder term.

There is a really good discussion of this in Chapter 9 of Concrete Mathematics, by Graham, Knuth and Patashnik. You might particularly enjoy the fourth example in Chapter 9.6, where they use Euler Maclaurin summation to show that, for any $N$, we have $$\sum_{k=-\infty}^{\infty} e^{-k^2/t} = \sqrt{\pi t} + O(t^{-N}) \ \mbox{as}\ t \to \infty.$$

Many examples of Euler-Maclaurin summation are actually harmonic sums and can be treated by Mellin transform methods.

In the present case put $$S(x) = \sum_{k\ge 1} e^{-x^2 k^2}$$ with so that we are interested in $S(1/\sqrt{t})$ as $t\rightarrow\infty.$

This sum can be evaluated by inverting its Mellin transform.

Recall the harmonic sum identity $$\mathfrak{M}\left(\sum_{k\ge 1} \lambda_k g(\mu_k x);s\right) = \left(\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} \right) g^*(s)$$ where $g^*(s)$ is the Mellin transform of $g(x).$

In the present case we have $$\lambda_k = 1, \quad \mu_k = k \quad \text{and} \quad g(x) = e^{-x^2}.$$

We need the Mellin transform $g^*(s)$ of $g(x)$ which is $$\int_0^\infty e^{-x^2} x^{s-1} dx.$$ Use the substitution $x^2 = u$ so that $2x \; dx = du$ to get $$\int_0^\infty e^{-u} u^{1/2s-1/2} \frac{1/2 \; du}{\sqrt{u}} = \frac{1}{2} \int_0^\infty e^{-u} u^{1/2s-1} du = \frac{1}{2} \Gamma(s/2).$$

The fundamental strip of this Mellin transform is $\langle 0,\infty\rangle.$

It follows that the Mellin transform $Q(s)$ of the harmonic sum $S(x)$ is given by

$$Q(s) = \frac{1}{2} \Gamma(s/2) \zeta(s) \quad\text{because}\quad \sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} = \sum_{k\ge 1} \frac{1}{k^s} = \zeta(s)$$ for $\Re(s) > 1.$

The Mellin inversion integral for this transform is $$\frac{1}{2\pi i} \int_{3/2-i\infty}^{3/2+i\infty} Q(s)/x^s ds$$ which we evaluate by shifting it to the left for an expansion about zero (recall that as $t\rightarrow\infty$ we have $1/\sqrt{t}\rightarrow 0.$)

Observe that the poles are at $s=1$ from the zeta function term and at the non-positive even integers from the gamma function term. However all of the latter except the one at zero are canceled by the trivial zeros of the zeta function term, leaving only the poles at $s=0$ and $s=1.$

For these two we have $$\mathrm{Res}\left(Q(s)/x^s; s=1\right) = \frac{1}{2} \Gamma(1/2)\frac{1}{x} = \frac{\sqrt{\pi}}{2x}$$ and $$\mathrm{Res}\left(Q(s)/x^s; s=0\right) = \frac{1}{2} \times 2 \times -\frac{1}{2} = -\frac{1}{2}.$$

It follows that as $t\rightarrow\infty$ we have $$S(1/\sqrt{t}) \sim \frac{1}{2} \sqrt{\pi t} - \frac{1}{2}$$ and in particular $$2S(1/\sqrt{t})+1 = \sum_{k=-\infty}^\infty e^{-k^2/t} \sim \sqrt{\pi t}.$$

As for the error term if we have shifted the integral to the line $\Re(s) = -q/2$ with $q>1$ and $q$ odd we have for the norm of the zeta function term on the line $-q/2+iv$ the bound $|v|^{1/2+q/2}$ and the gamma function term decays exponentially in $v$ along vertical lines and in $q$ at the values at $-q/2$, so the error term decays exponentially also.