The credits of this should go to @rchilton1980. I only want to validate the hypothesis of dispersion as the result of the response of a spike to the finite difference operator.
$\bf{ Basic \; Theory: }$
Physical dispersion accounts for the phenomena of waves moving in a way
that different frequencies or wavenumbers travel at different velocities.
We can make an analogy with a group of cliclist traveling along a flat road
at constant speed. Once they hit the mountain it could happen that the heaviest
cliclists will ride slower and the group will spread (or disperse) along the road. If we
think of a wavelet as the supperposition (synthesis) of many frequencies (and wavenumbers)
it could happen that the velocity of each frequency is different. However we
need to be more precise on the meaning of velocity here.
We define two velocities:
$\bf{Phase \; velocity} $
\begin{equation}
v_f= \frac{\omega}{k}
\end{equation}
$\bf { Group \; Velocity: }$
\begin{equation}
v_g= \frac{\partial \omega}{\partial k}
\end{equation}
The previous definitions indicate, in a way, that the angular frequency $\omega$
should be a function of the wavenumber $k$. That is, $\omega=\omega(k)$.
To better understand this let us consider the 1D wave equation with constant
wavespeed $c$.
\begin{eqnarray*}
\frac{\partial^2 u}{\partial x^2} -\frac{1}{c^2} \frac{\partial^2 u}{\partial t^2} = 0.
\end{eqnarray*}
If we take a double Fourier transform (from time $t$ to frequency $\omega$, and
from space $x$ to wavenumber $k$) we find
\begin{eqnarray*}
\left ( -k^2 + \frac{\omega^2}{c^2} \right ) U(k,\omega) = 0,
\end{eqnarray*}
where $U(k,\omega)$ is the Fourier transform of $u(x,t)$.
We get, from here
\begin{eqnarray*}
\omega = \pm c k
\end{eqnarray*}
The "$+$" sign indicates a wave traveling to the right, while
The ``$-$'' sign indicates a wave traveling to the left.
From now on, we will ignore the sign on this discussion. Then
\begin{eqnarray*}
\frac{\omega}{k} = \frac{\partial \omega}{\partial k} = c
\end{eqnarray*}
So we find that both phase and group velocity coincide with the wavespeed $c$.
So there is no physical dispersion.
However, when we find a finite difference approximation we could introduce
dispersion (numerical dispersion). That is, both phase and group velocities
could change with frequency (or wave number).
Let us consider the following finite difference scheme for the wave equation
without source.
\begin{eqnarray}
w_{i j+1} = \gamma^2 w_{i-1 j} + 2(1-\gamma^2) w_{ij} +
\gamma^2 w_{i+1 j} - w_{i j-1},
\end{eqnarray}
where $\gamma$ is the CFL number $c \Delta t/\Delta x$.
We take the Fourier transform of this expression, and using the shift property of the Fourier
transform to find
\begin{eqnarray*}
W_{ij} \left ( \mathrm{e}^{\mathrm{i} \omega \Delta t} +
\mathrm{e}^{-\mathrm{i} \omega \Delta t} \right ) = \left (
\gamma^2 \mathrm{e}^{-\mathrm{i} k \Delta x}
+ 2(1 - \gamma^2) W_{ij} + \gamma^2 \mathrm{e}^{\mathrm{i} k \Delta x} W_{ij} \right )
\end{eqnarray*}
where $W_{ij}$ is the Fourier transform of $w_{ij}$. That is, we find
\begin{eqnarray*}
\left [ 2 \cos \omega \Delta t - 2 \gamma^2 \cos k \Delta x
+ 2(1 - \gamma^2) \right ] W_{ij} = 0.
\end{eqnarray*}
We isolate $\omega$ from here so that
\begin{eqnarray*}
\omega \Delta t = \arccos \left ( \gamma^2 \cos k \Delta x - (1 - \gamma^2) \right ),
\end{eqnarray*}
That is,
\begin{eqnarray*}
v_f = \frac{\omega}{k} = \frac{1}{k \Delta t} \arccos
\left ( \gamma^2 \cos k \Delta x - (1 - \gamma^2) \right ),
\end{eqnarray*}
Clearly $v_f$ changes with wavenumber $k$. That is there is numerical dispersion.
The group velocity is given by
\begin{eqnarray*}
\frac{\partial \omega}{\partial k} = -\frac{1}{\Delta t}
\frac{ -\gamma^2 \sin k \Delta x}
{ \sqrt{1 - \gamma^2 \cos k \Delta x - (1 - \gamma^2) }}.
\end{eqnarray*}
Again, here group velocity $v_g$ changes with $k$
There is a special case given by $\gamma=1$. Here
\begin{eqnarray*}
w \Delta t = \arccos(\cos k \Delta x)= k \Delta x
\end{eqnarray*}
or
\begin{eqnarray*}
\frac{\omega}{k} = \frac{\Delta x}{\Delta t} \\
\frac{\partial \omega}{\partial k} = \frac{\Delta x}{\Delta t} \\
\end{eqnarray*}
and since $\gamma=c \Delta t/\Delta x =1$, then $c=\Delta x/\Delta t$ and
the phase and group velocities coincide with $c$. So there is no dispersion.
If $\gamma < 1$ then we can expect dispersion as shown in the original question.
Here I include the plots after a few fractions of second of propagation.
