I found the following identity on Wikipedia, and I am having a difficult time proving it.
For $m,n\in\Bbb N$,
$$I(m,n):=\int_0^1B_n(x)B_m(x)\mathrm{d}x=(-1)^{n-1}\frac{m!n!}{(m+n)!}b_{n+m}$$
Where $B_n(x)$ is the $n$-th Bernoulli polynomial, and $b_n=B_n(0)$ is the $n$-th Bernoulli number. Notably, when one uses $x!=\Gamma(x+1)$ on the RHS and then uses the Beta function, we arrive at
$$I(m,n)=(-1)^{n-1}b_{n+m}\int_0^1t^n(1-t)^m\mathrm{d}t$$
Here's what I've tried.
Because $B_n(x)$ satisfies $$B_n'(x)=nB_{n-1}(x)$$ We can integrate by parts with $\mathrm{d}v=B_n(x)\mathrm{d}x$: $$I(m,n)=\frac1{n+1}B_{n+1}(x)B_m(x)\bigg|_0^1-\frac m{n+1}\int_0^1B_{n+1}(x)B_{m-1}(x)\mathrm{d}x$$ $$I(m,n)=\frac1{n+1}\bigg(B_{n+1}(1)B_m(1)-b_{n+1}b_m\bigg)-\frac{m}{n+1}I(m-1,n+1)$$ Which is a start, but I don't know how to proceed. I would think that the integral should be easy given the fact that Bernoulli polynomials have so many identities, yet I can't seem to find the right one to use. Could you either point me in the right direction, or show me a complete proof? Thanks.