Considering an SDE of the form
$$dX_t = \mu(X_t, t)dt + \sigma(X_t, t)dW_t ,$$
where $W_t$ is a Wiener process, is there a set of necessary and sufficient conditions on the structure of the functions $\mu$ and $\sigma$ such that $X_t$ is multivariate normal (found perhaps by solving the Fokker-Planck equation)?
First, (to address some of the comments under the main post) it should be noted that we can easily formulate an SDE whose solution is not Gaussian at fixed times. A famous example is the geometric Brownian,
$$ X_t = \mu X_t dt + \sigma X_t dB_t $$
The solution to which is,
$$ X_t = X_0 \exp((\mu - \sigma^2/2)t + \sigma B_t) $$
as we can verify by a simple application of Ito's lemma. Letting $f(t, x) = \exp((\mu - \sigma^2/2)t + \sigma x)$
\begin{align*} dX_t &= f_t(t, B_t)dt + f_x(t, B_t) dB_t + 1/2f_{xx}(t, B_t)d[B_t,B_t] \\ &= ((\mu - \sigma^2/2)t dt + \sigma dB_t + \sigma^2/2 dt)X_t \\ &= \mu X_t dt + \sigma X_t dB_t \end{align*}
It is clear that, fixing $t$, we have $\log(X_t) \sim N( \log(X_0) + (\mu - \sigma^2/2)t, \sigma^2 t)$. So $X_t$ is not a GP.
Now, with that said, I don't know of a necessary condition of the type you are after. However a widely-known sufficient condition is that $\mu$ and $\sigma$ be constant in $X_t$. For details on a proof see this post, https://math.stackexchange.com/questions/4036081/why-is-int-0t-gsdb-s-sim-mathcaln-left0-int-0t-g2s-ds-right.
We can assume a Gaussian as solution and insert that into the Fokker Planck equations, to get the potential forms of $\mu$ and $\sigma$ that have the Gaussian as solution.
Let's describe the potential solutions as a Gaussian distribution with a mean $M(t)$ and variance $V(t)$ that can be time dependent:
$$g=g(x,t) = \frac{1}{\sqrt{2\pi V(t)}} e^{-\frac{(x-m(t))^2}{2V(t)}}$$
Where I write $g$ as short for $g(x,t)$
The Fokker Planck equation will have a time and place dependent diffusion term and an external force.
$$\begin{array}{rcl} g_t &=& - \overbrace{ \left[ \mu g\right]_{x}}^{\text{force term}} + \overbrace{ \left[Dg\right]_{xx}}^{\text{diffusion term}} \\ &=& g_x\mu + g\mu_x + D_{xx}g + 2 D_xg_x + Dg_{xx} \end{array} $$
where we use subscripts for partial derivatives and $\mu = \mu(x,t)$ and $D=D(x,t)$ are considered functions of time and space.
the partial derivatives of the Gaussian density $g$ (with potentially time varying mean and variance) are
$$\begin{array}{rcl} g_t/g &=& M_t \frac{x-M}{V} + V_t\frac{(x-M)^2-V}{2V^2} \\ g_{x}/g &=& - \frac{x-M}{V} \\ g_{xx}/g&=& - \frac{1}{V} +\left( \frac{x-M}{V}\right)^2 \\ \end{array}$$
which we can fill into the Fokker Planck equation
$$M_t(x-M) + 0.5 V_t\frac{(x-M)^2-V}{V} = (x-M) \mu + V (-\mu_x + D_{xx}) - 2 D_x (x-M) + D\frac{(x-M)^2-V}{V}$$
and after some rearrangement
$$0 = (x-M) (\mu-M_t) + V (-\mu_x + D_{xx}) - 2 D_x (x-M) + (D-0.5V_t) \frac{(x-M)^2-V}{V}$$
The time variations of variance $V_t$ and mean $M_t$ can be absorbed into the diffusion $D$ and drift $\mu$, and we can treat them as equal to zero without loss of generality. Also we can apply a transformation and scaling of the $x$ coordinate such that $M = 0$ and $V=1$.
This simplifies the equation to a stationary solution with $M=1$ and $V=1$
$$x\mu - \mu_x = D_{xx} - 2 x D_x + D (x^2-1)$$
At this time it might be better to go back to using the original equation
$$\left[Dg\right]_{xx} - \left[\mu g\right]_{x}= \left[(-xD+D_x-\mu)g\right]_{x} = 0$$
Note that if $h = c e^{x^2/2}$ then $\left[hg\right]_{x} = 0$
So for stationary solutions with zero mean and variance one we have the condition
$$\mu = -xD+D_x + c e^{x^2/2}$$
and we can get other solutions by scaling and translation of the $\mu$ and $D$ that satisfy this condition.
This condition is a necessary condition but not a sufficient condition. When we use $c \neq 0$ in the example below then we do not get a stationary solution. The reason is because it represents an unrealistic case where a constant flux is added. If the flux is constant everywhere then the change in time is zero. However, what does this constant flux mean? It means that at infinity there must be a (unrealistic) non-zero flux as well and the force term is blowing up exponentially in order to achieve this.
For example let's take $D = 2+x^2$ then $\mu = c \cdot e^{x^2/2} - x^3$
Then we could simulate the distribution with a random walk and it looks as following:
set.seed(1)
xs = seq(-5,5, 0.01)
drift = function(x) {
0*exp(x^2/2) - x^3 + 0*x
}
Diff = function(x) {
2+x^2
}
n = 10^6
x = rep(NA,n)
x[1] = 0
dt = 0.01
for (i in 2:n) {
b = rnorm(1,0,sqrt(dt))
x[i] = x[i-1] + drift(x[i-1])dt + sqrt(2Diff(x[i-1]))*b
}
hist(x, breaks = seq(-40,40,0.2), freq = 0, xlim = c(-5,5))
lines(xs,dnorm(xs))
