Theoretical background

Definitions of correlation functions

Positional correlation functions

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

The density of particles is denoted \(\rho(\boldsymbol{r},t)\) and defined as

\[\rho(\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<\rho(\boldsymbol{q},t)\rho(-\boldsymbol{q},0)\right> \quad= \quad \frac{1}{N}\sum _n ^N \sum _m ^N \left< \mathrm{exp} \left[i\boldsymbol{q}\cdot(\boldsymbol{r}_n(t)-\boldsymbol{r}_m(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 scattering function, i.e. \(n=m\), is defined as

\[F_s(\boldsymbol{q},t)=\frac{1}{N} \sum _n ^N \left < \mathrm{exp}\left[i\boldsymbol{q}\cdot(\boldsymbol{r}_n(t)-\boldsymbol{r}_n(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 structure factor is given by

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

whereas the dynamical 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 and thus the anharmonicity of the system. Computing \(S(\boldsymbol{q},\omega)\) from MD has the advantage of fully including anharmonic effects and allowing one to easily 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

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

where \(\boldsymbol{v}_i(t)\) is the velocity of particle \(i\) at time \(t\). The current correlation function \(J_{\alpha \beta}(\boldsymbol{q},t)\) is defined as

\[J_{\alpha \beta}(\boldsymbol{q},t) = \frac{1}{N} \left < j_\alpha(-\boldsymbol{q},0) j_\beta(-\boldsymbol{q},t) \right >\]

where \(\alpha\) and \(\beta\) refer to Cartesian coordinates. The current correlation can be split into a longitudinal and a transverse part, \(J_l(\boldsymbol{q},t)\) and \(J_t(\boldsymbol{q},t)\), respectivley. As in the case of \(F(\boldsymbol{q},t)\) the current correlation function can be Fourier transformed to obtain \(J_{\alpha \beta}(\boldsymbol{q},\omega)\). The longitudinal current correlation \(J_l(\boldsymbol{q},\omega)\) is related to the dynamical structure factor by

\[J_l(\boldsymbol{q},\omega) = \frac{\omega^2}{q^2} S(\boldsymbol{q},\omega),\]

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.

In multi-component systems one can introduce partial correlations functions, which enable one to separate 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 furtheremo frequently desirable to determine the above mentioned quantites along specific paths between high symmetry \(vec{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

If we assume that the particle density \(\rho(\boldsymbol{q},t)\) oscillates as a damped harmonic oscillator we obtain

\[\ddot{\rho}(\boldsymbol{q},t) + \Gamma \dot{\rho}(\boldsymbol{q},t) + \omega_0^2{\rho}(\boldsymbol{q},t)=f(\boldsymbol{q},t),\]

which yields

\[\frac{\mathrm{d}^2 }{\mathrm{d} t^2}F(\boldsymbol{q},t) + \Gamma \frac{\mathrm{d} }{\mathrm{d} t}F(\boldsymbol{q},t) + \omega_0^2 F(\boldsymbol{q},t) = 0\]

This equation in turn describes a damped harmonic oscillator with frequency \(\omega_0\) and damping factor \(\Gamma\). Assuming weak damping and \(\frac{\mathrm{d} }{\mathrm{d} t}F(\boldsymbol{q},t=0) = 0\) the solution can be found by Laplace transformation

\[F(\boldsymbol{q},t) = S(\boldsymbol{q}) \mathrm{e}^{-\frac{\Gamma}{2} |t|} \Big [\cos\Big(\sqrt{\omega_0^2-\frac{\Gamma^2}{4}}t \Big) + \frac{\Gamma}{2\sqrt{\omega_0^2-\frac{\Gamma^2}{4}}}\sin \Big(\sqrt{\omega_0^2-\frac{\Gamma^2}{4}}|t| \Big) \Big],\]

where the absolute value of the time is used to make the function even in time. While the structure factor \(S(\boldsymbol{q})\) is known for simplicity when fitting we can use

\[F(\boldsymbol{q},t) = A\mathrm{e}^{-B|t|}\Big(cos(\omega_et)+\frac{B}{\omega_e}\sin(\omega_e|t|)\Big),\]

where \(A\), \(B\) and \(\omega_e\) are fitting parameters. They are related to the original parameters by

\[\omega_0 = \sqrt{\omega_e^2 + \frac{\Gamma^2}{4}}\]
\[\Gamma = 2B\]
\[\tau = 2\pi\frac{2}{\Gamma},\]

where \(\tau\) is the lifetime of the oscillation with frequency \(\omega_0\). The factor of \(2\pi\) comes from the fact that \(\Gamma\) has units rad/s

One can also carry out the fitting in the frequency domain. The fitting function for the dynamical structure factor is given by the Fourier transform of the above function for the intermediate scattering function

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

Again for simplicity, \(S(\boldsymbol{q})\) can be replaced with a constant for the purpose of fitting.

Velocity correlation functions

The relation between the dynamical structure factor and the longitudinal current is known in the frequency domain and we obtain

\[C_l(\boldsymbol{q},\omega) = \frac{\omega^2}{|\boldsymbol{q}|^2}S(\boldsymbol{q}) \frac{2\Gamma\omega_0^2}{(\omega^2-\omega_0^2)^2 + \Gamma^2 \omega^2}\]

Fourier transforming this gives us \(C_l(\boldsymbol{q},t)\) as

\[C_l(\boldsymbol{q},t) = 0\]

These functions can be fitted in a similarly to the intermediate scattering function and the dynamical structure factor. 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.


This last part needs to be checked and finished.


Check all expressions for harmonic oscillators FTs in Göran