Background: In a question about the "sum of sums of $k$th powers of first natural numbers" someone asked -in principle- for the 2'nd iteration of the Faulhaber-problem. The Faulhaber's problem can be stated as:
Find polynomials $f_p(n)$ in $n$ for each $p$ such that $f_p(n)$ represents equivalently $$ S_p(n) = \sum_{k=1}^n k^p \tag 1$$ The second iteration is then assumed as to replace $k^p$ with the sums $S_p(k)$ $$ SS_p(n) = \sum_{k=1}^n S_p(k) \tag 2$$ In my answer I gave a solution, where I used the Faulhaber-polynomials organized in a matrix $G$ and that matrix been taken to the second power. The correctness of the solution can then nicely be seen using some cases of $n$ and the matrix to some handy size of, say $32 \times 32$ or so, which shows the results $SS_p(n)$ in the $p$th row of the result vector.
Current question: Having (2) now as proper second iteration I asked myself -just for curiousity- whether one can define also a "half-iterate" for this problem. But how can a "half iterate" of $S_p(n)$ be seen as somehow sum of $p$'th powers of $k$ at all? Surely such sums must have a more intricate structure than in formula (2) ? I have possibly an ansatz which interprets the "half iterate" $R_p(n)$ of the sum-of-like-powers $S_p(n)$ for small $n$ and $p$ in the interpretation as sum like $$\begin{array} {cccc} R_p(1) &= &1 \\ R_p(2) &= &1 \cdot 2^p &+ \frac 12 \cdot 1^p \\ R_p(3) &= &1 \cdot 3^p &+ \frac 12 \cdot 2^p& + \frac 38 \cdot 1^p \\ R_p(4) &= &1 \cdot 4^p &+ \frac 12 \cdot 3^p& + \frac 38 \cdot 2^p &+ \frac{10}{32} \cdot 1^p\\ R_p(5) &= &1 \cdot 5^p &+ \frac 12 \cdot 4^p& + \frac 38 \cdot 3^p &+ \frac{10}{32} \cdot 2^p&+ \frac{35}{128} \cdot 1^p\\ \vdots &=&\vdots&\vdots&\vdots&\vdots&\vdots &\ddots \end{array}$$ where the new fractional coefficients in each row seem to be $\binom {2j+1}{j} \cdot \frac 1{2^{2j+1}} $
I guessed that expressions with the help of a roughly approximate matrix-squareroot of $G$ which contain now coefficients for series instead of Faulhaber-style polynomials and evaluation of the numerical results (implying divergent summation of series with alternating signs). (See my ansatz in my own answer 1 below.)
Q: - can this interpretation of $R_p(n)$ as analogon of some half-iterate of the Faulhaber-formulae be justified?
- Or are there possibly more convincing ones?
Answer 2
There seems to be a general solution for this problem: for fractional "iterates" of any order, and possibly even for real "iterates". To show a solution for fractional "iterates" I'll simplify the formula in the question setting the exponent $p=1$ .
Let us assume an infinite set of coefficients $A=\{a_0,a_1,a_2,...\}$ to be determined, such that
$$ f(n) = a_0 \cdot n + a_1 \cdot (n-1) + ... + a_{n-1} \cdot 1 \tag 1$$ then to generalize this that let's introduce a notation for the "iteration" "height" $h$ as $$ f°^{h+1}(n) = a_0 \cdot f°^h(n) + a_1 \cdot f°^h (n-1) + ... + a_{n-1} \cdot f°^h(1) \tag 2 $$ with $f°^h(0)=0$ for all $h$. (From this follows also $f°^h(1)=1$ for all $h$.) Let's identify notationally $f°^1(n)=f(n)$ from now.
solution for "half-iterates"
Now, if we want to make the second "iterate" equal to the formula for the sum-of-consecutive-$n$ ($s(n)=n+(n-1)+...+1$) we can expand the formula for the second iteration symbolically using $$ \begin{array} {} f°^2(n) &=& a_0 \cdot f(n) + a_1 \cdot f (n-1) + ... + a_{n-1} \cdot f(1) \\ & = &a_0 \cdot (a_0 \cdot n + a_1 \cdot (n-1) + ... + a_{n-1} \cdot 1) \\ && + a_1 \cdot (a_0 \cdot (n-1) + a_1 \cdot (n-2) + ... + a_{n-2} \cdot 1) \\ && + a_2 \cdot (a_0 \cdot (n-2) + a_1 \cdot (n-3) + ... + a_{n-3} \cdot 1) \\ && \vdots \\ && + a_{n-1} \cdot (a_0 \cdot 1) \end{array} \tag 3 $$ The coefficients can now be determined iteratively: we assume $n=1$ first and get $$ \begin{array} {} f°^2(1) &=& a_0 \cdot f(1) + a_1 \cdot f (0) + ... \\ & = &a_0 \cdot 1 + a_1 \cdot 0 + a_2 \cdot 0 + ... \\ &=& a_0 \end{array} \tag {4.0}$$ From comparision with $s(1)=1$ we conclude $a_0=1$.
Now we use $n=2$ $$ \begin{array} {} f°^2(2) &=& 1\cdot f(2) + a_1 \cdot f (1) + 0 \\ & = & 1 \cdot (a_0 \cdot 2 + a_1 \cdot 1) + a_1 \cdot 1 \\ &=& 2 + 2 \cdot a_1 \\ \end{array} \tag {4.1a}$$ To equal this with $s(2)=3$ we solve $$ 3 = 2 + 2 a_1 \\ 1 = 2 a_1 \\ a_1 = 1/2 \tag {4.1b} $$ This procedure repeated gives us the set of coefficients $A=\{1,1/2,3/8,10/32...\}$ with as many known coefficients as needed which I had already guessed by the "square-root of matrix G" ansatz.
For this specific function $f$ we can even find an independent formula for the set of coefficients, namely $$ a_0=1 \qquad \qquad a_k = \binom{2k-1}{k} \cdot 2^{1-2k} \tag 5$$ which is suggested by the entry in OEIS.
So far we have used that algorithm with a simplified function $f(n)$ where we ignored exponents at $n, (n-1),...$ in $f°^1(n)=a_0 \cdot n + a_1 \cdot (n-1) + ... $. It is straightforward to show that the same coefficients $A$ occur, if we use similarly the generalized $s(n)$ (containing the same exponents) in the comparision-step.
general fractional "iterates"
Now what we've done to arrive at a set of coefficients $A$ to get some "half-iterate", we can analoguously determine such set for the "1/3-iterate" when we simply compare the values of $s(n)$ with $f°^3(n)$ - and the same can be made with any "$1/m$-iterate" for $m \in \mathbb N$ .
Here I show the beginning of the list of sets-of-coefficients $A_m$ for $m=1..8$
$\tag {6.1}$
That coefficients can be normalized by multiplying them by powers of $m^2$. This gives the following decomposition: $\tag {6.2}$
The sequences of numerators of the first couple of sets are also in the OEIS, for instance
alternative determination of coefficients $a_{m,k}$
That coefficients in row $m$ seem to fit polynomials in $1/m$ taken from the Stirlingnumbers $1$'st kind. Let $S1$ be the matrix of unsigned Stirling numbers $1$'st kind $S1 = [s_1(r,c)]_{r,c=0}^\infty$ such that $\tag {7.1}$
is the top left of $S1$ with row- and column-indexes beginning at zero, then the coefficients $a_{m,k}$ of the sets $A_m$ seem to be determined by
$$ a_{m,k} = \frac 1{m!} \cdot \sum_{c=0}^m {s1(m,c)\over m^c} \tag {7.2} $$
"real iterates"
In the comments at the OEIS-entries we find hints to generating-functions for those coefficients. Adapting that properly we find for the set $A_m$ the general generating function $$ \mathcal{gf}(A_m) = (1-x)^{- 1/m} \tag {7.3} $$ and because of that simple expression we might even generalize $m$ from natural numbers to real-numbers; for instance the set of coefficients $A_m$ such that "the $\pi$'th iterate" equals $s(n)$ are $A_\pi =\{1, 0.318310, 0.209816, 0.162139, 0.134507, 0.116169, 0.102970, 0.0929424,...\}$