Why does $\lim_{x \rightarrow 0} B(x,y)$ exist and how is it calculated? In evaluating integrals like (link to another example)
$$I=\int_0^1\frac{\log(x) \log^2(1-x)dx}{x}$$
one can make the substitution $x=\sin^2(\theta)$ to obtain
$$I=16\int_0^\frac{\pi}{2}\frac{\log(\sin(x)) \log^2(\cos(x)) \cos(x)dx}{\sin(x)}$$
which is the partial derivative of the beta function at $(x=0,y=1)$:
$$I=\lim_{x \rightarrow 0^+}\partial_y^2\partial_x B(x,1)$$
because
$$B(x,y)=2 \int_0^\frac{\pi}{2} \cos^{2x-1}(x)\sin^{2y-1}(x)dx$$
$$\partial^2_y\partial_xB(x,y)=16\int_0^\frac{\pi}{2} \log(\cos(x))^2\log(\sin(x)) \cos^{2x-1}(x)\sin^{2y-1}(x)dx.$$
My problem: $I$ is convergent, and taking the partial derivatives of $B(x,y)$ in the way above does yeild $2I$, however $B(x,y)$ is only defined for $\Re(x), \Re(y) >0$. So, how does one:


*

*Go about proving that the limit exists (directly using the partial derivatives of the beta function, without reference to $I$)? (Extra points for intuition, as I don't see how the partial derivative can exist at a place where the function doesn't).

*Find the limit in this case?

 A: [Initial apology: this is just a bunch of considerations too long to fit in a comment] Now we have a proof.
$$ I = -\int_{0}^{+\infty} x\cdot\log^2(1-e^{-x})\,dx = \int_{0}^{+\infty}\frac{x^2}{e^{x}-1}\log(1-e^{-x})\,dx,$$
$$ I = -\int_{0}^{+\infty}\frac{x^2}{e^{x}-1}\sum_{k=1}^{+\infty}\frac{e^{-kx}}{k}dx = -\int_{0}^{+\infty}x^2 e^{-x}\sum_{k=1}^{+\infty}H_k e^{-kx}dx,$$
$$ I = -\sum_{k=1}^{+\infty}H_k \int_{0}^{+\infty}x^2 e^{-(k+1)x}dx = -\sum_{k=1}^{+\infty}\frac{2H_k}{(k+1)^3}.$$
Interestingly, integrating by parts in another way we get:
$$ I = \int_{0}^{1}\operatorname{Li}_2(x)\left(\frac{\log(1-x)}{x}-\frac{\log x}{1-x}\right)dx = -\frac{1}{2}\operatorname{Li}_2^2(1)-\int_{0}^{1}\frac{\operatorname{Li}_2(x)\log x}{1-x}\,dx,$$
[This part is merely optional]
so the identity $I=-\frac{\zeta(4)}{2}$ can probably be proved through dilogarithm identities. By using the Taylor series of $\operatorname{Li}_2(x)$ and $\frac{\log(1-x)}{x}$ and the well-known identity $\int_{0}^{1}(1-x)^m x^n dx=\frac{m! n!}{(m+n+1)!}$ it is quite easy to notice that:
$$\int_{0}^{1}\frac{\operatorname{Li}_2(1-x)\log(1-x)}{x}\,dx = -\sum_{m,n=1}^{+\infty}\frac{m! n!}{m^2 n^2 (m+n)!},$$
and I would not be surprised if someone manage to put the RHS in $c\cdot\zeta(2)^2$ form though "reverse" creative telescoping or simple fractions decomposition. 
[End of the optional part]
For istance:
$$\int_{0}^{1}\frac{\operatorname{Li}_2(x)\log x}{1-x}dx = -\sum_{n=1}^{+\infty}\frac{1}{(n+1)^2}\sum_{m=1}^{n}\frac{1}{m^2}.\tag{1}$$
Finally, here is the trick:
$$\frac{\pi^4}{120}=\frac{\zeta(2)^2-\zeta(4)}{2} = \sum_{m>n}\frac{1}{m^2 n^2}$$
is exactly the opposite of the RHS in $(1)$. By collecting pieces we get:
$$ I = -\frac{\zeta(2)^2}{2}+\frac{\zeta(2)^2-\zeta(4)}{2} = -\frac{\zeta(4)}{2}$$
QED.
