First I will present the notation for higher integration, which I copied from Danya Rose
\begin{align} J_x^n(f(x)) = \int ...\int f(x) dx^n \end{align}
Here are the higher integrations of the principial branch of the Lambert W function.
\begin{align} J_x^1(W(x)) &=\frac{x}{W(x)}\left( W^2(x) -W(x) + 1\right) \\ J_x^2(W(x)) &=\frac{x^2}{W^2(x)}\left(\frac{1}{2} W^3(x) - \frac{3}{4}W^2(x) +\frac{3}{4}W(x) + \frac{1}{8}\right) \\ J_x^3(W(x)) &=\frac{x^3}{W^3(x)}\left(\frac{1}{6} W^4(x) -\frac{11}{36} W^3(x) + \frac{11}{36}W^2(x) +\frac{19}{216}W(x) + \frac{1}{81}\right) \\ J_x^4(W(x)) &=\frac{x^4}{W^4(x)}\left(\frac{1}{24} W^5(x) - \frac{25}{288} W^4(x) +\frac{25}{288} W^3(x) + \frac{115}{3456}W^2(x) +\frac{175}{20736}W(x) + \frac{1}{1024}\right) \end{align}
I checked the correctness of these expressions on Wolfram|Alpha. Constant terms are not important in this question, so I will not mention them. Then I made a hypothesis \begin{align} J^k_x(W(x)) = \frac{x^k}{W^k(x)} \mathcal{J}_{-k}(W(x)) \label{100} \end{align} where $\mathcal{J}_{-k}(W(x))$ is $\mathcal{J}$-polynomial.
$\mathcal{J}$-polynomial
Definition
I defined the $\mathcal{J}$-polynomial using $j$-coefficients as \begin{align} \mathcal{J}_{-n}(x) = \sum_{k=0}^{n+1}x^k j(n+1 - k, -n ) \end{align} for $n\geq 1$. I will immediately explain why $-n$ is there. The coefficients that define the $\mathcal{J}$-polynomial are in such an ugly form, because my stupidity defined the coefficients in an inelegant way. $\mathcal{J}$-polynomials also occur in expressions for higher derivates of the principal branch of the Lambert W function. So I will give a complete definition of this polynomial. \begin{align} \mathcal{J}_n(x) = \left\{ \begin{array}{ll} (-1)^{n-1}\sum_{k=0}^{n}x^k j( n - k-1, n) &\mbox{if} \ n\geq 1, n \in \mathbb{Z} \\ 1 & \mbox{if} \ n=0 \\ \sum_{k=0}^{|n|+1}x^k j(|n|+1 - k, n) & \mbox{if} \ n \leq -1, n \in \mathbb{Z} \end{array} \right. \end{align}
Coefficients
Here is table of coefficients j(k, n).
(I apologize for the poor quality)
Certain patterns can be noticed from this table
\begin{align}
j(1, n) &= - j(2, n) < 0 & \ \mbox{for} \ n \leq -1\\
j(0, n) &= \frac{1}{|n|!} & \ \mbox{for} \ n \leq -1\\
j(1, n) & = -\frac{H_{|n|}}{|n|!} & \ \mbox{for} \ n \leq -1\\
\end{align}
\begin{align}
j(0, n) &= (n-1)! & \ \mbox{for} \ n \geq 1\\
j(1, n) &= 2(n-1)(n-1)! & \ \mbox{for} \ n \geq 1\\
j(n, n) & =(n-1)^{(n-2)} & \ \mbox{for} \ n \geq 0\\
\mathcal{J}_n(x) &= \frac{(x+1)\mathcal{J}_n^{'}(x) - \mathcal{J}_{n+1}(x)}{1 - n x - 3n} & \ \mbox{for} \ n \geq 2 .
\end{align}
I did not discover the last formula, but found it in this article on page 340. An explicit formula for j-coefficients for $\mathcal{J}$-polynomials of higher derivatives is given in this
paper on page 121.
\begin{align}
j(n-k-1, n) =\sum_{m=0}^k \frac{1}{m!} \binom{2n-1}{k-m} \sum_{q=0}^m \binom{m}{q} (-1)^q (q+n)^{m+n-1}
\end{align}
for $n \geq 1$.
Explicit formula for the $\mathcal{J}$-polynomial
First I used the power expansion of the Lambert W function, which converges for $|x|<\frac{1}{e}$. \begin{align} W_0(z) = \sum_{n=1}^\infty \frac{(-n)^{n-1}}{n!}x^n \label{101} \end{align} then I took the nth integration of this expansion. \begin{align} J_x^n(W_0) = \sum_{k=0}^\infty \frac{(-k)^{k-1}}{k!} \frac{k!x^{k+n}}{(k+n)!} = \sum_{k=0}^\infty \frac{(-k)^{k-1}x^{k+n}}{(k+n)!} \end{align} I used the hypothesis I mentioned here \begin{align} J^k_x(W(x)) = \frac{x^k}{W^k(x)} \mathcal{J}_{-k}(W(x)) \end{align} I substituted into the equation and got formula \begin{align} \frac{x^n}{W^n(x)} \mathcal{J}_{-n}(W(x)) = \sum_{k=0}^\infty \frac{(-k)^{k-1}x^{k+n}}{(k+n)!} \end{align} so \begin{align} \mathcal{J}_{-n}(W(x)) = \frac{W^n(x)}{x^n}\sum_{k=0}^\infty \frac{(-k)^{k-1}x^{k+n}}{(k+n)!} \end{align} This unfortunately only holds for $|x|<\frac{1}{e}$.
Final formula
The partial formula for higher integration of the Lambert W function is
\begin{align} J^n_x(W(x)) = \sum_{k=0}^\infty \frac{(-k)^{k-1}x^{k+n}}{(k+n)!} \end{align} for $|x|<\frac{1}{e}$.
Questions
What is the formula for j-coefficients $j(k, n)$ with $n \leq -1$?
Is there a formula for the nth integration of the Lambert W function that holds to all real numbers?
Thank you for your help.
This is not an answer.
First remark
$$I_n=\int \big[W(x)\big]^n \,dx=(-1)^n \Big[\Gamma (n+1,-W(x))-\Gamma (n+2,-W(x))\Big]$$
Second remark
I shall use $W=W(x)$ to have more room.
$$J_n=\int J_{n-1} \,dx \qquad \text{with} \qquad J_1=\frac x{W}\big[W^2-W+1 \big]$$
Define $$K_n=\frac{W^n \,x^{-n}\, J_n-\,n^{-(n+1)}}{W}=\frac {P_n(W)}{a_n}$$ where tha $a_n$ for the sequence $$\{1,4,216,20736\}$$ and the first $P_n(W)$ are $$\left( \begin{array}{cc} n & P_n(W) \\ 1 & W-1 \\ 2 & 2 W^2-3 W+3 \\ 3 & 36 W^3-66 W^2+66 W+19 \\ 4 & 864 W^4-1800 W^3+1800 W^2+690 W+175 \end{array} \right)$$