Hankel determinants of harmonic numbers Let $H_n=\sum_{k=1}^n\frac 1 k$ be the $n$-th harmonic number with $H_0=0.$
Question: Is the following true?
$$\det\left(H_{i+j}\right)_{i,j=0}^n=(-1)^n \frac{2H_{n}}{n! \prod_{j=1}^n \binom{2j}{j} \binom{2j-1}{j}}.$$
Edit:
Comparing with the orthogonal polynomials whose moments are the numbers $\frac{1}{n+1}$ it suffices to show the following identity:
$$\sum_{j=0}^n (-1)^j\frac{\binom n j \binom{n+j} j}{\binom{2n} n} H_j \prod_{j=0}^{n-1}\frac{(j!)^3}{(n+j)!} = (-1)^n \frac{2H_n}{n! \prod_{j=1}^n \binom{2j}{j} \binom{2j-1}{j}}.$$
 A: I prove your identity $$\sum_{j=0}^n (-1)^j\binom{n}{j}\binom{n+j}{j} H_j= 2(-1)^n H_{n}$$
which you claim to imply the result.
The method is the same as here. 
At first, use $(-1)^k\binom{n+k}k=\binom{-n-1}k$. Then $$F(y):=\sum_k (-1)^k\binom{n}k\binom{n+k}ky^k=[x^n] (1+x)^n(1+xy)^{-n-1}.$$
Next, for any polynomial $F(y)=\sum c_ky^k$ we have 
$$
\sum c_kH_k=\int_0^1 \frac{F(y)-F(1)}{y-1}dy.
$$
Integration over $[0,1]$ and taking the coefficient of $x^n$ commute, thus we have to prove
$$
[x^n]\int_0^1\frac{(\frac{x+1}{1+xy})^n\cdot \frac1{1+xy}-\frac1{1+x}}{y-1}dy=2(-1)^nH_n.
$$
A natural change of variables here is $t=(1+x)/(1+xy)$, we get that our integral equals $$-\frac{1}{1+x}\int_1^{1+x}\frac{1-t^{n+1}}{t(1-t)}dt=\frac{-\log(1+x)+H_n}{1+x}-\sum_{i=1}^n\frac{(1+x)^{i-1}}i.$$
A coefficient of $x^n$ indeed equals $2(-1)^nH_n$.
A: As asked by Fedor Petrov I sketch the missing details.
If $a(n)$ is any sequence  with $a(0)=1$, such that all Hankel determinants $M_n=\det\left(a(i+j)\right)_{i,j=0}^n$ are $ \neq 0$, define a linear functional $L$ on the polynomials by $L(x^n)=a(n).$  Let $p_n(x)$ be the uniquely determined monic polynomials which are orthogonal with respect to $L.$  These polynomials are given by $$M_{n-1}p_n(x)= \det\left(r(i,j,x)\right)_{i,j=0}^n$$ with $r(i,j,x)=a(i+j)$ for $j<n$ and $r(i,n,x)=x^i.$
For $a(n)=\frac{1}{n+1}$ the corresponding polynomials are $p_n(x)=\sum_{j=0}^n (-1)^j\frac{\binom{n}{j}\binom{n+j}{j}}{\binom{2n}{n}}x^j.$ In this case we get $M_{n-1}=\prod_{j=0}^{n-1}\frac{(j!)^3}{(n+j)!}$ (This seems to be well known, cf. e.g. this preprint (4.2) for $a=b=q=1.$)
Now $\det\left(H_{i+j}\right)_{i,j=0}^n$  can be reduced by column operations to $\det\left(v(i,j)\right)_{i,j=0}^n$, where $v(i,0)=H_{i}$ and  $v(i,j)=\frac{1}{i+j}$ for $j>0$. This is the same as replacing $x^i$ in $r(i,n,x)$ by $H_{i}.$ Therefore we get the above identity.
A: We propose a proof (somewhat different from Fedor's) for the crucial relation
$$\sum_{j=0}^n (-1)^j\binom{n}{j}\binom{n+j}{j} H_j= 2(-1)^n H_{n}.\tag1$$
To this end, define the polynomials
$$P_n(x):=\sum_{j=0}^n (-1)^j\binom{n}{j}\binom{n+j}{j}\binom{x+j}j.$$
Zeilberger's algorithm returns the recurrence
$$(n+2)^2P_{n+2}(x)+(2n+3)(2x+1)P_{n+1}(x)-(n+1)^2P_n(x)=0.\tag2$$
Using the fact that $[x]\binom{x+j}j=H_j,\, [x]x\binom{x+j}j=1,\, P_{n+1}(0)=(-1)^{n+1}$ (see Remark below), induction on equation (1) and applied to (2) leads to:
$$(n+2)^2[x]P_{n+2}(x)+(2n+3)[2(-1)^{n+1}+2(-1)^{n+1}H_{n+1}]-(n+1)^22(-1)^nH_n=0.$$ 
A direct simplification shows
$$(n+2)^2[x]P_{n+2}(x)
=2(n+2)^2(-1)^{n+2}H_{n+2},$$
which completes the induction process and the proof.
Remark. The identity $(-1)^nP_n(0)=\sum_{j=0}^n (-1)^{n-j}\binom{n}{j}\binom{n+j}{j}=1$ is easily provable by the Wilf-Zeilberger methodology. See my answer here as a further illustration.
A: Identities involving harmonic numbers that are of interest for physicists, Utilitas Mathematica 83 (2010), 291-299, H. Prodinger.
This paper contains the identity (1) as well.
Now starts Johann Cigler's big birthday (in 5 minutes). Hereby, I will send my best regards on the occasion.
