I have a sum of the form
$$\sum_{n,m=-N}^N e^{-\alpha (n-m)^2}$$
where $α>0$ is some constant, and I don't mind if the limit $N\rightarrow\infty$ is taken. I know there is a possibility of exchanging the "variable of summation" by introducing some "relative coordinates" $\tau=n−m$ and $\eta=m+n$, but I don't know how to continue from there.
That particular sum will not be very well behaved, because you are summing many terms very close to one. Most of the mechanics with this type of sum are fairly easy; for example, you should definitely be able to do the conversion to relative coordinates (or the whole problem will be much too complicated). There is one tricky step, though, at the heart of the problem, and it is the sum $$ \sum_{k=-\infty}^\infty e^{-\alpha k^2}=\vartheta_3(0,e^{-\alpha}), $$ where $\vartheta_3(0,e^{-\alpha})$ is an auxiliary Jacobi theta function. Those are hard to deal with; the DLMF has a good toolbox on them but they are essentially impossible to reduce to anything 'nicer'.
One perspective on this function that I really like is what happens if you transform to the variable $\sigma=1/\sqrt{\alpha}$, in which case the sum above looks like two smoothly joined straight lines:
The slope on the asymptote is, weirdly enough, $\sqrt{\pi}$.
As in all things, though, it depends what you want do do with that expression. If you make the question more concrete then there may be better ways to help you.
I) OP mentions in a comment that his sum should be normalized with a factor $\frac{1}{2N+1}$. Hence OP is considering
$$S_N ~:=~\frac{1}{2N+1}\sum_{-N\leq k,\ell\leq N} e^{-\alpha (k+\ell)^2} ~\stackrel{(3)}{=}~ \frac{1}{2N+1}\sum_{n,m\in\mathbb{Z}}^{-2N\leq n\pm m\leq 2N} e^{-\alpha n^2}$$ $$\tag{1}~=~\sum_{-2N\leq n\leq 2N}\frac{2N-|n|+1}{2N+1} e^{-\alpha n^2}~=~\sum_{-2N\leq n\leq 2N}\left(1-\frac{|n|}{2N+1}\right) e^{-\alpha n^2}, $$
where we have defined
$$\tag{2} n~:=~k+\ell~\in~\mathbb{Z} \quad\text{and}\quad m~:= ~k-\ell~\in~\mathbb{Z},$$
or equivalently,
$$\tag{3} k~=~\frac{n+m}{2} \quad\text{and}\quad \ell~=~\frac{n-m}{2}. $$
II) Now let us analyze the sum (1) further.
For $\alpha \gg 1$, we only need sum the few first terms ($n$ small) in the sum (1).
For $N^{-2}\ll\alpha \ll 1$, we get
$$\tag{4} S_N ~\approx~\int_{-2N}^{2N} \! dx~\left(1-\frac{|x|}{2N+1}\right) e^{-\alpha x^2}~\approx~\sqrt{\frac{\pi}{\alpha}}-\frac{\alpha^{-1}}{2N+1}.$$
For $\alpha=0$, we get
$$\tag{5} S_{N}=\sum_{-2N\leq n\leq 2N}\left(1-\frac{|n|}{2N+1}\right) ~=~4N+1-\frac{2}{2N+1}\sum_{1\leq n\leq 2N}n~=~ 2N+1.$$
I might not know a crucial detail because in the other comments you guys seem to dicuss otherwise, but as far as I can tell you have only positive summands and $2N-1$ of them are of the form $\mathrm{e}^{-\alpha(k-k)}=1$. Hence your sum diverges as $2N-1$: $$\lim_{N\to\infty}\sum_{n,m=-N}^N e^{-\alpha (n-m)^2}=\mathrm{e}^{-\alpha(0-0)}+\mathrm{e}^{-\alpha(1-1)}+\mathrm{e}^{-\alpha(2-2)}+\dots=\infty$$
Emilio notes a "weird" $\sqrt{\pi}$ slope, but it's just the proportionality factor $1+f$ you get from the diagonal heuristic $f\approx\int \mathrm{e}^{-\alpha\, k^2}\mathrm{d}k=\sqrt{\frac{\pi}{\alpha}}$.
Just from counting squares, your sum should be $$\sum_{n,m=-N}^N \mathrm{e}^{-\alpha (n-m)^2}\approx (2N-1)\left(1+2\sum_{k=1}^{\infty}\mathrm{e}^{-\alpha\, k^2}\right).$$ The second bracket is "EllipticTheta[3, 0, E^-[Alpha]]" in Mathematica notation, and for $\alpha\approx 5$ you can consider it $1$.
