Finding $\sum_{n=2}^{\infty} (\zeta^k(n) - 1)$ for $k \in \mathbb{Z}_{\ge 0}$ I am interested in finding the value of $$\sum_{n=2}^{\infty} (P(\zeta(n)) - P(1))$$
where $P(x)$ is a polynomial in $x$. Solving this is analogous to solving $$f(k)=\sum_{n=2}^{\infty} (\zeta^k(n) - 1) $$ and then adding up the results. Trivially, $f(0) = 0$. It is also not hard to show that $f(1) = 1$, by rearranging the sums and rewriting as a sum of geometric series, such as in this answer. For calculating $f(2)$, I did as follows: $$\sum_{n=2}^{\infty}\left(\left(\sum_{m=1}^{\infty} \frac{1}{m^n}\right)^2-1\right)$$
Expanding the squared term yields $$\sum_{n=2}^{\infty} \left(\sum_{m=2}^{\infty} \frac{1}{m^{2n}} + \sum_{m=2}^{\infty}\sum_{p=1}^{m-1}\frac{2}{m^n p^n}\right)$$
This can be simplified to $$\frac{7}{4} - \zeta(2) + 2 \sum_{n=2}^{\infty} \sum_{m=2}^{\infty}\sum_{p=1}^{m-1} \frac{1}{m^n p^n}$$
By summing over $n$, this can be simplified to $$\frac{7}{4} - \zeta(2) + 2 \sum_{m=2}^{\infty}\sum_{p=1}^{m-1} \frac{1}{mp(mp-1)}$$
The summand can be broken up to get $$\frac{7}{4} - \zeta(2) + 2 \sum_{m=2}^{\infty}\sum_{p=1}^{m-1} \left(\frac{1}{mp-1} - \frac{1}{mp}\right) = \frac{7}{4} - \zeta(2) + 2 \sum_{m=2}^{\infty} \frac{H_{m-1-\frac{1}{m}} - H_{-\frac{1}{m}} - H_{m-1}}{m}$$
Unfortunately, I wasn't able to simplify that last sum.
My questions:

$1.$ What is a closed form for $f(k)$ (or at least for $f(2)$)?
$2.$ If a closed form does not exist, what is the power series representation?

Taking derivatives with respect to $x$, I get that the series representation for $f$ is $$f(x) = \sum_{n=1}^{\infty} \sum_{m=2}^{\infty} \frac{\ln(\zeta(m))^n}{n!} x^n$$
However evaluating $\sum_{m=2}^{\infty}\ln(\zeta(m))^n$ seems to be even more difficult considering $\zeta(m)$ is defined by a sum.
 A: I've worked on this question for quite some time now, and suspect Claude Leibovici is right in the sense that a closed form probably doesn't exist. In the answer to this MO question - in which I cite your work - Fedor Petrov obtains $$S: = \sum (\zeta(n)^2-1)= 1+
\int_0^1 \sum_{k=0}^\infty\frac{x^{k}}{1+x+\ldots+x^{k+1}}dx.$$
If we then proceed by interchanging the sum and integral, we obtain: \begin{align} S &= 1 + \sum_{k=0}^{\infty} \int_{0}^{1}\frac{x^{k}}{1+x+\ldots+x^{k+1}}dx \\
&= 1 - \sum_{k=0}^{\infty} \frac{H_{- \frac{1}{k+2}}}{k+2}  \\
&= 1- \sum_{m=2}^{\infty} \frac{H_{-\frac{1}{m}}}{m}.
\end{align}
Now, we can equate the two expressions for the sum, group like terms and find
$$\sum_{m=2}^{\infty} \frac{H_{-\frac{1}{m}}}{m} = \frac{3}{4} - \zeta(2) + 2 \sum_{m=2}^{\infty} \frac{H_{m - 1- \frac{1}{m}} - H_{m-1} }{m}.  $$
Let's focus on the expression on the left. We can proceed with the definition $$ H_{x} = \sum_{k=1}^{\infty} \binom{x}{k} \frac{(-1)^{k}}{k} , $$
and consider the fact that relates this expression to the unsigned Stirling numbers of the first kind:  $$\begin{align} (-1)^n {-m \choose n} = {n+m-1 \choose n} = \frac{1}{n!} \sum_{i=0}^n \left[ {n \atop i} \right] m^i , \end{align} $$ we can replace $m$ by $\frac{1}{m}$ and rearrange sums to obtain:
$$\begin{align} 
 \sum_{m=2}^{\infty}  \frac{H_{-\frac{1}{m}}}{m} &=  \sum_{m=2}^{\infty} \sum_{k=1}^{\infty} \sum_{i=0}^{k} \frac{1}{k k!} m^{-(i+1)} \left[ {k \atop i} \right] \\
&=  \sum_{k=1}^{\infty} \sum_{i=0}^{k} \frac{1}{k k!} \left[ {k \atop i} \right] \sum_{m=2}^{\infty} m^{-(i+1)} \\
&=  \sum_{k=1}^{\infty} \sum_{i=0}^{k} \frac{1}{k k!} \left[ {k \atop i} \right] (\zeta(i+1) - 1) \\ 
&=  \bigg{(} \sum_{k=1}^{\infty} \frac{1}{kk!} \bigg{(} \sum_{i=0}^{k}  \left[ {k \atop i} \right] \zeta(i+1) -  \sum_{i=0}^{k} \left[ {k \atop i} \right] \bigg{)}  \bigg{)} \\ 
&=  \bigg{(} \sum_{k=1}^{\infty} \frac{1}{kk!} \bigg{(} \sum_{i=1}^{k} \left[ {k \atop i} \right] \zeta(i+1) - k! \bigg{)}  \bigg{)}. 
\end{align}.$$
What I found interesting about this sum is that expressions for variations of the inner sum have been obtained. To be more precise: quite a bit is known about sums of the form $$\sum_{i=1}^{n}{n \brack i}\zeta(n-i) \tag{*}$$
For instance, if we define $$K_n(m) = \overbrace{\int_0^1 \dots \int_0^1}^{n-\mathrm{times}} \left(-\frac{\ln(1-x_1x_2\cdots x_n)}{x_1 x_2  \cdots x_n}\right)^m \mathrm{d}x_1\mathrm{d}x_2 \cdots \mathrm{d}x_n , $$ then $$K_1(m) = m\sum_{n=0}^{m-1}\left[ m-1 \atop n\right]\zeta(m+1-n) . $$ Other, related results can be found over here.
Unfortunately, little is known about sums of the form $$T_k := \sum_{i=1}^{k}  \left[ {k \atop i} \right] \zeta(i+1). \tag{**} $$ I've tried finding integral expressions for $T_{k}$ by collecting integral identities that amount to $T_{1}$ or something close to it, and then generalize these integrals by considering powers of the integrands. A subset of the integrals I considered can be found in the following answer of mine.
However, none of them have yielded the result I was looking for. After correspondence with a certain number theorist I believe such an expression indeed does not exist. The mathematician in question gave some good evidence this is most probably not a 'natural' sum to consider.
Interestingly, there are sums of products of zeta values that can be evaluated and are more 'natural' in a sense. For instance, we find in equation $(133)$ of the following page that $$\lim_{n \to \infty} \sum_{k=2}^{n} \zeta(k) \zeta(n-k) x^{k-1} = \frac{1}{x} - \psi_{0}(-x) - \gamma. $$ Here, $\psi_{0}(\cdot) $ is the digamma function. Similarly, in equation $(1.6)$ of p. 24 of the following article by Dilcher we find $$\sum_{j=1}^{n-1} \zeta(n) \zeta(2n - 2j) = \left(n+\frac{1}{2}\right) \zeta(2n), \qquad n \geq 2. $$
Such sums have analogues with sums of products of Bernoulli numbers, as indicated in Dilcher's article. It might be an interesting avenue for you to explore further.
