I know I can get the Hermite polynomials in a single variable with: HermiteH[n, x]
Now I need the bivariate Hermite polynomials. I thought about building them with this procedure:
basis = Flatten[TensorProduct[x^Range[0, n], y^Range[0, n]]]
orthbasis = Orthogonalize[basis, Integrate[Exp[-x^2 - y^2]*#1*#2, {x, -Infinity, Infinity}, {y, -Infinity, Infinity}] &]
Is this procedure correct? Or is there a better way with internal Mathematica functions like HermiteH, without involving Gram-Schmidt?