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 !
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.
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&lt;k \cdot \hat{F}&gt;$$ 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 $&lt;k \cdot \hat{F}&gt;$. 

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.