# Density of phonon states from dispersion relation

I have a dispersion relation $$\omega( \textbf{q} ) = \omega_0 \sqrt{ \sum_{j=1}^{D} \sin^2{\frac{q_ja}{2}}. }$$ Where D is the dimension D=1,2,3. And my excersise is to calculate (numericaly on a computer) the density of states. And then to compare this density with the classical Debye-model density. I have found a formula for the calculation. $$g(\omega) = \frac{1}{V} \int_S \frac{dS}{\vert\nabla_q(\omega(q))\vert}$$ Where S is the surface of constant frequency. My question is what exactly does this surface look like ? Especially how do I create it in the program ? I don´t know over what surface I am supposed to integrate is it something like $\omega +d\omega$ or am I supposed to solve the equation $\omega=const.$ first ? Thanks for any answers.

Calculating DOSs can be tricky business. It sounds like for your purposes, you don't have to be super accurate. So, I suggest using a quick and dirty way of calculating the DOS:

1. Choose a bunch (at least a million, but feel free to go wild) $\mathbf{q}$ points at random.
2. For each one, find $\omega\left(\mathbf{q}\right)$
3. Make a histogram of the $\omega$ values you calculated.

That histogram is (proportional to) your DOS.

How does this work? The DOS basically tells you "how many" states have a given energy (or really, how many states are in a small energy window). So, if there are more states at a given energy, then a randomly chosen state is also more likely to have that energy. In short, the DOS is proportional to the probability distribution of the energy for a randomly chosen $\mathbf{q}$. Making a histogram in the way described above is a way to calculate the probability distribution. You'll have to figure out the proportionality constant (or fudge it, I think the Debye model should produce the same results at small $\mathbf{q}$). All this follows from an alternate (but equivilent) definition of the DOS: $g\left(\omega\right) = \int \frac{d\mathbf{q}}{\left(2 \pi \right)^d} \delta\left(\omega - \omega\left(\mathbf{q}\right)\right)$.

Another option is to convert that surface integral into a volume integral using the divergence theorem. This can be nice because finding the surface itself is tricky, but finding if a given $\mathbf{q}$ is inside the volume is easy. Just check that $\omega\left(\mathbf{q}\right) < const$. Computers are much better at volume integrals than surface integrals. This method is tricky, altho it's not totally crazy in your case since you have an analytic form for your dispersion relation. You could take the divergence analytically then do the integral numerically. However, I'd go with the first method unless you need a really accurate answer. There are many more advanced methods out there too, but they're probably overkill for your situation.

(In answer to your question, the surface you'd want to integrate over is an equpotential, so you're right that it's defined by $\omega\left(\mathbf{q}\right)= const$. If you want to do the surface integral, you'd first need to find the surface and then integrate over it. I'd avoid that at all costs if I were you.)