# 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

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

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

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

whereas the dynamic structure factor \(S(\boldsymbol{q},\omega)\) is the Fourier transform of \(F(\boldsymbol{q},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

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

Now the correlation functions can be computed as

The longitudinal current can be related to the particle density as

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

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

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

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

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

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

In the time domain this function becomes

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

or in energy units as

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

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

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

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