Asymptotics of a Mathieu-Gaussian series

We consider the asymptotic expansion of the functional series $$S_{\mu,\gamma}(a;\lambda)=\sum_{n=1}^\infty \frac{n^\gamma e^{-\lambda n^2/a^2}}{(n^2+a^2)^\mu}$$ for real values of the parameters $\gamma$, $\lambda>0$ and $\mu\geq0$ as $|a|\to \infty$ in the sector $|\arg\,a|<\pi/4$. For general values of $\gamma$ the expansion is of algebraic type with terms involving the Riemann zeta function and a terminating confluent hypergeometric function. Of principal interest in this study is the case corresponding to even integer values of $\gamma$, where the algebraic-type expansion consists of a finite number of terms together with a contribution comprising an infinite sequence of increasingly subdominant exponentially small expansions. This situation is analogous to the well-known Poisson-Jacobi formula corresponding to the case $\mu=\gamma=0$. Numerical examples are provided to illustrate the accuracy of these expansions.