Let $\psi_n(x)=e^{-\frac{x^2}{2}} H_n(x)$ where $H_n(x)$ is the Hermite polynomial. We refer to $\phi_n(x)$ as the Hermite function.
We are interested in finding the following convolution. \begin{align} \exp(-b^2 x^2) \star \psi_n(a x)= \int_{-\infty}^\infty \exp(-b^2 (x-t)^2) \psi_n(a t) dt, \end{align} where $a$ and $b$ are some positive constants.
I tried using the Fourier technique here, but since both the Gaussian function and the Hermite function are of the same form in the Fourier domain, the problem doesn't really change much since we have to take the inverse Fourier of the product which maps back to the same convolution.
The problem does change in the Fourier domain, because the convolution becomes multiplication. Then you can multiply those two Gaussians (one from the Hermite function) to get a single Gaussian and scale $x$ appropriately to remove extra factors from the exponent. It will scale the terms in the polynomial as well, so it will not be Hermite anymore. However, you can still find the result of the inverse Fourier transform with some tricks involving differential operators.
More details below.
Assume the Fourier transform is defined as ${\cal F}\{f(x)\} = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty} f(x)e^{-ikx}dx$. Then
$$ \begin{aligned} \exp(-b^2 x^2) ⋆ \psi_n(ax) &= \sqrt{2\pi}{\cal F}^{-1}\{{\cal F}\{\exp(-b^2 x^2) \}{\cal F}\{\psi_n(ax) \}\}\\ &= \sqrt{2\pi}{\cal F}^{-1}\left\{\frac{1}{|b|\sqrt{2}}\exp\left(-\frac{k^2}{4 b^2} \right) \frac{(-i)^n}{|a|}\psi_n\left(\frac{k}{a}\right)\right\}\\ &= \sqrt{2\pi}{\cal F}^{-1}\left\{\frac{1}{|b|\sqrt{2}}\exp\left(-\frac{k^2}{4 b^2} - \frac{k^2}{2 a^2} \right) \frac{(-i)^n}{|a|}H_n\left(\frac{k}{a}\right)\right\}\\ &= (-i)^n\sqrt{\pi}\left|\frac{c}{ab}\right|{\cal F}^{-1}\left\{H_n\left(\frac{ck}{a}\right)\exp\left(-\frac{k^2}{2} \right)\right\}\left(cx\right),\\ \end{aligned} $$ where $c = \sqrt{\frac{1}{\frac{1}{2 b^2} + \frac{1}{a^2}}}$.
Now the scaled argument in $H_n$ presents some problems because $H_n\left(\frac{ck}{a}\right)\exp\left(-\frac{k^2}{2} \right)$ is not a Hermite function. To find its inverse Fourier transform, we can use the operator identity $$ \left(\frac{d}{dx} + 2 C x \right)^n = e^{-Cx^2}\frac{d^n}{dx^n}e^{Cx^2}$$ and its corollary $$ \left(\frac{d}{dx} + 2 C x \right)^n e^{-Dx^2} = (-\sqrt{D-C})^n H_n(\sqrt{D-C}x) e^{-Dx^2}.$$ After the inverse Fourier transform, $\left(\frac{d}{dk} + 2 C k \right)$ becomes $\left(-ix + 2 C \frac{d}{d(ix)} \right)$, and Gaussian transforms into another Gaussian, so
$${\cal F}^{-1}\left\{H_n(Ck) \exp\left(-\frac{k^2}{2}\right)\right\} = \left(i\sqrt{2C^2-1}\right)^n H_n\left(\frac{C}{\sqrt{2C^2-1}}x\right) \exp\left(-\frac{x^2}{2}\right).$$
Finally,
$$ \begin{aligned} \exp(-b^2 x^2) ⋆ \psi_n(ax) &= (-i)^n\sqrt{\pi}\left|\frac{c}{ab}\right|{\cal F}^{-1}\left\{H_n\left(\frac{ck}{a}\right)\exp\left(-\frac{k^2}{2} \right)\right\}\left(cx\right)\\ &= \sqrt{\pi}\sqrt{2\frac{c^2}{a^2}-1}^n\left|\frac{c}{ab}\right|H_n\left(\frac{c/a}{\sqrt{2\frac{c^2}{a^2}-1}} cx\right)\exp\left(-\frac{c^2x^2}{2} \right). \end{aligned} $$
Disclaimer: I haven't double-checked all this. Use at your own risk.