Why do these two sums define functions that are asymptotically so close in value?

Question

Let $$f_n(\epsilon)=(2n+1)e^{-n(n+1)\epsilon},\qquad Z_\text{odd}(\epsilon)=\sum_{p\ge0}f_{2p+1}(\epsilon),\qquad Z_\text{even}(\epsilon)=\sum_{p\ge0}f_{2p}(\epsilon).$$ (These functions are relevant to describing rotation of diatomic molecules in statistical physics.)

When $\epsilon\searrow0$, these functions are extremly close in value; numerically: $$Z_\text{even}(0.1)-Z_\text{odd}(0.1)\simeq 3\,10^{-9},\qquad Z_\text{even}(0.05)-Z_\text{odd}(0.05)\simeq 2\,10^{-19}.$$ Is there a reason for this? What is an asymptotic equivalent of the difference?

I tried to obtain a small $\epsilon$ expansion using the Euler-MacLaurin expansion and Mathematica, and I found that they both have the same expansion. Notice that Euler-MacLaurin involves derivatives of $f$ at 0 for $Z_\text{odd}$ and derivatives of $f$ at 1 for $Z_\text{even}$, so it is not obvious how this should lead to the same expansion.

[Paragraph added in an edit] In fact, to compare the expansions, it is easier to write $Z_\text{even}(\epsilon)-Z_\text{odd}(\epsilon)=2Z_\text{even}(\epsilon)-Z_\text{total}(\epsilon)$, where $Z_\text{total}(\epsilon)=\sum_{n\ge0} f_n(\epsilon)$. Then, applying Euler-MacLaurin to order $m$ to $Z_\text{even}$ and $Z_\text{total}$, one obtains $$Z_\text{even}(\epsilon)-Z_\text{odd}(\epsilon)=\frac{f_0(\epsilon)}2+\sum_{p=0}^m \frac{b_{2p}(2^{2p}-1)}{(2p)!}f^{(2p-1)}(\epsilon)+R_m(\epsilon).$$ The first terms in this summation are $\frac12$, $\frac14(-2+\epsilon)$, $\frac1{48}(-12\epsilon+12\epsilon^2-\epsilon^3)$, $\frac1{480}(-120\epsilon^2+180\epsilon^3-30\epsilon^4+\epsilon^5)$, $\frac{17}{80\,640}(-1680 \epsilon ^3+3360 \epsilon ^4-840 \epsilon ^5+56 \epsilon ^6-\epsilon ^7)$. Notice how each new term cancels the lowest power of $\epsilon$ remaining. Finally, the result is the remainder term $R_m(\epsilon)$ plus a polynomial of $\epsilon$ with powers going from $\epsilon^{m+1}$ to $\epsilon^{2m+1}$. I find it quite fascinating, but maybe the Euler-MacLaurin approach is not the good one here. [End of paragraph added in the edit]

A last remark: this result is not universal. If I pick the not-so-different function $\tilde f_n(\epsilon)=2n e^{-n^2\epsilon}$, then the corresponding $\tilde Z_\text{odd}$ and $\tilde Z_\text{even}$ differ by a quantity of order $\epsilon$.

So, what is so special about my choice of $f_n(\epsilon)$ ? Why are my two $Z$ functions so close? How can I compute the asymptotic equivalent of $Z_\text{even}(\epsilon)-Z_\text{odd}(\epsilon)$ as $\epsilon\searrow0$ ?

Answer 1

We have (for $\epsilon>0$) $$Z_\text{even}(\epsilon)-Z_\text{odd}(\epsilon)=\sum_{n\geqslant 0}(-1)^n f_n(\epsilon)=\frac12\sum_{n\in\mathbb{Z}}(-1)^n f_n(\epsilon)$$ since clearly $f_{-n-1}(\epsilon)=-f_n(\epsilon)$. Now apply the Poisson summation formula $$\sum_{n\in\mathbb{Z}} g(n)=\sum_{n\in\mathbb{Z}}\hat{g}(n),\qquad\hat{g}(y)=\int_{-\infty}^\infty g(x)e^{-2i\pi xy}\,dx$$ to $g(x)=(2x+1)e^{-x(x+1)\epsilon-i\pi x}$. To compute $\hat{g}(y)$, use the integral $$\int_{-\infty}^\infty te^{-at^2+2bt}\,dt=\frac{b\sqrt\pi}{a^{3/2}}e^{b^2/a}\qquad(a\in\mathbb{R}_{>0},b\in\mathbb{C})$$ (obtained via $t=z+b/a$ and rectangular contours). This way we get $$Z_\text{even}(\epsilon)-Z_\text{odd}(\epsilon)=\left(\frac\pi\epsilon\right)^{3/2}e^{\epsilon/4}\sum_{n\geqslant 0}(-1)^n(2n+1)e^{-(2n+1)^2\pi^2/(4\epsilon)}.$$ This is very small $\asymp(\pi/\epsilon)^{3/2}e^{-\pi^2/(4\epsilon)}$ when $\epsilon\to 0^+$.