The product of two Gaussian functions
$$ f(x)=\frac{1}{\sqrt{2 \pi} \sigma_{f}} \exp\left(-\frac{x^{2}}{2 \sigma_{f}^{2}} \right) \quad \text { and } \quad g(x-y)=\frac{1}{\sqrt{2 \pi} \sigma_{g}} \exp\left(-\frac{\left(x-y\right)^{2}}{2 \sigma_{g}^{2}}\right) $$
is another Gaussian function with standard deviation and mean
$$ \sigma_{f g}=\sqrt{\frac{\sigma_{f}^{2} \cdot \sigma_{g}^{2}}{\sigma_{f}^{2}+\sigma_{g}^{2}}} \qquad\qquad \mu_{f g}=\frac{y\cdot\sigma_{f}^{2}}{\sigma_{f}^{2}+\sigma_{g}^{2}} $$
rescaled by
$$ S(y) = \frac{1}{\sqrt{2 \pi\left(\sigma_{f}^{2}+\sigma_{g}^{2}\right)}} \exp \left[-\frac{y^{2}}{2\left(\sigma_{f}^{2}+\sigma_{g}^{2}\right)}\right] $$
Thus, the product is
$$ f(x) \cdot g(x-y)=\frac{S(y)}{\sqrt{2 \pi} \sigma_{f g}} \exp \left[-\frac{\left(x-\mu_{f g}\right)^{2}}{2 \sigma_{f g}^{2}}\right] $$ Namely, a rescaled Gaussian function. The scaling factor is equivalent to $$ S(y) = \int_{\mathbb{R}} f(x) \cdot g(x-y) dx = [f * g](y) $$
In other words, the scaling factor of the product of two Gaussian functions is exactly the result of the convolution of the two Gaussians. It doesn't seem to be a coincidence. Is there any deep reason why this is happening?
The convolution of the two Gaussians is $$\tag{1} \frac{1}{\sqrt{2 \pi\left(\sigma_{f}^{2}+\sigma_{g}^{2}\right)}} \exp \left[-\frac{\left(\color{red}{x\,-\,}\mu_{f}-\mu_{g}\right)^{2}}{2\left(\sigma_{f}^{2}+\sigma_{g}^{2}\right)}\right]. $$ You are saying that this convolution is equal to the product of the two PDFs when we set $\mu_f=\mu_g=0\,.$ Let's check \begin{align} f(x)g(x)&=\frac{1}{\sqrt{2 \pi\left(\sigma_{f}^{2}+\sigma_{g}^{2}\right)}} \exp \left[-\frac{x^2}{2\sigma_f^2}-\frac{x^2}{2\sigma_g^2} \right]\\ &=\frac{1}{\sqrt{2 \pi\left(\sigma_{f}^{2}+\sigma_{g}^{2}\right)}} \exp \left[-\frac{\sigma_g^2x^2+\sigma_f^2x^2}{2\,\sigma_f^2\,\sigma_g^2}\right]\\ &= \underbrace{\frac{1}{\sqrt{2 \pi\left(\sigma_{f}^{2}+\sigma_{g}^{2}\right)}}}_{\text{ very good }} \exp \underbrace{\left[-\frac{(\sigma_g^2+\sigma_f^2)\,x^2}{2\,\sigma_f^2\,\sigma_g^2}\right]}_{(???)}\,.\tag{2} \end{align} To me this looks different from the convolution (1) even when we assume $\mu_f=\mu_g=0\,.$
Without further background I don't think that the product is very interesting. What is more interesting are
the product $f(x)\,g(\color{red}{y})$ which is -as we know- the PDF of the joint distribution of two independent Gaussians;
the integral $\int_{\mathbb R} \frac{1}{\sqrt{2\pi t}}\exp(-\frac{(x-z)^2}{2t})\frac{1}{\sqrt{2\pi s}}\exp(-\frac{(z-y)^2}{2s})\,dz=\frac{1}{\sqrt{2\pi}(t+s)}\exp(-\frac{(x-y)^2}{2(t+s)})$ which shows that the Brownian heat kernels are a convolution semigroup.