I really prefer the usual conventions for integer variables, so I will consider:
$$S_m=\sum_{n=1}^m n^d \exp(−n^2/\sigma^2)$$
We have a general inequality for a monotone decreasing function:
$$\int _{N}^{M+1}f(x)\,dx\leq \sum _{n=N}^{M}f(n)\leq f(N)+\int _{N}^{M}f(x)\,dx$$
Note that our function is not monotone for $d>0$, it has a maximum:
$$f(x)=x^d e^{-x^2/\sigma^2}$$
$$f'(x)=\left(d -\frac{2x^2}{\sigma^2}\right)x^{d-1} e^{-x^2/\sigma^2}$$
So we have: $$x_0= \sigma\sqrt{d/2}$$
$$n_0=\lceil \sigma\sqrt{d/2} \rceil$$
Thus, we need to write for $n_0>1$:
$$S_m=\sum_{n=1}^{n_0-1} n^d \exp(−n^2/\sigma^2)+\sum_{n=n_0}^m n^d \exp(−n^2/\sigma^2)$$
The latter sum has the function decreasing monotonely, and we can use the bound:
$$S_m \leq \sum_{n=1}^{n_0} n^d \exp(−n^2/\sigma^2)+\int_{n_0}^{m} x^d \exp(−x^2/\sigma^2) dx$$
Thus:
$$c \leq \sum_{n=1}^{n_0} n^d \exp(−n^2/\sigma^2)+\int_{n_0}^{m} x^d \exp(−x^2/\sigma^2) dx$$
$$n_0=\lceil \sigma\sqrt{d/2} \rceil$$
Mind, this is definitely not a simplification of the bound, as the integral can only be expressed in terms of incomplete Gamma function for the general $d$.
But if $m$ is very large, we can approximate the integral by $\int_{n_0}^\infty x^d \exp(−x^2/\sigma^2) dx$, and write:
$$c \leq \sum_{n=1}^{n_0} n^d \exp(−n^2/\sigma^2)+\int_{n_0}^\infty x^d \exp(−x^2/\sigma^2) dx$$
If, in addition, $n_0<1$ (it depends on the values of $d$ and $\sigma$), we have:
$$c \leq \int_0^\infty x^d \exp(−x^2/\sigma^2) dx=\frac{\sigma^{d+1}}{2} \Gamma \left( \frac{d+1}{2} \right)$$
This is actually a good approximation for quite a lot of combinations of $d,\sigma,m$, just be careful to make sure the conditions for the approximation are satisfied.
For example, for $d=1$, $\sigma=5$ and $m=25$ we have:
$$S_m=12.41633011 \\ \frac{\sigma^{d+1}}{2} \Gamma \left( \frac{d+1}{2} \right)=12.5$$
(It somehow works even better for larger $d$, but I'm too tired to figure out why. I leave this to the OP).
