A von Staudt-type formula for $\displaystyle{\sum_{z\in\mathbb{Z}_n[i]} z^k }$

In this paper we study the sum of powers in the Gaussian integers $\mathbf{G}_k(n):=\sum_{a,b \in [1,n]} (a+b i)^k$. We give an explicit formula for $\mathbf{G}_k(n) \pmod n$ in terms of the prime numbers $p \equiv 3 \pmod 4$ with $p \mid \mid n$ and $p-1 \mid k$, similar to the well known one due to von Staudt for $\sum_{i=1}^n i^k \pmod n$. We apply this formula to study the set of integers $n$ which divide $\mathbf{G}_n(n)$ and compute its asymptotic density with six exact digits: $0.971000\ldots$.