Let us define the following recursive function involving the sum of divisors function $\sigma(n)$:
\begin{array}{ l } r(n,1)=\sigma(n) \\ r(n,2)=\sum_{d|n}r(d,1) \\ r(n,3)=\sum_{d|n}r(d,2) \\ \ldots \\ r(n,k)=\sum_{d|n}r(d,k-1) \\ \end{array}
If n is squarefree, then the following identiy seems to hold ($p$ are the prime factors of $n$ smaller than $n$):
$$r(n,k)=\prod_{p\mid n}(p+k)$$
An implementation in Mathematica verifies this empirically:
rn = Compile[{{n, _Integer}, {k, _Integer}},
Total[Nest[Catenate@*Divisors, {n}, k]]];
p[n_, k_] := Product[p + k, {p, Complement[Divisors[n], {1, n}]}];
Print[rn[10, 5]];
Print[p[10, 5]];
(* Output: 70 *)
My question: Does there exist a theorem for this identity or (if not) how can we proof this identity?
$n\mapsto r(n,k)$ is the Dirichlet convolution of $k$ copies of $\sigma$, which means it is multiplicative. Therefore, for all $n$ square-free, $$ r(n,k) = \prod_{p\mid n}r(p,k). $$ Now, we see that for all $k\ge 2$, we have $r(1,k)=\sum_{d\mid 1}r(1,k-1)=r(1,k-1)$ and it follows that $r(1,k)=1$. Next, for all $p$ prime and all $k\ge 2$ $$ r(p,k) = \sum_{d\mid p}r(d,k-1) = r(p,k-1) + r(1,k-1) = r(p,k-1) + 1, $$ from which we conclude that $r(p,k)=r(p,1)+k-1=p+k$. Thus, for all square-free $n$, $$ r(n,k) = \prod_{p\mid n}(p+k). $$