Where $P_a$ is a prime number and $\sigma_k(n)$ is divisor function
We encountered these two identities:
$$\sigma_k(P_a)+\sigma_k(P^2_a)=(P_a^k+1)^2+1\tag1$$
$$\sigma_k(P_a^2)+\sigma_k(P^3_a)=(P_a^k+1)^3-(P_a^k+1)^2+(P_a^k+1)+1\tag2$$
We would like to know if this identity can be generalise? $$\sigma_k(P_a^j)+\sigma_k(P^{j+1}_a)=F(j),\tag3$$ $j\ge1$
Since $$\sigma_k(P_a^j)=\frac{P_a^{k(j+1)}-1}{P_a^{k}-1}$$ You can generalise $$\sigma_k(P_a^j)+\sigma_k(P^{j+1}_a)=\frac{(P_a^k+1)P_a^{k(j+1)}-2}{P_a^{k}-1}$$
Now if you set $b=P_a^{k}+1$ you get $\frac{b(b-1)^{(j+1)}-2}{b-2}$.
LHS -> Expand $b(b-1)^{(j+1)}$ and you'll get binomial coefficients with alternate signs.
What you did will work for $j=1$ and $j=2$ but not for $j>2$ (at least not that way)
RHS -> $(b-2)[b^{(j+1)}\pm b^{(j)} \pm b^{(j-1)}...\pm b]+(b-2)$
Note: My guess was that you were seeking a formula with powers of $P_a^{k}+1$ on the RHS. Don't know if it was your intention.
Edit:
In fact, you could: On the LHS you would have $\frac{b(b-1)^{(j+1)}-2}{(b-2)}$ with the alternate binomial coefficient of $(b-1)^{(j+1)}$, and on the RHS, you would have $1$ plus all terms from $b^{(j+1)}$ down to $b^2$ (if $j$ is odd) or $b$ (if $j$ is even) with coefficients based on the signed binomial coefficient of the LHS plus 2 times the coefficient of the previous term of the RHS.
e.g. with $j=5$
The signed binomial coefficients of the LHS are: $1,-6,15,-20,15,-6,1$
On RHS, coeficient of $b^6$ is $1$, coefficient of $b^5$ is $-6+2\cdot (1)=-4$, coefficient of $b^4$ is $15+2\cdot (-4)=7$, coefficient of $b^3$ is $-20+2\cdot (7)=-6$, coefficient of $b^2$ is $15+2\cdot (-6)=3$ and we stop at $b^2$ since $j$ is odd.
So you have $$\frac{b(b-1)^6-2}{(b-2)}=b^6-4\cdot b^5+7\cdot b^4-6\cdot b^3+3\cdot b^2+1$$ with $b=P_a^{k}+1$