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

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?

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$.
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.