Theoretical background

Definitions of correlation functions

Positional correlation functions

In this section the correlation functions computed by dynasor are introduced. For a more detailed description and derivations of all relations see [JPB80].

The density of particles 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 in terms of the Fourier transform of the particle density

\[F(\boldsymbol{q},t)=\frac{1}{N}\left<n(\boldsymbol{q},t)n(-\boldsymbol{q},0)\right> \quad= \quad \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.)

The self part of the intermediate scattering function, i.e. \(i=j\), is defined as

\[F_s(\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. diffusion. \(F_s(\boldsymbol{q},t)\) is commonly referred to as the incoherent 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{dt}\]

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 broadening of these peaks are related to the phonon lifetimes. Computing \(S(\boldsymbol{q},\omega)\) from MD 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 rather 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 and liquids

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_{AA}(\boldsymbol{q},t)\), \(F_{BB}(\boldsymbol{q},t)\), \(F_{AB}(\boldsymbol{q},t)\) and so on for all correlation functions introduced above. In some cases, instead of analyzing the partial functions directly, it is interesting to consider linear combinations of them. This will be discussed in more depth in the liquid sodium chloride example.

In solids, it is furthermore 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.

Analyzing correlation functions

The correlation functions can be convoluted with atomic form factors in order to obtain predictions for various types of neutron and X-ray scattering experiments. In addition, they contain crucial information about the structure and dynamics of the system. In this section, we provide a concise description of some of the information that can be obtained in this fashion.

Positional correlation functions

The correlation functions can be fitted to analytical expressions for damped harmonic oscillators in order to extract frequencies and phonon lifetimes. While many different analytical forms can be considered, here we use functions with a continous derivative at \(t=0\). For each \(\boldsymbol{q}\)-point the intermediate scattering function can be expressed as

\[\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)\,\, , \,\, \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)\,\, , \,\,\omega_0 < \frac{\Gamma}{2}\end{split}\]

where the fitting parameters are \(\omega_0\), \(\Gamma\) and \(A\). The eigenfrequency of the damped oscillator is \(\omega_e = \sqrt{\omega_0^2 - \frac{\Gamma^2}{4}}\). The fitting can also be carried out in the frequency domain, meaning 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}\]

This analysis can be extended to the current correlations by considering the relation between the longitudinal current 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) , \,\, \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)\,\, , \,\,\omega_0 < \frac{\Gamma}{2}\end{split}\]

The phonon lifetime, \(\tau\), is given by

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

or in energy as

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

Here \(\hbar\) is used due to the fact that \(\Gamma\) has units rad/s.

The transverse part of the current correlation is assumed to have the same functional form. However, sometimes the transverse correlation function is made up by two oscillations, not one.