Does this series over $SL_2(\mathbb{Z})$ converge? Consider the following series of positive numbers
$$
\sum_{\substack{a,b,c,d \in \mathbb{Z} \\ ad-bc=1}}\frac{1}{((a-b)^2 + (c-d)^2 + 1)^2}
$$
Does it converge? How to prove it?
 A: The answer to your literal question is "no", but for a reason which may clarify what question you actually intend to ask:
Since $\pmatrix{a & b \cr c & d}\pmatrix{1 \cr -1}=\pmatrix{a-b\cr c-d}$, the expression in the denominator is the fourth power of the length of the latter vector, plus $1$.
A key point is that isotropy/stabilizer subgroups in $SL_2(\mathbb Z)$ of non-zero vectors in $\mathbb Z^2$ are infinite. Here, $\pmatrix{1 & 0 \cr 1 & 1}\pmatrix{1 \cr -1}=\pmatrix{1\cr 0}$, and the isotropy subgroup of the latter vector is recognizable as $\pmatrix{1 & n\cr 0 & 1}$ for $n\in\mathbb Z$.
So at the very least you may want your sum to be not over $SL_2(\mathbb Z)$, but over the coset space $SL_2(\mathbb Z)/H$, where $H$ is the isotropy subgroup of $\pmatrix{1 \cr -1}$ in $SL_2(\mathbb Z)$.
In that case, your sum is a subsum of the sum of $4$th powers of lengths of non-zero vectors in $\mathbb Z^2$. (In particular, it is the sum over vectors $v=(x,y)$ with $x,y$ coprime, but we don't need this.) That is, $\sum_{x,y} {1\over (x^2+y^2)^2}$ with not both $x,y$ being $0$. Fooling around with polar coordinates and/or integral tests shows that this does converge.
EDIT: the revised version of the question is significantly different, even if the formulaic assertion may seem similar. Reformulating it to help explain its meaning, let $\varphi(z,w)=1/(z-\overline{w})^{2k}$. With $2k=2$ this would essentially be a Bergman kernel for the upper half-plane. For $2k>2$, this will wind up to be a reproducing kernel for holomorphic cuspforms (something like this appeared in Petersson 1940). That is,
$$
K(z,w) \;=\; \sum_{\gamma\in SL_2(\mathbb Z)} {1\over (cz+d)^{2k}}\varphi(\gamma z,w)
$$
has the property that (up to a constant)
$$ \int_{\Gamma\backslash\mathfrak H} K(z,w)\,f(w)\;(\Im w)^{2k} {|dw|\over (\Im w)^2} \;\;=\;\; f(z)
$$
where $\Gamma=SL_2(\mathbb Z)$. Summing over integer matrices with determinants $N$, instead, gives the $N$th Hecke operator. (This was also known long ago. S. Lang's "Intro to modular forms" has an appendix by D. Zagier using this to compute traces. In the first edition of the book, there was an error due to the fragility of the convergence of the trace... corrected in a later edition.)
A non-obvious, but conceptually robust and not-so-computational approach to proof of convergence is to "observe" that the given series is a subsum of the pullback of the Siegel-type holomorphic Eisenstein series on the Siegel upper half-space $\mathfrak H_2$ to the diagonally imbedded $\mathfrak H\times \mathfrak H$. (This was discussed in my paper, and H. Klingen's paper, in the 1984 Katata conference volume "Automorphic forms in several variables".) The convergence of those Eisenstein series was proven by H. Braun by 1939. Analogous reproducing kernels for holomorphic cuspforms exist for all the hermitian classical domains.
Of course, none of this is strictly necessary to understand convergence, but it does explain why we might care, and so on.
