$n^\text{th}$ integral of $\operatorname{li}(x)$ $\newcommand{\li}{\operatorname{li}}\newcommand{\Ei}{\operatorname{Ei}}$The $n^\text{th}$ order integrals (without constants) of the log integral $\li(x)$ show an obvious pattern and can most definitely be generalized.
\begin{align}
& \li(x) \\[10pt]
& -\Ei(2\ln(x))+x\li(x) \\[10pt]
& \frac{1}{2}(-2x\Ei(2\ln(x))+\Ei(3\ln(x))+x^2\li(x)) \\[10pt]
& \frac{1}{6}(-3x^2\Ei(2\ln(x))+3x\Ei(3\ln(x))-\Ei(4\ln(x))+x^3\li(x)) \\[10pt]
& \frac{1}{24}(-4x^3\Ei(2\ln(x))+6x^2\Ei(3\ln(x))-4x\Ei(4\ln(x))+\Ei(5\ln(x))+x^4\li(x))
\end{align}
What would the general form be for an abritrary $n^\text{th}$ integral?
 A: $\newcommand{\li}{\operatorname{li}}\newcommand{\Ei}{\operatorname{Ei}}$
Well I've finally figured it out.
The $n^\text{th}$ integral of $\li(x)$ is equal to:
$$\frac{1}{n!}\sum^n_{k=0}(-1)^k{{n\choose k}}x^{n-k}\Ei((k+1)\ln(x))$$
$$=\frac{1}{n!}\sum^n_{k=0}(-1)^k{{n\choose k}}x^{n-k}\li(x^{k+1})$$
or
$$\sum^n_{k=0}(-1)^k\frac{x^{n-k}}{(n-k)!k!}\li(x^{k+1})$$
A: Inspired by OP's original answer, I explored why the result is the way it is. First, by reverse-engineering and then by forward derivation.
We want to express this in the form $\frac{1}{n!}(x-T)^n \operatorname{li}(x)$ where $T$ is an operator with the property
$$T\operatorname{li}(x^n)=\operatorname{li}(x^{n+1}).$$
We know that
$$\int_0^x \frac{t^n{\,\rm d}t}{\ln t}\overset{u=t^{n+1}}{=}\int_0^{x^{n+1}}\frac{du}{\ln u}=\operatorname{li}(x^{n+1})$$
so to transform $\operatorname{li}(x)$ into $\operatorname{li}(x^{n+1})$ you must differentiate, multiply by $x^n$, and integrate again. Therefore
$$T^n=\int x^n \partial$$
and because $\partial\int = 1$, we have $T^n T^m = T^{n+m}$ and $T^0=1$ as it should be. Of course, $\int$ is now an operator that works on everything to the right.

Now we go forward cleanly.
Your claim is now
$$\underbrace{\int\int\cdots\int}_n \operatorname{li}(x) = \frac{1}{n!}(x-\int x \partial)^n \operatorname{li}(x)$$
First, observe that the derivative of a product $\partial (x f(x))=f(x)+x \partial f(x)$ can be written in the operator form as
$$\partial x = x \partial + 1$$
or, as a commutator
$$\partial x - x\partial = 1$$
Which is just saying that order of differentiation and multiplication is important - this exact commutator leads to Heisenberg's uncertainty principle.
Now integrate the commutator to get
$$\int = x - \int x \partial$$
This is the essence of integration by parts. At this point it's useful to verify that $x$ and $\int x \partial $ actually commute, so that order of multiplication in binomial formula doesn't matter.
So now we have
$$\left(\int\right)^n = (x-\int x \partial )^n = (x-T)^n$$
The only thing to recall is that $\int$ is now an operator that always integrates all the way to the $x$, no matter how many times you apply it, while the multiple integral is actually $\ldots\int_0^{x}dt\int_0^{t} f(t')dt'$. So, in $(t,t')$ space, the double integral integrates over a triangle, while $\left(\int\right)^2$ integrates over a square. In multiple dimensions, this goes over all permutations of order of integration, so you get $n!$-times the value (see time-ordering). You are left with a completely general formula
$$\int_0^x\ldots\int_0^{t''}\int_0^{t'} f(t)dt\,dt'\,\cdots\,dt'' =\frac{1}{n!}\left(\int\right)^n f=\frac{1}{n!}\left(x-\int x \partial\right)^n f$$
This is true for any $f$ and is a general partial integration formula for multiple integration. In our particular case, we are lucky, that we actually don't have to integrate anything once we see the rule
$$ T\operatorname{li}(x^n)=\operatorname{li}(x^{n+1}).$$
