An informal explanation of why the Apéry's constant appears in the calculation of the integral $\int_0^\infty e^{- \sqrt[3]{x} }\log^3(x)dx$ The second paragraph of section Integrals of the Wikipedia's article dedicated to the Euler–Mascheroni constant shows two definite integrals in which this constant appears, see it here. From those after I did some variations I saw that one can find definite integrals involving $\zeta(3)$. Of course these variations and the following different seem difficult and artificious indefinite integrals. I've curiosity about if you can explain me if one can predict (isn't required a rigurous explanation) when one of such integrals that I evoke has as one of their terms $\zeta(3)$.

Question. Can you explain/predict (isn't required a rigurous explanation) why, for instance, in the definite integral $$\int_0^\infty e^{- \sqrt[3]{x} }\log^3(x)dx$$ appears $\zeta(3)$ as a term of the solution but not in $$\int_0^\infty e^{- \sqrt[7]{x} }\log^2(x)dx?$$ Many thanks

I did the calculations with the help of Wolfram Alpha
int e^{-x^(1/3)}log^3(x)dx, from x=0 to infinite
and the corresponding code in Wolfram Language for the indefinite integral 
int e^{-x^(1/3)}log^3(x)dx

I don't know if this integrals were in the literature, if in your explanation you need literature refers it, and I try find those literature. But as I've said your explanation can be in this questions without rigor. I want to know if one can to predict more or less when $\zeta(3)$ appears as a term in the solution of  previous integrals.
 A: We have
$$ \int_0^{\infty} x^{a-1} e^{-x^{1/p}} \, dx = p\Gamma(ap), $$
where $\Gamma$ is the Gamma-function. Differentiating $k$ times gives
$$ \int_0^{\infty} x^{a-1} e^{-x^{1/p}} (\log{x})^k \, dx = \frac{d^k}{da^k}p\Gamma(ap). \\
\int_0^{\infty} e^{-x^{1/p}} (\log{x})^k \, dx = \left.\frac{d^k}{da^k}p\Gamma(ap) \right|_{a=1} \tag{1} $$
Now, the derivative of the Gamma-function is normally given by $\Gamma'(z) = \Gamma(z)\psi(z)$. Since $\Gamma(z+1)=z\Gamma(z)$, $\psi$ satisfies
$$ \psi(z+1) = \frac{1}{z}+\psi(z), $$
which suggests that the $k$th derivative of $\Gamma$ will involve sums of $1/n^k$. Indeed, since $\Gamma'(1)=-\gamma$, $\psi(1)=-\gamma$, and we therefore find that $\psi(m) = \sum_{r=1}^{m-1} \frac{1}{r} - \gamma$. One can show more generally that for $z \notin \{ 0,-1,-2,\dotsc \}$,
$$ \psi(z) = -\gamma + \sum_{n=0}^{\infty} \left( \frac{1}{n+1} - \frac{1}{n+z} \right). $$
Now, the $k$th derivative of $\Gamma$ will involve the first $k-1$ derivatives of this expression, evaluated at an integer. We can pass the derivative into the sum (because it converges uniformly to an analytic function), and find
$$ \psi^{(m)}(z) = (-1)^{m+1}m! \sum_{n=0}^{\infty} \frac{1}{(n+z)^{m+1}}. $$
So for $z$ an integer, this gives a finite sum subtracted from $\zeta(m+1)$. Therefore in general, the integral in (1) with $a$ a positive integer can be expressed as a finite sum of products of rationals, $\gamma$, and $\zeta(m)$ for $m \in \{ 2,3,\dotsc,k \}$.
