Deciphering Non-Gaussianity of Diffusion Based on the Evolution of Diffusivity

Non-Gaussianity indicates complex dynamics related to extreme events or significant outliers. However, the correlation between non-Gaussianity and the dynamics of heterogeneous environments in anomalous diffusion remains uncertain. Inspired by a recent study by Alexandre et al. [Phys. Rev. Lett. 130, 077101 (2023)], we demonstrate that non-Gaussianity can be deciphered through the spatiotemporal evolution of heterogeneity-dependent diffusivity. Using diffusion experiments in a linear temperature field and Brownian dynamics simulations, we found that short- and long-time non-Gaussianity can be predicted based on diffusivity distribution. Non-Gaussianity variation is determined by an effective P\'eclet number (a ratio of the varying rate of diffusivity to the diffusivity of diffusivity), which clarifies whether the tail distribution expands or contracts. The tail is more Gaussian than exponential over long times, with exceptions significantly dependent on the diffusivity distribution. Our findings shed light on heterogeneity mapping in complex environments using non-Gaussian statistics.

Because the heterogeneity of a complex environment causes the non-Gaussianity of diffusion, the central in-quiry becomes understanding the correlation between the variations in non-Gaussianity and the degree of heterogeneity.The diffusivity distribution and its variation can be applied to map the structural or dynamical heterogeneity of a complex environment based on the DifD model, because they contain physical properties of the heterogeneous environment, such as viscosity, permeability, or energy barriers.The difficulty of quantifying the evolution of diffusivity, particularly in biological environments [19], makes the DifD model phenomenological and raises doubts regarding its validity.To address this problem, among many recent efforts [20][21][22], an inspiring approach that drew an analogy with Taylor dispersion was proposed to mathematically link non-Gaussianity and heterogeneity [22].Nonetheless, controversies regarding variations in non-Gaussianity and tail distribution remain unresolved.
In this Letter, we used thermophoresis of nanoparticles (NPs, diameters d =500 and 1000 nm) in a microfluidic chip [Fig.1(a)] [23] with controlled temperature field to construct DifD of underlying probability space.We focused on clarifying the correlation between non-Gaussianity and DifD distribution p(D) by viewing the heterogeneous field as a spatiotemporal evolution of diffusivity.As the NPs move along the constant temperature gradient (x -axis) with a constant thermophoretic speed u T [Fig.1(b)], their diffusion, perpendicular to the temperature gradient (y-axis), experiences DifD determined by the temperature gradient, providing better controllability and quantifiability than that reported in recent studies [20][21][22][23][24]. Consequently, we can predict the short-time non-Gaussianity by the type and range of the DifD distribution, and show two long-time destinations depending on the presence of external field-driven migration.We found that the temporal variation of non-Gaussianity is determined by an effective Péclet number, which identifies the competition between the varying rate of diffusivity and the diffusivity of diffusivity [13].Unlike the majority perspective assuming an exponential tail, we use Laplace's approximation of Bayes' theorem to explain why a Gaussian distribution is more suitable.We established stable temperature gradients of ∇T = 4.1 × 10 4 and ∇T = 6.6 × 10 4 K/m (temperature differences of 8.1 and 13.2 K) respectively in the central microchannel of the microfluidic chip by setting flow rates for cold and hot water in both side channels (Fig. 1(a)-1(c), Supplemental Material [25]).The NPs moved toward the cold side along the x-axis with thermophoretic velocities u T =−0.41 and −0.66 µm/s, respectively [23] [Fig.1(d)].At the beginning of each experiment, particles located at x ∈[-80, 90] µm were chosen for tracking lasting ∆t=20 s.Observation was performed at the plane of z =15 µm to avoid wall-induced hydrodynamic drag, as the relative distance 2z/d was large.Particle trajectories with a spatial resolution of 80 nm were obtained from successive images captured at 20 fps [25].The NP's displacement in the x -axis is ∆x(t) = √ 2D x tW x +u T t+t∂D x /∂x ≈ √ 2D x tW x +u T t, where W x denotes an independent stochastic process with a mean of zero and standard deviation of one, D x (T ) is the local diffusivity depending on the temperature T(x), and ∂D x /∂x can be neglected as it is less than 1 nm/s.In the y-axis, the displacement follows ∆y(t) = 2D y tW y , which is determined by the DifD D y (t) varying with NP's thermophoretic migration along the temperature gradient in the x -axis.
The MSDs of 1000 nm particles under two temperature gradients are shown in Fig. 2(a), calculated using ⟨∆r 2 (t)⟩ = ⟨[r(t 0 + t) − r(t 0 )] 2 ⟩.Here, r represents x or y, and ⟨•⟩ denotes ensemble average.At short times, t ≈ 0.1 s, all MSDs display a linear tendency, as Brownian diffusion is dominant.When t > 2D x /u 2 T ∼ 2 s, the x -MSD ⟨∆x 2 (t)⟩ manifests a super-diffusive behavior with a slope near two, due to thermophoresis with speed u T .By contrast, the y-MSD ⟨∆y 2 (t)⟩ approximates a Fickian behavior.
To investigate the non-Gaussianity of this system, normalized DPDs at t=0.05 s were compared with the stan- where σ is the standard deviation.To quantify the amplification of the tailed DPD compared with the standard Gaussian distribution when ∆y/σ > 3, we show the tail ratio ϕ in the inset of Fig. 2(b).The maximum ratio was larger than two even though the motion in the y-direction was purely diffusive, manifesting a nonnegligible non-Gaussian tail beyond the error bars.The non-Gaussian parameter α(t) = (⟨∆r 4 (t)⟩/⟨∆r 2 (t)⟩ 2 ) − 3 [22] is introduced to quantitatively assess the above non-Gaussian behavior.For a Gaussian DPD, α=0, whereas positive and negative α values signify fat-tail leptokurtic and platykurtic distributions, respectively.Because the uncertainty of α is less than ±0.09 [25], the experimental data of α y suggest a slow decay from a positive α y ≈ 0.12 when t < 0.1 s to α y ≈ 0 when t ≈ 1 s.When t > 1 s, α y is still considered to be zero based on the error bars, which is further confirmed by numerical simulations.By contrast, α x starts from a small positive value, similar to α y , but decays more rapidly to a negative value when t > 1 s.Greater temperature gradient (∇T = 6.6 × 10 4 K/m) causes larger non-Gaussian parameter α y and fatter tail at short times.Experiments using 500 nm NPs show similar MSD and DPD results [25].
To tackle the issue of limited long-time statistics in the experiments, we used Brownian dynamics simulations [25] to complement the long-time evolution of non-Gaussianity.By assigning the same thermophoretic speed and uniform diffusivity distribution along the xaxis as the experiment, the simulation results [dashed and solid curves in Fig. 2(a)-2(b)] closely match the experimental data.The long-time α x (t) from the simulation collapses onto the extension of the experimental data when t ∼ 5 s.Next, primarily based on the simulation results, we explore the effect of the spatiotemporal variation of DifD on non-Gaussianity.The initial distribution of D r and its varying rate ∂⟨D r ⟩/∂t ∼ (∂⟨D r ⟩/∂T )∇T u T are controlled [25], and the diffusivity of diffusivity, defined as ∂⟨D 2 r ⟩/∂t, and the varying range of D r are monitored.
We first test the non-Gaussianity at ∇T = 3 × 10 5 K/m.After subjecting the particles to thermophoresis with = −3.0µm/s for ∆t=15 s, the diffusivity distributions when ensemble statistics was conducted were regulated by setting the ranges of x positions to [−95, 95] µm, [−60, 60] µm, and [−25, 25] µm, respectively, from the initial positions x ∈[−50, 95] µm, [−15, 60] µm, and [20,25] µm (SM [25]).A broader x range indicates a wider diffusivity distribution.In addition to the typical MSDs [Fig.2(c)] and DPDs [Fig.2(d)], similar to the experimental tendency, the variations in the non-Gaussian parameters in both directions are shown in Fig. 2(e)-2(f).The short-time α y (t) (solid symbols) starts from approximate constant values of 0.37, 0.11, and 0.03, following a gradual decay at t > 2D y /u 2 T ∼ 0.1 s.Interestingly, the non-Gaussianity shown by green symbols is similar to the experimental result with ∇T = 6.6 × 10 4 K/m [Fig.2(b)], because they have similar diffusivity ranges though different ∇T.For a wider diffusivity distribution, the value of α y (t) is larger and the tail is broader.The value of α x d (t)[Fig.2(f)], based on ∆x d = ∆x − u T t, is always almost the same as α y (t).Interestingly, α x (t) [empty symbols in Fig. 2(e)] rapidly turns negative and reaches −2 at t ∼ 10 s, despite shares the same beginning as α y (t).
The above results demonstrate the significant influence of diffusivity distribution on non-Gaussianity and fat-tailed DPD.Recalling a recent study suggesting an analogy of DifD motion with Taylor dispersion [22], we find that the asymptotic features of non-Gaussianity can be accurately predicted only if the distribution of D r is known.The mathematical derivation by Alexandre et al. [22] gives the fourth cumulant as: When t → 0, the non-Gaussian parameter can be approximately calculated by dividing the variance Var (D r ) over the square of expectation ⟨D r ⟩ 2 [22,25]: (1) Eq. ( 1) correlates α r (t → 0) with the diffusivity distribution p(D r ) in a probability space, which explains the same initial values α y =α x =α x d in Fig. 2(e)-2(f) because they share the same diffusivity.This relation could be applied to various systems with structural or dynamical heterogeneities by mapping the heterogeneity onto DifD.We modify distribution p(D r ) from uniform distribution, which is suitable for our experiments, using uniform sampling, to other commonly-used distributions.Surprisingly, we find that α y (t → 0) can be predicted solely based on the range ratio β(∆t) = D min /D max .We take a gamma distribution p(D y ) = D To verify Eq. ( 2), we assigned gamma distribution to D r in simulations and maintained ∇T and u T unchanged.The simulation results [Fig.2(g), ∆t=0.5 s] for different shape parameters, k=1, 2, and 4, were in good agreement with the prediction curve from Eq. ( 2).The result of gamma distributions can be extended to other exponential or power-law distributions and can be widely used to evaluate non-Gaussianity in complex scenarios.Besides, we derive the expressions of α r (t → 0) for uniform and normal distributions as , respectively, which perfectly match the simulation data in Fig. 2(h).By mapping to a standard normal distribution truncated at z ∈ The expression of α y (t → 0) when N < 3 is given in SM [25], which predicts a larger non-Gaussianity owing to a stronger dispersion of D y [pink curve in Fig. 2(h) for N =1].Note that α y (t → 0) in our simulation only varies slightly with the duration ∆t as the change of β(∆t) is tiny during slow evolution of diffusivity distribution, which has been manifested by the short-time plateau of α y (t) in Fig. 2(e)-2(f).The data (empty blue diamonds) are close to the theoretical curve of uniform distribution for ∆t=10 s when the diffusivity distribution is still approximately uniform.The experimental data acquired from the uniform samples were located near our prediction curve in Fig. 2(h) based on uniform distribution, suggesting the good predictive power of our approach up to ∆t=20 s.Furthermore, our results illustrate that non-Gaussianity in the same system can vary significantly depending on the sampling range of p(D r ).In particular, ergodic sampling with a minimal value of β can exhibit a much larger non-Gaussianity than truncated sampling with a larger β, as illustrated in Fig. 2(g-h).
We then consider the long-time limit of non-Gaussianity in Fig. 2(e-f).Although calculating the temporal evolution of p(D r ) is complicated for various heterogeneities, previous work [22] has predicted that α r (t) should decay with 1/t at long times.We indeed observed a 1/t decay to zero of α y when t ∼100 s [inset of Fig. 2(f)] in our simulation, when boundary confinement was absent.In realistic cases with confinement [21,22], 1/t decay was observed later when t ∼1000 s.This longtime scaling 1/t is assumed to be an intrinsic feature owing to diffusivity dispersion, which is independent on boundary confinement.
Distinct from the positive non-Gaussianity, the most evident feature of α x (t) is the asymptotic value α An intriguing inference is that the asymptotic value α x (t → ∞) = −2 will appear for any external field-driven migration with a large speed u 4 T t 2 ≫ ⟨D 2 x ⟩.Unlike the experiment which reported a rapid increase in non-Gaussianity at short times [21], our experimental and simulation results show a constantly decreasing α y (t).Next, we discuss the tendency of α y (t) by the derivative  space of D y , in analogy to traditional P e number in Taylor dispersion.As the varying rate can be approximated as T , in Fig. 3(a-b) we adjust the temperature gradient ∇T to investigate its influence on the variations of P e and α y .   in the simulation by changing the NP from thermophobic to thermophilic [25].Nonetheless, the increase in short-time non-Gaussianity in our simulation was much weaker than the rapid dynamics reported by Pastore et al. [21].We speculate that the difference was due to a sudden dispersion of ∂⟨D 2 r ⟩ ∂t produced by the opticalillumination speckle [21].The different non-Gaussian behaviors between the present system with slow DifD and the systems with rapid dynamics, such as glassy materials [10], can serve to detect specific short-time mechanisms.

Both the diffusivity of diffusivity
After clarifying the correlation between non-Gaussianity and distribution p(D), we further discuss how p(D) determines fat-tailed DPD, which is the most prominent feature of FnGD.Whether the tail of the DPD is exponential remains controversial [7,12,22,26], although a good exponential fit is observed in Fig. 2(d This approximation could help clarify the controversy regarding the tail as the term exp[− ∆y 2 4Dy(x * )t ] suggests a more Gaussian tail for large ∆y, which contradicts many existing results of exponential tails [6,7,11,21,26].As shown in the double-logarithmic plot of ln(G s √ 2π)vs.∆y/σ in Fig. 3(c), the slopes of gamma and normal distributions are approximately −2(Gaussian) for the tail ∆y/σ > 3. The function type of p(D y ) can influence the result of the integral by prefactor C in the Laplace approximation.For the uniform distribution of p(D y ), the slope of the tail deviates from −2 when ∆y/σ > 3 and becomes −1 near ∆y/σ = 6.This is consistent with the prediction of an exponential tail if ∂ln(p(D y ))/∂D y is constant, as reported by Wang et al. [12], using the steepest descent analysis for the integral G s ∼ exp[ln(p(D y )) − ∆y 2 4Dyt ]/ D y • dD y .Although the Gaussian tail for the gamma distribution is distinct from the prediction in [12,13], it is in accordance with Alexandre et al. [22], where a generic Gaussian tail was proposed.Additionally, a slight contraction of the tail turning more Gaussian with time, is observed for uniform and gamma distributions as the Gaussian term exp[− ∆y 2 4Dy(x * )t ] becomes more significant for large ∆y.The transition from a short-time exponential-like tail to a long-time Gaussian-like tail could clarify the controversy that similar systems may display contradictory tail shapes.As an exception, diffusion in strongly confined media [7], where statistical data cannot reach a sufficiently large ∆y, typically manifests an exponential tail.
Conclusion.-We investigated non-Gaussian diffusion in a thermophoretic system with controlled DifD and correlated non-Gaussianity with the spatiotemporal evolution of diffusivity distribution.
Short-time non-Gaussianity was predicted based on the range ratio β of diffusivity distribution.We demonstrated the predictive power using uniform, normal, and gamma distributions.The variation in non-Gaussianity at intermediate times was determined by the effective P e number, which characterizes the competition between the varying rate of diffusivity and the diffusivity of diffusivity.This explains why the tail of the DPD usually contracts, unless it is influenced by sudden dynamics that enhance the diffusivity of diffusivity.The long-time decay of non-Gaussianity followed a 1/t scaling in a purely diffusive process, whereas in the presence of particle migration, the non-Gaussianity eventually reached −2.Furthermore, using Laplace's approximation of Bayes' theorem, we depicted that the tail was more Gaussian than exponential.We showed the transition from a short-time exponentiallike tail to a long-time Gaussian-like tail to clarify the controversies in existing experiments.Our results shed light on heterogeneity mapping in complex environments based on non-Gaussian statistics, echoing the idea that Bayesian deep-learning deciphers the physics encoded in diffusion data [27].

FIG. 1 .
FIG. 1.(a) Experimental setup.A constant temperature gradient was created in the central microchannel of a microfluidic chip.(b) In x -direction (spanwise of microchannel), particle moves toward the cold side along the temperature gradient, and its displacement is determined by thermophoresis and diffusion.In y-direction (streamwise), the motion is diffusive with DifD.(c) Experimental temperature distribution in xdirection when ∇T = 4.1 × 10 4 K/m.Inset: simulation result of the uniform temperature field [25], z -axis is the height direction.(d) Trajectories of 1000 nm particles moving toward the cold side (left).

k− 1 y
exp(−D y )/Γ(k) as an example, which is truncated at D y ∈ [D min , D max ].By mapping it to a standard gamma distribution p(z) = z k−1 exp(−z)/Γ(k) truncated at z ∈ [0, N ]: (D y − D min )/(D max − D min ) = z/N , we obtain Var (D y ) = Var (z)(D max − D min ) 2 /N 2 and ⟨D y ⟩ = ⟨z⟩(D max − D min )/N + D min .For large N such as N =10, Var (z) ≈ k and ⟨z⟩ ≈ k, α y (t → 0) of the gamma distribution is: to as the diffusivity of diffusivity and the varying rate of diffusivity respectively.This competition defines an effective Péclet number P e = (⟨D y ⟩ ∂⟨Dy⟩ ∂t )/ ∂⟨D 2 y ⟩∂t , which helps assess whether α r (t) will increase or decrease with time.Here, ⟨D y ⟩, ∂⟨Dy⟩ ∂t , and∂⟨D 2y ⟩ ∂t can be seem as equivalent length, velocity, and diffusivity, respectively, in a
) are negative [inset of Fig. 3(a)], whereas the domination of the diffusivity of diffusivity, i.e., | ∂⟨D 2 y ⟩ ∂t | > | ∂⟨Dy⟩ 2 ∂t |, results in P e < 1 and )dD y .D y (T ) is determined by the temperature field.This integral can be approximated as G s (∆y, t) ≈ Cexp[− ∆y 2 4Dy(x * )t ] based on the Laplace approximation for large ∆y, where prefactor C depends on p(D y ), and x * is the position of maximum D y (x * ).