Time-Dependent Fractional Diffusion and Friction Functions for Anomalous Diffusion

The precise determination of diffusive properties is presented for a system described by the generalized Langevin equation. The time-dependent fractional diffusion function and the Green-Kubo relation as well as the generalized Stokes-Einstein formula, in the spirit of ensemble averages, are reconfigured. The effective friction function is introduced as a measure of the influence of frequency-dependent friction on the evolution of the system. This is applied to the generalized Debye model, from which self-oscillation emerges as indicative of ergodicity that breaks due to high finite-frequency cutoff. Moreover, several inconsistent conclusions that have appeared in the literature are revised.


INTRODUCTION
In the early 20th century, Brownian motion became the subject of a theoretical investigation by Einstein, Langevin, Smoluchowski, and others [1][2][3]. The well-known conclusion was that the mean squared displacement of a force-free particle grows linearly in time. Because the Brownian trajectories were relatively short, Jean Perrin in 1908 used an ensemble average over many particle traces to obtain meaningful statistics. A few years later, Ivar Nordlune conceived a method for recording much longer time series. This let him determine time average individual trajectories and thus avoid average ensembles of particles that were probably not identical [3]. The average can be understood either as an ensemble average over a large number of trajectories or as a temporal moving average over a very long time trajectory. Nevertheless, the time average diffusion coefficient might be a random variable different from that of the ensemble, albeit the measurement time is long [4][5][6]. Since the ergodicity is broken in this sense, the diffusion coefficient can only be calculated by using the ensemble average. In theory, one mostly considers ensemble average; in the experiment, however, it is sometimes the former that taken the place of time averages because only one realization of a process is recorded. For an ergodic system, the time average is the same as the ensemble average [7]. The ensemble averaging result can be regarded as a comparable criterion.
Superseding Einstein's theory of Brownian movement, anomalous diffusion is characterized by the variance of the position of a diffusing particle increasing with time in a power-law form [8][9][10][11]. From an analysis of experimental data, the variance of the position for a force-free particle may be written in scalar form, where 〈/〉 denotes the ensemble average and 0 < α ≤ 2. This setup allows for various sorts of different physical realizations. The expectation for the particular behaviors of L 1 (t) and L 2 (t) is: L 1 (t) > 0, L 1 (t → ∞) 1, L 2 (0) 0, and lim t → ∞ t −α L 2 (t) 0. There have been a number of attempts to calculate the coefficients of transport for apparently dissipative systems, the Fourier transform of the velocity autocorrelation function (VACF) being widely used in semi-phenomenological spectral studies of fluids, magnets, and other systems [12]. Furthermore, the diffusion constant D α and exponent α deduced from data analyses need to be addressed accurately. In the anomalous diffusion sense, the generalization of the two basic statistical-dynamics laws-the Green-Kubo and the Stokes-Einstein relations-regarding measurable approaches for arbitrary frequency-dependent friction also has importance. This aids our understanding of the characteristic behaviors of anomalous, diffusive, and nonergodic systems [7,13].
There are many formalisms that describe anomalous diffusion, and growing interest has gathered around using the generalized Langevin equation (GLE) [14][15][16] as a viable alternative for investigating anomalous diffusion. One of the key features of the GLE is that it contains an after-effect function, termed a memory function [17,18], which can be obtained from experiments or molecular dynamics simulations [19,20]. Previous works on anomalous diffusion using the non-Ohmic friction model without high-frequency decay modulation showed that the VACF varies with time governed by the Mittag-Leffler function [21][22][23][24][25]. Indeed, a little transitive friction information can be extracted, although one knows that the diffusion coefficient determined by standard approaches vanishes for subdiffusion and diverges for superdiffusion. Furthermore, the issue of long-term memory [26] is always of central importance in nonequilibrium statistical mechanics, but there are several grounds for doubt [27], e.g., the dynamical effect related to the sharp cutoff in the spectrum being widely used, that require clarification.
The purpose of this paper is to set the diffusion coefficient and scaling exponent for anomalous diffusion processes; in particular, we will introduce another quantity, the effective friction function related to non-Stokesian relaxation, to account for the after-effect of anomalous diffusion. A measurement of this function does not present major difficulties for experiments yet provides a novel characterization of the underlying evolutionary process. The following questions also arise: What forms do the generalized Green-Kubo and Stokes-Einstein relations take in regard to anomalous diffusion? Can the coefficient of diffusion and the exponent be measured accurately? How strong is the effect of the high finite-frequency cutoff on the dynamics? The answers will become clear in the present study.

RELATIONS BETWEEN VARIOUS TEMPORAL FUNCTIONS 2.1 Time-Dependent Fractional Diffusion Function
To extract the pre-time coefficient D α in Eq. 1 when the exponent α is known, we have to generalize the definition of the coefficient of diffusion [28] using L'Hôpital's rule. The "fractional diffusion function" is assumed to be time-dependent, i.e.,D α (t) where 0z α t denotes a fractional Riemann-Liouville derivative of order α [23][24][25]. Here,D α (t) is related to D α L 1 (t) in Eq. 1, which converges asymptotically to the constant D α . For normal diffusion with α 1, Eq. 2 reduces to the standard diffusion function defined as D(t) 1 2 d dt 〈Δx 2 (t)〉 [28,29]. The latter has a clear physical meaning representing the temporal expansion rate of the spatial distribution. We establish a relation between the above two diffusion functions: where Γ(β) is the gamma function, and setting α n − β + 1, where n 0, 1, 2, . . . is a non-negative integer and β ≥ 0 a real number. In calculations concerning subdiffusion and normal diffusion (0 < α ≤ 1), we choose n 1 and then; for superdiffusion and ballistic diffusion [30,31] (1 < α ≤ 2), we choose n 2 and then β 3 − α. To date, no systematic numerical results regarding the fractional diffusion function for any contexts have been provided. However, they are required to understand the experimental observations.

Generalization of Two Famous Statistical-Dynamics Relations
The motion of a particle driven by a Gaussian distributed noise ε(t) is described by the following GLE: where M 1 (t) is the memory function and U(x) the potential. In general, colored noise ε(t) has a vanishing mean and is not correlated with the initial velocity. The noise correlation and memory function satisfy the fluctuation-dissipation theorem, expressed as 〈ε(t)ε(t ′ )〉 mk B TM 1 (|t − t ′ |), where k B denotes Boltzmann constant and T the temperature. We emphasize that, under the physical requirements of noisy dynamics [28,32,33], a rather free choice of the memory kernel or of the noise autocorrelation function is possible [34]. When the external potential is absent, i.e., when dealing with free diffusive motion, the formal solution of Eq. 4 is obtained by means of the Laplace-transform technique [35][36][37][38] The sum of the first three terms represents the square of the average position, and the sum of the latter two terms is the position variance. Note that the position variance is independent of the initial preparation of the particle.
We have found that the time-dependent diffusion functions, calculated by taking the time-derivative of the position variance and the MSD, are different. If the particle is confined initially at the origin x(0) 0 and its velocity obeys the Maxwell equilibrium distribution with {v 2 (0)} k B T/m rather than rest, using Eq. 5 and Therefore, the diffusion function obtained by differentiating . This is specifically an expression of the time-dependent Kubo relation encountered in the literature [11]. From the time derivative of the variance of position, i.e., Eq. 5, we obtain a rigorous exact relation between the standard diffusion function and the two relaxation functions, i.e., The VACF is calculated numerically from a homogeneous integro-differential equation, Here, we report a generalized time-dependent Green-Kubo relation, This expression is valid for anomalous diffusion and nonergodic processes as well as for arbitrary initial preparations. Equation 7 is referred to as the generalized time-dependent Kubo relation. In particular, under the condition C v (t → ∞) 0, D(t) approaches a constant when the upper limit of the integral is set to infinity, and, under {v 2 (0)} k B T m , Eq. 7 reduces to the standard Kubo relation [28]. The standard Stokes-Einstein expression that has just been reviewed applies to steady motion only and should be regarded as a zero-frequency theory. It is not consistent with using a coefficient of friction derived assuming steady motion describing changes in velocity. For this reason, the generalization of the Stokes-Einstein formula for an arbitrary frequency is of interest. For a broad class of systems, relaxation into the stationary state is exponentially fast; however, this is not true for all physical systems. Some systems may possess a stationary process, but the relaxation toward this state may be slow, and the well-known Stokes-Einstein formula thus needs to be generalized more. For anomalous diffusion, with a stationary VACF, the limiting result of the ratio of the time-dependent diffusion function to the velocity relaxation time, we report a generalized Stokes-Einstein relation: Here 0C v (t ′ )dt ′ is called the non-Stokesian relaxation time, which may become zero or infinity, and b C v (t → ∞) (1 + M ′ 1 (0)) − 1 represents the nonergodic strength of the first type [39,40]; however, b 0 for ergodic processes. In fact, Eq. 8 involves a spectral result [19], m , which also differs from recent work on the generalized asymptotic Einstein relation [20]. The latter reported a scale-dependent asymptotical Einstein relation for anomalous diffusion, i.e., is the noise correlation strength [20]. Fortunately, the present result [Eq. 8] seems universal as long as the integral over the VACF is used to replace that over the memory function.

Typical Long-Range Memory
We now consider a generic noise spectral density (NSD) generated by a model for non-Ohmic friction [21,22], for which the memory function takes the form whereω denotes the reference frequency for friction c δ having the dimensions of viscosity for any δ and f (ω) the modulation function of the frequency. In addition, the constant c δ allows for GLE (4) to have the correct dimension as well as f (ω) is a frequency modulate function. When 0 < δ < 2, the diffusive exponent α δ in Eq. 1, one can safely set f (ω) 1 [21]; nevertheless δ > 2, α 2, a decay form for f (ω), must be addressed.
In Figure 1, we plot the time-dependent fractional diffusion function calculated numerically using Eq. 3 combined with Eqs 7,9. All values of quantities used here on in are stated in dimensionless form, (i.e. k B T 1, m 1, c δ 1, andω 1). The scale-dependent fractional derivative is sufficient to produce finite fractional coefficients of diffusion in the longtime limit as expected. Notably, an overshooting peak arises during intermediate time periods, implying a strongly temporal "diffusion rate". For sub-diffusive situations, because of thermal fluctuations, the diffusion function starts from zero, increases with time, and finally decays to zero at late times. Hence, it remarkably ensures the existence of a temporal diffusive maximum.

Self-Consistent Extraction
More importantly, we want to introduce the temporal effective friction function, which is extracted consistently from Eq. 7 by assuming Frontiers in Physics | www.frontiersin.org February 2021 | Volume 9 | Article 567161 Starting from experimental sampling of the velocities of tagged particles with injection velocity v(0), we express the two functions h(t) and H(t) in alternative forms, An effective friction function that depends on both the physical features and the underling process is then defined, The advantage of this approach is that the entire configurational details are reduced to the particle's velocity, which is easily measured. Theoretically, it is easier to discuss a double average, one over time and the other over the ensemble of different trajectories.
For example, for the Ohmic friction (δ 1) with f (ω) ω 2 /(ω 2 + ω 2 c ), the Ornstein-Uhlenbeck (OU) colored noise is reached as M 1 (t) c 1 ω c exp(−ω c t) (the noise correlation time is τ c ω −1 c ). Then the effective friction function can be obtained analytically, where A (ω c + z 1 )/(ω c + 2z 1 ), B (ω c + z 2 )/(ω c + 2z 2 ), Figure 2 shows the time-dependent effective friction function for various δ. It starts from infinity because the system does not yet receive any dissipative energy from its environment at the initial time. This function approaches a constant value for normal diffusion, i.e., the low-frequency Markovian friction strength   0). With observations, a decrease/increase in the effective friction function corresponds to superdiffusion/subdiffusion. In contrast to the usual understanding, the effective friction function emphasizes the process dependence rather than the static result. From Eq. 13, the steady value of c eff (t) increases with decreasing ω c . Indeed, friction does vanish, but there is still no thermal excitation when ω c → 0. This implies that the effective friction function reflects its diffusive nature in that two effects are combined, the memory damping force and the random force, in generalized Brownian motion.

Determination of Diffusive Exponent
Lastly, the thorny problem is the measurement of the power-law exponent α from observed data. Of note is the fact that if L 2 (t) in Eq. 1 is a non-stationary function of time, the local differential method should be invalid for this situation. Although the leading term in the variance of the position may exceed the L 2 (t) oscillation term, as long as the time derivative of this variance is performed, the oscillating behavior is revealed. Previously, the approach α lim t → ∞ t d ln〈Δx 2 (t)〉 dt was used to evaluate the power exponent [41]. However, what is not valid for the situation of diffusion accompanied by oscillation when the spectral density of driving noise is the cut at finite or infinite frequency. Here, we propose a whole approach using the slope of the logarithm of the coordinate at large times, where t 1 and t 2 are two observation times at which the system has evolved into the asymptotic state.

The Effect of High-Frequency Cutoff
In contrast to usual understanding, the notion of the effective friction function proposed here is dependent on underlying processes rather than the static result. It is concluded from Figure 2, (e.g., the comparison between solid and dashed lines when δ 1) and Eq. 13 that the steady value of c eff (t) increases with a decrease of ω c . In fact, the friction does vanish, but there is yet no thermal excitation when ω c → 0. This implies that the effective friction function reflects a diffusive feature, which combines two effects of the memory damping force and random force in generalized Brownian motion. Nevertheless, the first choice in Eq. 9 cannot be adopted when δ > 2, which results in that the noise correlation function corresponding to the memory function does not approach vanishing at long times. This violates Kubo's requirements of noisy dynamics [28]. However, the cutoff function f (ω) in Eq. 9 should be applied safely. It is known that the treatment of highfrequency cutoff has been applied widely as well as α δ when 0 < δ < 2 and α 2 when δ > 2 [21]. However, little attention was given to the influence of this cutoff on dynamics.
The Laplace transform of the memory function for the non-Ohmic friction model combining with f (ω) Θ(ω c − ω), where Θ(ω c − ω) is the Heaviside function, equaling unity when ω ≤ ω c and vanishing when ω > ω c , is given by where 2F 1 denotes the hypergometric function. In particular, when δ 1, 2 F 1(1, 1/2, 3/2, −x 2 ) π √ arctan(x)/x. Assuming that s s 1 + is 2 , where θ arctan[(2s 1 ω c )/(s 2 1 + s 2 2 − ω 2 c )]. Although M 1 (s) is a multi-value function on the complex plane, where exists at least a fair of pure complex roots (i.e., s ± is 2 ) for the characteristic equation: s + A arctan(ω c /s) 0, which is given by s 2 According to the Khinchin theorem, if lim t → ∞ C v (t) ≠ 0, ergodicity is broken [42,43]. Two situations arise depending on whether C v (t → ∞) lim s → 0 [s h(s)] b{v 2 (0)} ≠ 0 or whether it does not exist. We establish here a universal condition for ergodicity breaking that there exists a zero root or at least a pair of real complex roots for the characteristic equation: s + M 1 (s) 0, namely, and the residues do not equal zero at these poles. The second term on the l.h.s of Eq. 17 is the Laplace transform of the memory function. We report also an alternative NSD, which can induce nonergodicity of three types, i.e., where ρ 0 (ω) is a function of the frequency, 0 ≤ ω L < ω 1 < ω 2 < ω H < ∞. Eq. 18 can be used to physically describe the generalized Debye model [16], the longwavelength limit of acoustic phonon [44], the Bethe lattice [45], and the localized mode [46]. Applying the residue theorem [47,48], we obtain a general result for the velocity relaxation function, Frontiers in Physics | www.frontiersin.org February 2021 | Volume 9 | Article 567161 The quantity is equal to a constant less than unity as long as ρ 0 (ω) ∼ ω δ (δ > 1) at low frequencies when ω L 0 and ρ 0 (ω) ∼ ω δ (δ < 1) at high frequencies if ω L ≠ 0. The existence of time oscillations in the VACF is because s ± iλ 1 and s ± iλ 2 may be numerical solutions of Eq. 17, which requires ω 1 < λ 1 < ω 2 and ω H < λ 2 ; also, A 1 and A 2 are each twice the residue at its respective pole. The last two terms in Eq. 19 vanish in the long-time limit because Re s n < 0 for any n.
Clearly, the upper-limit of the summation in Eq. 19 may be infinity because M 1 (s) calculated using Eq. 18 is a multi-valued function on the complex plane. Let us perform an integration via parts in Eq. 17 to render explicit the dependence of M 1 (s) on s, For a complex variable s, we have arctan(ω/s) (2i) − 1 ln|(s + iω)/(s − iω)| + 1 2 (θ + 2nπ) with n 0, ± 1, ± 2, . . ., where θ is the angle of the complex variable (s + iω)/(s − iω) obtained by expressing it as an exponential function. Further, it is noticed that both the wellknown Mittag-Leffler function [22] and the series solution [43,44] cannot be used to demonstrate the self-oscillating result at large times but Eq. 19 is reasonable. Figure 3 shows the resulting VACF for various δ and two different modulation functions of the frequency. Evidently, there exist three types of nonergodicity where the VACF 1) approaches a positive constant, 2) oscillates with time around a plateau value, or 3) oscillates with time around zero; the simple conditions satisfied by the NSD are low-hindering, band-passing, and high-hindering or band-hindering. This forces the ensemble to become unstable, although the GLE used is initially stationary. The oscillation implies that the particle tends to reverse its direction of motion frequently relative to its former step. This classification of the nonergodic behavior may be regarded as a probe of a more precise estimation of the low and high frequencies distribution for thermal colored noise. In addition, from the microscopic Hamiltonian viewpoint, ergodicity breakage is caused by a localized mode with an isolated frequency from the continued phonon spectrum [46,49].
In Figure 4, we plot the time-dependent diffusion function for various situations and show the diffusion function oscillating with time. Because the diffusion function produced carries out a time derivative over the MSD, the whole approach [Eq. 14] rather than a local differential method must be used to determine the power-law exponent of diffusion. Many models may induce such a phenomenon, and several famous examples have been mentioned before. A somewhat surprising result for us is that the system associated with the Debye frequency-cutoff does not yield ergodic behavior.

SUMMARY
This work aims at furnishing a connection between velocity autocorrelation function (VACF) and frictional kernel in the generalized Langevin equation (GLE) framework for the diffusion dynamics. On the one hand, the VACF is a quantity that can be obtained from experimental or numerical data; on the other hand, the friction kernel incorporates the noise correlation properties. However, this relation is hard to achieve by analytical means, and we thus introduce an effective friction function that depends solely on the VACF and is supposed to inform on the retarded memory effects on the dynamics. The notion of effective  friction function has been observed as an echo of the non-Markovian Brownian motion, which is associated with the diffusion function and is measured through a tagged particle's average velocity. This provides interesting information regarding the evolution of kinetics. We have obtained alternative timedependent Green-Kubo and generalized Stokes-Einstein relations to universal situations connected with anomalous diffusion and nonergodic processes. The effect of finite or high-frequency cutoff as a facile method on the dynamics is numerically investigated. A self-oscillation phenomenon emerges as a manifestation of ergodicity breakdown. In particular, several inconsistent conclusions in the literature, e.g., the analysis requirement of susceptibility, the friction feature extracted by the generalized Fokker-Planck equation, and the power-law exponent obtained from data, have been clarified. We are also confident that the present results will serviceably impact complex dissipative systems.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.