 Solution

On any compact set avoiding the integers, writeFor all sufficiently large $|n|$, uniformly on that compact set this is $O(|n{|}^{-2})$. The summands in the higher series are $O(|n{|}^{-k})$ for $k\ge 2$. The Weierstrass M-test therefore gives locally uniform convergence away from the integer poles and permits termwise differentiation. Near an integer $m$, isolate its singular summand; all the others converge uniformly on a small neighborhood. Thus the first function has a simple pole of residue one at each integer, and the $k$th function has a pole of order $k$ there with principal part $(z-m{)}^{-k}$. There are no other singularities, so all are meromorphic functions.

Pairing positive and negative indices is legitimate because the regularized first series converges absolutely. It becomesHere is a direct derivation of its cotangent partial-fraction expansion. Fix $z\notin \mathbb{Z}$ and integrate $F(\zeta )=\pi \cot (\pi \zeta )/({\zeta }^2-z^2)$ around the square with real and imaginary coordinates $\pm R$, $R=N+1/2$, taking $N$ large enough to enclose $\pm z$. On the vertical sides $|\cot \pi \zeta |\le 1$, and on the horizontal sides it is at most $\coth (\pi R)$. Thus the cotangent is uniformly bounded and $|{\zeta }^2-z^2|\ge R^2-|z{|}^2$, so the contour integral is $O(R^{-1})\to 0$.

The residues at integers are $1/(n^2-z^2)$; the two residues at $z$ and $-z$ sum to $\pi \cot (\pi z)/z$. The residue theorem and absolute convergence of the integer residue sum consequently giveUsing the exponential definitions of sine and cosine now yields, for $w=e^{2\pi iz}$,The same locally uniform convergence gives ${\epsilon }_k^{\prime }=-k{\epsilon }_{k+1}$, including $k=1$. Hence if ${\epsilon }_k(z)={\phi }_k(e^{2\pi iz})$, the chain rule givesThis is a rational function whenever ${\phi }_k$ is, proving the requested descent for every $k$ by induction, withIn particular periodicity is a consequence of the explicit descent, rather than an assumption about the summation order.

Finally, for $|z|<1$ expand each paired summand as a geometric series:On $|z|\le r<1$ the double series converges absolutely, since its absolute sum is at most $2r{\sum }_{n\ge 1}{n}^{-2}/(1-r^2)$. Interchanging the sums and using the Riemann zeta function gives the Laurent seriesThus the $z^{-1}$ coefficient is one; the positive odd coefficients are $-2\zeta (2j+2)$, and all other coefficients vanish.