Fourier Coefficients of Eisenstein Series $G_{2k}(\tau)=\sum\limits_{(m,n)\in\Bbb Z^2\setminus \{(0,0)\}}\frac{1}{(m+n\tau)^{2k}}$. 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$?
 A: Although this formula has been known for probably 200+ years, the computation is not at all obvious, I think, although with enough time a persistent person might stumble upon it.
It's not a short computation, either, and I hesitate to write it out here. Also, you will find it in any introductory notes on (holomorphic, elliptic) modular forms, such as in my http://www-users.math.umn.edu/~garrett/m/mfms/notes_2015-16/10_level_one.pdf  There are many other on-line and not-on-line sources for this, too. I'd wager that googling "Fourier coefficients of holomorphic Eisenstein series" will give more hits than you need.
Again, in my opinion, this computation is the kind of thing that one sees other people do, thereby learning from them, rather than just figuring it out for oneself. (I first saw it aeons ago in R. Gunning's little orange Princeton book "intro to modular forms".) It does use some mildly clever ideas, such as "winding-up" and "unwinding" integrals, and some "cancellation lemmas", ... which might not be obvious if one hasn't had that experience... as well as more-standard evalution of integrals by residues. :)
A: Apostol in his Modular functions and Dirichlet series in Number Theory gives a very simple proof based on the partial fraction expansion of cotangent function.
We have $$\pi\cot\pi \tau=\frac{1}{\tau}+\sum_{m\in\mathbb {Z}, m\neq 0}\left(\frac{1}{\tau+m}-\frac{1}{m}\right)\tag{1}$$ If $\tau$ has a positive imaginary part and $q=\exp(2\pi i\tau) $ then $|q|<1$ and we have $$\pi\cot\pi\tau=-\pi i\left(1+2\sum_{r=1}^{\infty}q^r\right)$$ so that we have $$\frac{1}{\tau}+\sum_{m\in\mathbb {Z}, m\neq 0}\left(\frac{1}{\tau+m}-\frac{1}{m}\right)=-\pi i\left(1+2\sum_{r=1}^{\infty} q^r\right)$$ Differentiating the above with respect to $\tau$ we get $$-\sum_{m\in\mathbb {Z}} \frac{1}{(\tau+m)^2}=-(2\pi i) \sum_{r=1}^{\infty}rq^{r-1}\frac{dq}{d\tau}=-(2\pi i) ^2\sum_{r=1}^{\infty} r\exp(2\pi i r\tau) \tag{2}$$ Differentiating this as many times as needed we can get expression for the sum $\sum(m+\tau) ^{-2k}$. The job is complete if we replace $\tau$ by $n\tau$ and sum over $n$.
Thus for example $$\sum_{m\in\mathbb {Z}} \frac{1}{(m+n\tau)^{4}}=\frac{(2\pi i) ^4}{3!}\sum_{r=1}^{\infty} r^3\exp(2\pi irn\tau) $$ and more generally we have $$\sum_{m\in\mathbb {Z}} \frac{1}{(m+n\tau)^{2k}}=\frac{(2\pi i) ^{2k}}{(2k-1)!}\sum_{r=1}^{\infty} r^{2k-1}\exp(2\pi irn\tau) =\frac{(2\pi i) ^{2k}}{(2k-1)!}\sum_{r=1}^{\infty} r^{2k-1}q^{nr}\tag{3}$$ Now we have
\begin{align}
G_{2k}(\tau)&=\sum_{m,n=-\infty,(m,n)\neq (0,0)}^{\infty} \frac{1}{(m+n\tau)^{2k}}\notag\\
&=\sum_{m\in\mathbb {Z}, m\neq 0}\frac{1}{m^{2k}}+\sum_{n=1}^{\infty} \sum_{m\in\mathbb {Z}} \left(\frac{1}{(m+n\tau)^{2k}}+\frac{1}{(m-n\tau)^{2k}}\right)\notag\\
&=2\zeta(2k)+2\sum_{n=1}^{\infty} \sum_{m\in\mathbb {Z}} \frac{1}{(m+n\tau)^{2k}}\notag\\
&=2\zeta(2k)+2\cdot\frac{(2\pi i) ^{2k}}{(2k-1)!}\sum_{n=1}^{\infty} \sum_{r=1}^{\infty} r^{2k-1}q^{nr}\notag\\
&=2\left(\zeta(2k)+(-1)^k\frac{(2\pi)^{2k}}{(2k-1)!}\sum_{r=1}^{\infty} \frac{r^{2k-1}q^r}{1-q^r}\right)\notag\\
&=2\left(\zeta(2k)+(-1)^k\frac{(2\pi)^{2k}}{(2k-1)!}\sum_{n=1}^{\infty} \sigma_{2k-1}(n)q^n\right)\notag
\end{align}

Elliptic function theory is difficult but not too difficult if one gets hold of good books. Essentially it requires a basic knowledge of calculus and a very deep skill in algebraic manipulation. At least thats how people like Legendre, Abel, Gauss, Jacobi and Ramanujan developed the subject. Then somehow Liouville and Weierstrass arrived on the scene and ditched algebraic manipulation altogether and exclusively used complex analysis making the topic widely inaccessible. 
