Theoretical background

Below follows a brief overview of the theory underpinning the correlation functions supported by dynasor as well as their analysis. For a more detailed description and derivations of the correlations functions see for example [Boon and Yip, 1980].

Dynamic structure factor

The particle density is denoted \(n(\boldsymbol{r},t)\) and defined as

\[n(\boldsymbol{r},t) = \sum _i ^N \delta (\boldsymbol{r} - \boldsymbol{r}_i(t)),\]

where \(N\) is the number of particles and \(\boldsymbol{r}_i(t)\) is the position of particle \(i\) at time \(t\). The intermediate scattering function \(F(\boldsymbol{q},t)\) is defined as the Fourier transform of the auto-correlation of the particle density

\[F(\boldsymbol{q},t)=\frac{1}{N}\left<n(\boldsymbol{q},t)n(-\boldsymbol{q},0)\right> \,= \, \frac{1}{N}\sum _i ^N \sum _j ^N \left< \mathrm{exp} \left[i\boldsymbol{q}\cdot(\boldsymbol{r}_i(t)-\boldsymbol{r}_j(0))\right] \right>,\]

where \(\boldsymbol{q}\) is a wave vector. The brackets denote an ensemble average, which in the case of dynasor is replaced by a time average. (The system under study should be ergodic in order for the time average to be a suitable substitute for the ensemble average.) Here, \(F(\boldsymbol{q},t)\) is commonly referred to as the coherent part.

The incoherent intermediate scattering function, i.e., \(i=j\), is defined as

\[F_\mathrm{incoh}(\boldsymbol{q},t)=\frac{1}{N} \sum _i ^N \left < \mathrm{exp}\left[i\boldsymbol{q}\cdot(\boldsymbol{r}_i(t)-\boldsymbol{r}_i(0))\right] \right >\]

and describes single particle motion, i.e., diffusive motion. \(F_\mathrm{incoh}(\boldsymbol{q},t)\) is commonly referred to as the incoherent scattering function, but sometimes also referred to as the self-part of the intermediate scattering function.

The static structure factor is given by

\[S(\boldsymbol{q}) = F(\boldsymbol{q},0),\]

whereas the dynamic structure factor \(S(\boldsymbol{q},\omega)\) is the Fourier transform of \(F(\boldsymbol{q},t)\)

\[S(\boldsymbol{q},\omega) = \int _{-\infty} ^\infty F(\boldsymbol{q},t) \, e^{-iwt} \mathrm{d}t,\]

where \(\omega\) is the angular frequency. \(S(\boldsymbol{q},\omega)\) exhibits peaks in the \((\boldsymbol{q},\omega)\) plane corresponding to the dispersion of lattice vibrations (phonons). The widths of these peaks are related to the phonon lifetimes. Computing \(S(\boldsymbol{q},\omega)\) via molecular dynamics simulation has the advantage of fully including anharmonic effects and allowing one to study the temperature dependence of phonon dispersions.

Velocity correlation functions

It is often convenient to also consider current correlations based on the particle velocities. The current density \(\boldsymbol{j}(\boldsymbol{r},t)\) is given by

\[\begin{split}\boldsymbol{j}(\boldsymbol{r},t) &= \sum_i^N \boldsymbol{v}_i(t) \, \delta (\boldsymbol{r} - \boldsymbol{r}_i(t)) \\ \boldsymbol{j}(\boldsymbol{q},t) &= \sum_i^N \boldsymbol{v}_i(t) \, \mathrm{e}^{\mathrm{i}\boldsymbol{q} \cdot \boldsymbol{r}_i(t)},\end{split}\]

where \(\boldsymbol{v}_i(t)\) is the velocity of particle \(i\) at time \(t\). This can be split into a longitudinal and transverse part as

\[\begin{split}\boldsymbol{j}_L(\boldsymbol{q},t) &= \sum_i ^N(\boldsymbol{v_i}(t) \cdot\hat{\boldsymbol{q}}) \, \hat{\boldsymbol{q}} \, \mathrm{e}^{\mathrm{i}\boldsymbol{q} \cdot \boldsymbol{r}_i(t)} \\ \boldsymbol{j}_T(\boldsymbol{q},t) &= \sum_i ^N \left[\boldsymbol{v_i}(t) - (\boldsymbol{v_i}(t) \cdot \hat{\boldsymbol{q}}) \, \hat{\boldsymbol{q}}\right] \, \mathrm{e}^{\mathrm{i}\boldsymbol{q} \cdot \boldsymbol{r}_i(t)}.\end{split}\]

Now the correlation functions can be computed as

\[\begin{split}C_L(\boldsymbol{q},t) &= \frac{1}{N}\left<\boldsymbol{j}_L(\boldsymbol{q},t)\cdot\boldsymbol{j}_L(-\boldsymbol{q},0)\right> \\ C_T(\boldsymbol{q},t) &= \frac{1}{N}\left<\boldsymbol{j}_T(\boldsymbol{q},t)\cdot\boldsymbol{j}_T(-\boldsymbol{q},0)\right>.\end{split}\]

The longitudinal current can be related to the particle density as

\[\begin{split}\frac{\partial }{\partial t}n(\boldsymbol{q}, t) &= \mathrm{i} \boldsymbol{q}\cdot \boldsymbol{j}(\boldsymbol{q}, t) \\ \omega^2 S(\boldsymbol{q},\omega) &= q^2 C_L(\boldsymbol{q},\omega),\end{split}\]

which means some features can be easier to resolve in one function than the other. The current correlation functions can be thought of as spatially-dependent generalization of the velocity correlation function and are closely related to the phonon spectral energy density.

Multi-component systems

In multi-component systems one can introduce partial correlations functions, which enables separation of the contributions by groups of particles, e.g., by type or site symmetry. For example, in the case of a binary system, one can define \(F_\mathrm{AA}(\boldsymbol{q},t)\), \(F_\mathrm{BB}(\boldsymbol{q},t)\), \(F_\mathrm{AB}(\boldsymbol{q},t)\) and so on for all correlation functions introduced above. In dynasor we define the partial correlaion functions as

\[F_\mathrm{AA}(\boldsymbol{q},t) = \frac{1}{N}\sum_{i\in A}^{N_\mathrm{A}} \sum_{j\in A}^{N_\mathrm{A}} \left< \mathrm{exp} \left[i\boldsymbol{q}\cdot(\boldsymbol{r}_i(t)-\boldsymbol{r}_j(0))\right] \right>\]

where the sums runs over all A atoms. Here, \(N_\mathrm{A}\) refers to number of atoms of type A, \(N_\mathrm{B}\) refers to number of atoms of type B, and \(N = N_\mathrm{A} + N_\mathrm{B}\) refers to the total number of atoms in the system. For the cross-term, \(F_\mathrm{AB}(\boldsymbol{q},t)\), we get

\[\begin{split}F_\mathrm{AB}(\boldsymbol{q},t) = & \frac{1}{N}\sum_{i\in A}^{N_\mathrm{A}} \sum_{j\in B}^{N_\mathrm{B}} \left< \mathrm{exp} \left[i\boldsymbol{q}\cdot(\boldsymbol{r}_i(t)-\boldsymbol{r}_j(0))\right] \right> + \frac{1}{N}\sum_{i\in B}^{N_\mathrm{B}} \sum_{j\in A}^{N_\mathrm{A}} \left< \mathrm{exp} \left[i\boldsymbol{q}\cdot(\boldsymbol{r}_i(t)-\boldsymbol{r}_j(0))\right] \right> \\=& \frac{2}{N}\sum_{i\in A}^{N_\mathrm{A}} \sum_{j\in B}^{N_\mathrm{B}} \left< \mathrm{exp} \left[i\boldsymbol{q}\cdot(\boldsymbol{r}_i(t)-\boldsymbol{r}_j(0))\right] \right>\end{split}\]

Here, the factor two in the final expression is due to considering both the A-B and B-A terms. These normalization choices means the total intermediate scattering function \(F(\boldsymbol{q},t)\) (as defined above) is given by

\[F(\boldsymbol{q},t) = F_\mathrm{AA}(\boldsymbol{q},t) + F_\mathrm{AB}(\boldsymbol{q},t) + F_\mathrm{BB}(\boldsymbol{q},t)\]

In some cases, instead of analyzing the partial functions directly, it is interesting to consider linear combinations of them. For example weighting them with relevant atomic form factors, masses or charges etc.

Relation to scattering experiments

The correlation functions can be convoluted with atomic form factors and cross sections in order to obtain predictions for various types of neutron and X-ray scattering experiments.

In solids, it is frequently desirable to determine the above mentioned quantites along specific paths between high-symmetry \(\boldsymbol{q}\)-points. In isotropic samples, such as for example liquids, it is on the other hand usually preferable to compute these functions with respect to \(q=|\boldsymbol{q}|\), i.e., a spherical average over wave vectors.

Damped harmonic oscillator model

The correlation functions can be fitted to analytical expressions for damped harmonic oscillators in order to extract phonon frequencies and phonon lifetimes. While many different analytical forms can be considered for the damped harmonic oscillator, here we use the assumption that the correlation function of the position of the damped harmonic oscillator has a zero derivative at \(t=0\). The position correlation function (corresponding to the intermediate scattering function) of a damped harmonic oscillator is given by

\[\begin{split}F(t) &= A \mathrm{e}^{-\Gamma t/2} \big(\cos{ \omega_e t} + \frac{\Gamma}{2\omega_e}\sin{ \omega_e t}\big), \quad \omega_0 > \frac{\Gamma}{2} \\ F(t) &= A \mathrm{e}^{-\Gamma t/2} \big(\cosh{|\omega_e| t} + \frac{\Gamma}{2|\omega_e|}\sinh{|\omega_e| t}\big), \quad \omega_0 < \frac{\Gamma}{2}\end{split}\]

where the fitting parameters are \(\omega_0\), \(\Gamma\) and \(A\). The first case, \(\omega_0 > \Gamma/2\) corresponds to an underdamped oscillator, and the second one to an overdamped oscillator. The eigenfrequency of the damped oscillator is \(\omega_e = \sqrt{\omega_0^2 - \Gamma^2/4}\). The fitting can also be carried out in the frequency domain by fitting the dynamic structure factor to the following analytical function

\[a(\omega) = A\frac{2\Gamma \omega_0^2}{(\omega^2 - \omega_0^2)^2 + (\Gamma\omega)^2}.\]

The intermediate scattering function and dynamic structure factor can thus be analyzed by fitting them to the expressions above in order to extract frequency \(\omega_0\) and damping \(\Gamma\).

This analysis can be extended to the velocity correlation for the damped harmonic oscillator (corresponding to the current correlations) by considering the relation between the longitudinal current correlation and the dynamic structure factor, which gives the following solutions

\[b(\omega) = B\frac{2\Gamma \omega^2}{(\omega^2 - \omega_0^2)^2 + (\Gamma\omega)^2}.\]

In the time domain this function becomes

\[\begin{split}b(t) &= B \mathrm{e}^{-\Gamma t/2} \big(\cos{ \omega_e t} - \frac{\Gamma}{2\omega_e}\sin{ \omega_e t} \big), \quad \omega_0 > \frac{\Gamma}{2} \\ b(t) &= B \mathrm{e}^{-\Gamma t/2} \big(\cosh{|\omega_e| t} - \frac{\Gamma}{2|\omega_e|}\sinh{|\omega_e| t} \big), \quad \omega_0 < \frac{\Gamma}{2}\end{split}\]

Lastly, the damping \(\Gamma\) is related to the inverse phonon lifetime \(\tau\) according to

\[\tau = \frac{2}{\Gamma}\]

or in energy units as

\[\tau = \frac{2\hbar}{\Gamma},\]

where \(\hbar\) appears since \(\Gamma\) has units of rad/s.

If multiple modes are present in a correlation function they should be modeled as a sum of damped harmonic oscillators.

Spectral energy density method

The spectral energy density (SED) is closely related to the phonon dispersion and thus also to the current correlations above. It is a measure of how the kinetic energy is partitioned over different wavelengths and frequencies in the system. The theoretical background for SED is quite simple and can be derived in several ways. The total, integrated kinetic energy of the infinite system can be written as

\[E = \int dt \sum_n \sum_i \frac{1}{2} m_i v_i(n, t)^2.\]

Here \(n\) is the cell index and \(i\) is the basis index, which together uniquely identify each atom in an infinite crystal. By using Parseval’s theorem for both the time integration and space summation the above is transformed to inverse spatial \(q\)-space and inverse temporal \(\omega\)-space

\[E = \int d\omega \sum_q \sum_i \frac{1}{2} m_i v_i(q, \omega)^2.\]

This is motivated by the Wiener–Khinchin theorem and the SED is simply defined as the quantity

\[\text{SED}(q, \omega) = \sum_i \frac{1}{2} m_i v_i(q, \omega)^2\]

For a more in-depth discussion see, e.g., [Thomas et al., 2010].