On $\int_0^{\pi} \operatorname{Si}^n(x) \ \mathrm{d}x $ The Sine Integral, $\operatorname{Si}(x)$, is defined as $$\operatorname{Si}(x)=\int_0^x \frac{\sin t}{t} \ \mathrm{d}t $$
It follows that $$\int_0^{\pi} \operatorname{Si}(x) \ \mathrm{d}x = \bigg[x \operatorname{Si}(x) \bigg]_0^{\pi}-\int_0^{\pi} \sin x \ \mathrm{d}x = \pi \operatorname{Si} (\pi)-2$$ and further $$\int_0^{\pi} \operatorname{Si}^2 (x) \ \mathrm{d}x = \bigg[ x \operatorname{Si}^2(x) \bigg]_0^{\pi} -2\int_0^{\pi} \operatorname{Si} (x) \sin x \ \mathrm{d}x \\ = \pi \operatorname{Si}^2(\pi) -2 \left( \big[ -\cos x \operatorname{Si} (x) \big]_0^{\pi}+\int_0^{\pi} \frac{\sin x \cos x}{x} \ \mathrm{d}x \right)\\ = \pi \operatorname{Si}^2(\pi)-2 \operatorname{Si} (\pi)-\int_{2x=0}^{2\pi}\frac{\sin 2x}{2x} \ \mathrm{d}(2x)\\ = \pi\operatorname{Si}^2(\pi)-2 \operatorname{Si} (\pi)-\operatorname{Si} (2\pi)$$
Following a similar procedure on $\int_0^{\pi} \operatorname{Si}^3(x) \ \mathrm{d}x$ doesn’t quite work, because it requires evaluating the integral $\int_0^{\pi} \operatorname{Si}^2(x) \sin x \ \mathrm{d}x$, or equivalently, $\int_0^{\pi} \frac{\operatorname{Si}(x) \sin 2x}{x} \ \mathrm{d}x$.
This leads me to ask, can (and how can) $\int_0^{\pi} \operatorname{Si}^3(x) \ \mathrm{d}x$ be expressed in terms of elementary or special functions? How about the general case, i.e.  $$\color{purple}{\int_0^{\pi} \operatorname{Si}^n(x) \ \mathrm{d}x} \; ? $$

Here is an observation that seems to hint a negative answer to OP's question.

Considering that $\pi$ does not play a particular role in the special values for $\operatorname{Si}(\cdot)$, the question would be almost equivalent to investigating the antiderivatives of $\operatorname{Si}(x)^k$'s. Now let us write
$$F(x) = \int_{x}^{\infty} \frac{\sin t}{t} \, \mathrm{d}t = \frac{\pi}{2} - \operatorname{Si}(x)$$
for the complementary version of $\operatorname{Si}(\cdot)$. In light of the binomial theorem, it is clear that $\int\operatorname{Si}(x)^n\,\mathrm{d}x$ can be computed if all of $\int F(x)^k\,\mathrm{d}x$ for $k=0,\dots,n$ can be computed. So we will shift our focus to finding the antiderivatives of $F(x)^k$'s.
Now employing OP's technique, we find that
\begin{align*}
\int F(x) \, \mathrm{d}x
&= x F(x) + \int \sin x \, \mathrm{d}x \\
&= x F(x) - \cos x + C
\end{align*}
and
\begin{align*}
\int F(x)^2 \, \mathrm{d}x
&= x F(x)^2 + 2\int F(x) \sin x \, \mathrm{d}x \\
&= x F(x)^2 - 2F(x)\cos x - \int \frac{2\sin x\cos x}{x} \, \mathrm{d}x \\
&= x F(x)^2 - 2F(x)\cos x + F(2x) + C.
\end{align*}
These computations, together with the asymptotic formula for $F(x)$, reveal the closed-forms for the following integrals:
$$ I_1 = \int_{0}^{\infty} F(x) \, \mathrm{d}x = 1 \qquad\text{and}\qquad I_2 = \int_{0}^{\infty} F(x)^2 \, \mathrm{d}x = \frac{\pi}{2}. $$
Similarly, if $\int F(x)^3 \, \mathrm{d}x$ admits a closed-form expression involving $F(\cdot)$, then we can expect that the integral
$$ I_3 = \int_{0}^{\infty} F(x)^3 \, \mathrm{d}x $$
will admit an elementary closed-form (or at least assume a simpler formula than $\int_{0}^{\pi}\operatorname{Si}(x)^3\,\mathrm{d}x$). However, we can prove that
$$I_3 = \frac{\pi^2}{4} - \frac{3}{2}\log^2 2 - \frac{3}{4}\operatorname{Li}_2\left(\frac{1}{4}\right) = 1.545982100082988\dots $$
holds, where $\operatorname{Li}_2(\cdot)$ is the dilogarithm. Since this expression seems non-elementary, I am skeptical about the ideal that $\operatorname{Si}(x)^3$ has an elementary antiderivative.

Addendum: Proof of the formula for $I_3$. We find that
\begin{align*}
\int F(x)^3 \, \mathrm{d}x
&= x F(x)^3 + 3\int F(x)^2 \sin x \, \mathrm{d}x \\
&= x F(x)^3 - 3F(x)^2\cos x - 6\int \frac{F(x)\sin x\cos x}{x} \, \mathrm{d}x \\
&= x F(x)^3 - 3F(x)^2\cos x - 3\int \frac{F(x)\sin(2x)}{x} \, \mathrm{d}x.
\end{align*}
From this, we get
$$ I_3 = \frac{3\pi^2}{4} - 3\int_{0}^{\infty} \frac{F(x)\sin(2x)}{x} \, \mathrm{d}x. $$
To compute this, we adopt the Feynman's trick. More precisely, we consider the parametrized integral
$$ I_3(s) = \frac{3\pi^2}{4} - 3\int_{0}^{\infty} \frac{F(x)\sin(sx)}{x} \, \mathrm{d}x. $$
Then
\begin{align*}
I_3'(s)
&= -3\int_{0}^{\infty} F(x)\cos(sx) \, \mathrm{d}x \\
&= -\frac{3}{s}\int_{0}^{\infty} \frac{\sin x \sin(sx)}{x} \, \mathrm{d}x \\
&= -\frac{3}{2s}\int_{0}^{\infty} \frac{\cos((s-1)x) - \cos((1+s)x)}{x} \, \mathrm{d}x
\end{align*}
The last integral can be computed using the Frullani's integral, yielding
\begin{align*}
I_3'(s)
&= -3\int_{0}^{\infty} F(x)\cos(sx) \, \mathrm{d}x \\
&= -\frac{3}{s}\int_{0}^{\infty} \frac{\sin x \sin(sx)}{x} \, \mathrm{d}x \\
&= \frac{3}{2s}(\log\left|s-1\right| - \log(s+1)).
\end{align*}
So it follows that
\begin{align*}
I_3 = I_3(2)
&= I_3(0) + \int_{0}^{2} I_3'(s) \, \mathrm{d}s \\
&= \frac{3\pi^2}{4}+ \int_{0}^{2} \frac{3}{2s}(\log\left|s-1\right| - \log(s+1)) \, \mathrm{d}s.
\end{align*}
By computing the last integral using the dilogarithm, we can complete the proof of the formula.
