Spurred on by this question, I decided to investigate a more generalised form:
\begin{equation}
I_{m,n} = \int_0^{\infty} \ln^m(x)\sin\left(x^n\right)\:dx
\end{equation}
Where $n,m \in \mathbb{N}$
I have formed a solution in terms of the Gamma Function but I'm unsure whether it can be expressed in terms of other Non-Elementary and/or Elementary Functions. Also very interested to see other approaches (Real + Complex Analysis).
To solve, we first observe that:
\begin{equation}
I_{n,k} = \lim_{\phi\rightarrow 0^+} \frac{d^m}{d\phi^m}\int_0^\infty x^\phi \sin\left(x^n\right)\:dx
\end{equation}
Here let:
\begin{equation}
J_{n}(\phi) = \int_0^\infty x^\phi \sin\left(x^n\right)\:dx
\end{equation}
We observe that we first must solve $J_{n,k}(\phi)$. To achieve we employ Feynman's Trick coupled with Laplace Transforms. This is allowable as the integrand conforms with both Fubini's Theorem and the Dominated Convergence Theorem. Here we introduce:
\begin{equation}
H_{n}(t,\phi) = \int_0^\infty x^\phi \sin\left(tx^n\right)\:dx
\end{equation}
Where
\begin{equation}
J_{n}(\phi) = \lim_{t\rightarrow 1^+} H_{n}(t,\phi)
\end{equation}
We now take the Laplace Transform of $H_{n}(t,\phi)$ with respect to $t$:
\begin{align}
\mathscr{L}_t\left[H_{n}(t,\phi)\right] = \int_0^\infty x^\phi \mathscr{L}_t\left[\sin\left(tx^n\right)\right]\:dx = \int_0^\infty x^\phi \frac{x^n}{s^2 + x^{2n}} \:dx = \int_0^\infty \frac{x^{\phi + n}}{s^2 + x^{2n}} \:dx
\end{align}
Thankfully (and as I address here) this integral can be evaluated easily:
\begin{align}
\mathscr{L}_t\left[H_{n}(t,\phi)\right] = \int_0^\infty \frac{x^{\phi + n}}{s^2 + x^{2n}} \:dx = \frac{1}{2n} \cdot \left(s^2\right)^{ \frac{\phi + n + 1}{2n} – 1}\cdot B\left(1 – \frac{\phi + n + 1}{2n}, \frac{\phi + n + 1}{2n} \right)
\end{align}
Using the relationship between the Beta Function and the Gamma Function:
\begin{equation}
\mathscr{L}_t\left[H_{n}(t,\phi)\right] = \frac{1}{2n} s^{ \frac{\phi + n + 1}{n} – 2}\Gamma\left(1 – \frac{\phi + n + 1}{2n} \right)\Gamma\left(\frac{\phi + n + 1}{2n} \right)
\end{equation}
We now resolve $H_{n}(t, \phi)$ by taking the Inverse Laplace Transform:
\begin{align}
H_{n}(t,\phi)&=\mathscr{L}_s^{-1}\left[ \frac{1}{2n} s^{ \frac{\phi + n + 1}{n} – 2}\Gamma\left(1 – \frac{\phi + n + 1}{2n} \right)\Gamma\left(\frac{\phi + n + 1}{2n} \right)\right]\\
& = \frac{1}{2n} \cdot \frac{1}{\Gamma\left(2 – \frac{\phi + n + 1}{n}\right)t^{-\left(\frac{\phi + n + 1}{n} – 2 + 1\right)} } \cdot \Gamma\left(1 – \frac{\phi + n + 1}{2n} \right)\Gamma\left(\frac{\phi + n + 1}{2n} \right)
\end{align}
We can now solve $J_n(\phi)$:
\begin{equation}
J_{n}(\phi) = \lim_{t\rightarrow 1^+} H_{n}(t,\phi) = \frac{\Gamma\left(1 – \frac{\phi + n + 1}{2n} \right)\Gamma\left(\frac{\phi + n + 1}{2n} \right) }{2n\:\Gamma\left(2 – \frac{\phi + n + 1}{n}\right) }
\end{equation}
And finally we have
\begin{equation}
I_{m,n} = \int_0^{\infty} \ln^m(x)\sin\left(x^n\right)\:dx = \lim_{\phi\rightarrow 0^+} \frac{d^m}{d\phi^m} \left[\frac{\Gamma\left(1 – \frac{\phi + n + 1}{2n} \right)\Gamma\left(\frac{\phi + n + 1}{2n} \right) }{2n\:\Gamma\left(2 – \frac{\phi + n + 1}{n}\right) } \right]
\end{equation}
For example, using the example as linked above we have $m = 2$, $n = 2$:
\begin{equation}
I_{2,2} = \int_0^{\infty} \ln^2(x)\sin\left(x^2\right)\:dx = \lim_{\phi\rightarrow 0^+} \frac{d^2}{d\phi^2} \left[\frac{\Gamma\left(1 – \frac{\phi + 2 + 1}{2\cdot 2} \right)\Gamma\left(\frac{\phi + 2 + 1}{2\cdot2} \right) }{2n\:\Gamma\left(2 – \frac{\phi + 2 + 1}{2}\right) } \right]
\end{equation}
I was too lazy to do it by hand, but evaluated through WolframAlpha we observe that:
\begin{equation}
I_{2,2} = \int_0^{\infty} \ln^2(x)\sin\left(x^2\right)\:dx = \frac{1}{32}\sqrt{\frac{\pi}{2}}(2\gamma-\pi+4\ln2)^2
\end{equation}
As required
Best Answer
A method relying on the Mellin transform of the sine:
For $s>1$ and $z \in \mathbb{C}$ with $-s < \operatorname{Re} (z) < s$ we have $$ f_s (z) \equiv \int \limits_0^\infty x^{z-1} \sin(x^s) \, \mathrm{d} x = \frac{1}{s} \int \limits_0^\infty t^{\frac{z}{s}-1} \sin(t) \, \mathrm{d} t = \frac{1}{s} \mathcal{M}(\sin) \left(\frac{z}{s}\right) = \frac{1}{s} \sin \left(\frac{\pi z}{2s}\right) \operatorname{\Gamma} \left(\frac{z}{s}\right) $$ with the limit $f_s(0) = \frac{\pi}{2s}$. For $m \in \mathbb{N_0}$ and $s>1$ this implies \begin{align} I_{m,s} &\equiv \int \limits_0^\infty \ln^m(x) \sin(x^s) \, \mathrm{d} x = f_s^{(m)}(1) = \frac{1}{s} \frac{\mathrm{d}^m}{\mathrm{d} z^m} \left[\sin \left(\frac{\pi z}{2s}\right) \operatorname{\Gamma} \left(\frac{z}{s}\right)\right] \Bigg\vert_{z=1} \\ &= \frac{1}{s^{m+1}} \frac{\mathrm{d}^m}{\mathrm{d} x^m} \left[\sin \left(\frac{\pi }{2}x\right) \operatorname{\Gamma} \left(x\right)\right] \Bigg\vert_{x= 1/s} = \frac{1}{s^{m+1}} \sum \limits_{k=0}^m {m \choose k} \left(\frac{\pi}{2}\right)^k \sin^{(k)} \left(\frac{\pi}{2s}\right) \operatorname{\Gamma}^{(m-k)} \left(\frac{1}{s}\right) \, , \end{align} where the last step follows from the general Leibniz rule. The derivatives of the gamma function can be expressed in terms of polygamma functions using Faà di Bruno's formula, but otherwise that is probably as elementary as it gets.
In the special case $s=2$ we can use $$ \sin^{(k)} \left(\frac{\pi}{4}\right) = \frac{1}{\sqrt{2}} (-1)^{\lfloor k/2 \rfloor} $$ for $k \in \mathbb{N}_0$ and the values (obtained from the Legendre duplication formula) \begin{align} \operatorname{\Gamma} \left(\frac{1}{2}\right) &= \sqrt{\pi} \, , \\ \operatorname{\psi}^{(0)} \left(\frac{1}{2}\right) &= - \gamma - 2 \ln(2) \, , \\ \operatorname{\psi}^{(n)} \left(\frac{1}{2}\right) &= (-1)^{n-1} n! (2^{n+1}-1) \zeta(n+1) \, , \, n \in \mathbb{N} \, , \end{align} to simplify the final result. The complexity of Faà di Bruno's formula prevents us from finding a reasonably nice general expression for $(I_{m,2})_{m \in \mathbb{N}_0}$ , but at least we know that these integrals can be written in terms of $\pi$, $\ln(2)$, $\gamma$ and zeta values.