To get the dispersion relation, we put back in the time dependence of the velocities to get the function $\vec{v}_i(t)$. From this we form the velocity distribution $\vec{v}(\vec{r},t) = \sum_i \vec{v}_i(t) \delta(\vec{r} - \vec{r}_i)$. Now we can again do a spatial fourier transform to get $\vec{v}(\vec{k},t)$. We expect an oscillation to with wavevector $\vec{k}$ to osciallate at some frequency $\omega$, to find this frequency $\omega$, we take the temportal fourier transform of $\vec{v}(\vec{k},t)$ to get $\vec{v}(\vec{k},\omega)$. Now if we want to know the frequency corresponding to a given wavevector $\vec{k}_0$, then we just look at the frequency dependence of $\vec{v}(\vec{k}_0,\omega)$, and the $\omega$ for which this is non-zero give the frequency of the oscillation. So this procedure from getting an $\omega$ from a $k$ is a one way of getting the dispersion relation.

What do we expect to see? Let's start at $k=0$ and go up (the dispersion relation should have the symmetry $\omega(-k)=\omega(k)$). We expect $\omega(0) = 0$, and we expect that $\omega(k)$ should have a finite derivate as you increase $k$ from $0$. (This derivative is the speed of sound for this material.) Now we expect $\omega$ to keep increasing (albeit nonlinearly and even discontinuously at one point) until it gets to some maximum. This maximum occurs when the velocity oscillation of sodium is $pi$ out of phase with the velocity oscillation of chlorine. The half-wavelength of the velocity oscillation is $a$, so the wavelength must be $2a$ and the wavenumber must be $k_M=\pi/a$.

So to get acoustic and optical phonons from the $\omega{k}$ in your example you need to find this second form of the "fundamental zone" in reciprocal space, which is conventionally called the "first brillouin zone", and you need to do the appropriate folding of $\omega(\vec{k})$ into this zone. This is probably a little trickier in 3d.

To get the dispersion relation, we put back in the time dependence of the velocities to get the function $\vec{v}_i(t)$. From this we form the velocity distribution $\vec{v}(\vec{r},t) = \sum_i \vec{v}_i(t) \delta(\vec{r} - \vec{r}_i)$. Now we can again do a spatial fourier transform to get $\vec{v}(\vec{k},t)$. We expect an oscillation to with wavevector $\vec{k}$ to oscillate at some frequency $\omega$, to find this frequency $\omega$, we take the temporal fourier transform of $\vec{v}(\vec{k},t)$ to get $\vec{v}(\vec{k},\omega)$. Now if we want to know the frequency corresponding to a given wavevector $\vec{k}_0$, then we just look at the frequency dependence of $\vec{v}(\vec{k}_0,\omega)$, and the $\omega$ for which this is non-zero give the frequency of the oscillation. So this procedure from getting an $\omega$ from a $k$ is a one way of getting the dispersion relation.

What do we expect to see? Let's start at $k=0$ and go up (the dispersion relation should have the symmetry $\omega(-k)=\omega(k)$). We expect $\omega(0) = 0$, and we expect that $\omega(k)$ should have a finite derivative as you increase $k$ from $0$. (This derivative is the speed of sound for this material.) Now we expect $\omega$ to keep increasing (albeit nonlinearly and even discontinuously at one point) until it gets to some maximum. This maximum occurs when the velocity oscillation of sodium is $pi$ out of phase with the velocity oscillation of chlorine. The half-wavelength of the velocity oscillation is $a$, so the wavelength must be $2a$ and the wavenumber must be $k_M=\pi/a$.

So to get acoustic and optical phonons from the $\omega(k)$ in your example you need to find this second form of the "fundamental zone" in reciprocal space, which is conventionally called the "first brillouin zone", and you need to do the appropriate folding of $\omega(\vec{k})$ into this zone. This is probably a little trickier in 3d.

Edit: How to get acoustic and optical bands.

I think the problem is when you do things this way, you don't put an upper limit on $\vec{k}$. Typically people restrict $\vec{k}$ to the first brillouin zone by means of defining a projection from the first brillouin zone by means of projection. Let me describe what I mean by an example.

Suppose we have a 1D system in 1D space. Suppose this system is a chain of ions whose species alternates between sodium and chlorine. Suppose the spacing between two adjacent ions is always $a$. Suppose the ions interact with a harmonic potential interaction: $U(x_\mathrm{Cl}, x_\mathrm{Na}) = \frac{1}{2}k(x_\mathrm{Cl}- x_\mathrm{Na}-a)^2 $, where $k$ is a spring constant.

Now let us find the dispersion relation. As before, we look at velocity as a function of space and time. We do a spatial fourier transform to pick out velocity oscillations of a particular wavevector, and then we do a temporal fourier transform to get the frequency corresponding to each wavevector.

What do we expect to see? Let's start at $k=0$ and go up (the dispersion relation should have the symmetry $\omega(-k)=\omega(k)$). We expect $\omega(0) = 0$, and we expect that $\omega(k)$ should have a finite derivate as you increase $k$ from $0$. (This derivative is the speed of sound for this material.) Now we expect $\omega$ to keep increasing (albeit nonlinearly and even discontinuously at one point) until it gets to some maximum. This maximum occurs when the velocity oscillation of sodium is $pi$ out of phase with the velocity oscillation of chlorine. The half-wavelength of the velocity oscillation is $a$, so the wavelength must be $2a$ and the wavenumber must be $k_M=\pi/a$.

What happens if we increase $k$ past $k_M=\pi/a$? Well due to aliasing effects, a velocity oscillation with wave number $k_M + \Delta k$ gives the same particle velocities as a velocity oscillation of $k_M - \Delta k$ (since $\mathrm{Re}(\exp[i(k_Mna \pm \Delta k na)])=\mathrm{Re}(\exp[i(\pi n \pm \Delta k na)])=(-1)^n\mathrm{Re}(\exp[\pm i \Delta k na)])$, which is independent of the choice of plus vs minus). This tells us that $\omega(k_M + k) = \omega(k_M - k)$. Now we know what $\omega(k)$ is for $k< 2k_M$. To go beyond this we note that $\omega(k+2k_M) = \omega(k)$.

From the equations $\omega(-k)=\omega(k)$, $\omega(k_M + k) = \omega(k_M - k)$, and $\omega(k+2k_M) = \omega(k)$. It is clear that the represent the whole dispersion relation we need only consider $\omega(k)$ in the region $0<k<k_M$. This is one way of picking a "fundamental zone" of reciprocal space, and concentrating on that zone.

Now I still haven't gotten to acoustic and optical phonons. These come up when you pick the "fundamental zone" differently. What people usually done is to pick the fundamental zone to be half the size of the previously described one. That is, it is the region $0<k<k_M/2$. So that the information about the dispersion relation for $k>k_M/2$ doesn't get lost, we need to make two dispersion relations in the region $0<k<k_M/2$. One is the acoustic phonon dispersion relation defined by $\omega_a(k) = \omega(k)$ and the other is the optical phonon dispersion relation given by $\omega_o(k) = \omega(k_M-k)$, so that the optical phonon dispersion relation is the right half of the full dispersion relation folded over $k_M/2$.

So to get acoustic and optical phonons from the $\omega{k}$ in your example you need to find this second form of the "fundamental zone" in reciprocal space, which is conventionally called the "first brillouin zone", and you need to do the appropriate folding of $\omega(\vec{k})$ into this zone. This is probably a little trickier in 3d.

The impression I got form the equation is that you model your system as being composed of a collection of point particles. The index $i$ indexes this set of point particles (as opposed to indexing the components of a vector).

One could create a spatial velocity distribution by saying $\vec{v}(\vec{r})=\sum_i \vec{v}_i \delta(\vec{r} - \vec{r}_i)$, where $\vec{v}_i$ is the velocity of the $i$th particle and $\vec{r}_i$ is the position of the $i$th particle. Then the spatial fourier transform of this $\vec{v}(\vec{r})$ is $\vec{v}(\vec{k})$ or, as you have written it, $\vec{v}_\vec{k}$.

That being said, it is unclear how to get a dispersion relation from this.

The impression I got form the equation is that you model your system as being composed of a collection of point particles. The index $i$ indexes this set of point particles (as opposed to indexing the components of a vector).

One could create a spatial velocity distribution by saying $\vec{v}(\vec{r})=\sum_i \vec{v}_i \delta(\vec{r} - \vec{r}_i)$, where $\vec{v}_i$ is the velocity of the $i$th particle and $\vec{r}_i$ is the position of the $i$th particle. Then the spatial fourier transform of this $\vec{v}(\vec{r})$ is $\vec{v}(\vec{k})$ or, as you have written it, $\vec{v}_\vec{k}$.

To get the dispersion relation, we put back in the time dependence of the velocities to get the function $\vec{v}_i(t)$. From this we form the velocity distribution $\vec{v}(\vec{r},t) = \sum_i \vec{v}_i(t) \delta(\vec{r} - \vec{r}_i)$. Now we can again do a spatial fourier transform to get $\vec{v}(\vec{k},t)$. We expect an oscillation to with wavevector $\vec{k}$ to osciallate at some frequency $\omega$, to find this frequency $\omega$, we take the temportal fourier transform of $\vec{v}(\vec{k},t)$ to get $\vec{v}(\vec{k},\omega)$. Now if we want to know the frequency corresponding to a given wavevector $\vec{k}_0$, then we just look at the frequency dependence of $\vec{v}(\vec{k}_0,\omega)$, and the $\omega$ for which this is non-zero give the frequency of the oscillation. So this procedure from getting an $\omega$ from a $k$ is a one way of getting the dispersion relation.

There is one more subtlety which is that $\vec{v}(\vec{k},\omega)$ has three components. Thus at a given $\vec{k}$ and $\omega$, the component of $\vec{v}(\vec{k},\omega)$ parallel to $\vec{k}$ may be non-zero but the perpendicular component may be zero. This tells you that longitudinal and transverse waves oscillate at different frequencies.