The Apéry's constant and the Airy function

You may refer to the entry 9.10.17 of DLMF, which describes the Mellin transform of the Airy function $\operatorname{Ai}$:

$$ \int_{0}^{\infty} t^{\alpha-1}\operatorname{Ai}(t)\,\mathrm{d}t = \frac{\Gamma(\alpha)}{3^{(\alpha+2)/3}\Gamma\big(\frac{\alpha+2}{3}\big)}. \tag{9.10.17}$$

Differentiating both sides with respect to $\alpha$ 3-times and plugging $\alpha = 1$ gives the answer to your integral in terms of polygamma functions with further simplifications available.

---

Sketch of proof of $\text{(9.10.17)}$. We begin with the integral representation

$$ \operatorname{Ai}(x) = \frac{1}{\pi} \int_{0}^{\infty} \cos\left(\frac{t^3}{3} + xt \right) \, \mathrm{d}t. $$

Taking Mellin transform and switching the order of integration,

\begin{align*} \int_{0}^{\infty} x^{\alpha-1}\operatorname{Ai}(x)\,\mathrm{d}x &= \frac{1}{\pi} \operatorname{Re}\bigg[ \int_{0}^{\infty}\int_{0}^{\infty} x^{\alpha-1} e^{\frac{1}{3}it^3 + ixt} \, \mathrm{d}x\mathrm{d}t \bigg] \\ &= \frac{1}{\pi} \operatorname{Re}\bigg[ \int_{0}^{\infty} \frac{\Gamma(\alpha)}{(-it)^{\alpha}} e^{\frac{1}{3}it^3} \, \mathrm{d}t \bigg] \\ &= \frac{1}{\pi} \operatorname{Re}\bigg[ \Gamma(\alpha) i^{\alpha} \frac{\Gamma\big(\frac{1-\alpha}{3}\big)}{3\cdot(-i/3)^{\frac{1-\alpha}{3}}} \bigg] \\ &= \frac{1}{\pi} \operatorname{Re}\bigg[ \Gamma(\alpha) i^{\frac{2\alpha+1}{3}} \frac{\Gamma\big(\frac{1-\alpha}{3}\big)}{3^{\frac{\alpha+2}{3}}} \bigg] \\ &= \frac{\Gamma(\alpha)\Gamma\big(\frac{1-\alpha}{3}\big)}{\pi 3^{\frac{\alpha+2}{3}}}\cos\big( \tfrac{1}{6}\pi + \tfrac{1}{3}\alpha\pi \big). \end{align*}

Now applying the Euler's reflection formula $\Gamma(s)\Gamma(1-s) = \pi \csc(\pi s)$ with $s = \frac{1-\alpha}{3}$ yields the right-hand side of $\text{(9.10.17)}$.

A little note. Another way to derive the first formula from Sangchul Lee's answer.

We use the fact that:

$$\mathrm{Ai}''(x)=x \mathrm{Ai}(x)$$

Let's consider the following integral:

$$I_n(y)=\int_0^\infty \mathrm{Ai}(y x) x^n dx$$

$$I'_n(y)=\int_0^\infty \mathrm{Ai}'(y x) x^{n+1} dx$$

$$I''_n(y)=\int_0^\infty \mathrm{Ai}(y x) x^{n+3} dx=I_{n+3}(y)$$

Now we can rewrite it:

$$I_n(y)= \frac{1}{y^{n+1}}\int_0^\infty \mathrm{Ai}(x) x^n dx=\frac{I_n(1)}{y^{n+1}}$$

$$I''_n(y)= \frac{(n+1)(n+2) I_n(1)}{y^{n+3}}=I_{n+3}(y)$$

So we have obtained an explicit recurrence relation for $y=1$:

$$I_{n+3}=(n+1)(n+2) I_n \tag{1}$$

We know that:

$$I_0=\frac{1}{3}$$ (a well known integral)

$$I_1=-\mathrm{Ai}'(0)=\frac{1}{3^{1/3} \Gamma(1/3)}$$

$$I_2=\int_0^\infty \mathrm{Ai}(x) x^2 dx= \int_0^\infty \mathrm{Ai}''(x) x dx=\mathrm{Ai}'(x) x \big|_0^\infty-\int_0^\infty \mathrm{Ai}'(x) dx=\mathrm{Ai}(0)=\frac{1}{3^{2/3} \Gamma(2/3)}$$

So, we now know the general recurrence relation for $I_n$ and three initial conditions needed to find any term.

After we derive the general expression for $I_n$, it's easy enough to use continuity of all the functions involved to generalize for non-integer parameters and obtain the desired Mellin transform.

Using the various provided answers, I present to you the craziest integral I've ever seen. Behold:

$$\int_0^\infty \mathrm{Ai}(x)\ln^3x\ dx=-\frac1{81}\left[(4\pi^2+\beta^2)\beta+52\zeta(3)\right]$$ where $\beta=2\gamma+\ln3$.