calculus - Closed-form of $sum_{k=1}^{infty }left(psi_1(k)right)^n$

Inspired by answers to this question, for which $n$ values could we specify a closed-form of

$$S(n)=\sum_{k=1}^{\infty }\left(\psi_1(k)\right)^n\,?$$

Here $\psi_1$ is the trigamma function, and $n\geq2$ is an integer.

I think for small $n$ values there is a closed-form. From the answers above we know that $S(2)=3\zeta(3)$. I think there is a generalization of Olivier Oloa's approach using techniques from Pedro Freitas' paper. Note that there is a known error in the paper. Furthermore I think robjohn's answer has a generalization using results from the paper by Philippe Flajolet and Bruno Salvy.

Question. Is there a closed-form of $S(n)$ for $2\leq n \leq 7\,?$

Answer

Using
$$\psi_1(n) = \sum_{k\geq n}\frac1{k^2},$$
write
$$S(3) = \sum_{k\geq1}\sum_{i,j,l\geq k}\frac{1}{i^2j^2l^2}.$$
The sum over $i,j,l\geq k$ decomposes into the sums
$$\{i=j=l\geq k\} + 3\{i=j, l\geq k ; i\neq l\} + \{i,j,l\geq k; i\neq j, i\neq l, j\neq l\}.$$
Hence
$$S(3) = \sum_{i\geq k\geq 1}\frac{1}{i^6} + 3\sum_{i\geq k, j\geq k, i\neq j}\frac{1}{i^4j^2} + \sum_{i,j,l\geq k, i\neq j, j\neq l, i\neq l} \frac{1}{i^2j^2l^2}.$$
Doing the sum over $k$ first, and reordering the summation appropriately gives

$$\begin{aligned} S(3) &= \sum_{i\geq 1}\frac1{i^5} + 3\sum_{i>l\geq1}\frac{1}{i^4l} + 3\sum_{i>l\geq1}\frac{1}{i^2l^3} +6\sum_{i>j>l\geq1}\frac{1}{i^2j^2l} \\&= \zeta(5) + 3\zeta(4,1) + 3\zeta(2,3) + 6\zeta(2,2,1) \end{aligned}$$
in terms of the multivariate zeta function.

Simplifying, this analysis gives:
$$S(2) = 3\zeta(3).$$
$$S(3) = 9 \zeta (3) \zeta (2)-\tfrac{25}{2} \zeta (5).$$
$$S(4) = 10 \zeta (5) \zeta (2)+51 \zeta (3) \zeta (4)- \tfrac{301}{4}\zeta (7)$$

$$S(5) = -\tfrac{1505}{4} \zeta (7) \zeta (2)+\tfrac{125}{2} \zeta (5) \zeta (4)+\tfrac{835}{4} \zeta (3) \zeta (6)+\tfrac{11791}{36} \zeta (9)-10 \zeta (3)^3.$$
Here I found $S(4)$ and $S(5)$ numerically without proof. For $S(m)$ for $m\geq 6$ there is probably no closed form except in terms of irreducible multivariate zeta values.

For example, if you allow closed forms in terms of double Euler sums, then
$$S(6) = \tfrac{820230}{901} \zeta (3) s_h(2,6)+\tfrac{192}{901} s_h(8,3)+\tfrac{10}{901} s_h(9,2)-\tfrac{1953509}{2703} \zeta (9) \zeta (2)-\tfrac{762298}{901} \zeta (3)^3 \zeta (2)-\tfrac{19487373}{7208} \zeta (7) \zeta (4)+\tfrac{48455}{34} \zeta (5) \zeta (6)-\tfrac{8977165}{1272} \zeta (3) \zeta (8)+\tfrac{1790142}{901} \zeta (3)^2 \zeta (5)+\tfrac{12188899}{3604} \zeta (11)$$
$$S(7) = \tfrac{6143922}{1007} \zeta (3) \zeta (2) s_h(2,6)+\tfrac{29361277}{2014} \zeta (5) s_h(2,6)-\tfrac{2696}{1007} \zeta (3) s_h(7,3)+\tfrac{3138}{1007} \zeta (3) s_h(8,2)+\tfrac{273}{1007} \zeta (2) s_h(8,3)-\tfrac{84}{1007} \zeta (2) s_h(9,2)+\tfrac{1046978}{1007} s_h(5,8)-\tfrac{805894}{1007} s_h(6,7)-\tfrac{1796398}{1007} s_h(7,6)-\tfrac{1530997}{1007} s_h(8,5)-\tfrac{720964}{1007} s_h(9,4)-\tfrac{144069}{1007} s_h(10,3)+\tfrac{12}{1007} s_h(11,2)-\tfrac{4525536}{1007} \zeta (11) \zeta (2)-\tfrac{11096255}{2014} \zeta (3)^2 \zeta (5) \zeta (2)-\tfrac{14173495}{2014} \zeta (9) \zeta (4)-\tfrac{19952365}{1007} \zeta (3)^3 \zeta (4)+\tfrac{8234835}{152} \zeta (7) \zeta (6)+\tfrac{339826375}{2544} \zeta (5) \zeta (8)-\tfrac{9006237625}{32224} \zeta (3) \zeta (10)+\tfrac{15302200}{1007} \zeta (3) \zeta (5)^2+\tfrac{42533615}{2014} \zeta (13)-\tfrac{39044583}{1007} \zeta (3)^2 \zeta (7)$$