The heterogeneous energy landscape expression of KWW relaxation

Here we show a heterogeneous energy landscape approach to describing the Kohlrausch-Williams-Watts (KWW) relaxation function. For a homogeneous dynamic process, the distribution of free energy landscape is first proposed, revealing the significance of rugged fluctuations. In view of the heterogeneous relaxation given in two dynamic phases and the transmission coefficient in a rate process, we obtain a general characteristic relaxation time distribution equation for the KWW function in a closed, analytic form. Analyses of numerical computation show excellent accuracy, both in time and frequency domains, in the convergent performance of the heterogeneous energy landscape expression and shunning the catastrophic truncations reported in the previous work. The stretched exponential β, closely associated to temperature and apparent correlation with one dynamic phase, reveals a threshold value of 1/2 defining different behavior of the probability density functions. Our work may contribute, for example, to in-depth comprehension of the dynamic mechanism of glass transition, which cannot be provided by existing approaches.

The famous Kohlrausch-Williams-Watts (KWW) relaxation function or the stretched exponential relaxation function is an important observation in complex systems from the intricate behavior of liquids and glasses, the folding of proteins, to the structure and dynamics of atomic and molecular clusters, describing well the phenomena of important time-dependent dynamic processes [1][2][3][4][5][6][7][8][9][10][11] . The ubiquitous character of the KWW relaxation has shown irreversibility on the atomic, molecular or electronic scale and the dynamic nature of irreversible processes can be scrutinized in the context of the H-theorem to equilibrium, with the glassy state highlighting the limiting non-equilibrium behavior 1 . The dynamics of protein conformational changes clearly follows the KWW relaxation modes 2 and geometric frustration can happen once lattice structure averts simultaneous minimization of local interaction energies 3 . KWW related slow dynamics and internal stress relaxation in bundled cytoskeletal network is essential for the mechanical properties of living cells 4 , in contrary to the stretched relaxation of flux-freezing breakdown in high-conductivity magnetohydrodynamic turbulence 5 . Most often, phenomena of the KWW relaxation are typical of glass forming liquids and other complex fluids and have been extensively investigated in such a context 1,10,12 .
Since there is no obvious mathematical means to analytically transform the function f (t) in spite of its simple form, so a proper resolution and understanding of the function imperatively relies on its relaxation time spectra, which is still evading due to the complexity of the function and non-closed analytic approaches used in the previous research. Nevertheless, attempts have been made to explicate the stretched exponential behaviour as a linear superposition of simple exponential decay 13,14 , t 0 taking τ KWW = 1 for brevity. Eq. 2 is an inhomogeneous Fredholm equation of the first kind, in which the problem is to get the function ρ(τ), provided the continuous kernel function τ − / e t and the function f (t) 18 . ρ(τ) plays the role of the distribution of relaxation times as the probability density function of the relaxation modes. The solution of ρ(τ) can be computed from the series expansion 14 . However, problems of oscillation and deviation arise due to truncations from calculating the series expansion in the non-closed form 13,19,20 . We shall use an alternative distribution, the modulus function τ ( ) G 17,18 , defined in the way of There is a simple relation between ρ τ ( ) and τ ( ) . The study of the KWW relaxation is turned into the computation of τ Evidently, an accurate inverse transformation of the KWW function in a closed form is of importance in applications, particularly relevant to processing experimental data 13,14,19,20 , but it needs tremendous efforts. From the viewpoint of the dynamic free energy distributions and heterogeneity of relaxation as well as the characteristics of a rate process, here we present a heterogeneous energy landscape scheme to obtain the relaxation time distributions of the dynamic modes of a dynamic process which is dependent on the stretching exponential. In this way, we put the stretched exponential function on a solid physical basis, resolving the dilemma that in spite of the widespread success in describing relaxation data, the function is by and large viewed as an expedient phenomenological approach short of fundamental significance.

Results and Discussion
The concept of energy landscapes has been well explored in separate disciplines 8,10 . The spectrum of the KWW relaxation times implies a distribution of the free energies associating with the corresponding relaxation modes. In order to get such a distribution for a homogeneous process, we consider a global free energy random variable, ∆F in the reduced form of the free energy ∆G relative to the thermal energy k T (k B the Boltzmann constant and T the temperature), of a system as the sum over infinite many energy random variables (fluctuations) around its mean value from a constant random energy variable, ∆F 0 , at different levels of stochastic cascading with exponential distributions 21,22 . Suppose ∆F n (n = 1, 2, 3, …, m, m→ + ∞) are those independently, nonidentically distributed random variables, with the exponential distribution of (∆ ) = . With the zero mean value and the limited magnitude of the standard deviation, ∆F n represents a fluctuating contribution to the global energy quantity of the system. The roughness strength of ∆F n may be quantified through its standard deviation. As n increases, the measure of σ (∆ ) F n shows a harmonic-like dwindle.
We are interested in the limiting probability density distribution of the global energy variable ∆F defined as ∆ = ∆ + ∑ ∆ = F F F 0 n 1 m n as m→ + ∞. In the equation, the global energy variable ∆F has the same expectation, set to µ, as that of the constant random energy variable ∆F 0 since (∆ ) = (∆ ) E F E F 0 . Moreover, ∆F has a finite standard deviation square of σ = ∑ is divergent when m→ + ∞, so ∆F spreads over the domain (− ∞, + ∞). By some mathematical manipulation 23,24 , we are able to formulate the general probability density distribution function of the global free energy quantity ∆F as is verified as the probability density function of the global free energy distribution for the homogeneous process with the three parameters α, ε, and µ.
In reality, relaxation is a rate process and the characteristic relaxation time is related to the corresponding free energy by the Arrhenius equation 25,26 . As a result, the general probability density distribution function of the global free energy quantity in Eq. 4 is converted to the relaxation time spectrum by the expression of τ α ε τ Furthermore, the realization of relaxation may pass a transient state during the rate process, which can go forward to a relaxed state or move back to the initial state without relaxing [25][26][27] . The rates of the forward and backward relaxation are probably correlated to the duration of dwelling on the transient state, which can be characterized by the corresponding relaxation time. Hence, the forward and backward transmission rate is assumed to have We turn to consider the fact that the heterogeneous dynamics in glasses and other complex systems is attributed to the transitory coexistence of two dynamical sub-processes (phases) characterized by a fast and a slow relaxation rate in general 10,28 . In this scenario, the two sub-processes or dynamic phases contribute to the total relaxation, probably separating during relaxation and mixing afterwards. Therefore, the modulus function τ ( ) G is composed of such two heterogeneous dynamic phases, , implicitly validating the expression of τ ( ) G 14 . For a given β, search is attempted (refer to Methods), based on the numerical data from Eq. 3 compared to Eq. 5, to find an optimized set of the parameters in the parameter space, as summarized in Table 1 for τ ( ) G , with the correlation coefficient reported to be 1. Fig. 1 shows the outcomes from the analyses of the numerical computation for β between 0.05 and 0.95. The modulus function τ ( ) G in Fig. 1a gives a strong dependence on the stretched parameter β. For the same β, the function reveals the monotonic trend of initial increasing, attaining the maximum and then decline. Moreover, a tighter distribution is found for a larger value of β than 1/2 versus the more spread distribution of β less than 1/2. The evaluation proves the accurate consensus between the numerical calculation by Eq. 3 and the derived results based on Eq. 5. In the supplemental Figure S1 we show the validity of τ ( ) G over a broader range of relaxation time. The analyses substantiate the proposition of dynamic phase coexistence in the KWW relaxation course 10 . Reviewing the parameters obtained, the transmission coefficient κ has a different weight for different values of β, less independence on the relaxation time τ for small β but bigger reliance large β, hinting a threshold value of β = 1/2.
The probability density distribution ρ τ ( ) is more revealing in the exposition of the heterogeneous dynamic behavior of the KWW relaxation. The computation results are summarized in and Table 2 for ρ τ ( ), with the correlation coefficient recorded to be 1 (refer to Methods for the detailed computation procedure). Furthermore, the integration of ρ τ ( ) over τ is automatically normalized, confirming the property of the probability density function and unambiguously demonstrating the self-consistence and effectiveness of our approach including the accuracy of the numerical calculation and the validity of the equations derived. Fig. 1b shows the probability density function ρ τ ( ) dependent on the stretched exponential β. The dissimilar heterogeneous behavior is evidently manifested with a larger β than 1/2 which shows a phenomenon of an initial decrease followed by an increase and then drop off after reaching the maximum, in contrast to the monotonic decline of the distributions with a smaller β than 1/2.
The behavior of the probability density function ρ τ ( ) as a function of the stretched exponential becomes more distinctive if we plot the data in the log-log scale, as shown in Fig. 2a. In the figure, the numerical data calculated from Eq. 3 and the derived results based on Eq. 5 coincide over a broader range of relaxation time, showing the rationality of the approach adopted in this work in a closed, analytic form of the relaxation time spectra of the KWW relaxation. As already manifested in Fig. 1, ρ τ ( ) or τ ( ) G becomes more and more peaked around τ KWW when β approaches 1. This limiting behavior turns out to more distinguishing by re-plotting the data in the normal coordinates, as demonstrated in Fig. 2b for the probability density function ρ τ ( ). Fig. 3 gives the decomposed dynamic phases of the probability density distribution ρ τ ( ) for several representative β values. The value of β = 1/2 has a defining property, of which the two dynamic phases merge to have the same behavior. On the basis of analyzing the parameters as acquired in Table 2  function ρ τ ( ) switches from the scenario that the probability densities of the two dynamic phases share analogous monotonic decrease with the relaxation time for β below 1/2 to the observation that the two dynamic phases present a more complicated pattern for β above 1/2.
In Fig. 4, we provide detailed decomposition analyses of heterogeneity of the probability density function ρ τ ( ) for β = 0.3, β = 0.8 and β = 0.95. In general, two dynamic phases mix for small β, but for large β the major phase dominates while the minor phases diminishes. The limiting feature is evidently manifested in the variation of the curves from β = 0.3 to β = 0.95, that is, the bimodal feature is rapidly diminishing as β approaches 1, with the fast growing magnitude of the major phase against the quick weakening contribution of the minor phase. Fig. 5 reports the verification of accuracy in the reverse computation outcomes of the KWW relaxation function ( ) f t applying the formulation of , using the parameters from Fig. 1 for β between 0.05 and 0.95 versus the theoretical curves (refer to Methods). The precise performance of the assessment is clearly exposed in the consistency of the computed data with the analytic results.
In order to analyze the KWW function in the domain of frequency, a Fourier transform is needed to explain dynamic susceptibilities and scattering experiments from the perspective of linear response theory 13,14,19,20 . Absent of analytical expression for the transform, nevertheless, previous numerical methods suffer from problems originating from approximations and truncations which yield undesired oscillations 13,14,19,20 . Our approach shuns the cutoff effects and scrubs out oscillations. The results of the Fourier transformation using the derived parameters from Fig. 1 are presented in Fig. 6 (refer to Methods). The susceptibilities, real part ε′ (Fig. 6a) and imaginary part ε′′ (Fig. 6b) as well as the loss tangents δ tan (Fig. 6c) demonstrate the relevant properties of well-defined smoothness with respect to the frequency domain and strong dependence on the stretching parameter β. The Cole-Cole plots in Fig. 6d illustrate the susceptible relation of the relaxation, indicating a robust β dependence. for β values between 0.05 and 0.95. The results manifest a strong dependence on β, and for the same β, τ ( ) G monotonically increases to attain a peak value and then decreases. b, Semi-log plots of the probability density function ρ τ ( ) for β values between 0.05 and 0.95. The outcomes reveal quite different, strong dependence on β, which divides the β values in two ranges split by β = 1/2. For the same β below 1/2, ρ τ ( ) shows the behavior of monotonic decrease, in contrary to the observation that for the same β above 1/2, ρ τ ( ) gives a rapid initial decrease, then increases to attain a peak value and then decreases. Associated with the KWW relaxation is one important issue in condensed matter physics concerning glass transition of glass-forming materials, sharing the characteristics of free energy landscape, non-equilibrium, and heterogeneity 10 Table 2. Derived parameters of the function ρ (τ). Note: (a) Log-log plots of the probability density function ρ τ ( ) for β values between 0.05 and 0.95. The outcomes reveal quite different, strong dependence on β, which divides the β values in two ranges split by β = 1/2. For the same β below 1/2, ρ τ ( ) shows the behavior of monotonic decrease, in contrary to the observation that for the same β above 1/2, ρ τ ( ) gives a rapid initial decrease, then increases to attain a peak value and then decreases. (b) Linear plots of the probability density function ρ τ ( ) for β values between 0.05 and 0.95. For a small β, the two dynamic phases mix, but when β approaches 1, the major phase exclusively dominates with the minor phase disappearing, revealing the limiting behavior of ρ τ ( ) as β approaches 1.
value, corresponding to the wide-spread relaxation time distribution (Fig. 1). The relation between β and temperature is interesting and has been examined by numerical simulations or experiments 31 , but a direct connection is still elusive. Indeed, we have tried to follow the direction to work out such a correlation between β and temperature but it requires more efforts to reach a conclusive result. Nonetheless, it may be constructive to point out that bimodal or bimodal like distributions are observed, for example, in treating the dynamic order-disorder transition in atomistic models of structural glass formers 32 . The coexistence of the bimodal order parameter distributions is clearly related to the ordered and disordered phases. In our work, the bimodal like shape is observed in the density distribution of the relaxation time. A correlation could exist between the two, but it is recognizable that more work is of necessity to establish such a direct association between the dynamic order-disorder transition and the KWW relaxation.

Methods
This work reports a closed, analytic expression, Eq. 5, for describing the relaxation time probability density distribution function which is numerically calculated according to Eq. 3. As described below, the parameters in Eq. 5 are derived from the fit to the numerical data of Eq. 3. In this work, the fact that the two data sets coincide proves our approach, namely, Eq. 5 can excellently describe the relaxation time probability density distribution of the KWW relaxation. No other equations like Eq. 5 have been reported yet. In other words, to our knowledge, no equation other than Eq. 5 has been reported up to now to satisfactorily describe the numerical data from Eq. 3. We have performed the calculation with the help of the Mathematica and Origin software packages.
Based on Eq. 3 or the expression of the numerical data of the modulus function τ ( ) G or the probability density function ρ τ ( ) of the relaxation time distributions of the KWW relation were obtained via Mathematica for a fixed β value as a function of the relaxation time. Specifically, each data point of τ was computed up to 10 6 terms.
Then, we used Eq. 5 or the expression of , via the Origin program to conduct nonlinear regression of the data from the above numerical computation for a given β. Search was repeated until an optimized set of the parameters in the parameter space was found, with the correlation coefficient reported to be 1. The results are summarized in Table 1 for τ ( ) G and Table 2 for ρ τ ( ). Subsequently, the integration of ρ τ ( ) over τ was calculated, using the parameters recorded in Table 2. It is found that ρ τ ( ) is normalized for all β values discussed in this work.  1 k (in symbol) and the calculated results based on Eq. 5 in the shadowed regions. The limiting behavior is revealed from the peaking when β approaches 1. a, β = 0.3. b, β = 0.8. c, β = 0.95. The bimodal feature is rapidly diminishing as β approaches 1, with the fast growing magnitude of the major phase against the quick weakening contribution of the minor phase, unveiling the limiting behavior of ρ τ ( ) as β approaches 1. The reverse computation of the KWW relaxation function ( ) f t applied the formulation of ( )

Conclusions
We have shown a heterogeneous energy landscape approach to describing the Kohlrausch-Williams-Watts (KWW) relaxation function in a closed, analytic form, which is effective both in time and frequency domains. The equations obtained ascribe the heterogeneous dynamics of the KWW relaxation to the transitory coexistence of two dynamic phases as well as the characteristics of a rate process. The relaxation time probability density distribution acquired in this way changes upon varying the stretched exponential and, in particular, it is found that β = 1/2 marks a crossover from a small β regime to a large β regime. Our work significantly advances the mechanism of the KWW relaxation which cannot be provided by existing schemes and offers physical insights into the dynamic processes of glass transition and other complex phenomena.