How naturally can functions defined by parametric integrals be interpolated from $\mathbb N$ to $\mathbb R^+$?

Question

It is well known that several definite integrals $I_n$ containing a parameter $n\in\mathbb N$ can be expressed recursively (e.g. doing integration by parts) in terms of $I_{n-1} $ or $I_{n-2} $, and thus written as some expression containing factorials, e.g. the gamma function itself $$\int_0^\infty t^{n}e^{-t}dt=n\int_0^\infty t^{n-1}e^{-t}dt=\cdots=n!$$

Or take the formula $I_n:=\int_0^\pi \sin^nx\;dx=\dfrac{n-1}nI_{n-2}$, allowing to obtain $$\int_0^\pi \sin^{2n}x\;dx=\dfrac{(2n)!}{2^{2n}n!^2}\pi\ \ \text{ and }\ \int_0^\pi \sin^{2n+1}x\;dx=\dfrac{2^{2n+1}n!^2}{(2n+1)!}$$ (which BTW easily yields the Wallis product).

As long as for such integrals the LHS is also defined for non-integer (say all positive real) $n$,

is it "automatically" guaranteed that replacing on the RHS $n!$ by $\Gamma(n+1)$ yields a valid formula?

We all know that the interpolation of the gamma function between the factorials to $\mathbb R$ is not unique but that the Bohr-Mollerup theorem assures a unique extension when adding the mild (?) condition of log-convexity. The problem: the LHS expressions above don't "know" anything about log-convexity...

We can ask a similar question (even though there is no simple recursion formula in this case) about the Riemann zeta function with, for $n>1$, $$\int_{0}^{\infty} \frac{t^n}{e^t - 1} \; \frac{dt}{t}=\zeta(n) \; \Gamma(n) $$ or, closely related, this interpolation of the Bernoulli numbers $$ 4n\int_{0}^{\infty} \frac{t^{2n}}{e^{2\pi t}-1} \frac{dt}{t}=4n\frac{2^{2n-1}}{2^{2n-1}-1}\int_{0}^{\infty} \frac{t^{2n}}{e^{2\pi t}+1} \frac{dt}{t}=(-1)^{n+1}B_{2n}$$ or this one of the Euler numbers $$ \int_{0}^{\infty} \frac{t^{2n}}{\cosh\frac{\pi t}2}\; dt =(-1)^{n}E_{2n}.$$ The two last ones are somewhat intriguing because for half-integers $n$, the integrals are obviously not $0$, unlike the odd Bernoulli and Euler numbers. For Euler numbers, this can of course be explained by the fact that the "entire" integral from $-\infty$ to $\infty$ does vanish in this case (odd function), but for the Bernoulli numbers the situation is different.

What is going on here?

Answer 1

I don't have a full answer to this question, but I do have a rather good semi-answer. Considering Carlson's theorem there is something that can be said about this.

Now Carlson's theorem is a rather niche result that doesn't pop up as much as I think it should but it has some rather unwieldy results. To be short it gives an identity theorem for natural numbers. The result can be stated a few ways, so I'll choose a modest version of it which is a consequence of Ramanujan's Master theorem (Carlson's theorem is a stronger version of this).

Let $z = x+iy$ for $\Re(z) > 0$

Let $\phi$ be holomorphic. If $|\phi(z)| < C e^{\kappa|y| + \rho|x|}$ for $C,\kappa, \rho \in \mathbb{R}^+$ where $\kappa < \pi$. If $\phi\Big{|}_{\mathbb{N}} = 0$ then $\phi = 0$.

This helps us in your current situation through the following

If $\phi(z) = \int_0^\pi \sin^{2z}(x)\,dx - \pi\frac{\Gamma(2z+1)}{2^{2z}\Gamma(z+1)^2}$ is bounded as above, then since it equals zero on the naturals, they must equal. Successfully interpolating your values. Now this isn't the case here, but very often it does pop up; and it quickly solves a result like this. For your Bernoulli number identity it should work through some manoeuvring. Similarly with the $\Gamma\zeta$ case, however one should divide by the Gamma function first and then apply the identity. It also works with the Euler numbers, in this case.