2.1. Diffusion as Random Walk
Self-diffusion usually describes the displacement of a test particle immersed in a medium with no external gradients. A canonical example is the Brownian motion, representing a random motion of macroscopic particles suspended in a liquid or a gas. Here, we are interested in atomic scales and, hence, consider displacements of a labeled atom in a fluid of unlabeled, but otherwise identical, atoms. If this motion can be considered a random walk process, then the diffusion coefficient in three spatial dimensions can be defined as [
2]
where
r is an actual (variable) length of the random walk,
is the time scale, and we focus on sufficiently long times (
). Consider first an ideal gas as an appropriate example. The atoms move freely between pairwise collisions. If the distribution of free paths between collisions follows the
scaling, then
and
, where
is the mean free path [
2]. Combining this with the relation for the average atom velocity
, we recover the elementary kinetic formula for the diffusion coefficient of an ideal gas
The dynamical picture is very different in liquids and this simple consideration clearly does not apply. The very concept of random walk, however, remains relevant, although characteristic length and time scales associated with a random walk process in liquids are very different from those in gases.
Below, the model proposed by Zwanzig [
1] to describe relations between the self-diffusion and shear viscosity coefficients of liquids, is discussed in some detail. In doing so, we naturally repeat some arguments and formulas from Zwanzig’s original work and later publications (for instance, from a recent paper by the present author [
3]). The emphasis is, however, not on the Stokes–Einstein relation per se, but rather on the possibility of presenting self-diffusion as a random walk process, and on defining the associated length and time scales. The emerging picture represents an alternative to the often assumed hop** mechanism of diffusion in the liquid state.
Zwanzig’s approach is based on the assumption that atoms in liquids exhibit solid-like oscillations about temporary equilibrium positions corresponding to a local minimum on the system’s potential energy surface [
2,
4]. These positions do not form a regular lattice like in crystalline solids. They are also not fixed, and change (or drift) with time (this is why liquids can flow), but on much longer time scales. Local configurations of atoms are preserved for some time until a fluctuation in the kinetic energy allows rearranging the positions of some of the atoms towards a new local minimum in the multidimensional potential energy surface. The waiting time distribution of the rearrangements scales as
, where
is a lifetime. Atomic motions after the rearrangements are uncorrelated with motions before rearrangements [
1].
Within this ansatz, a simplest reasonable approximation for the velocity autocorrelation function of an atom
j is
corresponding to a time dependence of a damped harmonic oscillator. Here,
T is the temperature in energy units,
m is the atomic mass, and
is an effective vibrational frequency. The self-diffusion coefficient
D is given by the Green–Kubo formula
Zwanzig then assumed that vibrational frequencies
are related to the collective mode spectrum and performs averaging over collective modes. After the evaluation of the time integral, this yields
where the summation runs over
normal mode frequencies. The dynamical picture used by Zwanzig makes sense only if the waiting time
is much longer than the inverse characteristic frequency of the solid-like oscillations. In this case, we can rewrite Equation (
5) as
where the conventional definition of averaging,
has been used.
Equation (
6) allows for a simple physical interpretation. It represents a diffusion coefficient for a random walk process, Equation (
1). The length scale of this process is identified as
which is twice the mean-square displacement of an atom from its local equilibrium position due to solid-like vibrations [
5]. The coefficient of two appears, because the initial atom position is not at the local equilibrium, but randomly distributed with the same properties as the final one (after the waiting time
). The characteristic time scale of the random walk process is just the waiting time
. Moreover, this waiting time should be associated with the Maxwellian shear relaxation time [
2,
6]
where
is the shear viscosity coefficient,
is the infinite frequency (instantaneous) shear modulus,
n is the density, and
is the transverse sound velocity.
Thus, self-diffusion in the liquid state can be viewed as a random walk due to atomic vibrations around temporary equilibrium positions over time scales associated with rearrangements of these equilibrium positions. In this paradigm, consecutive changes of temporary equilibrium positions (jumps of liquid configurations between two neighboring local minima of the multidimensional potential energy surface in Zwanzig’s terminology) are relatively small, much smaller than the vibrational amplitude. Hop** events with displacement amplitudes of the order of interatomic separation may be present, but they are relatively rare and do not contribute to the diffusion process. This picture is very different from the widely accepted hop** mechanism of self-diffusion in liquids. Previously, the concept of random walk was suggested in the context of molecular and atomic motion in water and liquid argon [
7]. Here, we provide a more quantitative basis for this treatment.
Substituting Equation (
8) into Equation (
6), we obtain a relation between the self-diffusion and viscosity coefficients in the form of the Stokes–Einstein (SE) relation,
where
is the mean interatomic separation and
is the SE coefficient.
Formula (
9) particularly emphasizes the relation between the liquid transport and collective mode properties. Since the exact distribution of frequencies is generally not available, Zwanzig originally used a Debye approximation, characterized by one longitudinal and two transverse modes with acoustic dispersion. The sum over frequencies can be converted to an integral over
k using the standard procedure
, where
V is the volume. This yields
where the cutoff
is chosen to provide
n modes in each branch of the spectrum. This ensures that the averaging procedure applied to a quantity that does not depend on
k does not change its value. Substituting
and
into Equation (
10) we arrive at
This essentially coincides with Zwanzig’s original result, except he expressed the SE coefficient in terms of the longitudinal and shear viscosity
. The equivalence was pointed out in Reference [
6]. Note that since the sound velocity ratio
is confined in the range from 0 to
, the coefficient
can vary only between ≃0.13 and ≃0.18 [
1,
6]. Possible relations between the viscosity and thermal conductivity coefficients of dense fluids that can complement the SE relations of Equations (
9) and (
11) have been discussed recently [
8].
An important time scale of a liquid state is a structure relaxation time. This can be defined as an average time it takes an atom to move the average interatomic distance
(sometimes it is referred to as the Frenkel relaxation time [
9,
10,
11]). Taking into account diffusive atomic motions, we can write
. From Equation (
1), we immediately get
This implies that
. The time scale ratio
has a maximum at melting conditions, where, according to the Lindemann melting criterion
[
5,
12]. This picture is consistent with the results from numerical simulations (see, e.g., Figure 3 from Reference [
11]). Thus, there is a huge separation between the structure relaxation and individual atom dynamical relaxation time scales.
2.2. Relation to Collective Modes Properties
Despite the simplifications involved, the predictive power of Zwanzig’s model is quite impressive. Although the model does not allow making independent theoretical predictions of viscosity and self-diffusion coefficients, its prediction of the product, in the form of Equation (
9), is highly accurate in some vicinity of the liquid–solid phase transition of many simple liquids [
6,
13,
14]. Moreover, the coefficient
can be correlated with the potential softness (via the ratio of the sound velocities), as the model predicts [
6]. Some of the assumptions, such as the effect of the waiting time distribution, were critically examined in Reference [
15]. In particular, it was demonstrated that the SE relation of the form (
9) is not obeyed if the distribution of waiting times is not exponential. In this section, we address another interesting question: how sensitive is the value of
to the assumptions about liquid collective mode properties?
To be specific, we consider a model one-component plasma (OCP) system. The OCP fluid is chosen for the following three main reasons: (i) vibrational (caging) motion is most pronounced due to extremely soft and long-ranged character of the interaction potential [
16,
17]; (ii) Zwanzig’s original derivation is not directly applicable to the OCP case, because the longitudinal mode is not acoustic (but plasmon) and, thus, it is a good opportunity to examine how the model should be modified in this case; (iii) collective modes in the OCP system are well studied and understood (for example, simple analytical expressions for the long-wavelength dispersion relations are available, see
Appendix A).
The OCP model is an idealized system of mobile point charges immersed in a neutralizing fixed background of opposite charge (e.g., ions in the immobile background of electrons or vice versa) [
18,
19,
20,
21,
22,
23,
24]. From the fundamental point of view, OCP is characterized by a very soft and long-ranged Coulomb interaction potential,
, where
q is the electric charge. The particle–particle correlations and thermodynamic properties of the OCP are characterized by a single dimensionless coupling parameter
, where
is the Wigner–Seitz radius. At
, the OCP is strongly coupled, and this is where it exhibits properties characteristic of a fluid phase (a body centered cubic phase becomes thermodynamically stable at
, as the comparison of fluid and solid Helmholtz free energies predicts [
22,
25,
26]). Dynamical scales of the OCP are usually expressed by the plasma frequency
. For example, the Einstein frequency is
. The transverse sound velocity at strong coupling is
[
27].
From extensive molecular dynamics simulations, it is known that the SE relation is satisfied to a very high accuracy in a strongly coupled OCP fluid with
at
[
14,
28,
29].
Figure 1 demonstrates that
approaches the strongly coupled asymptote already at
. Note that the OCP value of the SE coefficient is not truly universal, but rather representative for soft long-ranged pairwise interactions, in which case the transverse-to-longitudinal sound velocity ratio is small [see Equation (
11)]. For example, the same value (≃0.14) is reached in weakly screened Coulomb (Yukawa) fluids, while for Lennard-Jones fluids it increases to
and further to
in hard-sphere fluids [
14].
Now, we examine the sensitivity of the theoretical value of the SE coefficient
to concrete assumptions about the collective excitation spectrum. We start with the simplest approximation that all atoms are oscillating with the same Einstein frequency
(known as the Einstein model in the solid state physics). This approximation results in
, which is too low compared to the actual value from MD simulations (see
Figure 1).
As a next level of approximation a Debye-like vibrational density of states (VDOS), is assumed (averaging is performed using a standard definition ). Using the cutoff Debye frequency and requesting that we arrive at . This yields , which is somewhat better, but still considerably smaller than the actual result.
The most accurate theoretical estimate would be obtained if the exact VDOS were known. However, this is not the case. Nevertheless, accurate knowledge of the real dispersion relations for the longitudinal and transverse modes can be already quite useful. We make use of simple expressions based on the quasi-localized charge approximation (QLCA) [
30] combined with the excluded cavity model for the radial distribution function [
27]. The corresponding expressions for
and
are provided in the
Appendix A. Substituting these in Equation (
10), we have obtained
and
. This is very close to the exact result from MD simulations, as expected. Note that the exact result
is reproduced by construction.
The last demonstration uses a heuristic VDOS of the form
which reproduces the Debye model at low frequencies and implements the Gaussian cutoff at high
. This form was inspired by the observation that the functional form
can fit the numerically obtained VDOS of Lennard-Jones liquids reasonably well [
31]. We just substituted the linear scaling at low frequencies with the quadratic one to make the integral converging. This is clearly not a valid physical argument, but we use it here merely for illustrative purposes. The two normalization conditions yield
Application of this VDOS results in , which almost coincides with the exact MD result. Thus, implementation of the Gaussian cutoff to the Debye-like VDOS improves the situation considerably.
It should be noted that very long wavelengths and low frequency parts of the spectra are not relevant for the present consideration, because dynamics on time scales shorter than the relaxation time
is considered. However, since
needs to be satisfied, this corresponds to only a small part of the entire spectrum, and we therefore included low frequencies for simplicity, similar to what Zwanzig did originally [
1]. This also allows us to disregard the effects associated with the
k-gap in the dispersion relation of the transverse mode, an important property of liquid dynamics [
32,
33,
34,
35,
36,
37].