I've came across an identity once (I don't remember where) concerning convolutions of Gaussian measures. If I'm not mistaken, this identity was \begin{eqnarray} (\mu_{C}*f)(y) = \exp\bigg{[}\frac{1}{2}\bigg{(}\frac{\partial}{\partial x}, C\frac{\partial}{\partial x}\bigg{)}\bigg{]}f(x)\bigg{|}_{x=y} \tag{1}\label{1} \end{eqnarray} where $C=(C_{ij})$ is a $n\times n$ Gaussian covariance, i.e. a positive-definite matrix and $$\bigg{(}\frac{\partial}{\partial x},C\frac{\partial}{\partial x}\bigg{)} := \sum_{i,j}\frac{\partial}{\partial x_{i}}C_{ij}\frac{\partial}{\partial x_{j}}$$ is a differential operator.
I recently discovered that D. Brydges states a very similar result in his notes on functional integrals (page 34, exercise 4.3.3). He mentions that this identity is related to Wick's theorem for Gaussian measures. Note, however, that Brydges states (\ref{1}) in the case $f$ is a polynomial.
I checked some references I know to look for this version of Wick's Theorem but I have found nothing and I'd like to know more about this identity. If $f$ is assumed to be a polynomial, it seems to me that the natural way to prove it is to consider the Fourier transform of $\mu_{C}$ in order to obtain its moments. But does this identity hold in more general cases? If this is the case, how to address it and prove it? Also, any reference is welcome here. Thanks!