Does $\operatorname{div}\left(\nabla G +xG\right)=0\Longleftrightarrow \nabla G +xG=0$? Let $G$ be a smooth function defined on $\textbf{R}^d$. What are the assumptions I should use to assume that
$$\operatorname{div}\left(\nabla G(x) +xG(x)\right)=0 \quad (\forall x\in \textbf{R}^d)$$
implies that $G$ is a Gaussian? (Several answers are possible I guess...) 
Many thanks !
 A: I'm going to use an argument along the lines of what @hOff proposed.  In addition to the smoothness of $G$ I am going to suppose that it decays rapidly at infinity so that all of the following calculations are justified.  For example, we could assume that $G$ is Schwartz class.
As the first order of business we expand your equation as 
$$
\Delta G(x) + x \cdot \nabla G(x) + d G(x) =0.
$$
To attack this problem we're going to use the Fourier transform, which we define as
$$
  \hat{f}(\xi) = \int_{\mathbb{R}^d} f(x) e^{-2\pi i x\cdot \xi} dx.
$$
We have the following useful identities:
$$
\widehat{\Delta f}(\xi) = -4 \pi^2 |\xi|^2 \hat{f}(\xi) \\
\widehat{x \cdot \nabla f}(\xi) = -\xi \cdot \nabla \hat{f}(x) - d \hat{f}(\xi).
$$
These both follow by integrating by parts in the formula for the Fourier transform.  For example, 
$$
\widehat{x \cdot \nabla f}(\xi) = \sum_{j=1}^d \int_{\mathbb{R}^d} x_j \partial_j f(x) e^{-2\pi i x\cdot \xi} dx \\
=  -\sum_{j=1}^d \int_{\mathbb{R}^d}  f(x) \partial_j( x_j e^{-2\pi i x\cdot \xi})  dx  \\
= -\sum_{j=1}^d \int_{\mathbb{R}^d} f(x) [x_j(-2\pi i \xi_j) + 1] e^{-2\pi i x\cdot \xi} dx  \\
= - d \hat{f}(\xi) -  \sum_{j=1}^d \xi_j \int_{\mathbb{R}^d} f(x)(-2\pi i x_j)e^{-2\pi i x\cdot \xi} dx \\
=  - d \hat{f}(\xi) -  \sum_{j=1}^d \xi_j \partial_{\xi_j} \int_{\mathbb{R}^d} f(x)e^{-2\pi i x\cdot \xi} dx \\
= - d \hat{f}(\xi)  - \xi \cdot \nabla \hat{f}(\xi).
$$
Now we apply the Fourier transform to our equation for $G$ to find that 
$$
0= - 4\pi^2 |\xi|^2 \hat{G}(\xi) -\xi \cdot \nabla \hat{G}(\xi) - d \hat{G}(\xi) +  d \hat{G}(\xi) \\
= - 4\pi^2 |\xi|^2 \hat{G}(\xi) -\xi \cdot \nabla \hat{G}(\xi).
$$
This is a first order PDE that we can solve using the method of characteristics.  Suppose for the moment that $\xi_0 \in \mathbb{R}^d$ is such that $|\xi_0|=1$.  Consider the function $g(t) = \hat{G}(e^t \xi_0)$.  We compute 
$$
\frac{d}{dt} g(t) = \nabla \hat{G}(e^t \xi_0) \cdot e^t \xi_0 = -4\pi^2 |e^t \xi_0|^2 \hat{G}(e^t \xi_0) = -4\pi^2 e^{2t} g(t),
$$
which implies that (using standard ODE theory)
$$
g(t) = g(0) \exp\left(-2\pi^2(e^{2t} -1) \right)
$$
and hence that 
$$
\hat{G}(e^t \xi_0) = \hat{G}(\xi_0) \exp\left(-2\pi^2(e^{2t} -1) \right).
$$
Now, for any $\xi \neq 0$ we write $\xi = e^t \xi_0$ for $\xi_0 = \xi/|\xi|$ and $|\xi| = e^t$.  Plugging in above then shows that 
$$
\hat{G}(\xi) = \hat{G}\left(\frac{\xi}{|\xi|} \right)\exp\left(-2\pi^2(|\xi|^2 -1) \right).
$$
Since we have assumed that $G$ is integrable, we know that $\hat{G}$ is continuous, and so we have derived the most general form for $\hat{G}$, namely
$$
\hat{G}(\xi) = K\left(\frac{\xi}{|\xi|} \right) e^{-2\pi^2 |\xi|^2}
$$
for some continuous function $K$.  However, the exponential term is unity at $\xi=0$, so the continuity of $\hat{G}$ requires that actually $K$ is constant.
Hence $\hat{G}(\xi) = Ke^{-2\pi^2 |\xi|^2}$, which then implies that
$$
G(x) = \frac{K}{(2\pi)^{d/2}} e^{-|x|^2/2}.
$$
Let's now go back to your question of what assumptions have to be made on $G$.  To make the above analysis valid we need that $G$, $\Delta G$, and $x\cdot \nabla G$ all decay fast enough to justify applying the Fourier transform.  For example, we can assume that
$$
G, \Delta G, x \cdot \nabla G \in L^1(\mathbb{R}^d).
$$
A last remark: You should also be able to prove that $G$ is Gaussian if and only if $G$ is radial.  One direction is trivial.  For the other you assume that $G$ is radial and then rewrite the PDE as a second order ODE in $r$ and the resulting solution should be a Gaussian.
A: First, let me introduce the function $f : \mathbb{R}^n \rightarrow R | f(x) = \exp( \frac{1}{2} x \cdot x)$. Since it is nowhere zero and $\frac{1}{f}\nabla f = x$, you can rewrite $\nabla G + x G = 0$ as $\nabla G + \frac{G}{f}\nabla f = 0$, which is equivalent to :
$$ \nabla ( fG ) = 0 $$
Now, your second order equation reads $\Delta (fG) = \operatorname{div}\nabla(fG) = 0 $. Fortunately, you asked the question in $\mathbb{R}^n$, which has trivial cohomology so that Poincaré Lemma applies (see this Wikipedia page for more info). Therefore there exists a vector field $A$ defined on $\mathbb{R}^n$ such that :
$$ \nabla (fG) = \operatorname{curl} A$$ 
This means that $A$ and $fG$ are members of the same Helmoltz decomposition.
Basically, what you are asking is that $\operatorname{curl} A$ be zero. Again using Pointcaré lemma, this means there exists a function $\phi$ such that :
$$ A = \nabla \phi$$
Hope this helps !
A: For the functions from Schwartz functional space $S_{1,1}$ where $S_{\alpha,\beta} = |x^{\alpha}\frac{\partial^{\beta}}{\partial x^{\beta}}G| \lt \infty$ where norm means existence if integral over $R^d$ we can use multidimensional Fourier transform.
This end up in: $$0 = \digamma[\nabla(\nabla G + xG)] = i<k \cdot \hat{F}>$$ where $\hat{F} = ik\hat{G} +i\nabla\hat{G}$ and $\hat{G} = \digamma[G]$
Now, $\hat{F}$ vector is parallel to $k$, i.e. we do not have any perpendicular components that vanish in scalar product $<k \cdot \hat{F}>$. 
It is known that Schwartz space is dense in $L^p$, so we can approximate function from $L^p$ by the sequence of function from $S$ to converge weakly by $L^p$ norm or a.e. approximation.
