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 weights, \(w_i\), i.e.,

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

where the weights could for example be masses, charges, or even atomic form factors (scattering lengths) to allow for direct comparison to structure factors measured in various neutron and X-ray scattering experiments. For some probes, e.g., neutrons, the scattering lenghts are in general constant but differ for the coherent and incoherent correlation functions. X-ray form factors, on the other hand, are \(q\)-dependent and thus take the form

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

Conveniently, since the weights only depend on the atom type, the weighting can be done as a post-processing step using partial correlation functions. For a binary (A-B) system, the total weighted correlation function is

\[F_\mathrm{coh}(\boldsymbol{q},t)=f^{AA}(q) F^{AA}_\mathrm{coh}(\boldsymbol{q},t) + f^{AB}(q) F^{AB}_\mathrm{coh} (\boldsymbol{q},t) + f^{BB}(q) F^{BB}_\mathrm{coh}(\boldsymbol{q},t).\]

Furthermore, it is frequently desirable to determine the above mentioned quantites along specific paths between high-symmetry \(\boldsymbol{q}\)-points when working with solids. 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.

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. Formally the kinetic action \(W\) of the system can be written as

\[W = \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

\[W = \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].

Damped harmonic oscillator model

The correlation functions can be fitted to analytical expressions for a damped harmonic oscillator (DHO) in order to extract phonon frequencies and lifetimes. The governing differential equation for the displacement (or coordinate) of a DHO, \(x(t)\), is

\[\ddot{x}(t) + \Gamma \dot{x}(t) + \omega_0^2 x(t) = f(t).\]

Here, \(\omega_0\) is the natural angular frequency and \(\Gamma\) is the damping constant in units of reciprocal time. Note that we work in units where the mass \(m=1\). Below follow the analytical solutions for the equations of motion of a DHO under the assumption that the stochastic force, \(f(t)\), is white-noise.

The DHO can be either underdamped or overdamped. In the underdamped case \(\omega_0\tau > 1\), where \(\tau = \frac{2}{\Gamma}\) is commonly referred to as the lifetime of the DHO. The position autocorrelation function (ACF) in the underdamped case is given by

\[F(t) = A \mathrm{e}^{- t/\tau} \left ( \cos{ \omega_e t} + \frac{1}{\tau\omega_e}\sin{ \omega_e t}\right ),\]

where the angular eigenfrequency of the DHO is \(\omega_e = \sqrt{\omega_0^2 - 1/\tau^2}\) and \(A\) is a constant amplitude \(A=F(t=0)\). The functional form is only valid for time lag \(t\geq 0\) and should be interpreted as even, so that \(t\rightarrow -t\) for negative times. Alternatively replace \(t\) with \(|t|\) in the expressions for the ACFs.

In the overdamped case, \(\omega_0 \tau < 1\), the autocorrelation function is given by

\[F(t) = \frac{A}{\tau_\text{L} - \tau_\text{S}} \left( \tau_\text{L} \text{e}^{-t / \tau_\text{L}} - \tau_\text{S} \text{e}^{-t / \tau_\text{S}} \right),\]

where \(\tau_\text{S}\) and \(\tau_\text{L}\) correspond to the short and long relaxation times, respectively, defined as

\[\tau_\text{S,L} = \frac{\tau}{1 \pm \sqrt{1-(\omega_0\tau)^2}}.\]

The corresponding underdamped and overdamped velocity ACFs are related to the position ACFs via \(F_v(t) = -\ddot{F}(t)\):

\[\begin{split}F_v(t) &= A \omega_0^2 \mathrm{e}^{- t/\tau} \left ( \cos \omega_e t - \frac{1}{\tau\omega_e}\sin \omega_e t\right) \quad (\text{underdamped}) \\ F_v(t) &= \frac{A}{\tau_\text{L}-\tau_\text{S}} \left( \frac{1}{\tau_\text{S}} \text{e}^{-t / \tau_\text{S}} - \frac{1}{\tau_\text{L}} \text{e}^{-t / \tau_\text{L}} \right ) \quad (\text{overdamped})\end{split}\]

The analysis can also be carried out in the frequency domain by computing the power spectral density (PSD) functions, which are Fourier transforms of the corresponding ACFs. The position PSD of a DHO is given by

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

which exhibits a peak at \(\omega_\text{max} = \sqrt{\omega_0^2 - \Gamma^2/2}\).

The PSD of the velocity, \(S_v(\omega) = \omega^2 S(\omega)\), is

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

The velocity PSD has a peak at \(\omega_\text{max} = \omega_0\). The amplitude \(A\) can be related to the variance (strength) \(\sigma^2=\langle f^2 \rangle\) of the noise by \(2\Gamma\omega_0^2 A=\sigma^2\). Taking the strength to be equal to \(\sigma^2=2\Gamma k_B T\) the mean squared displacement becomes \(\langle x^2 \rangle=\frac{k_B T}{\omega_0^2}\) and mean squared velocity is given by \(\langle v^2 \rangle=k_B T\).

The position and velocity ACF and PSD solutions for the DHO with natural angular frequency \(\omega_0 = 1\) for a few different values of damping constant \(\Gamma\). The red line corresponds to the critical damping. For the position ACF a critically damped oscillator goes to zero without any oscillation. When the damping is larger than \(\sqrt{2}\) the peak for the corresponding PSD fixates at zero.