# Phonon density of states from velocity autocorrelation function

I'm using molecular dynamics and I autocorrelate the velocities and Fourier transform them to obtain the phonon density of states (DOS). I have many doubts about this:

1. The definition of DOS is: number of states per interval of energy at each energy level that ARE AVAILABLE to be occupied. If this is ok, then what I obtain autocorrelating velocities shouldn't be the DOS. I think that what I obtain is the population of the states, because with the velocities that the system have you should obtain the phonons that the system have, not the phonons that the system could have. I think it is something like the product of the probability distribution function and the DOS, but not just the DOS. Am I wrong? And if what I think is correct, the probability distribution function is the Bose-Einstein distribution (because the phonons are quantum particles) or the Maxwell-Boltzmann distribution (because of the classic behaviour of the atoms in my simulations)?

2. In all places that I read about this they say that the system is in equilibrium. If I do a non-equilibrium molecular dynamic (NEMD) simulation, can't I correlate the velocities and obtain information about the phonons? Can I take a small part of the system, assume that this small part is in local equilibrium, correlate the velocities of the atoms in this part and obtain information about the phonons? What would I obtain?

I rencently stumbled on a more or less similar question. Actually the Fourier transform of the velocity autocorrelation does not gives you the phonon DOS but rather the phonon population of your system. In other words it gives you normal modes.

The actual calculation of the DOS can be found, for example, in the thesis of Dr Hugo Ruiz, gracefully linked by Prof. G. Naumis in this discussion: https://www.researchgate.net/post/How_do_I_calculate_Phonon_Density_of_states_from_VACF

As Prof. Naumis says, the actual calculation is rare in the litterature and it could clarify some question to develop it here.

Unfortunately it is in spanish but I can give here a quick explaination of his calculations.

A. Computation of the phonon DOS

Given the Fourier transform of the velocities from the space of time ($t$) to the space of frequencies ($\omega$):

$$\mathbf{v}_n(\omega)=\int_{-\infty}^{\infty}\mathbf{v}_n(t)e^{i\omega{}t}dt$$

where the subscript $n$ corresponds to the $n^{th}$ atom and $i$ is the imaginary unit. From here one can get the spectrum of potential kinetic energy:

$$|\mathbf{v}_n(\omega)|^2 = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \mathbf{v}_n(t') \mathbf{v}_n(t) e^{i\omega{}(t-t')}dtdt'$$

In a stationnary oscillatory state, the coordinate $r_n(t)$ can be rewritten as a function of the normal modes of oscillation:

$$r_{nj}(t) = \sum_s Q_{snj} e^{-i\omega_{sj} t}$$

with the $\omega_s$ the normal frequencies of oscillation and j the direction of the 3D space, $Q_{snj}$ being the mean coordinate of the particle. Hence one can express velocities as the derivative regarding time:

$$v_{nj}(t) = \sum_s Q_{snj}(-i\omega_{sj}) e^{-i\omega_{sj} t}$$

Using this expression in the second integral, one can deduce:

$$\sum_{i=1}^N |\mathbf{v}_n(\omega_s)|^2 = \sum_s\sum_{(n*j)=1}^{3N} \int_{-\infty}^{\infty} |Q_{snj}|^2\omega_s^2 e^{i(\omega+\omega_s)t''}dt''$$

where $t''=t-t'$.

At thermal equilibrium, given the equipartition of the energy, one has directly $\sum_j |Q_{snj}|^2\omega_{sj}^2=3k_BT$. As the normale modes are supposed to be dominant, the integral of the complex exponential can be reduced to a Dirac distribution such as $\rho(\omega) = \sum_s \delta(\omega+\omega_s)$ is the phonon DOS.

Hence: $$\sum_{i=1}^N |\mathbf{v}_n(\omega_s)|^2 = 3Nk_BT\sum_s\delta(\omega+\omega_s)$$ $$\downarrow$$ $$\rho(\omega) = \frac{\sum_{i=1}^N |\mathbf{v}_n(\omega_s)|^2}{3Nk_BT}$$

Taking $t'=0$ for the origin of the time, this final expression reduces to $$\rho(\omega) = \frac{\int_{-\infty}^{\infty}\sum_{i=1}^N <\mathbf{v}_i(t)\cdot\mathbf{v}_i(0)>}{3Nk_BT}$$ linking it to the autocorrelation. If you just Fourier transform the autocorrelation of the velocities, you get the normal modes. You have to had them (i.e integrate) to get the actual DOS. More informations on the assumptions are in the thesis.

B. Concerning phonons

Phonon actually follow a Bose-Einstein distribution as they can be created and annihilated by energy fluctuations just like bosons particles (photons for example) independently of a classical/quantum behaviour of the system. As such they follow a Bose-Einstein distribution with a chemical potential equal to zero:

$$<n_i> = \frac{1}{e^{\frac{\hbar\omega_i}{k_BT}-1}}$$

Furthermore concerning NEMD it depends on the kind of constraint you impose on your system. As you can see in the derivation of the DOS, equipartition of the energy and thermal equilibrium are important assumptions. The main problem in NEMD will be the fluxes, not only of matter (which can appen even in solids and I don't know what kind of system you are studying) but also of heat which can end up causing some flux of momentum. The problem with fluxes is that they break the isotropy of your system and you can end up with local non equipartition of the kinetic energy, breaking one of the assumptions of the derivation. To assume local equilibrium you have to be sure that the input and output of energy and/or matter compensate in every directions for a long enough time and on a large enough region to do your statistics. This is far to be obvious even (or especially) in stationnary out of equilibrium states.

I think a better way to determine the normal modes of higher energy is to simulate your system at a higher temperature in the canonical ensemble with a good thermostat (Nosé-Hoover Chain if you can).

I hope this helps. Any comment is welcome.