Let $h_0, h_1, \dots$ be the classical univariate Hermite polynomials, renormalized to have constant norm. Is $$x\mapsto\prod_{j=1}^n h_{l_j}(x_j), \quad l_j\in \mathbb N$$ a complete orthogonal system for $L^2(\mathbb R^n)$, as in the univariate case? The orthogonality is simple, but I still have some doubts on the completeness.

In one dimension, you have $$ h_k(t)e^{-π t^2}=e^{π t^2}(\frac{d}{dt})^k\bigl(e^{-2π t^2}\bigr), $$ and $h_k$ is easily proven to be with degree $k$. The completeness question amounts to proving $L^2=\overline{\text{span}\{t^ke^{-π t^2}\}}=E$. Taking $u\in E^\perp$, you obtain $$ \int u(t)\sum_k\frac{t^ki^k\xi ^k}{k!}e^{-π t^2} dt=0, $$ and by injectivity of the Fourier transform $u(t)e^{-π t^2}=0$, so that $u=0$ (note that this argument works for $u$ tempered distribution). In several variables, you get by tensorization $$ E=\overline{\text{span}\{x^\alpha e^{-π \vert x\vert^2}\}}, $$ and the same Fourier argument works.

In addition to @Bazin's good answer, it is also possible to understand this as an example of a more general phenomenon. Namely, at least up to normalizing constants, $h_k(x)\cdot e^{-x^2/2}$ are eigenfunctions for a Schrodinger operator $-\Delta+|x|^2$. The key point is that (the Friedrichs self-adjoint extension of) this operator provably has compact resolvent, so gives a Hilbert-space basis of eigenfunctions (that is, no continuous spectrum). In one dimension, the Schrodinger operator factors as $(i{d\over dx}-ix)(i{d\over dx}+ix)+1$, and the "lowest" eigenfunction is obtained by solving $(i{d\over dx}+ix)u=0$. The raising and lowering operators $i{d\over dx}\pm ix$ map eigenfunctions to other eigenfunctions, and one shows that all eigenfunctions are obtained by repeated applications of raising operators $i{\partial\over \partial x_j}+ix_j$ to the lowest eigenfunction $e^{-|x|^2/2}$.

That is, the completeness arises from compactness of the resolvent (as in Sturm-Liouville problems) and then explicit identification of all eigenfunctions.

The answer of Bazin is, of course, impeccable but it is perhaps worth noting that the result you are after is a special case of an abstract one on tensor products of bases in Banach spaces (the pioneering article was by Gelbaum and de la Madrid---Pac. J. Math 11 (1961) 1281-1286--- but there have been many further ones in the meantime). This works for many different tensor products---in your case you would be using the version for the Hilbert space product of two such spaces.