Infinite Series $‎\sum\limits_{n=2}^{\infty}\frac{\zeta(n)}{k^n}$ ‎If $f\left(z \right)=\sum_{n=2}^{\infty}a_{n}z^n$ and $\sum_{n=2}^{\infty}|a_n|$ converges then‎,
$$\sum_{n=1}^{\infty}f\left(\frac{1}{n}\right)=\sum_{n=2}^{\infty}a_n\zeta\left(n\right)‎.$$
‎Since if we set $C:=\sum_{m=2}^{\infty}|a_m|<\infty$‎, ‎then‎
$$\sum_{n=1}^{\infty}\sum_{m=2}^{\infty}|a_m\frac{1}{n^m}|\leq\sum_{n=1}^{\infty}\sum_{m=2}^{\infty}|a_m|\frac{1}{n^2}\leq C\sum_{n=1}^{\infty}\frac{1}{n^2}<\infty‎$$
and ‎by Cauchy's double series theorem‎, ‎we can switch the order of summation‎:
$$\sum_{n=1}^{\infty}f\left(\frac{1}{n}\right)=\sum_{n=1}^{\infty}\sum_{m=2}^{\infty}a_m\frac{1}{n^m}=\sum_{m=2}^{\infty}a_m\sum_{n=1}^{\infty}\frac{1}{n^m}=\sum_{n=2}^{\infty}a_n\zeta\left(n\right)‎.$$
This shows that $‎\sum_{n=2}^{\infty}\frac{\zeta(n)}{k^n}=\sum_{n=1}^{\infty}\frac{1}{kn(kn-1)}$.
My Questions:
1) It's obvious that $\sum_{n=1}^{\infty}\frac{1}{2n(2n-1)}=\log(2)$, but how can I evaluate $\sum_{n=1}^{\infty}\frac{1}{3n(3n-1)}$?
2) Is there another method to evaluate $\sum_{n=2}^{\infty}\frac{\zeta(n)}{k^n}$? 
 A: Answering you first question.
Don't know whether it's simpler but still (since you'll have to deal with hypergeometric functions).
$$ \sum_{n=1}^{\infty}\frac{1}{3n(3n-1)}=\frac{1}{3}\sum_{n=1}^{\infty}\frac{1}{3n-1}\int_0^1x^{n-1} \mathrm dx=\frac{1}{3}\int_0^1\sum_{n=1}^{\infty}\frac{x^{n-1}}{3n-1}\mathrm dx$$
$$\sum_{n=1}^{\infty}\frac{x^{n-1}}{3n-1}=\frac{1}{2} \, _2F_1\left(\frac{2}{3},1;\frac{5}{3};x\right)$$
And $$\frac{1}{6} \int_0^1 \, _2F_1\left(\frac{2}{3},1;\frac{5}{3};x\right) \, \mathrm dx=\frac{1}{6} \left(3\log(3) -\frac{\pi }{\sqrt{3}}\right)$$
So $$\sum_{n=1}^{\infty}\frac{1}{3n(3n-1)}=\frac{1}{6} \left(3\log(3) -\frac{\pi }{\sqrt{3}}\right)$$
A: This is an aswer for every $k$, using your method. For $k=2$ and $3$ it fits with other answers, so it's probably correct :) It is completely elementary, in the sense that it uses just the Taylor expansion of $\log(1-x)$ and the fact that the sum $\sum_{\alpha^k=1}\alpha^n$ ($\alpha$ runs over the $k$-th roots of $1$) is $k$ if $k$ divides $n$ and $0$ otherwise.
There are some $\log$'s of complex numbers. Those numbers have always non-negative real part, for the Arg we take the angle between $-\pi/2$ and $\pi/2$, so that it fits with the power series for $\log(1-x)$. 
$$\sum_n x^{kn}/kn=-\log(1-x^k)/k=-\frac{1}{k}\sum_{\alpha^k=1}\log(1-\alpha x)$$
$$\sum_{\alpha^k=1}\sum_{m=1}^\infty\alpha(\alpha x)^m/m=k\sum_{n=1}^{\infty}x^{kn-1}/(kn-1)$$
but also
$$\sum_{\alpha^k=1}\sum_{m=1}^\infty\alpha(\alpha x)^m/m=-\sum_{\alpha^k=1}\alpha\log(1-\alpha x).$$
We thus have
$$\sum_n x^{kn}/(kn(kn-1))=\sum_n x^{kn}(\frac{1}{kn-1}-\frac{1}{kn})=$$
$$=\frac{1}{k}\sum_{\alpha^k=1}(1-x\alpha)\log(1-\alpha x).$$
We have to take the limit $x\to 1$. The $\alpha=1$ term disappears, so we get
$$\sum_n\frac{1}{kn(kn-1)}=\frac{1}{k}\sum_{\alpha^k=1,\alpha\neq1}(1-\alpha)\log(1-\alpha)=$$
($\alpha=\exp(2\pi i m/k)$)
$$=\frac{1}{k}\sum_{m=1}^{k-1}((1-\cos\frac{2\pi m}{k})-i\sin\frac{2\pi m}{k})(\log(2\sin\frac{\pi m}{k})+\pi i(m/k-1/2))=$$
$$=\frac{1}{k}\sum_{m=1}^{k-1}(1-\cos\frac{2\pi m}{k})\log(2\sin\frac{\pi m}{k})+\pi(m/k-1/2)\sin\frac{2\pi m}{k}.$$
A: In the intro above the question is asked how one might go about computing $$\sum_{n=1}^\infty \frac{1}{3n(3n-1)}.$$
I can contribute an asymptotic expansion and an infinite series to this discussion, shown and proved below. In addition to posting this answer I am also asking two questions I. Can you verify this expansion using a different proof technique? II. Might there even be a reference to the below formula somewhere?
This is the proof, which is procedural, with the result at the end. Suppose we seek to evaluate
$$T(q) = \sum_{n\ge 1} \frac{1}{qn(qn-1)} =
\sum_{n=1}^p \frac{1}{qn(qn-1)} + \sum_{n\ge 1} \frac{1}{(qn+pq)(qn+pq-1)} \\=
\sum_{n=1}^p \frac{1}{qn(qn-1)} + \frac{1}{q^2}\sum_{n\ge 1} \frac{1}{(n+p)(n+(pq-1)/q)}$$
with $p, q\ge 2$ two integers. 
The sum term, call it $S(x)$, is harmonic  and may be evaluated by inverting its Mellin transform.
Recall the Mellin transform identity for harmonic sums with base function $g(x)$, which is
$$\mathfrak{M}\left(\sum_{k\ge 1}\lambda_k g(\mu_k x); s\right) =
\left(\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s}\right) g^*(s)$$
where $g^*(s)$ is the Mellin transform of $g(x).$
In the present case we have
$$\lambda_k = 1, \quad
\mu_k = k \quad \text{and} \quad
g(x) = \frac{1}{(x+p)(x+(pq-1)/q)}.$$
The Mellin transform of $g(x)$ is
$$\int_0^\infty  \frac{1}{(x+p)(x+(pq-1)/q)} x^{s-1} dx.$$
We chose $p\ge 2$ so that the base function does not have a pole at $x=-1$, which would cause the expansion about infinity that we eventually obtain not to converge. (As we evaluate the harmonic sum at $x=1$ this would be situated right on the boundary between the two disks about zero and about infinity, causing an evaluation at $x=1$ to fail.)
We evaluate the Mellin transform next. Use a circular contour to get
$$ g^*(s) (1 - e^{2\pi i (s-1)}) \\=
2\pi i \left(
\operatorname{Res}\left(\frac{x^{s-1}}{(x+p)(x+\frac{pq-1}{q})}; s=-p\right) +
\operatorname{Res}\left(\frac{x^{s-1}}{(x+p)(x+\frac{pq-1}{q})}; s=-\frac{pq-1}{q}\right)
\right) \\=
2\pi i \left(-q(-p)^{s-1} + q(-(pq-1)/q)^{s-1}\right)
= 2\pi i e^{\pi i (s-1)} q \left(((pq-1)/q)^{s-1}-p^{s-1}\right).$$
This implies that
$$g^*(s) = 2\pi i  q \frac{-e^{\pi i s}}{1 - e^{2\pi i s}}
\left(\left(\frac{pq-1}{q}\right)^{s-1} - p^{s-1}\right)\\=
-2\pi i q \frac{1}{e^{-\pi i s} - e^{\pi i s}}
\left(\left(\frac{pq-1}{q}\right)^{s-1} - p^{s-1}\right) \\=
q\pi\frac{2i}{e^{\pi i s} - e^{-\pi i s}}
\left(\left(\frac{pq-1}{q}\right)^{s-1} - p^{s-1}\right) =
q\frac{\pi}{\sin(\pi s)} \left(\left(\frac{pq-1}{q}\right)^{s-1} - p^{s-1}\right).$$
It follows that the Mellin transform of the sum term for $S(x)$ is given by
$$ Q(s) =
q \frac{\pi}{\sin(\pi s)}
\left(\left(\frac{pq-1}{q}\right)^{s-1} - p^{s-1}\right) \zeta(s)
\quad\text{because}\quad \sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} = \zeta(s).$$
Finally we invert the Mellin transform with the inversion integral
$$\frac{1}{2\pi i}\int_{3/2-i\infty}^{3/2+i\infty} Q(s)/x^s ds$$
to obtain an expansion about infinity of $S(x)$ starting with the pole at $s=2$ and getting
$$ S(x) \sim
- q \sum_{m\ge 2} (-1)^m
\left(\left(\frac{pq-1}{q}\right)^{m-1} - p^{m-1}\right) \frac{\zeta(m)}{x^m}.$$
Setting $x=1$ and substituting into the original equation we finally have
$$ T(q) = \sum_{n=1}^p \frac{1}{qn(qn-1)} 
- \frac{1}{q} \sum_{m\ge 2} (-1)^m
\left(\left(\frac{pq-1}{q}\right)^{m-1} - p^{m-1}\right) \zeta(m).$$
Addendum Wed Apr  2 20:14:50 CEST 2014. The following MSE link shows how to work with divergent Mellin transforms where there is a pole on the positive real axis (even ones).
A: $\newcommand{\+}{^{\dagger}}
 \newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle}
 \newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack}
 \newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,}
 \newcommand{\dd}{{\rm d}}
 \newcommand{\down}{\downarrow}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,{\rm e}^{#1}\,}
 \newcommand{\fermi}{\,{\rm f}}
 \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}
 \newcommand{\half}{{1 \over 2}}
 \newcommand{\ic}{{\rm i}}
 \newcommand{\iff}{\Longleftrightarrow}
 \newcommand{\imp}{\Longrightarrow}
 \newcommand{\isdiv}{\,\left.\right\vert\,}
 \newcommand{\ket}[1]{\left\vert #1\right\rangle}
 \newcommand{\ol}[1]{\overline{#1}}
 \newcommand{\pars}[1]{\left(\, #1 \,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\pp}{{\cal P}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,}
 \newcommand{\sech}{\,{\rm sech}}
 \newcommand{\sgn}{\,{\rm sgn}}
 \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}}
 \newcommand{\ul}[1]{\underline{#1}}
 \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}
 \newcommand{\wt}[1]{\widetilde{#1}}$
\begin{align}
&\color{#c00000}{\sum_{n = 2}^{\infty}{\zeta\pars{n} \over k^{n}}}
=\sum_{n = 2}^{\infty}{1 \over k^{n}}
\sum_{\ell = 0}^{\infty}{1 \over \pars{\ell + 1}^{n}}
=\sum_{\ell = 0}^{\infty}\sum_{n = 2}^{\infty}{1 \over \bracks{k\pars{\ell + 1}}^{n}}
=\sum_{\ell = 0}^{\infty}
{\bracks{k\pars{\ell + 1}}^{-2} \over 1 - \bracks{k\pars{\ell + 1}}^{-1}}
\\[3mm]&={1 \over k^{2}}\sum_{\ell = 0}^{\infty}
{1 \over \pars{\ell + 1}\pars{\ell + 1 - 1/k}}
={1 \over k^{2}}\,{\Psi\pars{1} - \Psi\pars{1 - 1/k} \over 1 - \pars{1 - 1/k}}
\\[3mm]&=\color{#c00000}{%
{1 \over k}\bracks{\Psi\pars{1} - \Psi\pars{1 - {1 \over k}}}}
\end{align}
where $\ds{\Psi\pars{z}}$ is the Digamma Function.

$$
\color{#00f}{\large\sum_{n = 2}^{\infty}{\zeta\pars{n} \over k^{n}}
=-\,{1 \over k}\,\bracks{\gamma + \Psi\pars{1 - {1 \over k}}}}\,,\qquad\verts{k} > 1
$$
  since $\ds{\Psi\pars{1} = -\gamma.\quad}$ $\gamma$ is the
  Euler-Mascheroni Constant. 

$\large\tt\mbox{Just one of your questions !!!}$:
\begin{align}
\sum_{n = 1}^{\infty}{1 \over 3n\pars{3n - 1}}&
={1 \over 9}\sum_{n = 1}^{\infty}{1 \over n\pars{n - 1/3}}
={1 \over 9}\sum_{n = 0}^{\infty}{1 \over \pars{n + 1}\pars{n + 2/3}}
={1 \over 9}\,{\Psi\pars{1} - \Psi\pars{2/3} \over 1 - 2/3}
\\[3mm]&={1 \over 3}\,\bracks{\Psi\pars{1} - \Psi\pars{{2 \over 3}}}
\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\pars{1}
\end{align}
Also,

$$
\Psi\pars{{2 \over 3}} = -\gamma + {\root{3} \over 6}\,\pi - {3 \over 2}\,\ln\pars{3}
$$ 
  By replacing in $\pars{1}$, we'll find:

$$
\color{#00f}{\large\sum_{n = 1}^{\infty}{1 \over 3n\pars{3n - 1}}
={1 \over 18}\bracks{9\ln\pars{3} - \root{3}\pi}} \approx 0.2470
$$
A: Extended Harmonic Numbers
Normally, we think of Harmonic Numbers as
$$
H_n=\sum_{k=1}^n\frac1k\tag{1}
$$
However, an alternate definition is often useful:
$$
H_n=\sum_{k=1}^\infty\left(\frac1k-\frac1{k+n}\right)\tag{2}
$$
For integer $n\ge1$, it is not too difficult to see that the two definitions agree. However, $(2)$ is easily extendible to all $n\in\mathbb{R}$ (actually, to all $n\in\mathbb{C}$). We can say some things about $H_n$ for some $n\in\mathbb{Q}\setminus\mathbb{Z}$.
Note that for $m,n\in\mathbb{Z}$,
$$
\begin{align}
H_{mn}-H_n
&=\sum_{k=1}^\infty\left(\frac1k-\frac1{k+mn}\right)-H_n\\
&=\sum_{k=1}^\infty\sum_{j=0}^{m-1}\left(\frac1{km-j}-\frac1{km-j+mn}\right)-H_n\\
&=\frac1m\sum_{j=0}^{m-1}\sum_{k=1}^\infty\left(\left(\frac1k-\frac1{k-j/m+n}\right)-\left(\frac1k-\frac1{k-j/m}\right)\right)-H_n\\
&=\frac1m\sum_{j=0}^{m-1}\left(\left(H_{n-j/m}-H_n\right)-H_{-j/m}\right)\tag{3}
\end{align}
$$
Since $H_0=0$ and $H_n=\log(n)+\gamma+O\left(\frac1n\right)$, where $\gamma$ is the Euler-Mascheroni Constant, if we let $n\to\infty$ in $(3)$, we get
$$
\sum_{j=1}^{m-1}H_{-j/m}=-m\log(m)\tag{4}
$$
Using identity $(7)$ from this answer,
$$
\begin{align}
\pi\cot(\pi z)
&=\sum_{k\in\mathbb{Z}}\frac1{k+z}\\
&=\sum_{k=1}^\infty\left(\frac1{k-1+z}-\frac1{k-z}\right)\\
&=\sum_{k=1}^\infty\left(\frac1k-\frac1{k-z}\right)-\sum_{k=1}^\infty\left(\frac1k-\frac1{k+z-1}\right)\\
&=H_{-z}-H_{z-1}\tag{5}
\end{align}
$$
which implies
$$
H_{-j/m}-H_{-(m-j)/m}=\pi\cot\left(\frac{\pi j}{m}\right)\tag{6}
$$

Using $(4)$ and $(6)$ for $m=3$ yields
$$
H_{-1/3}+H_{-2/3}=-3\log(3)\tag{7}
$$
and
$$
H_{-1/3}-H_{-2/3}=\pi\cot\left(\frac\pi3\right)\tag{8}
$$
Averaging $(7)$ and $(8)$ yields
$$
H_{-1/3}=-\frac32\log(3)+\frac\pi{2\sqrt3}\tag{9}
$$
Finally,
$$
\begin{align}
\sum_{n=1}^\infty\frac1{3n(3n-1)}
&=-\frac13\sum_{n=1}^\infty\left(\frac1n-\frac1{n-1/3}\right)\\
&=-\frac13H_{-1/3}\\[6pt]
&=\frac12\log(3)-\frac\pi{6\sqrt3}\tag{10}
\end{align}
$$

Values for Future Reference
Using $(4)$ and $(6)$, we can also compute
$$
\begin{align}
H_{-1/4}&=\pi/2-3\log(2)\\
H_{-1/2}&=-2\log(2)\\
H_{-3/4}&=-\pi/2-3\log(2)
\end{align}\tag{11}
$$
A: Answer to Second Question
We remark that
$$\Gamma(z)\zeta(z)=\int_0^\infty \frac{u^{z-1}}{e^u-1}du\tag{1}$$
$$\sum_{n=2}^\infty \frac{t^n}{(n-1)!}=(e^t-1)t \tag{2}$$
$$\psi_0(s+1)=-\gamma+\int_0^1 \frac{1-x^s}{1-x}dx \tag{3}$$ 
where $\psi_0(s)$ is Digamma Function and $\gamma$ is the Euler's Constant.
Then
$$\begin{align*}
\sum_{n=2}^\infty \frac{\zeta(n)}{k^n} &= \sum_{n=2}^\infty \frac{1}{k^n \Gamma(n)}\int_0^\infty \frac{u^{n-1}}{e^u-1}du\\
&=\int_0^\infty\frac{1}{u(e^u-1)}\left( \sum_{n=2}^\infty \frac{1}{(n-1)!}\left(\frac{u}{k}\right)^n\right)du \\
&= \frac{1}{k}\int_0^\infty \frac{e^{\frac{u}{k}}-1}{e^u-1}du \tag{4}
\end{align*}$$
Substituting $t=e^{-u}$, we get
$$
\begin{align*}
\sum_{n=2}^\infty \frac{\zeta(n)}{k^n}&=-\frac{1}{k}\int_0^1 \frac{1-t^{-1/k}}{1-t}dt \\
&= \frac{-1}{k}\left\{ \gamma+\psi_0 \left(1-\frac{1}{k} \right)\right\} \tag{5}
\end{align*}
$$

If $k(\geq 2)$ is an integer, equation $(5)$ can be further simplified using Gauss' Digamma Theorem.
$$\sum_{n=2}^\infty \frac{\zeta(n)}{k^n}=-\frac{\pi}{2k}\cot \left( \frac{\pi}{k}\right)+\frac{\log k}{k}-\frac{1}{k}\sum_{m=1}^{k-1}\cos \left(\frac{2\pi m}{k} \right) \log \left( 2\sin \frac{\pi m}{2}\right)$$
