$L$-function, easiest way to see the following sum?

What is the easiest way to see that$$\sum_{(m, n) \in \mathbb{Z}^2 \setminus \{0, 0\}} (m^2 + n^2)^{-s} = 4\zeta(s)L(s, \chi)?$$Here $\chi$ is the homomorphism $(\mathbb{Z}/4\mathbb{Z})^\times \to \mathbb{C}^\times$ which sends $3$ mod $4$ to $-1$.

Since $\mathbb{Z}[i]$ is an Euclidean domain, it is a UFD domain. So we have that the numbers that are the sum of two squares are a semigroup, due to the Lagrange identity: $$ (a^2+b^2)(c^2+d^2) = (ac+bd)^2+(ad-bc)^2 \tag{1}$$ that is equivalent to the norm on $\mathbb{Z}[i]$ being multiplicative.

That also gives that the representation function: $$ r(n) = \#\{(a,b)\in\mathbb{Z}^2 : a^2+b^2 = n \} \tag{2}$$ is a constant times a multiplicative function. We have, indeed: $$ r(n) = 4\left(\chi_4 * 1\right)(n) = 4\sum_{d\mid n}\chi_4(n) \tag{3}$$ where $*$ is the Dirichlet convolution and $\chi_4$ is the non-principal Dirichlet character $\!\!\pmod{4}$.

Now back to our series. Assuming $\text{Re}(s)>1$, absolute convergence allows us to rearrange the series as: $$ S=\sum_{(m,n)\in\mathbb{Z}^2\setminus(0,0)}\frac{1}{(m^2+n^2)^s} = \sum_{n\geq 1}\frac{r(n)}{n^s} = 4\sum_{n\geq 1}\frac{(\chi_4*1)(n)}{n^s}\tag{4}$$ so by Dirichlet convolution we have: $$ S = 4 \sum_{n\geq 1}\frac{1}{n^s}\sum_{n\geq 1}\frac{\chi_4(n)}{n^s} = 4\zeta(s) L(\chi_4,s)\tag{5}$$ as wanted.

Footnote: $(3)$ may be proved also by manipulating Lambert series and exploiting the Jacobi triple product, or through modular forms. Anyway, I believe that the algebraic approach is way easier to follow. If we replace $m^2+n^2$ by $m^2+Dn^2$ and $m^2+Dn^2$ is the only reduced binary quadratic form of discriminant $-4D$ (aka $h(-4D)=1$, class number one), then $r(n)$ is still a constant times a multiplicative function and $S$ is still a constant times the product of two $L$-functions.