A "half-iterate" of the Faulhaber-problem of "summing like powers"?

Question

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 1

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}$

          a_0   a_1    a_2  ...
     A_1:  1     1       1        1           1              1                 1                   1                       1                        1
     A_2:  1   1/2     3/8     5/16      35/128         63/256          231/1024            429/2048              6435/32768              12155/65536
     A_3:  1   1/3     2/9    14/81      35/243         91/729          728/6561          1976/19683              5434/59049           135850/1594323
     A_4:  1   1/4    5/32   15/128    195/2048       663/8192        4641/65536        16575/262144          480675/8388608         1762475/33554432
     A_5:  1   1/5    3/25   11/125      44/625      924/15625        4004/78125        17732/390625           79794/1953125           363506/9765625
     A_6:  1   1/6    7/72  91/1296  1729/31104    8645/186624    267995/6718464    1416545/40310784     60911435/1934917632  2984660315/104485552128
     A_7:  1   1/7    4/49   20/343    110/2401      638/16807       3828/117649      164604/5764801        1028775/40353607        6515575/282475249
     A_8:  1   1/8   9/128  51/1024  1275/32768    8415/262144    115005/4194304     805035/33554432     45886995/2147483648    331406075/17179869184
     A_9:  1   1/9    5/81  95/2187   665/19683    4921/177147    113183/4782969     889295/43046721       7114360/387420489    519348280/31381059609
    A_10:  1  1/10  11/200  77/2000  2387/80000  97867/4000000  1663739/80000000  14498297/800000000  1029379087/64000000000  9264411783/640000000000

That coefficients can be normalized by multiplying them by powers of $m^2$. This gives the following decomposition: $\tag {6.2}$

 m  A_m  a_0  a_1    a_2      a_3        a_4         a_5            a_6
---+----+-----------------------------------------------------------------
 1  A_1:  1  1/1^1   1/1^3    1/1^5      1/1^7       1/1^9         1/1^11
 2  A_2:  1  1/2^1   3/2^3   10/2^5     35/2^7     126/2^9       462/2^11
 3  A_3:  1  1/3^1   6/3^3   42/3^5    315/3^7    2457/3^9     19656/3^11
 4  A_4:  1  1/4^1  10/4^3  120/4^5   1560/4^7   21216/4^9    297024/4^11
 5  A_5:  1  1/5^1  15/5^3  275/5^5   5500/5^7  115500/5^9   2502500/5^11
 6  A_6:  1  1/6^1  21/6^3  546/6^5  15561/6^7  466830/6^9  14471730/6^11

The sequences of numerators of the first couple of sets are also in the OEIS, for instance

  entry in oeis            first few coeffs  short comment              contributor
 http://oeis.org/A000012   1,1,1,1
 http://oeis.org/A001700   1,3,10,35,126
 http://oeis.org/A034171   1,6,42,315,2457 "related to triple factorials"   W.D.Lang 
 http://oeis.org/A034255   1,10,120,1560   "related to quartic factorials"  W.D.Lang
 http://oeis.org/A034687   1,15,275,5500   "related to quintic factorials"  W.D.Lang
 http://oeis.org/A034789   1,21,546,15561  "related to sextic factorials"   W.D.Lang       

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}$

 1   .   .   .   .  .
 0   1   .   .   .  .
 0   1   1   .   .  .
 0   2   3   1   .  .  = top-left segment of S1
 0   6  11   6   1  .
 0  24  50  35  10  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,...\}$

Answer 2

Answer 1
Here I show how I tried to get a solution which lead to my guesses for interpretation of $R_p(n)$ in my question.

Using approximate Matrix-logarithm/Matrix-exponential as well as alternatively the Newton-squareroot-algorithm applied to the matrix of Faulhaber-polynomials G I got the following approximate matrix-squareroot H where (if I could use infinite size for matrixes) $H^2 = G$ . I tested increasing matrixsizes $16 \times 16$,$32 \times 32$ and $64 \times 64$ - which seemed to converge to some stable entries.

Top-left segment of $G$ ("Faulhaber-matrix")

  .      1      .    .     .    .
  .    1/2    1/2    .     .    .
  .    1/6    1/2  1/3     .    .
  .      .    1/4  1/2   1/4    .
  .  -1/30      .  1/3   1/2  1/5
  .      .  -1/12    .  5/12  1/2             

In the rows we see the coefficients of the Faulhaber-polynomials for $n$ with exponents beginning at $0$ from left hand. So for instance the polynomial expression for $S_2(n)$ is $f_2(n) = 1/6 n +1/2 n^2 + 1/3 n^3$ . As the $f_p(n)$ are polynomials (with finitely many coefficients) we'll get exact solutions for the $S_p(n)$ sums.
The top-left segment of the matrix $G^2$ was shown in my answer in the linked question, showing the polynomials $ff_p(n)$ for calculation of the "sums-of-sums" $SS_p(n)$
With a vector-definition $V(n) = [1,n,n^2,n^3,n^4,...]$ (of appropriate size) and an according vectordefinition $S(n)$ we can write
$$ G \cdot V(n) = V(1)+V(2)+V(3)+...+V(n) = S(n) $$ and the second iterate
$$ G \cdot S(n) = V(1)+(V(1)+V(2))+(V(1)+V(2)+V(3))+...+S(n) = SS(n) $$ but also simply using the second power of $G$ : $$ G^2 \cdot V(n) = V(1)+(V(1)+V(2))+(V(1)+V(2)+V(3))+...+S(n) = SS(n) $$

---

---

Now what I try is analoguously to define a matrix-squareroot of $G$, say $H$ with $H^2=G$ (in infinite size), such that $$ H^2 \cdot V(n) = S(n)$$ and analoguously to the previous $$ H \cdot V(n) = R(n) \\ H \cdot R(n) = S(n) \\ $$ which would then mean, the matrix $H$ provides us with one form of half-iterate.

The top-left of H after application of the Newton-squareroot algorithm a couple of times looks like

  1/16384     1.99871447249    -2.75904493629     5.11632425994    -9.75077175984     18.2463308729
        .    0.666666666667    0.409137092579   -0.111165751636   0.0554471279037  -0.0324380613176
        .    0.133333333333    0.615160751858    0.305076523661  -0.0778431753098   0.0378700882927
        .  -0.0190476190472    0.216881797347    0.600594596146    0.240852973131  -0.0562906965015
        .  -0.0190476190481  -0.0370864566058    0.292351894073    0.596033988324    0.197374698454
        .  0.00865800865845  -0.0524065889929  -0.0599548589806    0.365583952038    0.594744375947
    ...      ...                     

where the first couple of digits seem to stabilize as well with iterating the Newton-algorithm as well as increasing the size. (Note that the top-left small fraction vanishes with more iterations and larger matrix size)

The rows provide now coefficients for series (instead of polynomials) and also the resulting ranges-of-convergence are/seem to be small (or possibly even zero). Noticing, that the signs seem to strictly alternate when going to higher column-indexes we can apply Euler-summation to arrive at approximate results.

Using $H \cdot V(2) = R(2)$ I got for the $R_p(2)$ (which are in the $p$'th row of $R(2)$) the approximations

R_0(2)=  1.49999999825   ~ 2^0 +1/2*1^0
R_1(2)=  2.49999999970   ~ 2^1 +1/2*1^1
R_2(2)=  4.49999999014   ~ 2^2 +1/2*1^2
R_3(2)=  8.49999979057   ~ 2^3 +1/2*1^3
R_4(2)=  16.4999967231   ~ 2^4 +1/2*1^4
R_5(2)=  32.4999594569   ~ 2^5 +1/2*1^5
...

Re-inserted the guessed exact values in the rightmost column of the above table as $H \cdot R(2) = S(2) $ gives the following approximations

S_0(2)=  1.99999999825  ~ 2^0 + 1^0
S_1(2)=  2.99999999965  ~ 2^1 + 1^1
S_2(2)=  4.99999998926  ~ 2^2 + 1^2
S_3(2)=  8.99999978099  ~ 2^3 + 1^3
S_4(2)=  16.9999966461  ~ 2^4 + 1^4
S_5(2)=  32.9999589653  ~ 2^5 + 1^5

which confirm nicely the hypothese, that indeed $H$ gives Faulhaber-analogue series (when H is thought of infinite size) for the half-iterate of the sum-of-like-powers problem. (The actual computations for larger $n$ need of course higher orders for the Euler-summation, I could reproduce meaningful approximations up to $n=6$ so far... )

---

The sum-expressions of the $R_p(n)$ in my question are then guesses from deconstruction of the empirical numerical values for small $n$ and $p$.