Self-Diffusion in Simple Liquids as a Random Walk Process

Abstract

It is demonstrated that self-diffusion in dense liquids can be considered a random walk process; its characteristic length and time scales are identified. This represents an alternative to the often assumed hopping mechanism of diffusion in the liquid state. The approach is illustrated using the one-component plasma model.

1. Introduction

About 40 years ago, Robert Zwanzig published an influential paper on the relation between self-diffusion and viscosity of liquids (Stokes–Einstein relation) [1]. The purpose of the present paper is to demonstrate that the dynamical picture behind Zwanzig’s result is equivalent to a random walk process, with well defined length and time scales. It is also demonstrated that a theoretical prediction for the numerical factor relating the self-diffusion and viscosity coefficients, in the form of the Stokes–Einstein relation, is quite sensitive to concrete assumptions about the liquid collective mode spectrum. The results provide a consistent picture of the diffusion mechanism in dense liquids with soft isotropic pairwise interactions.

2. Results

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]

$$D=\frac{1}{6}\frac{\langle r^2\rangle }{\tau },$$

(1)

where r is an actual (variable) length of the random walk, $\tau$ is the time scale, and we focus on sufficiently long times ($t\gg \tau$). 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 $e^{-r/\lambda }/\lambda$ scaling, then $\langle r\rangle =\lambda$ and $\langle r^2\rangle =2{\lambda }^2$, where $\lambda$ is the mean free path [2]. Combining this with the relation for the average atom velocity $\langle v\rangle =\lambda /\tau$, we recover the elementary kinetic formula for the diffusion coefficient of an ideal gas

$$D=\frac{1}{3}\langle v\rangle \lambda .$$

(2)

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 hopping 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 $exp(-t/\tau )/\tau$, where $\tau$ 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

$$Z_j\left(t\right)\simeq \left(\frac{T}{m}\right)cos\left({\omega }_{j}t\right)e^{-t/\tau },$$

(3)

corresponding to a time dependence of a damped harmonic oscillator. Here, T is the temperature in energy units, m is the atomic mass, and ${\omega }_j$ is an effective vibrational frequency. The self-diffusion coefficient D is given by the Green–Kubo formula

$$D=\frac{1}{N}{\int }_0^{\infty }\sum _j{Z}_j\left(t\right)dt.$$

(4)

Zwanzig then assumed that vibrational frequencies ${\omega }_j$ are related to the collective mode spectrum and performs averaging over collective modes. After the evaluation of the time integral, this yields

$$D=\frac{T}{3mN}\sum _{\mathbf{k}}\frac{\tau }{1+{\omega }_{\mathbf{k}}^2{\tau }^2},$$

(5)

where the summation runs over $3N$ normal mode frequencies. The dynamical picture used by Zwanzig makes sense only if the waiting time $\tau$ is much longer than the inverse characteristic frequency of the solid-like oscillations. In this case, we can rewrite Equation (5) as

$$D=\frac{T}{m\tau }〈\frac{1}{{\omega }^2}〉,$$

(6)

where the conventional definition of averaging, $\langle {\omega }^{-2}\rangle =(1/3N){\sum }_{\mathrm{k}}{\omega }_{\mathbf{k}}^{-2}$ 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

$$\langle r^2\rangle =\frac{6T}{m}〈\frac{1}{{\omega }^2}〉,$$

(7)

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 $\tau$). The characteristic time scale of the random walk process is just the waiting time $\tau$. Moreover, this waiting time should be associated with the Maxwellian shear relaxation time [2,6]

$${\tau }_{\mathrm{M}}=\frac{\eta }{G_{\infty }}=\frac{\eta }{mn{c}_t^2},$$

(8)

where $\eta$ is the shear viscosity coefficient, $G_{\infty }$ is the infinite frequency (instantaneous) shear modulus, n is the density, and $c_t$ 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. Hopping 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 hopping 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,

$$D\eta \left(\frac{\mathrm{\Delta }}{T}\right)=\frac{c_t^2}{{\mathrm{\Delta }}^2}〈\frac{1}{{\omega }^2}〉={\alpha }_{\mathrm{SE}},$$

(9)

where $\mathrm{\Delta }=n^{-1/3}$ is the mean interatomic separation and ${\alpha }_{\mathrm{SE}}$ 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 ${\sum }_{\mathbf{k}}\to V\int d\mathbf{k}/{\left(2\pi \right)}^3$, where V is the volume. This yields

$$〈\frac{1}{{\omega }^2}〉=\frac{1}{6{\pi }^{2}n}{\int }_0^{k_{\max }}{k}^{2}dk\left(\frac{1}{{\omega }_l^2}+\frac{2}{{\omega }_t^2}\right),$$

(10)

where the cutoff $k_{\max }={\left(6{\pi }^{2}n\right)}^{1/3}$ 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 ${\omega }_l=c_{l}k$ and ${\omega }_t=c_{t}k$ into Equation (10) we arrive at

$${\alpha }_{\mathrm{SE}}=\frac{2}{{\left(6{\pi }^2\right)}^{2/3}}\left(1+\frac{c_t^2}{2c_l^2}\right)\simeq 0.13\left(1+\frac{c_t^2}{2c_l^2}\right).$$

(11)

This essentially coincides with Zwanzig’s original result, except he expressed the SE coefficient in terms of the longitudinal and shear viscosity ${\alpha }_{\mathrm{SE}}\simeq 0.13(1+\eta /2{\eta }_l)$. The equivalence was pointed out in Reference [6]. Note that since the sound velocity ratio $c_t/c_l$ is confined in the range from 0 to $\sqrt{3}/2$, the coefficient ${\alpha }_{\mathrm{SE}}$ 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 $\mathrm{\Delta }$ (sometimes it is referred to as the Frenkel relaxation time [9,10,11]). Taking into account diffusive atomic motions, we can write ${\tau }_{\mathrm{R}}={\mathrm{\Delta }}^2/6D$. From Equation (1), we immediately get

$${\tau }_{\mathrm{R}}=\frac{{\mathrm{\Delta }}^2}{\langle r^2\rangle }{\tau }_{\mathrm{M}}.$$

(12)

This implies that ${\tau }_{\mathrm{R}}/{\tau }_{\mathrm{M}}\gg 1$. The time scale ratio ${\tau }_{\mathrm{R}}/{\tau }_{\mathrm{M}}$ has a maximum at melting conditions, where, according to the Lindemann melting criterion ${\mathrm{\Delta }}^2/\langle r^2\rangle \sim 100$ [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 ${\alpha }_{\mathrm{SE}}$ 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 ${\alpha }_{\mathrm{SE}}$ 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, $\varphi \left(r\right)=q^2/r$, where q is the electric charge. The particle–particle correlations and thermodynamic properties of the OCP are characterized by a single dimensionless coupling parameter $\mathsf{\Gamma }=q^2/aT$, where $a={(4\pi n/3)}^{-1/3}$ is the Wigner–Seitz radius. At $\mathsf{\Gamma }\gtrsim 1$, 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 $\mathsf{\Gamma }\gtrsim 174$, 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 ${\omega }_{\mathrm{p}}=\sqrt{4\pi q^{2}n/m}$. For example, the Einstein frequency is ${\mathrm{\Omega }}_{\mathrm{E}}^2\equiv \langle {\omega }^2\rangle ={\omega }_{\mathrm{p}}^2/3$. The transverse sound velocity at strong coupling is $c_t^2=(3/100\pi ){(4\pi /3)}^{1/3}{\omega }_{\mathrm{p}}^2{\mathrm{\Delta }}^2\simeq 0.015{\omega }_{\mathrm{p}}^2{\mathrm{\Delta }}^2$ [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 ${\alpha }_{\mathrm{SE}}\simeq 0.14\pm 0.01$ at $\mathsf{\Gamma }\gtrsim 50$ [14,28,29]. Figure 1 demonstrates that ${\alpha }_{\mathrm{SE}}$ approaches the strongly coupled asymptote already at $\mathsf{\Gamma }\simeq 10$. 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 ${\alpha }_{\mathrm{SE}}\simeq 0.15$ and further to ${\alpha }_{\mathrm{SE}}\simeq 0.17$ in hard-sphere fluids [14].

Now, we examine the sensitivity of the theoretical value of the SE coefficient ${\alpha }_{\mathrm{SE}}$ to concrete assumptions about the collective excitation spectrum. We start with the simplest approximation that all atoms are oscillating with the same Einstein frequency ${\mathrm{\Omega }}_{\mathrm{E}}$ (known as the Einstein model in the solid state physics). This approximation results in ${\alpha }_{\mathrm{SE}}\simeq 0.046$, 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), $g\left(\omega \right)\propto {\omega }^2$ is assumed (averaging is performed using a standard definition $\langle {\omega }^{\ell }\rangle =\int {\omega }^{\ell }g\left(\omega \right)d\omega /(\int g\left(\omega \right)d\omega )$). Using the cutoff Debye frequency ${\omega }_{\mathrm{D}}$ and requesting that $\langle {\omega }^2\rangle ={\mathrm{\Omega }}_{\mathrm{E}}^2$ we arrive at $\langle {\omega }^{-2}\rangle =9/5{\mathrm{\Omega }}_{\mathrm{E}}^2$. This yields ${\alpha }_{\mathrm{SE}}\simeq 0.083$, 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 ${\omega }_l\left(k\right)$ and ${\omega }_t\left(k\right)$ are provided in the Appendix A. Substituting these in Equation (10), we have obtained $\langle {\omega }_{\mathrm{p}}^2/{\omega }^2\rangle \simeq 9.76$ and ${\alpha }_{\mathrm{SE}}\simeq 0.150$. This is very close to the exact result from MD simulations, as expected. Note that the exact result $\langle {\omega }^2/{\omega }_{\mathrm{p}}^2\rangle =1/3$ is reproduced by construction.

The last demonstration uses a heuristic VDOS of the form

$$g\left(\omega \right)=\mathcal{A}{\omega }^{2}exp(-\mathcal{B}{\omega }^2),$$

(13)

which reproduces the Debye model at low frequencies and implements the Gaussian cutoff at high $\omega$. This form was inspired by the observation that the functional form $g\left(\omega \right)=2\alpha \omega e^{-\alpha {\omega }^2}$ 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

$$\mathcal{A}=\frac{4}{\sqrt{\pi }}{\left(\frac{3}{2{\mathrm{\Omega }}_{\mathrm{E}}^2}\right)}^{3/2},\mathcal{B}=\frac{3}{2{\mathrm{\Omega }}_{\mathrm{E}}^2}.$$

(14)

Application of this VDOS results in ${\alpha }_{\mathrm{SE}}\simeq 0.139$, 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 ${\tau }_{\mathrm{M}}$ is considered. However, since $\omega {\tau }_{\mathrm{M}}\gg 1$ 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].

3. Discussion and Conclusions

While transport phenomena in gaseous and solid phases can be well described at the quantitative level, transport in liquids is still much less understood, even at the qualitative level. Here, we have demonstrated that self-diffusion in dense liquids can be described as a random walk process with well defined time and length scales. The length scale is related to the amplitude of solid-like vibrations around local temporary equilibrium positions. The time scale is set by the Maxwellian shear relaxation time. This dynamical picture results in the Stokes–Einstein relation between the coefficients of self-diffusion and viscosity, which is satisfied in many simple liquids. Importantly, the hoping mechanism of atomic diffusion in liquids is irrelevant in this picture of microscopic atomic dynamics.

The dynamical picture involved requires that the atomic motion be dominated by fast solid-like oscillations around the local equilibrium positions. This limits the model applicability to regions on the phase diagram located not too far from the liquid–solid phase transition (high densities and low temperatures). Additionally, it applies to sufficiently soft interaction potentials with pronounced oscillation dynamics. In the hard sphere interaction limit, this model is clearly inadequate (although SE relation is still satisfied even in this limit [14]).

Finally, we have demonstrated that a theoretically obtained numerical factor in the SE relation is sensitive to concrete assumptions about the liquid collective modes properties. This highlights the necessity of accurate knowledge of the vibrational density of states and dispersion relations in liquids.