Okay, the estimate actually holds. I had a "moment" yesterday when I first read gerw's argument, which is actually subtly flawed. Let me quickly illustrate:

To estimate $\| u^2\|_{H^1}$ you need to estimate $\|u^2\|_2$ and $\|\nabla (u^2)\|_2$. The first term is bounded by $\|u\|_4^2$ as claimed, but the second term is not $\|\nabla u\|_4^2$! What you have is in fact $\|u \nabla u\|_2 \leq \|u \|_4 \|\nabla u\|_4^2$. This behaves a bit better than $\|\nabla u\|_4^2$ since the "total number of derivatives is fewer".

By examining the binomial formula, we have that

$$ \| \nabla^k(vw)\|_{2} \lesssim \sum_{j = 0}^k \| \nabla^j v \nabla^{k-j}w\|_2 $$

Only one of the two terms can receive more than half the derivatives. We put that term in $L^2$ and the other in $L^\infty$, and we get very roughly

$$ \|\nabla^k(vw)\|_2 \lesssim \| v\|_{W^{\alpha,\infty}}\| w\|_{H^{k,2}} + \| w\|_{W^{\alpha,\infty}}\| v\|_{H^{k,2}} $$

where $\alpha = \lfloor k/2\rfloor$. Sobolev embedding gives

$$ W^{\alpha + n/2 + 1,2} \to W^{\alpha,\infty} $$

So for $k >n + 2$ we have shown that $W^{k,2}$ is an algebra.

The bound $k > n+2$ is not sharp. If you do better interpolation between the terms you can get it all the way down to almost the boundary for Sobolev embedding to $L^\infty$ to hold.

Basically, you need to apply the product rule to distribute the derivatives on the two terms and estimate each of the two terms in the optimal Sobolev spaces given by the Sobolev inequalities. I leave the details as an exercise for you to check.