Suppose $\tau\in\Bbb C$ and $\Im(\tau)>0$. Also, let $k\in\Bbb Z_{>2}$, and $A=\Bbb Z^2\setminus\{(0,0)\}$. Then the Eisenstein series $G_{2k}(\tau)$ is given by $$G_{2k}(\tau)=\sum_{(m,n)\in A}\frac{1}{(m+n\tau)^{2k}}.$$ It is then evident that $G_{2k}(\tau+1)=G_{2k}(\tau)$ for all $\tau$. Thus we may write $G_{2k}$ as a Fourier series $$G_{2k}(\tau)=\sum_{n\ge0}g_nq^n,$$ where $q=e^{2i\pi \tau}$. Supposedly, we can find these coefficients $g_n$ by computing the integral $$g_n=\int_0^1 e^{-2i\pi n\tau}G_{2k}(\tau)d\tau=\sum_{(u,v)\in A}\int_0^1 \frac{e^{-2\pi n\tau}}{(u+v\tau)^{2k}}d\tau.$$ Apparently, these coefficients have an explicit formula: $$\begin{align} g_0&=2\zeta(2k)\\ g_n&= (-1)^k\frac{2^{2k+1}\pi^{2k}}{(2k-1)!}\sigma_{2k-1}(n),\tag{1} \end{align}$$ but I have no clue how to prove this.
I thought about writing $(u+v\tau)^{-2k}$ as a power series $(u+v\tau)^{-2k}=\sum_{l\ge0}\alpha_l\tau^l$, but that seems overly complicated. This may be simplified by writing $$G_{2k}(\tau)=\sum_{n\ge1}\frac{1}{(a_n+b_n\tau)^{2k}},$$ where $(a_n,b_n)=($$\text{A174344}$$(n),$$\text{ A274923}$$(n))$.
Take note that I'm assuming the partial sums $\sum_{n=1}^{N}(a_n+b_n\tau)^{-2k}$ converge uniformly to $G_{2k}(\tau)$ as $N\to\infty$, hence I interchanged the sum and integral. Please correct me if this is not the case.
So, what's left is the integral $$j_n(a,b)=\int_0^1\frac{e^{-2i\pi n\tau}}{(a+b\tau)^{2k}}d\tau.$$ I'm fairly certain, since $\tau$ is a complex variable, that this integral will be taken over some path in the complex plane $\gamma$ that starts at $0$ and ends at $1$. This (potentially, but I'm not sure) being the case, I do not know which path to pick.
Could I have some help proving $(1)$? Is there a better approach than the evaluation of $j_n$?