Does this interpretation $\phi'(x)=-x\phi(x)$ of the normal distribution have any significance? For the standard normal distribution $\phi(x)$, we can see that $\phi'(x)=-x\phi(x)$. Put differently, $\frac{\mathrm{d}\ln(\phi(x))}{\mathrm{d} x}= -x $. I see this as the fall in the value of the function being proportional to the distance from the mean, which is 0. In a way, this is like characterizing light tails, as the farther from the mean, the steeper is the fall in the value of the function.
Is the above interpretation sound? Is this the basis behind the ubiquity of the normal distribution?
 A: I like to think of it in a similar way but with slightly different differential equations. (edit: below I managed to make it also intuitive for $\phi'(x) = -x \phi(x)$)
Case: heat equation
$$\frac{d}{dt} \left[ \frac{1}{\sqrt{t}} \phi(x/\sqrt{t}) \right]= 0.5 \frac{d^2}{dx^2} \left[ \frac{1}{\sqrt{t}} \phi(x/\sqrt{t}) \right]$$
This is how I actually like to see the normal distribution and the central limit theorem. The sum of variables which is in the limit like a diffusion process (which follows a differential equation)
For a random walk where the steps are some variable $X_i$ with unit variance and zero mean then you get that the distribution function is approaching a scaled normal distribution
$$P \left( \sum_{i=1}^t X_i = x\right) \approx \frac{1}{\sqrt{t}} \phi \left(x/\sqrt{t}\right) = f(x,t)$$
The change of this function can be seen as a bit simultaneously as the differential function for Brownian motion and is like a wave spreading out in time following the equation:
$$\frac{\partial}{\partial t} f(x,t) = \frac{\partial^2}{\partial x^2} f(x,t)$$
See also https://en.m.wikipedia.org/wiki/Normal_distribution#Exact_normality

Case: $\mathbf{\phi(x) +x\phi'(x)+ \phi''(x) =0}$
Now when we divide by $t$
$${P\left(\frac{\sum_{i=1}^t X_i }{t}= x\right) }\approx \phi\left({x}\right) = g(x,t)$$
and we could describe it as some scaled random walk like: if $Z_t = Z_{t-1} + X_t$ then for $Y_t = Z_t/\sqrt{t}$
$$Y_t -Y_{t-1} = - Y_{t-1} \frac{\sqrt{t-1}-\sqrt{t}}{\sqrt{t}} + X_t \frac{1}{\sqrt{t}}$$
which you could see as shrinking the current value by a factor and adding a variable $X_t$ scaled with some factor determined by $t$. Then in the limit the change of $\frac{\partial}{\partial t} g(x,t)$ should balance these two processes
$$ \frac{\partial g(x,t)}{\partial t} = 
\underbrace{ a  g(x,t) }_{\substack{
\text{shrinking}\\
\text{moves values up}
}}
+\underbrace{
 \overbrace{x}^{
 \substack{
\llap{\text{further }}
\rlap{\text{away the}} \\
\llap{\text{shrinking }}
\rlap{\text{is stronger}}\\ \,
}}
\frac{\partial g(x,t)}{\partial x}}_{
\substack{
\text{shrinking} \\
\text{shifts/squeezes the function}
}} + 
\underbrace{ c \frac{\partial^2 g(x,t) }{\partial x^2 }}_{\text{diffusion}}=0$$
And the normal function is the function that such that the derivative of time (which is now expressed in terms of derivatives of space only) is zero.
So in this way we can relate intuitively $\phi(x) +x\phi'(x)+ \phi''(x) =0$ to a diffusion process with shrinking where the function remains unchanged.
Case: $\mathbf{x \phi(x) + \phi'(x) =0}$
This case is relatively similar to the one above.

*

*The shrinking of the formula can be related to a flux that is pulling the density to the inside. The flux is the product of the amount of mass (which is moving) and the speed of the mass (which relates to the distance)
$$ \text{flux}_\text{shrinking} =  -x \phi(x)$$


*The diffusion can be related to a flux that is related to the slope of the function. If at some point there is more density to one direction than the other then the convolution/diffusion will cause some density to flow donwards on the slope.
$$ \text{flux}_\text{diffusion} =  -\phi'(x)$$
When these two fluxes are opposite then there is no net flux and the function remains stable. So that is how you could view the relationship $\phi'(x) = -x\phi(x)$
I have made a computation to make the above intuitive idea more clear. In the computation I compute 1000 points according to some random distribution. And then I transform each point by scaling it with a factor $(1-c)$ and I am adding a centered Bernoulli variable to it with a factor $\sqrt{2c-c^2}$. This transformation will turn over time the distribution into a stable situation where the effect of the scaling is equal to the effect of the addition of the Bernoulli variable.

Below I have made two sketches for the intuition behind the terms in the differential equation.
(It is not a rigid derivation and one should go from the difference equations to differential and take the limit to very all linearizing the function, and also one could generalize the variable that is added, and represents the diffusion, which is now just a Bernoulli distributed variable. But I guess that in this way it is more intuitive and captures the essence more clearly)


# to plot points in the distribution
histpoints <- function(x, min, max) {
  counts <- rep(0, length(min:max))
  y <- rep(0,length(x))
  for (i in 1:length(x)) {
    bin <- ceiling(x[i]-min)
    counts[bin] <- counts[bin]+1
    y[i] <- counts[bin]
  }
  points(x,y, pch = 21, col = 1, bg = 1, cex = 0.4)
  counts
}

# transforming the points by
#  - scaling/shrinking
#  - and adding a Bernoulli variable
convertpoints <- function(x,c,var = 25) {
  x <- x * (1-c)  # scaling
  x <- x + sqrt(2*c-c^2) * (-1+2*rbinom(length(x), size = 1, prob = 0.5))*sqrt(var) # adding noise term
  return(x)
}

# make 2000 points according to some funny distribution
set.seed(1)
start <- seq(-20,20,0.01)
x <- sample(start, 1000, replace = TRUE, prob = 20+start^2-(20^-2+20^-3)*start^4)

# plot initial histogram
layout(matrix(1:8,4))
par(mar=c(3,1,2,1))
hist(x, breaks = c(-40:40), xlim=c(-25,25), ylim = c(0,80), main = "begin", xlab = "", yaxt = "n", ylab = "", xaxt = "n")
bins <- histpoints(x,-30,30)

for (j in 1:7) {
  for (i in 1:(100)) {
    x <- convertpoints(x,0.003)
  }
  
  #plot transformed
  hist(x, breaks = c(-40:40), 
       xlim=c(-25,25), ylim = c(0,80), main = paste0("after ",j*100," transforms"), xlab = "", yaxt = "n", ylab = "", xaxt = "n")
  bins <- histpoints(x,-30,30)
}

A: That differential equation is how Gauss arrived at the normal distribution in 1809.
Gauss wanted to rationalize the choice of the average as an estimator of a location parameter.
He imposed the following conditions for the distribution of errors:


*

*The density function $\phi(x)$ is differentiable.

*$\phi(-x) = \phi(x)$.

*$\phi(x)$ is maximum at $x=0$.

*Given multiple measurements of the same quantity corrupted by additive i.i.d. errors, the most likely value of the quantity is the average of the measurements.


From these conditions, he obtained the differential equation $\phi'(x) = -hx \phi(x)$ from which the normal pdf follows (recognizing $h$ as the precision parameter).  You can find the full derivation (in modern notation) in "The Evolution of the Normal Distribution" by Saul Stahl.
