For a paper I was working on recently I needed to find the value of the following sum:

$$S(n,k) = \sum_{i_1 = 1}^n \sum_{i_2 = i_1+1}^n \cdots \sum_{i_k=i_{k-1}+1}^n \frac{1}{i_1 i_2 \cdots i_k}.$$

I found a couple of references (by Adamchik and Cheon and El-Mikkaway) that have an expression for $S(n,k)$ as a polynomial containing generalized harmonic numbers $H_n^{(r)}$, where $$H_n^{(r)} = \sum_{j=1}^n \frac{1}{j^r}.$$ For example, $$S(n,2) = \frac{1}{2}\left(H_n^2 - H^{(2)}_n \right),$$ $$S(n,3) = \frac{1}{6}\left(H_n^3 - 3H_n H^{(2)}_n + 2 H_n^{(3)}\right),$$ $$S(n,4) = \frac{1}{24}\left(H_n^4 - 6 H^2_n H_n^{(2)} + 3 (H_n^{(2)})^2 + 8 H_n H_n^{(3)} - 6 H_n^{(4)}\right).$$

Neither of these papers considers the corresponding polynomial sequence of indeterminates (the polynomials before substituting in the generalized harmonic numbers), though. Calculations for small values of $n$ indicate that these are the cycle index polynomials of the symmetric groups, with the sign pattern such that each factor of $H_n^{(r)}$ contributes a $+1$ if $r$ is odd and a $-1$ if $r$ is even.

Could someone give a proof of this, particularly one with a combinatorial flavor that gives some real insight into why the cycle index polynomials of the symmetric groups show up here (assuming that they do)?

(For the record, I don't need this answered for my paper. I just want to know for my own sake.)

As a side note, the papers also give the extremely (and, to me, surprisingly) simple expression $$S(n,k) = \frac{1}{n!} \left[ n+1 \atop k+1 \right],$$ where $\left[ n \atop k \right]$ is an unsigned Stirling number of the first kind.