It's a well known fact that $$\sum_{k=2}^{\infty} \zeta(k) x^k=-x \psi(1-x)-x\gamma \space (|x|<1) $$ but I didn't meet yet a version for squared Riemann zeta function $$\sum_{k=2}^{\infty}\zeta(k)^2 x^k$$
Do you know such a generating function? If yes, what is this one and how to derive it?
On
There are really no miracles here, the best we can do is to use a double integral:
$$\zeta(k+1)^2=\frac{1}{(k!)^2}\int _0^1\int _0^1\frac{(\log(u)\log(v))^k}{(1-u)(1-v)}\,du\,dv$$
which stems from one of the polylogarithm formulae valid for the positive integers (source is Wikipedia).
Now, let's set the generating function to the below form for convenience:
$$\sum _{k=1}^{\infty} x^{2k}\zeta(k+1)^2$$
A series with squared factorials involves the Bessel function, and the one which is useful for our case has order 0:
$$J_0(x)-1=\sum_{k=1}^{\infty}\frac{(-1)^k}{(k!)^2}\left(\frac{x}{2}\right)^{2k}$$
So, after we put it all together, we find the below not so simple formula, which should be valid for $|x|<$1:
$$\sum _{k=1}^{\infty} x^{2k}\zeta(k+1)^2=\int _0^{1}\int _0^{1}\frac{-1+J_0\left(2i x\sqrt{\log(u)\log(v)}\right)}{(1-u)(1-v)}\,du\,dv$$
I've tested it numerically and it's matching.
It's possible that this formula reduces to a simpler, more human-like, form if the zeta at the odd integers are not part of the series.
Here's my crack at it...
Recall $$\zeta^2(s)=\sum_{k\ge1}\frac{\sigma_0(k)}{k^s},\ \quad\ \psi(s+1)+\gamma=\int_0^1\frac{1-x^s}{1-x}\ dx,\qquad\sigma_k(n)=\sum_{d\mid n}d^k$$ Now\begin{align}\sum_{k\ge2}\zeta^2(k)x^k&=\sum_{k\ge2}x^k\sum_{n\ge1}\frac{\sigma_0(n)}{n^k}=\sum_{n\ge1}\sigma_0(n)\frac{x^2/n^2}{1-x/n}\\ &=\sum_{n\ge1}x\,\sigma_0(n)\bigg(\frac1{n-x}-\frac1n\bigg)=x\sum_{n\ge1}\sigma_0(n)\int_0^xy^{n-x-1}\,dy\\ &=x\int_0^xy^{-x-1}\sum_{n\ge1}\frac{y^n}{1-y^n}\,dy=x\int_0^{x^n}\sum_{n\ge1}\frac{z^{(n-x-1)/n}}{n(1-z)}z^{1/n-1}\,dz\tag{$z=y^n$}\\ &=x\int_0^{x^n}\sum_{n\ge1}\frac{z^{-x/n}-1+1}{n(1-z)}\,dz=x\sum_{n\ge1}\frac1n\bigg[-\ln(1-x^n)-\upsilon(x^n,-x/n)\bigg]\end{align} with $\upsilon(a,b):=\int_0^a\frac{z^b-1}{z-1}\,dz\quad$(I used a Lambert series for the third line.) I don't know where to go from here, but with Gauss's digamma theorem the summand becomes elementary for rational $x$.
[Edit] I found a simpler formulation using a clever trick:\begin{align}\sum_{n\ge1}\sigma_0(n)\frac{x^2}{n(n-x)}&=x^2\sum_{k,n\ge1}\frac1{kn(kn-x)}=x^2\sum_{k\ge1}\frac1{k^2}\sum_{n\ge1}\frac1{n(n-x/k)}\\ &=\fbox{$x\sum_{k\ge1}\frac1k\bigg(-\gamma-\psi\Big(1-\frac xk\Big)\bigg)$}\end{align} and also derived a functional equation for $\,\upsilon(a,b)\;$that may allow the original series to be analytically continued:\begin{align}\upsilon(a,b)&=\Big[\int_0^1+\int_1^a\Big]\frac{z^b-1}{z-1}\,dz=\gamma+\psi(b+1)+\int_{1/a}^1\frac{du}{u^2}\cdot\frac{u^{-b}-1}{1/u-1}\tag{$u=1/z$}\\ &=\gamma+\psi(b+1)+\int_{1/a}^1\cdot\frac{u^{-b}-u+u-1}{u(1-u)}\,du\\ &=\gamma+\psi(b+1)-\ln\frac1{1/a}+\bigg[\int_0^1-\int_0^{1/a}\bigg]\frac{u^{-b-1}-1}{1-u}\,du\\ &=\gamma-\ln a+\psi(b+1)-\big(\gamma+\psi(-b)\big)+\upsilon(1/a,-b-1)\\ &=-\pi\cot b\pi-\ln a+\upsilon(1/a,-b-1) \end{align}(the last line follows from Digamma's reflection formula.)