Fractional Time Fluctuations in Viscoelasticity: A Comparative Study of Correlations and Elastic Moduli

We calculate the transverse velocity fluctuations correlation function of a linear and homogeneous viscoelastic liquid by using a generalized Langevin equation (GLE) approach. We consider a long-ranged (power-law) viscoelastic memory and a noise with a long-range (power-law) auto-correlation. We first evaluate the transverse velocity fluctuations correlation function for conventional time derivatives C^NF(k→,t) and then introduce time fractional derivatives in their equations of motion and calculate the corresponding fractional correlation function. We find that the magnitude of the fractional correlation C^F(k→,t) is always lower than the non-fractional one and decays more rapidly. The relationship between the fractional loss modulus GF″(ω) and C^F(k→,t) is also calculated analytically. The difference between the values of G″(ω) for two specific viscoelastic fluids is quantified. Our model calculation shows that the fractional effects on this measurable quantity may be three times as large as compared with its non-fractional value. The fact that the dynamic shear modulus is related to the light scattering spectrum suggests that the measurement of this property might be used as a suitable test to assess the effects of temporal fractional derivatives on a measurable property. Finally, we summarize the main results of our approach and emphasize that the eventual validity of our model calculations can only come from experimentation.


Introduction
The prominent role that time correlation functions have played in the description of non-equilibrium properties of fluids stems from their close connection with their transport coefficients, a relation that can be obtained from time-dependent correlation functions of suitable fluxes [1]. In a simple fluid in near-to-equilibrium states the central limit theorem (CLT) applies and a well-defined separation between the time scales associated with the macroscopic transport processes and the microscopic processes giving rise to them exists. As a consequence, the stochastic dynamics of the thermal fluctuations around equilibrium can be characterized as a Gaussian-Markovian process [2][3][4][5]. However, the presence of long noise-correlation or long time-memory in the time evolution equations for the fluctuations may destroy this separation, and the usual description of fluctuations in terms of the conventional Langevin equations may no longer be adequate [6][7][8][9][10][11]. This situation may occur in a large variety of relaxation processes in complex systems like viscoelastic fluids, glassy materials, synthetic polymers or biopolymers, all of which have in common that their relaxation functions are non-exponential, due to the large number of highly coupled elementary units responsible for the relaxation. As a consequence, the CLT is not applicable and the requirement of a high cooperation between these elements leads to slower decays of the relaxation of fluctuations which are often modeled by empirical rheological power-laws [12,13]. One systematic way of dealing with these types of memory effects is to replace the first order time derivative in the conventional hydrodynamic (transport) equations by a fractional time derivative which is interpreted as the memory or the after-effect of the underlying stochastic process [14]. These effects on correlation functions have been studied for some complex fluids [9,14,15].
The purpose of this work is to study and compare the effect of time fractional derivatives on the correlation function of the transverse velocity fluctuations of a (homogeneous) viscoelastic fluid by using a GLE. We consider a long-ranged (power-law) viscoelastic memory and a noise with a long-range (power-law) auto-correlation. More specifically, we first evaluate the transverse velocity fluctuations correlation function for conventional time derivatives and then introduce time fractional derivatives in the equation of motion of this correlation to calculate the corresponding fractional correlation function [16]. Since for finite frequencies the imaginary part of the dynamic shear modulus (loss modulus), G (ω), of the dynamic viscosity η(ω) can be expressed in terms of the time-correlation function of the transverse velocity fluctuations, we compare the moduli for the non-fractional and fractional cases and find that the fractional modulus may be three times larger than the non-fractional one.
To this end the article is organized as follows. In Section 2 we set up the GLE for the dynamics of internal fluctuations of a viscoelastic fluid. Then an analytic exact expression for the one-time non-fractional correlation function for transverse velocity fluctuations (NF) is derived for power-law viscoelasticity. In Sections 3 and 4, we introduce both, noise and fractional time derivatives into the GLE and the fractional temporal transverse velocity fluctuations correlation function, where → k is the wavevector, is calculated analytically from them. We find that its magnitude is always lower than the non-fractional one and decays more rapidly. The relationship between G (ω) and C F → k , t is calculated analytically. The difference between the values of G (ω) for two specific viscoelastic fluids for accessible ranges of frequencies to our calculations is quantified. We find that the fractional effects on this measurable modulus may be as large as~300% when compared with its non-fractional values. The fact that the dynamic shear modulus is related to the light scattering spectrum suggests that the measure of this property might be used as a suitable test to assess the effects of temporal fractional derivatives on a measurable property. Finally, in Section 5 we summarize the main results of our approach and emphasizes that the eventual validity of our model calculations can only come from experiments.

Model Formulation
The deformation of spatially homogeneous viscoelastic liquids near equilibrium is described by the linear response theory [17]. In this regime, the most general constitutive equation for the linear stress-strain relation is of the form [18], The homogeneous character of these fluids comes from the assumption that the two independent, scalar moduli, the shear G(t) and the bulk (compressional) K(t), may only depend on time. Here, → r is the position vector → r = (x, y, z), σ ij → r , t is the symmetric stress tensor, p → r , t is the pressure and . γ ij → r , t is the rate of strain tensor, where the upper dot denotes the time derivative and Einstein's summation convention for repeated indices is implied.
Consistency with linear response requires to linearize the hydrodynamic equations in the small deviations with respect to a reference equilibrium state identified by the subscript 0. This yields where we have taken into account Galilei invariance for δv i . The complete system of linearized hydrodynamic equations for δp and → v then turns out to be [9,14,19], These equations are further simplified by choosing the direction of the z-axis as the longitudinal component and by separating where i = x, y, identifies the two transverse components. In what follows we shall only consider one of them which will be denoted as v

Hydrodynamic Fluctuations
One of the simplest formulations to describe fluctuations in fluids near equilibrium is the GLE with additive Gaussian random forces [4]. According to this approach, the stochastic dynamic of v is described by where the additive fluctuating force f (t) is defined as a Gaussian, stationary, stochastic process with zero mean, f (t) = 0 and with a, so far, arbitrary (long range) correlation C(t) subject to the condition Here, the angular brackets denote an average over both, the realizations of the noise and over an equilibrium ensemble of initial conditions. The rationale and an experimental evidence for assuming a Gaussian noise are as follows. It is well known that this noise describes the fluctuations around any equilibrium state and are dealt with in statistical mechanics with a variety of standard methods [5,20]. But when this state is (slightly) perturbed by changing the initial constraints in such a way that the state of the fluid remains within the linear response regime, the system relaxes to a new equilibrium state. The dynamics of the relaxation of the fluctuation is described by the GLE, Equation (6), and it is adequate to model them by a Gaussian process. This assumption has been experimentally shown to be appropriate for other complex systems involving the motion of tracers suspended in a fluid of swimming microorganisms [21]. In this system, the displacement of the tracers has a self-similar probability density function (pdf) with a Gaussian spatial effect which can be modeled by using a fractional diffusing equation [22,23]. On this basis, it seems reasonable to consider a Gaussian noise with a long-time correlation in Equation (6).
We define the combined Fourier-Laplace transform of an arbitrary field A → r , t by where s = iω is the Laplace variable. In what follows, the caret Â will denote its Laplace or Fourier transform with respect to one of its variables, whereas the tilde A will indicate a transform with respect to both. Thus, from Equation (6) whereĜ(s) is the Laplace transform of the so far arbitrary shear modulus G(t).

Power Law Viscoelasticity
There are many classes of materials in which the stress relaxation following a step strain is not close to an exponential, but is best represented by a power law in time, G(t) ≈ t −β . Examples of such materials-power law materials-include physically crossed-linked polymers, soft glassy materials and hydrogels. Non-exponential stress relaxation in the time domain also implies power law behavior in the viscoelastic storage modulus, G (ω), and in the loss modulus, G (ω), measured in the frequency domain by using small-amplitude oscillatory shear deformations. This broad spectral response is indicative of the wide range of distinctive relaxation processes available to the microstructural elements that compose the material, and there is no single relaxation time [24]. Let us assume that the viscoelasticity of the fluid is well represented by a long-range power-law rheological equation of state, i.e., where G 0 denotes the zero-frequency shear viscosity. Then where Γ(x) is the Gamma function. The parameter λ, measures the degree of viscoelasticity of the flow field; a low λ implies a weakly elastic flow field, whereas a large λ indicates an exceedingly elastic one.

Transverse Velocity Correlation
The single time auto-correlation function of the transverse velocity fluctuations iŝ where v 0 ≡v * → k and the asterisk denotes complex conjugation. The notation indicates the following.
Take a certain real constant value v 0 at t = 0, calculate the average over the realizations of the noise f (t) conditional on the given v 0 . Multiply it by v 0 and average this product over the values of v 0 , as they occur in the initial equilibrium distribution. Therefore, from Equations (10) and (13) it follows that Consistently with representing the viscoelasticity by Equation (11), we assume that the auto-correlation of the noise, Equation (7), is also a long-range power-law, Then, for given v 0 Equation (10) yields The inverse Laplace transform of Equation (16) is well defined and is given by [25] where E q (−bt q ) is the Mittag-Leffler function and R q,ν (a, t) denotes the special function defined by the infinite series Note that R q,ν reduces to the exponential function e at when q = 1 and ν = 0 (see Equation (40) in Reference [25]), i.e., R 1,0 (a, t) = e at . Therefore, the conditional average is From Equations (13) and (19), it then follows that the normalized non-fractional (NF) transverse velocity correlation function for power-law viscoelasticity is given bŷ The explicit dependence of q ≡ 2 − λ and b ≡ ρ −1 0 G 0 Γ(1 − λ)k 2 on λ shows that the dynamics ofĈ NF is indeed affected by the viscoelasticity of the fluid.
To examine the behavior of Equation (20) quantitatively, we consider two specific viscoelastic fluids, namely, silicon oil (S 2 ) and a solution of 0.02% separan MG500+2% water in glucose (E 1 ). Chhabra et al. [26] have studied the rheological properties of these fluids and according to their shear stress and normal stress data, S 2 would be classified as a weakly elastic fluid with a small λ, whereas E 1 is exceedingly elastic and has a large λ. Figure 1 shows the plot ofĈ NF ( → k , t) as a function of time t, as given by Equation (20), for the small values λ = 0.03, 0.06, which would correspond to a fluid like S 2 . This figure shows that as λ increases and the viscoelasticity decreases, the amplitude and the range of the correlation also decrease. Actually, the same behavior is observed in Figure 2 for E1, as λ increases from λ = 0.3 to λ = 0.4. However, in this case ˆ( , )  NF C k t decays three orders of magnitude faster than for S2.

Time Fractional Derivatives
In Equation (6), we now replace ∂/∂t by a Caputo left handed fractional time derivative (LHCD) 0+ D μ defined by [28][29][30], This figure shows that as λ increases and the viscoelasticity decreases, the amplitude and the range of the correlation also decrease. Actually, the same behavior is observed in Figure 2 for E 1 , as λ increases from λ = 0.3 to λ = 0.4. However, in this caseĈ NF ( → k , t) decays three orders of magnitude faster than for S 2 . This figure shows that as λ increases and the viscoelasticity decreases, the amplitude and the range of the correlation also decrease. Actually, the same behavior is observed in Figure 2

Time Fractional Derivatives
In Equation (6)

Time Fractional Derivatives
In Equation (6), we now replace ∂/∂t by a Caputo left handed fractional time derivative (LHCD) D µ 0+ defined by [28][29][30], where µ is the order of the derivative, m is an integer such that m − 1 < µ < m and v (m) ≡ ∂ m v/∂t m . The order µ of the fractional derivative should be chosen within the interval 0< µ <1 (and m = 1), because in this way, the integral in Equation (21) takes into account the contribution of the past values of this first order non-fractional derivative, and Equation (6) becomes After ensemble averaging and Fourier transforming with respect to x, this equation reads Since the Laplace transform of a Caputo derivative is given by [31] L and since for the power-law viscoelasticityĜ(s) is given by Equation (12), from Equation (23) with By using Equation (17) to invert Equation (25), we finally arrive at the following fractional time transverse velocity correlation function for power-law viscoelasticitŷ which is also given by a Mittag-Leffler function with different parameters. This correlation function is plotted in Figure 3 for S 2 and for the same parameter values used in Figure 1. where μ is the order of the derivative, m is an integer such that The order μ of the fractional derivative should be chosen within the interval 0< μ <1 (and m = 1), because in this way, the integral in Equation (21) takes into account the contribution of the past values of this first order non-fractional derivative, and Equation (6) becomes After ensemble averaging and Fourier transforming with respect to x, this equation reads Since the Laplace transform of a Caputo derivative is given by [31] ( ) (24) and since for the power-law viscoelasticity ( ) G s is given by Equation (12), from Equation (23) By using Equation (17) to invert Equation (25), we finally arrive at the following fractional time transverse velocity correlation function for power-law viscoelasticity ( ) (27) which is also given by a Mittag-Leffler function with different parameters. This correlation function is plotted in Figure 3 for S2 and for the same parameter values used in Figure 1.  The behavior ofĈ F ( → k , t) for E 1 is shown in Figure 4.
Entropy 2018, 20, 28 8 of 13 The behavior of ˆ( , )  F C k t for E1 is shown in Figure 4.

Dynamic Shear Modulus
If Equation (14) is compared with the corresponding one for a Newtonian fluid, namely [32], ( )  (28) where ηs is the shear viscosity, one can see that ( ) G s plays the role of a dynamic shear viscosity. It is useful to define the dynamic shear modulus as where its real ′ G and imaginary parts ′′ G are given, respectively, by the sine-Fourier and cosine-

Dynamic Shear Modulus
If Equation (14) is compared with the corresponding one for a Newtonian fluid, namely [32], where η s is the shear viscosity, one can see thatĜ(s) plays the role of a dynamic shear viscosity. It is useful to define the dynamic shear modulus as where its real G and imaginary parts G are given, respectively, by the sine-Fourier and cosine-Fourier transforms shown below. The real and imaginary parts of the dynamic viscosity η(ω) are related to G (ω) and G (ω) by η (ω) = ω −1 G (ω) and η (ω) = ω −1 G (ω), respectively. Note that in the limit of vanishing frequency, ω→0, G (ω) → 0 and η (ω) → 0, whereas G (ω) → 1 and η (ω) → η s . For finite frequencies G (ω) or η (ω) can be expressed in terms of the time-correlation function of the transverse velocity by setting s = iω and by equating the imaginary parts of the equation to obtain [18] where A → k , ω and B → k , ω are given, respectively, by It should be emphasized that Equation (30) provides a method of calculating G (ω) if the time-correlation function of the transverse velocity can be calculated from a model, or if it can be measured by some experimental technique [18]. By inserting Equations (20) and (27) into Equations (31) and (32), and by substituting the result into Equation (30) we arrive at a complicated but analytic expressions for G NF (ω) and G F (ω), which we do not write down because it is not necessary for our purpose. For S 2 this yields the solid curve for G NF (ω) and the dotted curve for G F (ω) in Figure 5.
It should be emphasized that Equation (30)   G ω . A feature of our model that shows that fractional effects in this modulus is a large effect that might be measurable. A similar behavior for these moduli is obtained for E1 as shown in Figure 6. However, the difference between the fractional and non-fractional results can be better quantified if we plot the ratio R ( ) ( ) ( ) ; , 1 ; , ; , 1 of the fractional to the non-fractional moduli, as shown in Figures 7 and 8, for S2 and E1, respectively. For the frequency interval considered this figure shows that G F (ω) is always larger than G NF (ω). A feature of our model that shows that fractional effects in this modulus is a large effect that might be measurable. A similar behavior for these moduli is obtained for E 1 as shown in Figure 6. However, the difference between the fractional and non-fractional results can be better quantified if we plot the ratio R of the fractional to the non-fractional moduli, as shown in Figures 7 and 8, for S 2 and E 1 , respectively. Entropy 2018, 20, 28 10 of 13 The curve in Figure 7 shows These results indicate that the relative changes in R, are not a small effect and might be measurable. The same behavior of R is observed for S2 in Figure 8.  The curve in Figure 7 shows R(ω; λ = 0.3, µ = 0.95) for E 1 ; R is larger than one for the frequency intervals considered. Note, for instance, that R has a maximum value of about R~3.4 at ω = 8.3 × 10 7 s −1 and a minimum value of R~1.63 for ω = 1 × 10 5 s −1 . This means that fractional fluctuations have a large effect on the transverse velocity fluctuations correlation than on the velocity correlation. These results indicate that the relative changes in R, are not a small effect and might be measurable. The same behavior of R is observed for S 2 in Figure 8. The curve in Figure 7 shows These results indicate that the relative changes in R, are not a small effect and might be measurable. The same behavior of R is observed for S2 in Figure 8.

Discussion
In this work, we have analyzed the effects produced on the transverse velocity fluctuations correlation function ˆ( , )  NF C k t , due to the presence of fractional temporal derivatives in the dynamics of these fluctuations. More specifically, we considered the case where the liquid possesses a longrange power viscoelasticity given by Equation (11) [29] or Grünwald-Letnikov [31]. The reason for choosing LHCD in this work lies in the physical initial conditions necessary to obtain a particular solution of Equation (22). When the LHCD is used in (22) to obtain its Laplace transform, the initial value ( ) ,0 v k and the initial value of the integer-order derivatives should be known, but this is precisely the type of initial condition that can be specified or controlled in the relaxation process described by Equation (22) where n is an integer such that 1 − < < n n α , would have to be known; however, for real physical systems like the viscoelastic liquid considered in this work, they are unknown. Thus, the RL fractional time derivative was discarded and the physically proper choice was the LHCD [33]. It should be remarked

Discussion
In this work, we have analyzed the effects produced on the transverse velocity fluctuations correlation functionĈ NF ( → k , t), due to the presence of fractional temporal derivatives in the dynamics of these fluctuations. More specifically, we considered the case where the liquid possesses a long-range power viscoelasticity given by Equation (11) and the transport equations have long-range correlated stochastic terms. The interplay between the fluctuations arising from these two features motivated the introduction of fractional time derivatives into the hydrodynamic equations, since in this way it is possible to take into account the fact that the transverse velocity fluctuations evolve on vastly different time scales.
The most important results of our analysis are the analytic expressions for the different transverse velocity correlation functions,Ĉ NF ( → k , t) andĈ F ( → k , t) for the long-range power viscoelasticity and given, respectively, by Equations (20) and (27), both expressed in terms of Mittag-Leffler functions. Since for finite frequencies, G (ω) can be expressed in terms of the time-correlation function of the transverse velocity, these fractional effects on the correlations do affect the frequency behavior of this measurable property of the system. To our knowledge, this issue has been scarcely considered in the literature of viscoelastic fluids.
Another issue that deserves further comments is the justification of the choice of the Caputo's left-handed time derivative (LHCD), in spite of the fact that there exist other fractional time derivatives, such as those of Riemann-Liouville (RL) [29] or Grünwald-Letnikov [31]. The reason for choosing LHCD in this work lies in the physical initial conditions necessary to obtain a particular solution of Equation (22). When the LHCD is used in (22) to obtain its Laplace transform, the initial valuev(k, 0) and the initial value of the integer-order derivativesv (m) (k, t = 0), 1 ≤ m ≤ α, should be known, but this is precisely the type of initial condition that can be specified or controlled in the relaxation process described by Equation (22); thus, it is consistent to use the LHCD in this equation. If we had used the RL time fractional derivative to obtain the Laplace transform (25), the initial values of the fractional derivatives D µ 0+ v → r , t with β = α − k − 1, k = 0, 1, . . . , n − 1, where n is an integer such that n − 1 < α < n, would have to be known; however, for real physical systems like the viscoelastic liquid considered in this work, they are unknown. Thus, the RL fractional time derivative was discarded and the physically proper choice was the LHCD [33]. It should be remarked that for relaxation processes quite different from the study of the time evolution of fluctuations in viscoelasticity, other time fractional derivatives, different from that of Caputo, have been proposed. For example, in Ref. [34], a fractional anomalous-growth equation employing a Riemann-Liouville derivative has been introduced to describe nucleation-and-growth processes in superconducting materials like weakly viscoelastic cuprate systems where a relationship between the order of the fractional derivative and the fractal dimension of an anomalous random walk process was found. Another example is the use of Weyl derivatives to describe nucleation-and-growth processes in model lipid membranes [35]. In this latter example, a Caputo derivative of an order centering around~1/2 is proposed to explain the plausible multitude of time scale of relevance. This form of introducing fractional derivatives has been also explored and its limitations have been discussed for a binary mixture in [33].
Finally, it should be emphasized that we are not aware of any specific experimental results to compare with the predictions of our model, and therefore, it is not possible to conclude from our analysis if this fractional behavior of various correlation functions and elastic moduli may be measurable; this is an open issue that remains to be assessed.