# Constructing the density of states of multiple independent harmonic oscillators

I have a system of $$N$$ uncoupled 1D quantum harmonic oscillators, each with its own frequency $$\omega_i$$. The density of states for a single quantum harmonic oscillator shall be defined as

$$\rho(E) = \frac{dn}{dE}$$ It is trivial in the 1-D case to obtain $$n$$ as function of energy by rearranging the eigenenergy equation, $$E(n) = \hbar \omega(n+1/2) \leftrightarrow n(E) = \frac{E}{\hbar \omega } - \frac{1}{2\hbar \omega}$$ The density for a single oscillator is then $$\rho(E) = \frac{dn}{dE}=\frac{1}{\hbar \omega}$$ So far so good. My problem is now to derive the density of a system with the eigenenergy defined as $$E(n_1,n_2,\dots, n_N) = \sum^N_{i=1}\hbar \omega_i(n_i+1/2)$$

Is there a way to obtain/construct the density from the 1D solution for the N-dimensional system ? There seems no straightforward way to define a total $$n$$ and to apply the derivative. Or is there a way ?

The definition of the density-of-states is $$\rho(E)=\frac{dN(E)}{dE},$$ where $$N(E)$$ is the number of states with energies less than $$E$$.

Note, that as the equation involves a derivative, it is not defined for non-continuous spectra, such as that of a Harmonic oscillator. However, it can be still treated in terms of generalized functions, such as Heaviside step-function, $$\theta(x)$$,a nd delta-function, $$\delta(x)$$.

The number of oscillator states with energies less than $$E$$ is: $$N(E) =\sum_{n=0}^{+\infty}\theta(E-E_n) = \sum_{n=0}^{+\infty}\theta\left(E-\hbar\omega_0(n+\frac{1}{2})\right)$$ The resulting density-of_states is then $$\rho(E) =\sum_{n=0}^{+\infty}\delta(E-E_n) = \sum_{n=0}^{+\infty}\delta\left(E-\hbar\omega_0(n+\frac{1}{2})\right).$$

With the correct formulas the generalization to the case of multiple oscillators with frequencies $$\omega_k$$ is trivial: $$N(E) =\sum_{n=0}^{+\infty}\sum_k\theta(E-E_n^{(k)}) = \sum_{n=0}^{+\infty}\sum_k\theta\left(E-\hbar\omega_k(n+\frac{1}{2})\right),$$ and similarly for the density-of-states. If the oscillators spectrum is continuous, the summation over $$k$$ becomes an integral, which is easily evaluated using the properties of the delta-function, resulting in a continuous density-of-states.

Continuous approximation for a harmonic oscillator
If we are interested in energy scales much greater than the oscillator energy level spacing, we can approximate the sum in the oscillator dos by an integral, with $$\epsilon=\hbar\omega_0 n$$, $$d\epsilon=\hbar\omega_0$$: $$\rho(E) =\sum_{n=0}^{+\infty}\delta\left(E-\hbar\omega_0(n+\frac{1}{2})\right)\frac{d\epsilon}{\hbar\omega_0}= \frac{1}{\hbar\omega_0}\int_{0}^{+\infty}\delta\left(E-\epsilon -{\hbar\omega_0}{2})\right)d\epsilon=\frac{1}{\hbar\omega_0}.$$ Note that we are essentially in the classical limit here.

• Is there also an approximation to the final expression with the double sums ? I would like to avoid the "counting" if possible. I assumed there would be some sort of approximate function that can be obtained by integration, as often done in statistical mechanics. May 12, 2021 at 10:11
• It depends on the nature of your oscillators - if they form a continuous spectrum, you will get a continuous function. Also, if you are interested in energy scales bigger than $\omega_0$, you can approximate the sum by an integral. But these are two different things (two different sums). May 12, 2021 at 10:21
• I guess i should post a second question. My system are quantum oscillators. So your answer without approximations seems appropriate and correct. But my actual problem which caused this question is numerically and i need something that i can evaluate. At the moment i have 66 oscillators and that is not a particular large number, just a test system. Doing the sum explicitly is not practical. But i realize now that this was not part of the question as i asked it. May 12, 2021 at 10:32