This is related to the probability density of the quantum harmonic oscillator. (The Hermite polynomials are the eigenstates of that problem.) For large $k$ the density $$P_k(x)=\frac{1}{\sqrt{2\pi}k!} [h_k(x)]^2 \, e^{-x^2/2},$$ normalized to unity, has the approximation $$p_k(x)=\frac{1}{\pi}(4k-x^2)^{-1/2},\;\;|x|<2\sqrt{k}.$$ The boundaries $\pm 2\sqrt{k}$ are the classical turning points of the harmonic oscillator. In that context $E=4k$ is the energy of the particle and $V(x)=x^2$ is the harmonic potential in which it moves. The large-$k$ approximation is known as the WKB approximation. More generally, for an arbitrary potential $V(x)$ the density varies as $[E-V(x)]^{-1/2}$.
Alternatively, one can use the asymptotic expansion of the Hermite polynomials, $$e^{-x^2/4}h_k(x) \propto \cos \left(x \sqrt{ k}- \pi k/2 \right)\left(4k-x^2\right)^{-1/4}.$$ Squaring and averaging over the rapidly oscillating $\cos^2$ gives $p_k(x)$.
The plot compares $P_k$ and $p_k$ for $k=100$. The function $p_k$ is a smoothed version of $P_k$, suitable for averaging quantities that vary little on the scale of $1/\sqrt{k}$. For example, one can check by explicit computation that the average of $x^n$, $n\in\mathbb{N}$, agrees to leading order in $k$, $$\lim_{k\rightarrow\infty} k^{-n}\int P_k(x)x^{2n}\,dx=\frac{2^{2n}\,\Gamma(n+1/2)}{\sqrt{\pi}\,\Gamma(n+1)}=k^{-n}\int p_k(x) x^{2n}\,dx.$$