Recently, interesting problems related to sums over the zeta function reappeared, see e.g. [1] or [2].
In approaching these problems it is useful to study the generating sums. The generating sum of the zeta function is well known (see formula $(3)$ in [1]):
$$g_1(z) :=\sum_{n=2}^{\infty}\zeta(n)z^n = - z H_{-z}\tag{1}$$
For ease of reference and because it is similar to what follows I repeat the derivation here:
$$\begin{align} \sum_{n=2}^{\infty}\zeta(n)z^n & =\sum _{n=2}^{\infty } z^n \left(\sum _{k=1}^{\infty } \frac{1}{k^n}\right)=\sum _{k=1}^{\infty } \left(\sum _{n=2}^{\infty } \frac{z^n}{k^n}\right)\\ & =\sum _{k=1}^{\infty } \frac{z^2}{k (k-z)}\\ &=-z \sum _{k=1}^{\infty } \left(\frac{1}{k}-\frac{1}{k-z}\right)=-z H_{-z} \end{align}\tag{1a}$$
Here $H$ is the harmonic number.
It seems natural to ask for the generating sum of the square of the zeta function:
$$g_2(z) = \sum_{n=2}^{\infty}\zeta(n)^2 z^n\tag{2}$$
My (unsuccessful) attempts are
(i) Harmonic number form
Using $\zeta(n) = \sum_{k \ge 1} k^{-n}$ gives, using $(1)$, recursively
$$\begin{align}g_{2,1}(z) & =\sum _{n=2}^{\infty } \sum _{k=1}^{\infty } \zeta (n) \left(\frac{z}{k}\right)^n=\sum _{k=1}^{\infty } \sum _{n=2}^{\infty } \zeta (n) \left(\frac{z}{k}\right)^n\\ & =\sum _{k=1}^{\infty } g_1\left(\frac{z}{k}\right)=\sum_{k=1}^{\infty }\left( -\frac{z}{k}\right) H_{-\frac{z}{k}}\end{align}\tag{3}$$
(ii) Bessel form
Using
$$\zeta (n)=\frac{1}{\Gamma (n)}\int_0^{\infty } \frac{t^{n-1}}{e^t-1} dt\tag{4}$$
the kernel is
$$\frac{1}{\Gamma (n)^2} \int _0^{\infty }\int _0^{\infty }\frac{t^{n-1} u^{n-1} z^n}{\left(e^t-1\right) \left(e^u-1\right)}dudt\tag{5}$$
and the $n$-sum becomes
$$\begin{align}& \sum _{n=2}^{\infty } \frac{z^n \left(t^{n-1} u^{n-1}\right)}{\Gamma (n)^2 \left(e^t-1\right) \left(e^u-1\right)}\\ & =\frac{z \left(I_0\left(2 \sqrt{t} \sqrt{u} \sqrt{z}\right)-1\right)}{\left(e^t-1\right) \left(e^u-1\right)} \end{align}\tag{6}$$
so that the generating sum has a representation of a double integral over a Bessel function:
$$g_{2,2}(z) = \int _0^{\infty }\int _0^{\infty }\frac{z \left(I_0\left(2 \sqrt{t} \sqrt{u} \sqrt{z}\right)-1\right)}{\left(e^t-1\right) \left(e^u-1\right)}dudt\tag{8}$$
Here I am stuck.
Questions
- Can you do better and find $g_2(z)$?
- Or, alternatively, can you find $g_2(\frac{1}{2})$ or $g_2(\frac{1}{4})$?
This is not an answer, since it neither provides a closed expression nor does it prove that there is none, but it is a continuation of reformulatinons which lead to perhaps interesting results.
§1. The divisor function $d(n)$ appears naturally
User C-RAM, in a comment, pointed out thet the square of the zeta function can be identified as the product of two Dirichlet series
$$\begin{align} (f*g)(s)& = \sum_{j\ge 1} \frac{f(j)}{j^s}\sum_{k\ge 1} \frac{g(k)}{k^s} \\ =& \sum_{j,k\ge 1}f(j) g(k) \frac{1}{(j k)^s}\overset{j k \to n}=\sum_{n \ge 1, k \ge 1\land k|n}f(n) g(\frac{n}{k}) \frac{1}{n^s}\end{align}\tag{s1} $$
specifically, for $g(k)=1$,
$$\begin{align} (f*1)(s)& =\sum_{n \ge 1, k \ge 1\land k|n}f(n) (1) \frac{1}{n^s}=\sum_{n \ge 1}\frac{f(n)}{n^s} \sum_{k \ge 1\land k|n}(1) =\sum_{n \ge 1}\frac{f(n)}{n^s} d(n)\end{align}\tag{s2} $$
where $d(n) = \sum_{k|n}1$ is the number of values of $k\ge 1$ which divide $n$, i.e. the number of divisors of $n$. Example $n=12$ has the 6 divisors $(1,2,3,4,6,12)$ hence $d(12)=6$.
The first few terms are $d(1..20) = \{1,2,2,3,2,4,2,4,3,4,2,6,2,4,4,5,2,6,2,6\}$
The sequence is listed in https://oeis.org/A000005. where also the generating function $\lambda(x)$ (called Lambert series) is given
$$\lambda(x)= \sum_{k\ge 1} \frac{x^k}{1-x^k}=\sum_{k\ge 1}x^k \sum_{j\ge0} x^{k j}=\sum_{k\ge 1}\sum_{j\ge 1} x^{k j}\overset{j k\to n}=\sum_{n\ge 1}\sum_{j|n} x^n=\sum_{n\ge 1}d(n) x^n\tag{s3}$$
§2. The Mittag-Leffler expansion of $g_m=\sum_{n \ge 2}\zeta(n)^m z^n$
From $s2$ we have $\zeta(n)^2 = (1*1)(n) = \sum_{k\ge 1}\frac{ d(k)}{k^n}$
This leads to another representation of $g_2$ called Mittag-Leffler expansion
$$g_2(z):=\sum_{n \ge 2}\zeta(n)^2 z^n =\sum_{n \ge 2} z^n \sum_{k \ge 1}\frac {d(k)}{k^n}=\sum_{k\ge 1}d(k) \sum_{n\ge 2}(z/k)^n = \sum_{k\ge 1} \frac{z^2 d(k)}{k(k-z)}\tag{s4}$$
For the $m$-th power of $\zeta$ we obtain similarly
$$g_m(z) := \sum_{n \ge 2}\zeta(n)^m z^n=\sum_{k\ge 1} \frac{z^2 d(m,k)}{k(k-z)}\tag{s5}$$
where, recursively,
$$d(m,k) = \left\{ \begin{array}{ll} 1 & m=1 \\ \sum_{j|k}d(m-1,j) & m\ge2 \\ \end{array} \right. \tag{s6}$$
§3. Attempt to calculate $g_2(1/2)$ using the Mittag-Leffler expansion
From $(s4)$ we have
$$g_2(1/2) = \frac{1}{4}\sum_{k\ge 1} \frac{d(k)}{k(k-\frac{1}{2})}=\frac{1}{2}\sum_{k\ge 1} \frac{d(k)}{k(2k-1)}\tag{s7}$$
Now it would be nice to have an explicit formula for $d(n)$. But we are in an even better position: we have its generating function $\lambda(x)$.
And we just need to generate the denominator which can be easily done observing that
$$\frac{1}{k(2k-1)} = \int_0^1 \frac{1}{p^2}\int_0^{p^2} \frac{x^k}{x} \, dx \, dp\tag{s8}$$
Hence we have to calculate the summand of $\lambda$, i.e.
$$\int_0^1 \frac{1}{p^2}\int_0^{p^2} \frac{x^k}{x(1-x^k)} \, dx \, dp\tag{s9}$$
The $x$-integral gives
$$\int_0^{p^2} \frac{x^k}{x \left(1-x^k\right)} \, dx=-\frac{\log \left(1-p^{2 k}\right)}{k}\tag{s10}$$
the remaning $p$-integral is then ... suspense ...
$$g_{2,3}=\int_0^1 -\frac{\log \left(1-p^{2 k}\right)}{k p^2} \, dp = -\frac{1}{k} H_{-\frac{1}{2 k}}\tag{s11}$$
and disappointment! We have (taking into account the temporarily suppressed factor $\frac{1}{2}$) just rediscovered the formula $(3)$ of the OP, hence $g_{2,3}=g_{2,1}$.
Little reward: at least one formula for the harmonic number popped up which was new to me and can be written as
$$H_{z} = -z \int_0^1 t^{z-1} \log (1-t) \, dt\tag{s12}$$
§4. Micellaneous
The function q-Polygamma and the relation to $\lambda(x)$
The q-digamma function is defined as
$$\psi_q(z) = -\log(1-q) + \log(q)\sum_{n\ge 0} \frac{q^{n+z}}{1-q^{n+z}}\tag{s13}$$
Hence we have
$$\lambda(x) = \frac{1}{\log(x)}\left(\psi_q(z=1) +\log(1-x) \right)\tag{s14}$$
3rd power of $\zeta$
$d(3,n)$ is listet in https://oeis.org/A007425 as A007425 d_3(n), or tau_3(n), the number of ordered factorizations of n as n = r s t.
G.f.: $\sum_{k>=1} d(k)*x^k/(1-x^k)$.
The generating function of $d(3,n)$ is given by
$$\sum_{n\ge 1} x^n d(3,n) = \sum_{k>=1} d(k) \frac{x^k}{1-x^k}$$