Hankel determinants of harmonic numbers

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

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.

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.