Permit me a brief introduction before I state the question, three questions in fact.
Inspired by this MSE link I computed the following harmonic sum identities: $$1/6\, \left( {H_{{n}}}^{(1)} \right) ^{3}-1/2\,{H_{{n}}}^{(2)}{H_{{n}}}^{(1)}+1/3\,{H_{{n}}}^{(3)} = \frac{1}{n!} \left[n+1\atop 4\right]$$ and $$1/24\, \left( {H_{{n}}}^{(1)} \right) ^{4}-1/4\,{H_{{n}}}^{(2)} \left( {H_{{n}}}^{(1)} \right) ^{2}+1/3\,{H_{{n}}}^{(3)}{H_{{n}}}^{(1)} \\+1/8\, \left( {H_{{n}}}^{(2)} \right) ^{2}-1/4\,{H_{{n}}}^{(4)} = \frac{1}{n!} \left[n+1\atop 5\right]$$ and finally $$\frac{1}{120}{ \left( {H_{{n}}}^{(1)} \right) ^{5}} -1/12\,{H_{{n}}}^{(2)} \left( {H_{{n}}}^{(1)} \right) ^{3}+1/6\,{H_{{n}}}^{(3)} \left( {H_{{n}}}^{(1)} \right) ^{2}+1/8\,{H_{{n}}}^{(1)} \left( {H_{{n}}}^{(2)} \right) ^{2}\\-1/4\,{H_{{n}}}^{(4)}{H_{{n}}}^{(1)} -1/6\,{H_{{n}}}^{(2)}{H_{{n}}}^{(3)}+1/5\,{H_{{n}}}^{(5)} = \frac{1}{n!} \left[n+1\atop 6\right].$$ The parameter here call it $q$ takes on the values $q=3,4$ and $q=5.$
I verified the above identities numerically and they seem to hold. The values of the coefficients indicate that the Polya Enumeration Theorem plays a role here, and indeed it is not difficult to present a general conjecture for the above pattern which is
$$\color{#A00}{\left.Z(P_q) \left(\sum_{k=1}^n \frac{1}{k^{s}} \right)\right|_{s=1} = \frac{1}{n!} \left[n+1\atop q+1 \right].}$$
We use $Z(P_q)$ to refer to the cycle index of the unlabeled set operator $\mathfrak{P}_{=q}$, which admits the following simple recursive definition:
$$Z(P_0) = 1 \quad\text{and}\quad Z(P_n) = \frac{1}{n} \sum_{l=1}^n (-1)^{l+1} a_l Z(P_{n-l}).$$
Recall that we must distinguish between multisets ($\mathfrak{M}$) and sets ($\mathfrak{P}$). Here are the cycle indices $Z(P_3), Z(P_4)$ and $Z(P_5):$ $$\begin{array}{|l|l|} \hline Z(P_3) & \frac{1}{6}\,{a_{{1}}}^{3}-1/2\,a_{{2}}a_{{1}}+1/3\,a_{{3}}\\ \hline Z(P_4) & \frac{1}{24}\,{a_{{1}}}^{4}-1/4\,a_{{2}}{a_{{1}}}^{2} +1/3\,a_{{3}}a_{{1}}+1/8\,{a_{{2}}}^{2}-1/4\,a_{{ 4}}\\ \hline Z(P_5) & {\frac {1}{120}}\,{a_{{1}}}^{5}-\frac{1}{12}\,a_{{2}}{a_{{1}}}^{3} +1/6\,a_{{3}}{a_{{1}}}^{2}+1/8\,a_{{ 1}}{a_{{2}}}^{2}-1/4\,a_{{4}}a_{{1}} -1/6\,a_{{2}}a_{{3}}+1/5\,a_{{5}}\\ \hline\end{array}$$
The substitution mechanism for PET here is $$a_l = \sum_{k=1}^n \frac{1}{k^{ls}}$$ and it should be immediate to see how the formulas were obtained from the cycle indices.
The proof of the conjecture now morphs into a simple one-liner. The LHS is the ordinary generating function of sets having cardinality $q$ chosen from the repertoire $\{1, 1/2^s, 1/3^s, \ldots, 1/n^s\}$ evaluated at $s=1.$ By inspection this same generating function is also given by
$$[z^q] \prod_{k=1}^n \left(1 + \frac{z}{k^s}\right).$$ Putting $s=1$ we obtain $$[z^q] \prod_{k=1}^n \left(1 + \frac{z}{k}\right) = \frac{1}{n!} [z^q] \prod_{k=1}^n \left(z+k\right) = \frac{1}{n!} [z^q] \frac{1}{z} \prod_{k=0}^n \left(z+k\right) \\= \frac{1}{n!} [z^{q+1}] \prod_{k=0}^n \left(z+k\right) = \frac{1}{n!} \left[ n+1\atop q+1 \right].$$
We have used the ordinary generating function of the Stirling numbers of the first kind in the last step of this calculation and obtained the RHS of the conjectured formula which completes the proof.
These are the three questions. First, I'd be interested in possible references for this formula, second, I'd like to know if perhaps it admits a more elementary proof, and third does this formula have computational uses e.g. it would appear to produce asymptotic expansions for certain Stirling numbers.
Thanks for reading.
Here is the Maple code that was used to verify the three formulas presented at the beginning. Even $q=25$ gives instant results. The combinatorics package is included for its Stirling number routines.
with(combinat);
pet_cycleind_set :=
proc(n)
local p, s;
option remember;
if n=0 then return 1; fi;
expand(1/n*add((-1)^(l-1)*a[l]*pet_cycleind_set(n-l), l=1..n));
end;
v :=
proc(n, q)
local sl, l;
option remember;
sl := [seq(a[l]=harmonic(n, l), l=1..q)];
subs(sl, pet_cycleind_set(q));
end;
Addendum. Here are the asymptotic expansions for the initial three identities: $$\frac{1}{n!}\left[n+1\atop 4\right] \sim 1/6\, \left( \ln \left( n \right) \right) ^{3}+1/2\, \left( \ln \left( n \right) \right) ^{2}\gamma+ \left( 1/2\,{\gamma}^{2}-1/12 \,{\pi }^{2} \right) \ln \left( n \right)\\ +1/6\,{\gamma}^{3}-1/12\, {\pi }^{2}\gamma+1/3\,\zeta \left( 3 \right)$$ and $$\frac{1}{n!}\left[n+1\atop 5\right] \sim 1/24\, \left( \ln \left( n \right) \right) ^{4}+1/6\,\gamma\, \left( \ln \left( n \right) \right) ^{3}\\+ \left( 1/4\,{\gamma}^{2 }-1/24\,{\pi }^{2} \right) \left( \ln \left( n \right) \right) ^{ 2}+ \left( 1/6\,{\gamma}^{3}-1/12\,{\pi }^{2}\gamma+1/3\,\zeta \left( 3 \right) \right) \ln \left( n \right) \\+1/24\,{\gamma}^{4} -1/24\,{\pi }^{2}{\gamma}^{2}+1/3\,\zeta \left( 3 \right) \gamma+{ \frac {{\pi }^{4}}{1440}}$$ and finally $$\frac{1}{n!}\left[n+1\atop 6\right] \sim {\frac { \left( \ln \left( n \right) \right) ^{5}}{120}}+1/24\, \gamma\, \left( \ln \left( n \right) \right) ^{4}+ \left( 1/12\,{ \gamma}^{2}-{\frac {{\pi }^{2}}{72}} \right) \left( \ln \left( n \right) \right) ^{3}\\+ \left( 1/12\,{\gamma}^{3}-1/24\,{\pi }^{2} \gamma+1/6\,\zeta \left( 3 \right) \right) \left( \ln \left( n \right) \right) ^{2}\\+ \left( 1/24\,{\gamma}^{4}-1/24\,{\pi }^{2}{ \gamma}^{2}+1/3\,\zeta \left( 3 \right) \gamma+{\frac {{\pi }^{4}}{ 1440}} \right) \ln \left( n \right)\\ +{\frac {{\gamma}^{5}}{120}}-{ \frac {{\pi }^{2}{\gamma}^{3}}{72}}+1/6\,\zeta \left( 3 \right) { \gamma}^{2}+{\frac {\gamma\,{\pi }^{4}}{1440}}-1/36\,{\pi }^{2} \zeta \left( 3 \right) +1/5\,\zeta \left( 5 \right).$$ The Maple code for these is as follows.
aex :=
proc(q)
local sl, l;
option remember;
sl := [a[1]=log(n)+gamma,
seq(a[l]=Zeta(l), l=2..q)];
subs(sl, pet_cycleind_set(q));
end;
