Transport dynamics of complex fluids

Transport Equation of Complex Fluids

Here, to address these issues, we present a transport equation that provides a quantitative description of thermal motion for various complex fluids. This equation can be derived by considering the continuum limit of a random walk model with a general sojourn time distribution, ψΓ(t), coupled to arbitrary hidden environmental variables Γ (Fig. 1); Γ designates the entire set of dynamical variables that affect transport dynamics in disordered fluids. For this model, the joint probability, p(m,Γ,t), that a random walker is located at the mth site and the environmental state is at Γ at time t can be written as (49, 59)p(m,Γ,t)=∑N=0∞p(m|N)PN(Γ,t).[1]In Eq. 1, p(m|N) denotes the conditional probability that the random walker is located at the mth site, given that it has undergone N jumps, which is well known: p(m|N)=(N!/m+!m−!)2−N with m±=(N±m)/2 (N≥|m|), and p(m|N)=0 (N<|m|) (60). On the other hand, PN(Γ,t) denotes the joint probability that the total number of jumps made by the random walker is N and the environmental state is at Γ at time t. PN(Γ,t) is the crucial factor determining the time dependence of p(m,Γ,t). Using a generalized version of Sung and Silbey’s (61) master equation, which provides a formally exact description of the time evolution of PN(Γ,t), we derive the following transport equation from Eq. 1 in the continuum limit (SI Appendix, Text S1):p˙^(r,Γ,s)=D^Γ(s)∇2p^(r,Γ,s)+L(Γ)p^(r,Γ,s).[2]Here, p˙^(r,Γ,s) and p^(r,Γ,s) denote the Laplace transform of ∂p(r,Γ,t)/∂t and p(r,Γ,t), respectively; p(r,Γ,t) denotes the joint probability density that a particle is located at position r and the hidden environment is at state Γ at time t. This joint probability density satisfies the following normalization condition: ∫dr∫dΓp(r,Γ,t)=1. Throughout this work, f^(s) denotes the Laplace transform of f(t), i.e., ∫0∞dte−stf(t). In Eq. 2, D^Γ(s) designates the diffusion kernel that is determined by the environmental state-dependent sojourn time distribution, ψΓ(t) of our random walk model; i.e., D^Γ(s)= limε→0(ε2/2d)sψ^Γ(s)/[1−ψ^Γ(s)] =limε→0ε2κ^Γ(s)/2d with ε and d denoting the lattice constant and the spatial dimension, respectively. κ^Γ(s)(≡sψ^Γ(s)/[1−ψ^Γ(s)]) denotes the jump rate kernel of the random walker, which is dependent on lattice constant ε; for the continuum limit description, we assume limε→0ε2κ^Γ(s) exists. In Eq. 2, L(Γ) designates a mathematical operator describing the dynamics of the hidden environmental variables Γ.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Schematic representation of our random walk model with environmental state-dependent dynamics. In our model, the sojourn time distribution ψΓ(t) of a random walker is dependent on environmental state variables, Γ. The probabilistic dynamics of this random walker model can be described by Eq. 2 in the continuum limit.

For the sake of generality, we do not assume a particular model of environmental state dynamics, nor do we assume a particular form of mathematical operator L(Γ) at this point. A correct mathematical form of L(Γ) is dependent on the environment surrounding the system in question; when environmental state dynamics are a non-Markov process, L(Γ) may be dependent on Laplace variable s. As demonstrated in this work, quantitative information about transport dynamics coupled to hidden environmental variables can be extracted from simultaneous analysis of the MSD and NGP time profiles. This information can then be used to construct a more explicit model of transport dynamics of complex fluidic systems.

Eq. 2 encompasses the CTRW model and the diffusing diffusivity models (45⇓⇓⇓⇓⇓–51) (see Discussion for more details). A further generalization of Eq. 2 for complex fluidic systems under a spatially heterogeneous external potential is presented in SI Appendix, Text S2.

Analytic Expressions of the Moments

From Eq. 2, we obtain the exact analytical expressions of the first two nonvanishing moments, 〈|r(t)−r(0)|2〉(≡Δ2(t)) and 〈|r(t)−r(0)|4〉(≡Δ4(t)), of the displacement distribution:Δ^2(s)=2ds2〈D^Γ(s)〉,[3a]Δ^4(s)=(1+2d)2sΔ^2(s)2[1+sC^D(s)].[3b]In obtaining Eqs. 3a and 3b, we assume that the hidden environment is initially in a stationary state, such as the equilibrium state or the nonequilibrium steady state (see SI Appendix, Text S3 for details). The bracket notation, 〈⋯〉, designates the average over the stationary initial distribution of the environmental state. We can extend Eqs. 3a and 3b to the case where the initial state of the hidden environment is a nonstationary state, as is the case for glass, for which the analytic expressions of Δ^2(s) and Δ^4(s) are more complicated than Eqs. 3a and 3b (SI Appendix, Text S4). However, transport dynamics of various complex fluids can be quantitatively explained by Eqs. 3a and 3b, as demonstrated in this work. Hereafter, we focus on quantitative analysis of transport dynamics of complex fluidic systems with use of Eqs. 3a and 3b. We leave a quantitative explanation of transport dynamics in glass to future research.

The mean diffusion kernel, 〈DΓ(t)〉, in Eq. 3a is nothing but the two-point VAF; i.e., 〈DΓ(t)〉=〈v(t)⋅v(0)〉/d with v(t) being the velocity vector. This can be seen by comparing Eq. 3a and the Laplace transform of the well-known relation, Δ2(t)=2∫0tdτ2∫0τ2dτ1〈v(τ2−τ1)⋅v(0)〉(62), exact as long as 〈v(τ2)⋅v(τ1)〉 =〈v(τ2−τ1)⋅v(0)〉. Knowing this and utilizing the Tauberian theorem, we obtain lims→∞s〈D^Γ(s)〉=limt→0〈v(t)⋅v(0)〉/d=〈|v|2〉/d =kBT/M, with kBT and M denoting thermal energy and the mass of our tracer particle, respectively. This means that 〈D^Γ(s)〉 is proportional to the mean-square velocity in the large-s limit; i.e., 〈D^Γ(s)〉≅s−1〈|v|2〉/d (s→∞). On the other hand, in the small-s limit, the value of 〈D^Γ(s)〉 approaches ∫0∞dt〈v(t)·v(0)〉/d, which is simply the diffusion constant, D¯, according to the Green–Kubo relation (63). Substituting the small (large)-s limit asymptotic behavior of 〈D^Γ(s)〉 into Eq. 3a, we recover the well-known asymptotic behavior of the MSD: d(kBT/M)t2 at short times and 2dD¯t at long times.

While the second moment, Δ2(t), is dependent only on 〈DΓ(t)〉, the fourth moment, Δ4(t), is dependent on the environment-coupled fluctuation of the diffusion kernel, DΓ(t). In Eq. 3b, C^D(s) is the Laplace transform of the diffusion kernel correlation (DKC) or TCF of the diffusion kernel fluctuation (see SI Appendix, Eq. S3-12 for the precise definition). At long times, where the MSD is linear in time, the diffusion kernel becomes the diffusion coefficient; i.e., D^Γ(s)≅D^Γ(0)(≡DΓ) so that CD(t) can be identified as the TCF of the diffusion coefficient fluctuation; i.e., CD(t)≅〈δD(t)δD(0)〉/〈D〉2 =ηD2ϕD(t), where ηq2 and ϕq(t)(q∈{v2,D}) designate the relative variance, 〈δq2〉/〈q〉2, and the normalized TCF of q, 〈δq(t)δq(0)〉/〈δq2〉, respectively. At short times, on the other hand, CD(t) can be identified as the TCF of squared speed v2(t)(≡|v(t)|2); i.e., CD(t)≅dηv22ϕv2(t)[=d〈δv2(t)δv2(0)〉/〈v2〉2] (SI Appendix, Text S5). Given that the initial speed distribution obeys the Maxwell–Boltzmann distribution, we obtain 〈v4〉=(1+2/d)〈v2〉2, and the initial value of CD(t) can then only be given by limt→0CD(t)=d((〈v4〉−〈v2〉2))/〈v2〉2 =2. We find this is true for supercooled water and dense hard disk fluids (Fig. 2D and SI Appendix, Fig. S1D).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Model-based quantitative analysis of the mean-square displacement (MSD) and non-Gaussian parameter (NGP) for the TIP4P/2005 water system. (A and B) Time profiles of the MSD and NGP at various temperatures. Open circles, simulation results; solid lines, the best fits of Eqs. 9 and 18 to the simulation results. (C) Mean diffusion kernel obtained from Eq. 10 with optimized parameter values. (D) Solid lines, diffusion kernel correlation, CD(t), extracted from the MSD and NGP data (Methods); dotted line, the mean-scaled TCF of squared speed of supercooled water at 193 K (SI Appendix, Fig. S11). In A, B, and D, the navy-blue triangles and the red squares represent the caging times, τc, and the NGP peak times, τng, respectively. In D, the yellow triangles represent the alpha relaxation times, τα (SI Appendix, Fig. S12). (E) Green circles, total disorder, limt→∞〈r2(t)〉α2(t)/σO2, with σO denoting an oxygen atom’s Lennard-Jones diameter, 3.1589 Å; yellow and cyan lines, intrinsic and extrinsic disorder (see text below Eq. 8); red circles, 4dC^D(0)/τD, where d, C^D(0), and τD(=σO2/D¯), respectively, denote the spatial dimension, the whole-time integration of CD(t), and the diffusion timescale. Tm and TW denote, respectively, the melting temperature (64) and the Widom line temperature (91, 92) (SI Appendix, Fig. S13).

In our theory, CD(t) is the essential dynamic quantity that characterizes environment-coupled mobility fluctuation. It serves as an ideal measure of mobility fluctuation for complex fluids exhibiting non-Fickian thermal motion, for which the TCF of the diffusion coefficient, the conventional measure of mobility fluctuation, is not well defined. As is shown below, there exists an exact relationship between CD(t) and the four-point and two-point VAFs, valid at all times and at any spatial dimension.

Using Eqs. 3a and 3b and the definition of the NGP, α2(t)[≡Δ4(t)/[(1+2/d)Δ2(t)2]−1], we can quantitatively explain the MSD and NGP of various complex fluids. From the MSD and NGP data, we can extract the time profiles of 〈DΓ(t)〉 and CD(t) either by assuming analytic functional forms for the MSD and NGP or without making any such assumption (Methods). In Fig. 2, we demonstrate our quantitative analysis of the molecular dynamics (MD) simulation results of the MSD and NGP for 4-point transferable intermolecular potential/2005 (TIP4P/2005) water (64), assuming specific analytic forms of the MSD and NGP but without assuming a particular model of the hidden environment or its influence on the diffusion kernel. As shown in this work, this information is useful in constructing an explicit model of transport dynamics of complex fluid systems; from this explicit model, we can predict or quantitatively understand the time dependence of the displacement distribution. Meanwhile, a model-free analysis of the MSD and NGP based on Eqs. 3a and 3b yields accurate and robust quantitative information about the time profiles of 〈DΓ(t)〉 and CD(t), but, on its own, is not physically interpretable (Fig. 3).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Model-free quantitative analysis of the mean-square displacement (MSD) and non-Gaussian parameter (NGP) for the TIP4P/2005 water system at 193 K. (A and B) Time profiles of the MSD and NGP: circles, simulation results; blue lines, best fits of Eqs. 9 and 18. (C) Mean diffusion kernel: circles, two-point VAF obtained from simulation; red line, second-order time derivative of MSD data; blue line, Eq. 10 with optimized parameters given in Table 1. (D) Time profile of diffusion kernel correlation extracted by analyzing the MSD and NGP data using (red line) model-free theory, (blue line) model-based theory, and (yellow dashed-dotted line) contribution of the unbound mode. (D, Inset) Time profile of CD(t) in linear timescale. The difference between CD(t) and the unbound-mode contribution results from the presence of the bound modes, accounted for by the second term on the R.H.S. of Eq. 12 (SI Appendix, Fig. S14).

Diffusion Kernel Correlation and Velocity Autocorrelation Functions

The diffusion kernel correlation, CD(t), is closely related to the two- and four-point VAFs through the BCF. It is known that the NGP, α2(t), or the fourth cumulant, X4(t)[=(1+2/d)Δ2(t)2α2(t) =Δ4(t)−(1+2/d)Δ2(t)2], of displacement is related to the BCF, β(t), by (54⇓–56)X4(t)=4!d∫0tdt1(t−t1)β(t1),[4]where β(t) is defined byβ(t1)=∫0t1dt2∫0t2dt3[〈vα(0)vα(t1)vα(t2)vα(t3)〉−〈vα(0)vα(t1)〉〈vα(t2)vα(t3)〉−〈vα(0)vα(t2)〉〈vα(t1)vα(t3)〉−〈vα(0)vα(t3)〉〈vα(t1)vα(t2)〉].[5]In Eq. 5, vα indicates a Cartesian component of velocity vector, v, and 〈vα(0)vα(t1)vα(t2)vα(t3)〉 and 〈vα(0)vα(t)〉[=(kBT/M)ϕvα(t)] denote the four- and two-point VAFs, respectively. According to Wick’s theorem, β(t) vanishes when vα(t) is a Gaussian process. Taking the Laplace transform of Eq. 4, we obtain X^4(s)=Δ^4(s)−(1+2/d)Δ^22(s)=4!dβ^(s)/s2, where Δ^22(s) denotes the Laplace transform of Δ2(t)2. Substituting Eqs. 3a and 3b into the latter equation, we obtain an exact relation of the DKC to the BCF and the normalized VAF (see SI Appendix, Text S6 for the derivation):C^D(s)〈D^Γ(s)〉2=∫0∞dte−st[32+dβ(t)+(kBTM∫0tdt′ϕvα(t′))2 +kB2T2M2∫0tdt1∫0t1dt2[ϕvα(t)−ϕvα(t−t1)]ϕvα(t2)].[6]Using Eq. 6 along with Eq. 5 and 〈D^Γ(s)〉=(kBT/M)ϕ^vα(s), we can calculate the time profile of CD(t) from the four- and two-point VAFs, as demonstrated in Fig. 4B for supercooled water at 193 K.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Microscopic measurement of the bound- and unbound-mode components of diffusion kernel correlation for the TIP4P/2005 water system at 193 K. (A) Burnett correlation function (BCF): red circles, β1(t) obtained from Eq. 5 and the four- and two-point VAFs obtained from simulation; blue circles, β2(t) obtained from β(t)=∂t2Χ4(t)/4!d and Χ4(t)=(1+2/d)Δ2(t)2α2(t). (Inset) The fourth cumulant, Χ4(t), of displacement. (B) Short-time diffusion kernel correlation, CD(t), obtained from (red circles) Eq. 6 and β1(t) and (blue circles) Eq. 6 and β2(t). (C) Long-time diffusion kernel correlation estimated by (squares) the mean-scaled TCF of diffusion coefficient fluctuation (Methods), (red dashed line) result of Eq. 14, and (Inset) time dependence of MSD scaled by 6t. The simulation results of the diffusion coefficient fluctuation are calculated using various bin times, ranging from 7 ns to 15 ns. In B and C, the black and yellow lines represent, respectively, the result of model-free theory and the contribution of the unbound mode to CD(t) shown in Fig. 3D. (D–G) Displacement distributions and representative time traces of three water molecules (D and E) at five different short times and (F and G) at three different long times. In E and G, the initial positions of all water molecules have been relocated to the center of the circle. In E, the dashed line represents a sphere centered at the initial position with radius Δ2(τc)1/2(≈0.6 Å).

Ergodicity and Long-Time Limit of Diffusion Kernel Correlation and Non-Gaussian Parameter

For fluidic systems showing terminal Fickian diffusion, CD(t) has the same long-time limit value as the NGP; i.e., CD(∞)=α2(∞), which can be shown by using Eqs. 3a and 3b and the definition of the NGP (SI Appendix, Text S7). This result indicates that the long-time limit NGP value, α2(∞), vanishes for ergodic fluid systems for which the TCF, CD(∞), of the diffusion coefficient vanishes in the long-time limit. However, for nonergodic systems with finite CD(∞), α2(∞) may not vanish, either. Therefore, α2(∞) can serve as an ergodicity measure for fluidic systems exhibiting long-time Fickian diffusion (SI Appendix, Fig. S2), which was noted by Odagaki (65) for the glass formation process. There exist transport systems with anomalous diffusion, Δ2(t)∝tv, and weak-ergodicity breaking. For such systems, α2(∞) deviates from CD(∞); even if CD(∞)=0, α2(∞) is finite and given by α2(∞)=vΓ(v)2/Γ(2v)−1 with Γ(z) denoting the Gamma function defined by Γ(z)=∫0∞dt tz−1e−t (SI Appendix, Text S7). This result was previously reported by Odagaki and Hiwatari (66) on the basis of the so-called coherent-medium approximation, which corresponds to a vanishingly small DKC; i.e., CD(t)=0. Finally, for nonergodic systems exhibiting anomalous diffusion, Δ2(t)∝tv, at long times, we find that the relationship between CD(∞) and α2(∞) deviates from the result for the weak-ergodicity breaking system; instead, α2(∞) is given by α2(∞)=vΓ(v)2[1+CD(∞)]/Γ(2v)−1 (SI Appendix, Text S7). These results suggest that the finite value of α2(∞) can serve an alternative measure of nonergodicity, which is constant in time unlike the ergodicity-breaking (EB) parameter proposed by He, Burov, Metzler, and Barkai (67).

Intrinsic and Extrinsic Disorder

For ergodic fluid systems, the long-time limit value of the product between the MSD and NGP serves as a measure of disorder strength. This disorder strength measure is decomposable into intrinsic and extrinsic disorder. To show this, consider the long-time asymptotic behavior of the MSD and NGP:Δ2(t)≅2dD¯t+Δc (t→∞)[7a]α2(t)≅2C^D(0)+Δc/dD¯t (t→∞).[7b]Eq. 7a can be obtained by substituting the Maclaurin series of 〈D^Γ(s)〉, 〈D^Γ(s)〉=(ε2/2d)〈κ^Γ(s)〉=(ε2/2d) × [〈κ^Γ(s)〉+〈κ^′Γ(s)〉s+⋯], into Eq. 3a and by taking the inverse Laplace transform of the resulting equation. We present the derivation of Eqs. 7a and 7b in SI Appendix, Text S7. In Eq. 7a, Δc is given by Δc=2∫0∞dt〈DΓ(t)〉t =2∫0∞dt〈v(t)⋅v(0)〉t and vanishes only when the VAF decays infinitely fast or only when velocity is white noise. In our random walk model, Δc emerges whenever the waiting-time distribution is a nonexponential function (SI Appendix, Text S8). On the other hand, in Eq. 7b, C^D(0)[=∫0∞dtCD(t)] emerges whenever the diffusion kernel or the waiting-time distribution is coupled to environmental variables.

From Eqs. 7a and 7b, we obtain the long-time limit value of the product of the MSD and NGP aslimt→∞〈r2(t)〉α2(t)=2[Δc+2dD¯C^D(0)].[8]We define the dimensionless disorder strength of complex fluids as limt→∞〈r2(t)〉α2(t)/σ2, with σ being the effective diameter of a tracer particle. Eq. 8 tells us that disorder strength has two different components, 2Δc/σ2 and 4dD¯C^D(0)/σ2, which originate from a finite relaxation time of the mean diffusion kernel and from environment-coupled fluctuation of the diffusion kernel, respectively. We designate the latter term extrinsic disorder, which is quite sensitive to temperature and density of the environment, as demonstrated in Fig. 2E and SI Appendix, Fig. S1E. On the other hand, we designate 2Δc/σ2 intrinsic disorder, because this term persists even when environment-coupled fluctuation in transport dynamics is negligible.

Intrinsic disorder is far less sensitive to the temperature and density of media than extrinsic disorder and can be easily estimated from Eq. 7a, the asymptotic long-time behavior of the MSD. Extrinsic disorder can be estimated by a direct numerical calculation of 4dD¯C^D(0)/σ2 or, more simply, by subtracting intrinsic disorder from total disorder strength (Fig. 2E); i.e., 4dD¯C^D(0)/σ2 =[limt→∞〈r2(t)〉α2(t)−2Δc]/σ2.

Disorder strength, defined in Eq. 8, is directly related to the whole-time integration of the BCF; that is, limt→∞〈r2(t)〉α2(t)/2 =Δc+2dD¯C^D(0) =6d2+dβ^(0)/D¯. This follows from the small-s limit of X^4(s)=Δ^4(s)−(1+2/d)Δ^22(s) =4!dβ^(s)/s2, obtained from Eq. 4, and from Eqs. 3a and 3b (SI Appendix, Text S6).

Model-Based Quantitative Analysis of the MSD and NGP

Intrinsic disorder causes the MSD to deviate from 2dD¯t, the prediction of the simple diffusion equation. We find that the following formula provides an excellent approximate description of the entire time range of the MSD for various disordered fluids (SI Appendix, Text S9):Δ2(t)=2dkBTMγ02c0(γ0t−1+e−γ0t)+2dkBTM∑i=1nciω0,i2[1−e−γit(cosh⁡ωit+γiωisinh⁡ωit)].[9]This equation represents the MSD of a bead in a Gaussian polymer, but also quantitatively explains the MSD of liquid water and dense hard disk fluids (Fig. 2A and SI Appendix, Fig. S1A).

The applicability of Eq. 9 to various disordered fluid systems implies a universality in the MSD time profile of disordered fluids, which is decomposable into one unbound-mode dynamic and multiple bound-mode dynamics, comparable to viscoelastic motion of a bead in a polymer network. At short times, a tracer molecule is trapped by the surrounding molecules. This bound state consists of multiple bound modes, each with their own characteristic frequencies. Meanwhile, at long times, a tracer molecule escapes the cage of the surrounding molecules and moves around in the media, repeatedly being caged and escaping the cage. The first term on the right-hand side (R.H.S.) of Eq. 9 accounts for the unbound mode, and the second term accounts for the bound modes. In Eq. 9, ci and γi designate the weight coefficient and relaxation rate of the ith mode (0≤i≤n). The weight coefficients are normalized by ∑i=0nci=1. ω0,i is the natural frequency of the ith bound mode and is related to ωi as ωi=γi2−ω0,i2.

At all times, temperatures, and densities investigated, Eq. 9 with only two bound modes (n=2) already provides a good quantitative explanation of the simulation results for the anomalous MSD of supercooled water and hard disk fluids (Fig. 2A and SI Appendix, Fig. S1A) and the experimental results for colloidal beads moving along lipid tubes (SI Appendix, Fig. S3) (36). According to Eq. 3a, the analytic expression of the mean diffusion kernel yielding the MSD given in Eq. 9 can be obtained by〈DΓ(t)〉=kBTMc0e−γ0t+kBTM∑i=1ncie−γit[cosh(ωit)−γiωisinh(ωit)].[10]Fig. 2C shows the mean diffusion kernel, or the VAF, calculated from Eq. 10 with parameter values optimized against MSD data from MD simulation shown in Fig. 2A for supercooled water.

The NGP is dependent not only on the mean transport dynamics, 〈D^Γ(s)〉, but also on fluctuation of transport dynamics, CD(t). For simple diffusion, we have 〈D^Γ(s)〉=D¯ and CD(t)=0, and Eqs. 3a and 3b yield Δ2(t)=2dD¯t and Δ4(t)=(1+2/d)Δ2(t)2 so that the NGP vanishes. However, whenever 〈DΓ(t)〉 is not constant and/or CD(t)≠0, the NGP does not vanish. We find that, for disordered fluid systems investigated in this work, the time profiles of the NGP cannot be quantitatively understood when we neglect fluctuation in the diffusion kernel or when we assume CD(t)=0 (SI Appendix, Fig. S4).

The NGP of disordered fluids is a nonmonotonic function of time with a single peak. According to our model, the NGP quadratically increases with time, α2(t)∝t2 at short times (SI Appendix, Text S10 and Eq. S10-9), but decreases with time, α2(t)∝t−1, at long times following Eq. 7b. As shown in Fig. 2A, it is only after the NGP peak time that Fickian diffusion emerges. These properties of the NGP are not specific to supercooled water but common across various disordered fluids (9, 12).

It was recently shown that diffusion coefficient fluctuation strongly correlates with string-like cooperative motion in dense fluids (57, 58), which is reportedly related to the NGP peak height (58). We find that the NGP peak height, α2(τng), serves as a measure of the relative variance of the diffusion coefficient for supercooled water. From the displacement distribution at the NGP peak time, we can extract the distribution of the diffusion coefficient using the method proposed in ref. 37. We find the relative variance, ηD2, of the extracted diffusion coefficient distribution has the same value as the NGP peak height (SI Appendix, Fig. S5D). This is not a coincidence. We find the NGP peak height has the same value as the relative variance of the diffusion coefficient at the Fickian diffusion onset time or the NGP peak time, τng (SI Appendix, Text S11). Both the NGP peak height and the NGP peak time increase with inverse temperature and density (Fig. 2B and SI Appendix, Fig. S1B).

Explicit Model for Diffusion Kernel

In the previous section, we demonstrated that the time profile of CD(t) can be extracted from the MSD and the NGP based on Eqs. 3a and 3b using Eq. 9 for the MSD. To achieve a physical understanding of the time profile of CD(t), we construct an explicit model of the environment-coupled diffusion kernel. Let us first consider the Laplace transform of Eq. 10, 〈D^Γ(s)〉=c0f^0(s)+∑i=1ncif^i(s), where f^0(s) and f^i(s) denote the diffusion kernels associated with the unbound-mode dynamics and the ith bound-mode dynamics, given by f^0(s)=(kBT/M)(s+γ0)−1 and f^i(s)=(kBT/M)s[(s+γi)2−ωi2]−1, respectively. We can extend this equation by assuming the weight coefficients {c0,c1,⋯,cn} are dependent on environmental state variables Γ, obtaining the following model of the diffusion kernel:D^Γ(s)=c0(Γ)f^0(s)+∑i=1nci(Γ)f^i(s).[11]From this model, we obtain the analytic expression for C^D(s) (SI Appendix, Text S12),C^D(s)=〈D^Γ(s)〉−2[γ02〈δD2〉ϕ^D(s)(s+γ0)2+∑i=0n∑j=0n′C^ij(s)f^i(s)f^j(s)],[12]where the prime notation signifies that the sum excludes the term with i=j=0, and Cij(t) denotes the time correlation between weight coefficients; i.e., Cij(t)=〈δci(t)δcj(0)〉. Noting that lims→0f^0(s)=kBT/Mγ0 and lims→0f^i>0(s)=0, we obtain D^Γ(s)≅c0(Γ)(kBT/Mγ0)(≡DΓ) from Eq. 11 in the small-s regime, s≪γi. Therefore, we can relate the weight coefficient TCF, 〈δc0(t)δc0(0)〉, of the unbound mode to the TCF of the diffusion coefficient fluctuation by 〈δc0(t)δc0(0)〉 =〈δD(t)δD(0)〉(Mγ0/kBT)2 at times longer than any element of {γi−1}. We find that the largest element of {γi−1} is γ0−1, whose value has order of 1 ps for the unbound mode for supercooled water (Table 1). At times longer than the NGP peak time, τng, which is greater than the largest velocity relaxation time, γ0−1 (Fig. 2 B and C), the bound-mode terms are negligible compared with unbound-mode term, so that the first term on the R.H.S. of Eq. 12 contributes the most to the relaxation of diffusion kernel fluctuation, leaving us with CD(t)≅ϕD(t)ηD2 (t>τng). This result indicates that the DKC becomes the TCF of the diffusion coefficient at long times. Using this result and recalling that ηD2≅α2(τng) (SI Appendix, Fig. S5D), we can then extract ϕD(t) from the time profile of CD(t)/α2(τng) at times longer than the NGP peak time. The long-time tail of CD(t) or ϕD(t) extracted from the MSD and NGP can be explained by an explicit model of the diffusion coefficient fluctuation described later in this work.

We note here that the whole-time integration of the diffusion kernel, CD(t), the determining factor of extrinsic disorder, is mostly contributed from the unbound-mode component (SI Appendix, Fig. S6). The bound-mode contribution to CD(t) has a comparable magnitude to the unbound-mode contribution, but has a negligibly smaller relaxation timescale than the unbound-mode contribution; consequently, the unbound-mode component makes the dominant contribution to C^D(0)=∫0∞dtCD(t). The major contributor to CD(t) is the unbound mode at long times but the bound modes at short times. For example, for supercooled water at 193 K, the bound-mode components of CD(t) are dominant at times shorter than 30 ps, and the unbound-mode component is dominant at times longer than 5 ns (Fig. 3D).

Model-Free Quantitative Analysis of the MSD and NGP

By analyzing the numerical data of the time-dependent MSD and NGP using Eqs. 3a and 3b, we can extract the time profiles of the mean diffusion kernel and the DKC without assuming a physical model, such as Eq. 9. In Fig. 3, we demonstrate this model-free analysis of the MSD and NGP data for supercooled water at 193 K (Methods).

According to our simulation, shown in Fig. 3A, the MSD of supercooled water exhibits oscillatory behavior with a slight bump and dip between 0.1 ps and 1 ps. This mysterious nonmonotonic oscillation in the MSD time profile of supercooled water was previously reported in the literature (5, 68, 69). We find the nonmonotonic MSD time dependence is unrelated to the finite-size effect and emerges not only under a constant temperature condition but also under the constant energy condition (SI Appendix, Fig. S7). The origin of the slight oscillation in the MSD time profile may be attributable to the intermolecular hydrogen-bond stretching vibration in supercooled water, which was previously identified in the quenched normal mode spectrum of the TIP4P/2005 water model at low temperatures (70) (SI Appendix, Fig. S8) and may also be the origin of the small oscillatory behavior in the NGP time profile between 0.1 ps and 1 ps (Fig. 3B). We find that these oscillations in the MSD and NGP time profiles are absent in liquid water above the melting temperature (Fig. 2A) and hard disk fluids (SI Appendix, Fig. S1A) at any density and cannot be accurately represented by Eqs. 9 and 18, used for the model-based analysis of the MSD and NGP in the previous section.

The mean diffusion kernel and the diffusion kernel correlation extracted from the model-free analysis of the MSD and NGP transiently display oscillatory behaviors at times around 0.1 ps, which are more complicated than the behavior of their counterparts extracted using Eqs. 9 and 18 of the MSD and NGP (Fig. 3 C and D). However, at times shorter than 0.01 ps or longer than 0.4 ps, the model-free analysis yields essentially the same results as the model-based analysis. As shown in Fig. 3D, Inset, both the model-free and model-based methods yield essentially the same result for the long-time DKC or the TCF of the diffusion coefficient fluctuation, quantitatively explainable by the unbound-mode component, or the first term on the R.H.S. of Eq. 12, only. Consequently, both methods yield the same value for the whole-time integration of the DKC, or extrinsic disorder, 4dD¯C^D(0)/σ2(=[limt→∞〈r2(t)〉α2(t)−2Δc]/σ2), and hence for intrinsic disorder as well, because total disorder can only be the sum of intrinsic and extrinsic disorder.

Microscopic Measurement of the Bound- and Unbound-Mode Components of Diffusion Kernel Correlation

To test the correctness of our results in the previous sections, we perform an alternative, microscopic measurement of the mean diffusion kernel and DKC using MD simulation and compare the results with those obtained in the previous sections. We then show that the bound-mode and unbound-mode transport dynamics, separately embodied in our transport model, clearly manifest, respectively, on the short-time and long-time dependence of the displacement distribution and the spatial volume spanned by the MD simulation trajectories for supercooled water at 193 K.

The mean diffusion kernel and DKC calculated from the MD simulation results of the two- and four-point TCF are found to be in excellent agreement with those extracted from the MSD and NGP. This is demonstrated for an example of supercooled water at 193 K in Figs. 3 and 4.

The short-time DKC, dominantly contributed from the bound-mode component, can be directly calculated using Eq. 6 and the direct MD simulation results of the VAFs; the BCF appearing in Eq. 6 can be calculated from its definition, Eq. 5. Note that the BCF can also be calculated from the MSD and NGP data, with use of the following relation: β(t)=∂t2X4(t)/4!d and X4(t)=(1+2/d)Δ2(t)2α2(t). These two methods produce similar, noisy BCF time profiles as shown in Fig. 4A for supercooled water at 193 K. When substituted into Eq. 6, these two BCF time profiles yield essentially the same results for the DKC, which is also in perfect agreement with the DKC extracted from our analysis of the MSD and NGP in the previous section (Fig. 4D).

We find that, at long times, the DKC is linearly proportional to the BCF (SI Appendix, Text S6):CD(t)≅32+dβ(t)/D¯2.[13]Recalling that the DKC becomes the TCF of the diffusion coefficient fluctuation, CD(t)≅ϕD(t)ηD2 =〈δD(t)δD(0)〉/D¯2 at long times, we obtain the long-time approximation of Eq. 13 as 〈δD(t)δD(0)〉≅3β(t)/(d+2). This result tells us that, for a one-dimensional system, the BCF is the same as the diffusion coefficient fluctuation, which was previously recognized by Nieuwenhuizen and Ernst (54) and others (71, 72) for a one-dimensional system of independent charged particles hopping on a lattice with static disorder. Our result here indicates that the long-time BCF has the same time profile as the long-time DKC multiplied by (d+2)/3, which is shown in Fig. 3D, Inset. It is not easy to calculate the long-time profile of the BCF directly from its definition, Eq. 5, because of the large computational expense involved in accurately estimating the multipoint VAFs from MD simulation and the 2D integral appearing in Eq. 5.

An independent estimation of the long-time DKC can be made by using MD simulations to measure the diffusion coefficient fluctuation along each time trace and calculating the TCF of the diffusion coefficient (Methods). This is because the DKC becomes the diffusion coefficient at long times, as mentioned earlier. In Fig. 4C, for supercooled water at 193 K, we show that ηD2ϕD(t) calculated from direct MD simulation is actually in good agreement with the long-time profile of CD(t) extracted from the MSD and NGP data. The long-time DKC is dominantly contributed from the unbound-mode component, or the first term on the R.H.S. of Eq. 12, as shown in Fig. 4C.

An alternative estimation of the long-time DKC can be made from the time profile of the NGP. The long-time DKC is directly related to the NGP and its time derivatives as follows:CD(t)≅12∂2∂t2[t2α2(t)]=α2(t)+2tα˙2(t)+2−1t2α¨(t).[14]This equation can be obtained by substituting the long-time expression of the fourth cumulant, X4(t)≅(1+2/d)(2dD¯t)2α2(t), into the well-known equation, β(t)=∂t2X4(t)/4!d, and then substituting the result into Eq. 13. Eq. 14 tells us that the NGP carries the complete information about the long-time relaxation of mobility fluctuation of complex fluids. As shown in Fig. 4C for supercooled water at 193 K, the long-time DKC estimated by Eq. 14 quantitatively agrees with the long-time DKC obtained from three other methods, namely, extraction from the MSD and NGP data, using Eq. 6 and MD simulation results of the VAF, and MD simulation of the diffusion coefficient fluctuation and calculation of its TCF.

Bound-mode (unbound-mode) transport dynamics are reflected on the time dependence of the displacement distribution at short (long) times. For supercooled water at 193 K, the displacement distribution broadens rapidly before 1 ps but after this does not greatly change for several picoseconds (Fig. 4D). This bound-mode feature of transport dynamics also manifests itself on the time-dependent volume spanned by the MD simulation trajectories of a water molecule. As shown in Fig. 4E, the spatial volume spanned by the simulation trajectories rapidly increases with time before 1 ps, but afterward, this trajectory volume tends to saturate to a certain critical value over several picoseconds while the trajectory length continues increasing with time. In contrast, at long-time scales, the displacement distribution and the trajectory volume exhibit unbound-mode dynamics, as demonstrated in Fig. 4 F and G.

Quantitative Explanation of Fickian Yet Non-Gaussian Displacement Distribution

Disordered fluids exhibit non-Gaussian diffusion; that is, the displacement distribution is non-Gaussian even at long times where the MSD linearly increases with time (9, 12, 36⇓–38). The displacement distribution of disordered fluids starts as Gaussian with variance given by r(Δt)−r(0)≅v(0)Δt at short times but deviates at any finite time. At the NGP peak time, where the deviation from Gaussian is greatest, the displacement distribution often appears similar to a Laplace distribution with an exponential tail (50). It is at this time that the non-Gaussian displacement distribution begins its relaxation to Gaussian and the MSD becomes linear in time. This phenomenon is widely observed across various disordered fluid systems (7, 9, 12), but has yet to be quantitatively explained.

To understand the time-dependent relaxation of the non-Gaussian displacement distribution in the Fickian diffusion regime, we need an explicit model of the diffusion coefficient fluctuation for the fluid system in question. In the literature, the diffusion coefficient is often modeled as D=A⁡exp(−βE), where A, β, and E are the entropic factor, inverse thermal energy, and activation energy, respectively (73). This model can be generalized by assuming that the entropic factor and activation energy are stochastic variables dependent on hidden environmental variables. In this work, we make this generalization and consider two exactly solvable models.

The first model assumes that the diffusion coefficient is given by DΓ=AΓ⁡exp(−βEΓ), where the fluctuation of EΓ around its mean value, 〈EΓ〉, is given by EΓ(t)−〈E〉=∑kbkΓk(t), where {bk} and AΓ are constants and {Γk} are stationary Gaussian Markov processes, also known as Ornstein–Uhlenbeck (OU) processes, satisfying 〈Γk(t)Γl(0)〉=δkl⁡exp(−λkt) (SI Appendix, Text S13) (74). For this model, exact analytic expressions of 〈D〉, ηD2, and ϕD(t)ηD2 can be obtained (Table 1 legend and SI Appendix, Eq. S13-3), where the adjustable parameters are optimized against the diffusion coefficient value, the NGP peak value, and the long-time DKC data (SI Appendix, Text S13). The optimized parameter values are presented in Table 1. Using the first model with optimized parameter values, we can now predict the time-dependent relaxation of the non-Gaussian displacement distribution in the Fickian diffusion regime (SI Appendix, Text S13). The prediction of this model is in excellent agreement with the MD simulation results for the displacement distribution of supercooled water, as shown in Fig. 5A.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Quantitative explanation of displacement distributions for supercooled water and colloidal beads on lipid tubes. (A) Displacement distributions of a water molecule along a Cartesian coordinate at three different times, τα(≅2.9ns), 5τα(≅14.4ns), and τEB(≅36.4ns). Circles, simulation results for TIP4P/2005 water at 193 K; lines, theoretical predictions of our first model. (B) Scaled displacement distribution of colloidal beads with diameter σ along a lipid tube at various times. Circles, experimental results reported in ref. 36; lines, theoretical results of our second model (SI Appendix, Texts S12 and S14).

In the second model, we model the diffusion coefficient as DΓ(t)=AΓ(t)exp(−βEΓ(t)) ≅〈D〉∑kakΓk2(t), where {ak} and {Γk(t)} are constants and OU processes, respectively. This is a generalization of the model in ref. 47 and yields an analytic expression of the displacement distribution (Table 2 legend and SI Appendix, Eq. S14-3). We find this expression provides an excellent quantitative explanation of the experimentally measured displacement distribution of colloidal beads diffusing on lipid tubes reported in ref. 36 (Fig. 5B). The optimized parameters of the second model are presented in Table 2. Eq. 12 with optimized parameter values allows us to calculate the time profiles of the DKC for the colloidal bead system (Table 2 legend and SI Appendix, Fig. S3).

The displacement distribution approaches a Gaussian distribution only after individual displacement trajectories become statistically equivalent. If individual displacement trajectories are statistically equivalent, the EB parameter proposed by He, Burov, Metzler, and Barkai (67) is linear in t/tmax, where t and tmax denote the time lag, or the interval over which the time-averaged MSD is calculated, and the maximum trajectory length, respectively (75). Otherwise, the EB parameter deviates from its linear dependence on t/tmax. At temperatures lower than 230 K, the EB parameter of supercooled water shows an anomalous power-law dependence on t/tmax at short times but resumes normal linear dependence on t/tmax at times longer than the characteristic time τEB (SI Appendix, Fig. S9). Deviation of the displacement distribution from Gaussian becomes negligible only at times much longer than τEB, as demonstrated in SI Appendix, Fig. S9B for supercooled water. On the other hand, CD(t) is negligibly small at the characteristic time τEB (SI Appendix, Fig. S9B). This is due to the fact that the long-time relaxation of the NGP is contributed not only from extrinsic disorder leading to the trajectory-to-trajectory variation in the transport dynamics, but also from intrinsic disorder, whose effects persist even for homogeneous systems with statistically equivalent displacement trajectories (Eq. 7b). This analysis shows that the long-time tail of CD(t) better characterizes the relaxation of diffusivity fluctuation than the NGP (76).