Convergence of a family of series
Edited on 2019 11/18.
This is an incorrect approach. Given $a,b,c$, there are infinitely many $\gamma\in SL_2(\mathbb{Z})$. For a correct approach, see fedja's answer.
Old answer
This is an incomplete answer using the Theory of Quadratic Forms. The constants $C$ appearing in this solution may not be the same. We begin with
Definition
Let $ax^2+bxy+cy^2$ be a quadratic form with $a, b, c\in \mathbb{Z}$ and $b^2-4ac=D$. We say that a quadratic form $a_1x_1^2+b_1x_1y_1+c_1y_1^2$ with $a_1,b_1,c_1\in\mathbb{Z}$ and $b_1^2-4a_1c_1=D$ is equivalent to $ax^2+bxy+cy^2$ if there is $\begin{bmatrix}p & q \\ r & s \end{bmatrix}\in\mathrm{SL}_2(\mathbb{Z})$ such that $$ a_1x_1^2+b_1x_1y_1+c_1y_1^2=a(px_1+qy_1)^2+b(px_1+qy_1)(rx_1+sy_1)+c(rx_1+sy_1)^2. $$
This can be written again with matrices $$ \begin{bmatrix}a_1 & \frac{b_1}2\\ \frac{b_1}2 & c_1 \end{bmatrix} = \begin{bmatrix}p & r \\ q & s \end{bmatrix} \begin{bmatrix}a & \frac{b}2\\ \frac{b}2 & c \end{bmatrix} \begin{bmatrix}p & q \\ r & s \end{bmatrix} $$ and $$ \begin{bmatrix}x\\y\end{bmatrix}=\begin{bmatrix}p & q \\ r & s \end{bmatrix}\begin{bmatrix}x_1\\y_1\end{bmatrix}. $$
The next one is the reducing step.
Lemma 1
Any integer quadratic form with discriminant $D$ is equivalent to $ax^2+bxy+cy^2$ with $|b|\leq |a|\leq |c|$.
The proof is a repeated application of the matrices in $\mathrm{SL}_2(\mathbb{Z})$. $$ \begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix}, \ \ \begin{bmatrix} 0 & -1 \\ 1 & 0\end{bmatrix}. $$
Corollary 1 The number of equivalence classes of integer quadratic forms with discriminant $D$ is finite.
If $b^2-4ac=D$ and $|b|\leq |a|\leq |c|$, then $|D|\geq 4|ac|-|b^2| \geq 3|ac|$ gives only finitely many choices for $a$ and $c$, subsequently for $b$.
Then it suffices to prove the following for $A=\begin{bmatrix}a& b/2\\ b/2&c\end{bmatrix}$ with $b^2-4ac=D$ and $|b|\leq |a|\leq |c|$, as the sum can be written as a finite sum of the following kind.
The sum is absolutely convergent for $k>1$: $$\sum_{\gamma\in \mathrm{SL}_2(\mathbb{Z})} \frac1{\left(| \gamma z - \alpha | |\gamma z - \alpha'|\right)^k|rz+s|^{2k}} \ \ \ (*) $$ where $\gamma z=\frac{pz+q}{rz+s}$, $\alpha=\frac{-b-\sqrt{D}}{2a}$, and $\alpha'=\frac{-b+\sqrt{D}}{2a}$.
Note that after the reduction, we have $ac<0$ and $|\alpha-\alpha'|>1$.
Let $\Delta_{\gamma}$ be the area of a triangle with vertices at $\alpha, \alpha',\gamma z$. Let $\theta_{\gamma}$ be the angle formed by the corner $\alpha, \gamma z, \alpha'$. Then $\Delta_{\gamma}=\frac12| \gamma z - \alpha | |\gamma z - \alpha'|\sin \theta_{\gamma}=\frac12 |\alpha-\alpha'|\Im(\gamma z)=\frac12|\alpha-\alpha'| \frac{\Im(z)}{|rz+s|^2}$.
Thus, it suffices to prove the convergence of $$ \sum_{\gamma\in \mathrm{SL}_2(\mathbb{Z})} (\sin \theta_{\gamma})^k $$
Consider $\Gamma=\mathrm{SL}_2(\mathbb{Z})$ and $H=\left\{\pm\begin{bmatrix} 1&2m_0m\\0&1\end{bmatrix}, m\in\mathbb{Z}\right\}$. Fix a fundamental domain $F_{\infty}$ for $H$ as $F_{\infty}=\{z| u-m_0 \le \Re z < u+m_0\}$ and $u=\frac{\alpha+\alpha'}2$. Here we set the integer $m_0$ large enough that $u-m_0<\min(\alpha,\alpha')<\max(\alpha,\alpha')<u+m_0$. For $\gamma z\notin F_{\infty}$, $\gamma z$ is sufficiently away from any of $\alpha$, $\alpha'$. Thus, we have $$ \sum_{(r,s)=1, \ \begin{bmatrix} *&*\\r&s\end{bmatrix} z \notin F_{\infty}}\frac1{\left(| \gamma z - \alpha | |\gamma z - \alpha'|\right)^k|rz+s|^{2k}}\ll \sum_{(r,s)=1} \frac1{|rz+s|^{2k}} $$ due to $\sum \frac 1{n^k}$ is convergent. Moreover, the sum on the right side is $$ \ll \sum_{(r,s)=1, \ rs\neq 0} \frac1{|rs|^k} $$ which is also convergent. The implied constants depend on $k,\alpha,\alpha',m_0,z$.
Thus, we consider the following sum $$ \sum_{\gamma\in H\backslash \Gamma} (\sin\theta_{\gamma})^k=\sum_{\gamma\in\Gamma, \ \gamma z \in F_{\infty}} (\sin \theta_{\gamma})^k. $$ We need a bound for the number of orbit points $\gamma z$ in the following figure.
It is actually rather straightforward. If $a=0$, then $b$ can take at most 2 values, none of them $0$, and we are left with the absolutely convergent sum $\sum_c |c+bi|^{-k}$. If $c=0$, we can run a similar argument. So it is enough to consider $a,b,c\ne 0$. In this case the roots $u,v$ of the corresponding quadratic polynomial $P(x)=ax^2+bx+c$ are real and lie at the distance $\sqrt D$ from $-\frac b{2a}$. Thus $|P(z)|=|a||z-u||z-v|$. This is never less than $|a||\Im z|^2$ and for $|b|>(6\sqrt D+4|\Re z|)|a|$ is not less than $|a||\frac b{4a}|^2\approx \frac{|b|^2}{|a|}$. Thus, $|P(z)|\gtrsim\max(|a|,\frac{|b|^2}{|a|})\ge |b|$ and we only need to consider $\sum_{a,b,c}|b|^{-k}$. However, for fixed $b$, the number of non-zero pairs $(a,c)$ with $0\ne b^2-D=4ac$ is at most the number of divisors of $b^2-D$, which grows slower than any positive power of $|b|$.