Esta es una pregunta interesante que puede ser abordado de muchas maneras, hay muchas probabilidades de que un buen pedazo de matemáticas va a salir de ella. Por ahora, me limitaré a seguir recopilando y ordenando las observaciones, hasta llegar a una respuesta completa.
Tenemos $H_x=\gamma+\psi(x+1)$$\int_{0}^{1}\psi(x+1)\,dx = \log\frac{\Gamma(2)}{\Gamma(1)}=0$, por lo que la integral es igual a
$\gamma^2+\int_{0}^{1}\psi(x+1)^2\,dx$. La función de $\psi(x+1)^2$ es positivo y convexo en $(0,1)$ y los valores de la $\psi$ función en puntos racionales en $(0,1)$ puede ser calculado en un modo explícito a través de Gauss Digamma Teorema, por lo tanto la evaluación numérica de la integral dada es bastante simple, a través de la regla de Simpson o enfoques similares.
En un buen barrio de el origen tenemos
$$ H_x = \zeta(2)x-\zeta(3)x^2+\zeta(4)x^3-\zeta(5)x^4+\ldots\tag{1} $$
por lo tanto
$$ \int_{0}^{1}H_x^2\,dx = \sum_{m,n\geq 2}\frac{(-1)^{m+n}}{m+n-1}\zeta(m)\zeta(n) = \sum_{j\geq 3}\frac{(-1)^{j+1}}{j}\sum_{k=2}^{j-1}\zeta(k)\,\zeta(j+1-k) \tag{2}$$
donde podemos recordar el teorema de Euler sobre $\sum_{n\geq 1}\frac{H_n}{n^q}$:
$$ \sum_{k=2}^{j-1}\zeta(k)\,\zeta(j+1-k) = (2+j)\,\zeta(j+1)-2\sum_{n\geq 1}\frac{H_n}{n^j}=j\,\zeta(j+1)-2\sum_{n\geq 1}\frac{H_{n-1}}{n^j}. \tag{3}$$
Este enfoque debería permitir convertir el original de la integral en una serie simple, ya que
$$ \sum_{j\geq 3}(-1)^{j+1}\zeta(j+1) \stackrel{\text{Abel reg.}}{=} 1-\zeta(2)+\zeta(3).$$
En particular, el problema se reduce a la aproximación de evaluación de la siguiente serie:
$$ \sum_{n\geq 1}\left[\frac{1-2n}{2n^2}+\log\left(1+\frac{1}{n}\right)\right]H_{n-1} \tag{4}$$
cuyo término general que sin embargo se comporta como $\frac{\log n}{n^3}$, lo bastante rápida convergencia.
Si aplicamos la sumación por partes, obtenemos un término general que es más sencillo, pero con un lento decaimiento hacia cero:
$$ \begin{eqnarray*}(4)&=&\lim_{N\to +\infty}\left[\left(-\gamma+\frac{\pi^2}{12}\right)H_{N-1}-\sum_{n=1}^{N-1}\frac{\frac{1}{2}H_n^{(2)}-H_n+\log(n+1)}{n}\right]\\&=&\frac{1}{2}\zeta(3)+\sum_{n\geq 1}\frac{H_n-\log(n+1)-\gamma}{n}\tag{5} \end{eqnarray*}$$
Ahora podemos emplear la forma asintótica de la serie de armónicos de números en orden a escribir $(5)$ en términos de números de Bernoulli, los valores de la Riemann $\zeta$ función y la serie
$$ \sum_{n\geq 1}\frac{\log(n+1)-\log(n)}{n}\stackrel{SBP}{=}\sum_{n\geq 1}\frac{\log(n+1)}{n(n+1)}=\int_{0}^{1}\frac{(1-x)\log(1-x)}{x\log x}\,dx \approx 1.25775 \tag{6}$$
que puede ser re-escrita en términos de Gregorio coeficientes o simplemente como $\sum_{m\geq 1}\frac{(-1)^{m+1}\zeta(m+1)}{m}$.
(Continúa)
