Original Research ARTICLE
Linear Viscoelastic Responses: The Prony Decomposition Naturally Leads Into the Caputo-Fabrizio Fractional Operator
- Department of Chemical Engineering, University of Chemical Technology and Metallurgy, Sofia, Bulgaria
The study addresses the physical background and modeling of linear viscoelastic response functions and their reasonable relationships to the Caputo-Fabrizio fractional operator via the Prony (Dirichlet series) series decomposition. The problem of interconversion with power-law and exponential (single and multi-term functions) has been discussed. Special attentions have been paid on the Prony series decomposition approach, the related interconversion problems and the expression of the viscoelastic constitutive equations in terms of Caputo-Fabrizio fractional operator.
Nothing is done in a vacuum; we must all stand on our forefathers to better ourselves and the world around us.
(Sir Isaac Newton, in his letter to rival Robert Hooke in 1676)
Science is built up of facts, as a house is built of stones; but an accumulation of facts is no more a science than a heap of stones is a house.
This article addresses the physical background of modeling of dissipative phenomena, precisely response functions such as stress-strain relationships in the framework of the linear viscoelasticity and their reasonable relationships to the Caputo-Fabrizio fractional operator  via the Prony (Dirichlet series) series decomposition .
The appearance of the new definition of fractional derivatives with non-singular kernels was provoked by needs to model dissipative transport processes in many new materials appearing in modern technologies . Recently the main achievements, especially results related to diffusion problems, were analyzed Hristov in  and we will avoid the thorough browsing and comments of published results (see also the rich list of reference in Hristov ). In the context of diffusion problems it was demonstrated that the Caputo-Fabrizio operator appears naturally in diffusion models [3, 4] when the flux-gradient relationship is expressed by a Jeffrey relaxation kernel and the Maxwell-Cattaneo concept of flux.
At the same time the Caputo-Fabrizio fractional operator was criticized [5–7] with points of view based on the classical fractional calculus with singular power-law memory kernel and with examples from the signal processing  and to some extent touching formal rheological relationship . Despite the mathematical exactness of the counterexamples they do not focus the attention on the physical basis leading to exponential memory kernels. As a matter of fact, we need a clear answer what really this new operator models and how it appears in the modeling of physical problems, despite the fact that some properties known from the classical fractional calculus, such as the index law  , are not satisfied.
1.1. Motivation of This Study
The analysis done here and the principle task are oriented to modeling of viscoelastic constitutive relationships in terms of Caputo-Fabrizio operator naturally appearing through Prony series decomposition  of stress relaxation functions and not obeying power-law behaviors. That is, we focus on viscoelastic materials with dynamics which cannot be modeled with the classical fractional derivatives of Riemann-Liouville and Caputo .
This introduction avoids huge citations of works on Caputo-Fabrizio operators since without understanding of the physical background and the logic of its appearance any analysis of models created by formalistic fractionalization (see the comments in Hristov et al. ) is unproductive. Moreover, to the author personal experience as editor many manuscripts devoted to applications of the Caputo-Fabrizio derivative in formalistically fractionalized existing models are directly rejected by the reviewers with motivations based of the opinions in the criticizing articles [5–7]. This situation resembles that in the Catch 22 movie without perspective for escape. The existing situation is a consequence of some main reasons: (1) The Caputo-Fabrizio operator does not hold some properties such as the semi-groups, which are existing with the classical fractional derivatives with singular kernels, the strange form of the associated fractional integral and these issues cast doubts when it is applied inadequately (in blind manner) to various functions in a manner known from the integer-order calculus. (2) The formalistic approach by simple replacements of integer or fractional order (with singular kernels) derivatives with the Caputo-Fabrizio operator in existing models without taking into account the physics behind. (3) Last but not least, the human factor of author's rivaling which is dividing the scientific society into competing groups rejecting the achievements of each other. All these elements of the current situation create a discouraging atmosphere and disbelieve that this new fractional operator really cannot model natural phenomena and actually stops the further research, generally among the young researchers highly aspiring publications of submitted manuscripts.
The deep physics behind the fractional operator with exponential kernel (and the author's long time experience in science) motivated this study and the efforts are oriented to show that the existing knowledge and models, as well as techniques of data treatment, in the framework of linear viscoelasticity, lead naturally to formulation of the Caputo-Fabrizio fractional operator. This is in the context of the Sir Isaac Newton quote at the beginning of the article: the steps ahead on the shoulder of existing facts and results on the road to creations of new information are natural ways and actually the exciting moments in the beautiful journey in the world of science.
1.2. Aim and Paper Organization
This article is organized as follows: section 3 presents briefly the main properties of the Caputo-Fabrizio operator that will be used further in the analysis of the viscoelastic problems. Section 3 addresses the constitutive equations of viscoelasticity based on the Boltzmann superposition principle  and fading memory approach incorporated the hereditary stress and strain integrals. The principle problems of the quality and choice of the response function are discussed and the main properties requited are outlined in section 3.3. The interconversion of relaxation and creep response functions is discussed in section 4: the linear and non-linear scale-invariant (power-law) (related to application of Riemann-Liouville and responses are analyzed. Section 5 focuses on decompositions of experimental response functions by Prony series leading to discrete spectra and the related interconversions by examples with single term, two-terms and multi-term responses of exponential type, as well as interconversion of relaxation and creep compliance expressed as Prony series. Section 6 demonstrates how the constitutive equations (based on the Maxwell model) can be expressed in terms of the Caputo-Fabrizio fractional operator. The subsection 6.3 demonstrates that approximation of the response function by Bessel functions (of Maxwell-like materials) of first kind and expressed as infinite Dirichlet series naturally leads to incorporation of the Caputo-Fabrizio operator in the constitutive equations. Section 7 demodulates briefly how constitutive relationship (following the idea of Bagley and Torvik) with two fractional operators of Caputo-Fabrizio (of different orders) can be formulated.
2. Caputo -Fabrizio Operator
The Caputo -Fabrizio operator is defined as 
In Equation (1) M(α) is a normalization function such that M(0) = M(1) = 1. This definition is of Caputo-type because there is a convolution of the derivative df(t)/dt. The explanations in  relate the development of Equation (1) to the classical Caputo derivative  by mechanistic replacements of the singular kernel (by a non-singular exponential kernel) and the normalization function (see the explanations in Hristov .
The integration by parts in Equation (1) results in an alternative form Caputo et al. 
Futher, the integer order differentiation of follows the rule 
Hereafter we will use the definition (Equation 6).
2.2. Laplace Transform
The Laplace transform of with a = 0 has the following Laplace transform LT given with p variable  taking into account the general rule of Laplace transform of a convolution, namely
2.3. Fractional Derivative of Elementary and Transcendental Functions
Linear function f(t) = Ct 
For f(t) = Ctβ and β > 0 the fractional derivative is
For α = 1 the expression (Equation 10) reduces to the classical (integer-order) result Cβtβ−1.
For α → 1 the second exponential term in the nominator of Equation (11) goes to zero and therefore .
For β < 0 we have
For α → 1 the second exponential term in the nominator of Equation (12) goes to zero and therefore .
2.4. Caputo-Fabrizio Fractional Operator: Determination of the Fractional Parameter
In the Caputo-Fabrizio operator there is formal ambiguity because the stretched time is multiplied by a dimensional factor α/(1 − α) which should have a dimension s−1. This contradicts the definition of the fractional order (parameter) because physically α is dimensionless. The answer given in [3, 13] resolved the problem by nondimesionalization of the exponential function by help of characteristic time scale of the relaxation process t0 (the maximum time of the experiment in the sense of the rheological tests discussed here) , namely
and relates it to data that can be really recovered from experimental data, such the relaxation and retardation times (see the sequel)
The relationship (Equation 14) says that for τ/t0 = 1 we get α = 1/2. Further, depending on the ratio τ/t0 we may have fractional orders roughly arranged in two groups : a) when 0 ≤ τ/t0 ≤ 1 we have fractional orders α ∈ [0.5, 1.0] and the relaxation time τ is shorter than the macroscopic process observation time t0, and b) 1 ≤ τ/t0 < ∞ the fractional orders are α∈(0, 0.5] since the relaxation time τ is larger than the macroscopic process observation time t0. Qualitatively, for τ/t0 < 1 the relaxations could be considered as fast (rapid) relaxations, while τ/t0>1 is related to slow relaxations (see in Hristov  the comments of numerical values of Prony decomposition and relaxation times)
Last to this point, but not least, as it was commented in Hristov , the ratio τ/t0 is not integer and expressing the memory kernel as where we get a fractional operator. Therefore, the memory kernel of the Caputo-Fabrizio operator is controlled by a non-integer parameter in the context of what is needed to say that this operator is fractional, despite the fact that it does not repeat exactly the properties of the Classical Riemann-Liouville and Caputo derivatives.
3. Constitutive Equations of Viscoelasticity: Fading Memory Approach
3.1. Boltzmann Superposition Principle and Fading Memory Concept
This definition is actually the Boltzmann linear superposition functional (Equation 16))
relating the present state of the flux to its history [9, 16–18] through the influence function (memory kernel) R(t, z) during the time interval defined by τ. In Equation (16) m and λ are transport coefficients (diffusivities) with real physical meanings as it will be demonstrated in the sequel.
The memory function could be unbounded and scale invariant such as R(t, τ) = t−μ with integration singularity at or bounded and not scale invariant such as R(t, τ) = e−t/τ (see the analysis in Hristov  and the comments further in this article).
3.2. Stress-Strain Viscoelasticity Response and Hereditary Integral Construction and Response Functions
The linear theory of elasticity  (chapter 1) considers the stress in a sheared solid body as a quantity proportional to the shear, while in the liquids the shearing stresses are proportional to the rate of shear . Most solid materials, for example polymers, compromising both effects are called viscoelastic. When a slab of solid material under a shearing motion caused by a step change in the stress load applied to it (Heaviside unit step function H0(t)) exhibits a strain (in one dimension)  (chapter 1)
With perfect elastic behavior of the body  σ(t) = σ0H0(t) for t > 0. On the contrary, in an ideal viscous fluid the stress is infinite and for t > 0 and the strain is ε(t) = (σ0/η)t, thus introducing the coefficient of viscosity η. Real materials do not shear with infinite speeds that is the reason of the concept of a finite relaxation time τ . Precisely, in solids the stresses attain finite values for long times. In contrast, in viscous fluids the stresses approach zero.
The task of this study addresses the functional representations of the viscoelastic material responses, that is : 1) the stress relaxation function R(k, t), that is the stress history due to a shear step of size ε, and 2) the creep function C(t) (shear history) due to unit stress σ applied. In the linear viscoelastic theory  the responses can be approximated as R(ε, t) = G(t)ε + O(ε3) and C(σ, t) = J(t)σ + O(σ3), thus defining the linear stress relaxation modulus G(t) and the linear creep compliance J(t). The common mechanistic models explaining the behavior of viscoelastic materials utilize linear springs and dashpots coming from the Maxwell interpretation [20–22] and modeling the pure elastic behavior and the pure viscous state, correspondingly.
Superposition of single-step material responses results in functional relationships of stress and strain in the linear viscoelasticity concept  incorporating the a time lag in G(t) and J(t) through hereditary stress (Equation 18) and creep integrals (Equation 19), namely
Both σ(t) a and ε(t) are causal functions and therefore the lower terminals in Equations (18, 19) are set at t = 0.
Applying the fading memory concept [9, 20] it is possible to relate the instantaneous responses G∞ and J∞ corresponding to equilibrium states (for long times) when the effects of memory terms (convolution integrals) fade out, namely
3.3. Response Function: General Properties and Requirements
For a linear isotropic viscoelastic body exerting uniaxial stress the Boltzmann superposition principle formulates the constitutive equation relating the strain responses by help of the following hereditary integral 
Coleman and Noll, in their article on foundation principles of linear viscoelasticity , point out that there is no universal approach in definition of unique relaxation function R(s) [the fading memory R(s) in Equation (15) its is equivalent to G(s) in Equation (22) . Despite this, two principle features are required to be obeyed by R(s):
1) R(s) is defined for 0 ≤ s < ∞ and R(s) > 0, and
2) R(s) decays monotonically to zero for large s , that is .
a) Generalized power-law memory function R(s) = sr(s+1)−μ of order r when r < μ. For r = 0 the memory function R(s) is unbounded (singular) at t0+ and this scale invariant kernel forms the constructions of the classical fractional integrals and derivatives .
b) Ordinary memory kernel , bounded at t0+ as exponential function e−βs , which is not invariant with respect to the time scale. The definition of Caputo-Fabrizio fractional operator utilizes this memory kernel.
The adequate mathematical structure in the construction of the relaxation (memory) function is the main problem that should be resolved in the modeling of viscoelastic responses. Generally, the relaxation function should adequately describe the natural process and therefore its structure should be tested (defined)by a data conversion algorithm (by fitting experimental data). Despite this intuitive and to greater extent logical approach the response function should satisfy some requirement defined in Garbarski  and Winter , and summarized in Table 1; with some author's comments related to application of fractional operators in the viscoelastic constitutive equations. Further, in the response of viscoelastic material to strain excitation, the stress relaxation spectrum (relaxation modulus) provides complete information concerning the time-dependent part of the material response. In the opposite situation, with materials undergoing stress excitations the retardation spectra (the compliances)  provide the required information. Hence, if the spectra are known it is possible to calculate the response of any excitation .
Besides, the thermodynamic theory of linear viscoelastic materials was developed extensively in the sixties of the last century . In this context, the second law of thermodynamics yields severe restrictions on the constitutive properties [32–34]. This problem is beyond the scope of the present analysis but we may say that when the constitutive equation (Equation 20) contains exponential kernel it is compatible to second law of thermodynamics .
The principle relaxation functions used and the applicability of the aforementioned requirements will be analyzed next.
4. Interconversion of Relaxation and Creep: Problem and Power-law Responses
4.1. Interconversion Problem
In the linear viscoelastic materials the relationships between the creep compliance J(t) and the relaxation response G(t) are expressed as convolution integrals in accordance with the Boltzmann superposition principle  (see Equations 20, 21)
With the Laplace transforms of Equations (20, 21) we get equivalent relationships in the p-space (p is the transform variable)
Hence, we may recast Equation (24) as
Equating the right-hand sides of Equation (25) we get an implicit fundamental relationship
with the constraints that
Explicit form of Equation (26) can be obtained if the analytical forms of either G(t) orJ(t) is known. Examples of interconversion related to power-law and exponential memory kernels are discussed in the sequel.
4.2. Scale-Invariant (Power-Law) Response
4.2.1. Physical Preliminaries Related to the Power-Law Response
Relatively short-time relaxation history modeled by time-dependent power-law can be observed in the time evolution of G(t) from zero to the end of the time of observation t0 : commonly at the beginning of the relaxation process (short times) and along the long tail (long times). The power law is scale-invariant. Materials exhibiting such behavior are called power-law viscoelastic materials [35–42] and can be easily detected by the linear behaviors of J(t) ~ tp and G(t) ~ t−p in logarithmic scales . Actually, we have to memorize that this type of relaxation was modeled in the article of Caputo (1967) where the construction of the Caputo derivative was conceived . The power-law functional relationships of the relaxation and the compliance, actually are the reasons to recognize the integral and derivatives of the classical fractional calculus  as adequate modeling tools since they are based on the same type of kernels [39, 42, 45], known also as weak singular kernels. As commented by Tschoegl  the scaling relationship G(t) ~ t−p is a fractional version of the Trouton law . Here, we have to mention that Metzler et al.  discussing relaxations of filled /polymers stated that non-exponential (non-Debay) relaxation implies memory (Sic!). We may consider this declarative opinion (the results of Glockle and Nonnemacher ) as adequate to the situation in the 90s of the last century, when fractional calculus was only related to power-law kernels .
Moreover, as a valuable analysis showing when the power-law is the adequate response function we refer to Chapter 5 of Findley et al.  where it is clearly demonstrated that for many materials the power-law response function tμ with μ < 0.5  is a good short-time approximation for materials such as plastics, metals and concrete. For short-time loading the creep of many different rigid plastics with sufficient accuracy can be presented as linear power-law , where μ is stress-independent and nearly temperature-independent, too. The same approach following from the linear approximation approach of Pipkin  is commented in 3.2 . In this case the corresponding retardation spectrum is Findley et al. 
The function (Equation 28) may be considered as continuum spectrum of material retardation times proportional to t−α . From this position, the step from Equation (28) toward modeling of the relaxation processes by the classical Riemann-Liouville derivative is straightforward. In this context, Bagley and Torvik  using (Equation 28) demonstrated that following constitutive equation holds
After application of the Leibniz rule to Equation (29) the result is
For ε(0) = 0, since no strain at t = 0 exists we get
Then, for G1 = G0 = const. the stress relaxation modulus G(t) and the creep compliance J(t) are
4.2.2. Example 1. Interconversion of Power-Law Relaxation and Creep: Linear Case
Here, we consider again the power-law response in the context of the response function interconversion and in the sense of the relationships (Equations 26, 27). With G(t) = At−μ , where A is a data fitting coefficient and 0 < μ < 1 , the Laplace transform in Equation (26) provides a power-law creep (Equation 33) compliance 
Following Hanyga, neither experimental nor theoretical reasons lead to the assumption that the creep rate function and the stress relaxation function are bounded and regular at t = 0 . The same point of view is valid for the unbounded kernels singular at t = 0 . In the seminal study of Boltzmann  the value of α = 1 was suggested. The condition 0 < α < 1 comes from the fact that for α > 1 we get an infinite propagation speed . Actually, the idea to apply fractional calculus, as is demonstrated by the model developed in the foundation studies of Scott-Blair [51–53] and Bagley and Torvik  comes from the experimental (empirical) findings of Nutting [54, 55]: precisely, for some viscoelastic materials the power-law relationship [ε(t)/Const] ~ σntα is satisfied. Consequently, we may obtain (Equation 32).
Since the power-law relaxation is not the main problem discussed in this work, at the end of this point, we refer to the study of Carillo and Giorgi  (section 4 in this chapter) where it is clearly demonstrated that the response of the material may be for short-time or as long tail (for long times) modeled by time-dependent power-law function. Here, we may repeat again (see the comments in Hristov , too) that, if the power-law is not exhibited by the material response functions then it is inadequate to apply the power-law memory kernels.
4.2.3. Example 2.Interconversion of Power-Law Relaxation and Creep: Non-linear Case
If the material exhibits a non-linear viscoelastic behavior, for instance, as commented by Lakes et al. , that is
The relaxation function G(t, ε) can be expressed as a sum of products, precisely sum of products of functions of time and functions of strain, that is Lakes and Vanderby  G(t, ε) = Gt(t)g(ε) . Assuming no interaction between the steps (immediate and delayed Heaviside steps) strains in the summation series and may write 
Consequently, the creep compliance is
and for the linear case we will obtain the relationship (Equation 26).
To have explicit form of these relationships, if we assume a power-law approximation in time  J(σ, t) = A(σ)tμ, that is , then we get
The nonlinear material exhibits a relaxation response which contains a sum of power-law terms, as given in equation (Equation 38) (see more comments in Lakes and Vanderby ).
Hereafter, non-linear behavior generally related to temperature-dependent and aging viscoelastic materials will not be discussed since the principle task of the present analysis is to show the correct origin and the physical background leading to the formulation of the Caputo-Fabrizio fractional operator.
5. Interconversion of Non Power-law Relaxation and Creep
5.1. Non Power-Law Response: Relaxation Curve Approximation
The selection of approximation function which would fit the experimental points depends on the type of materials and would be established by simple trial error process. In this context, two principle issues arise when experimental data for linear viscoelastic materials should be treated : (1) Appropriate relaxation curve approximation and (2) interconversion problem. The experimental data are commonly taken in the time-domain (or frequency-domain) tests. Now, we stress the attention on discrete spectrum approximation by Prony series  where the condition for the monotonicity of the functions (see the sequel) is satisfied.
5.2. Discrete Relaxation Spectrum by Prony Series Decompositions and Relevant Interconversions
Unlike the linear elasticity, where the material functions are related algebraically, the relationships in the linear viscoelasticity are time-dependent, as already was demonstrated in the previous section of the article. One known method is to find the unknown functions G(t) and J(t) by fitting the data points by a known function. However, when the focus is on the adequate construction of a fractional operator with a memory kernel obeying all mandatory conditions imposed on it, then the choice of approximating function is highly restricted. If the power-law is not the adequate choice then a series approximation is the more suitable solution. Since we address the construction of the Caputo-Fabrizio operator, thus it is natural to tackle data fitting by series of exponential terms (matching the function of the desired memory kernel) known as Prony series (also as Dirichlet series ).
5.2.1. Prony Series Decompositions
In Equations (39, 40) the parameters g∞ and gi are equilibrium (at large times) values and the relaxation moduli (stiffness), respectively, are constrained according to Brinson and Brinson 
The Prony series components have spectral strength gi and relaxation time τi. This is the so-called discrete Prony series representation known also as discrete relaxation spectrum [20, 65–69]. Besides, the series approximation (Equation 39) may be substituted in formulation of the materials law such as Maxwell, Kelvin Voigt, etc. The relaxation time τi associated with the ith element is related to the characteristic time of the spring-damper (dashpot) element and can be defined as the ratio of the viscosity over the elastic modulus, that is τi = gi/λi .
The first derivative of the Prony series (Equation 42) for t = 0 is finite, a fact irrespective of the number of terms used in the approximation.
As commented by Bradshaw and Brinson , an exact solution of interconversion problem is possible if the known function is approximated by Prony series where all coefficients are positive [using the forms (Equation 39, 40) ] and the basic interconversion relationship (Equation 26) is satisfied. In this case, by applying Laplace transform solutions, the resulting Prony series are analytically exact [70–73].
The Prony series approximations lead to the generalized Maxwell viscoelastic body (known also as Maxwell-Wiechert model) with Nϕ spring-dashpot elements in parallel. The Prony series decompositions are applicable to any viscoelastic models through their time-dependent shear and bulk moduli [74–76]. For a deep thermodynamic analysis of such type of materials termed viscoelastic solids of exponential type (VESET) we refer to works of Fabrizio et al. [32, 33].
A common step in approximation by the truncated exponential sums is the definition of preliminary stipulated decay rates βi , in a logarithmic scale that is , for Nϕ = 4 , for instance β1…4 = 100, 10, 1, 0.1 . This approach, is useful because the corresponding fractional parameters in the αi in this case can be easily calculated (see further Equation 65) as α1…4 = 0.009, 0.09, 0.5, 0.9, respectively. Alternatively, fixing the points (in normal or logarithmic scale) ti along the time axis we may define the corresponding βi = 1/τi.
The parameter estimation is important and the first step was done in the seminal work of Prony (1795)  and several algorithms have been developed among them: graphically by log−log plots [78, 79], least squares method [80–82], nonlinear optimization methods , genetic based algorithm , from the continuous relaxation time spectrum [68, 72], logarithm equidistant distribution of relaxation times (known as R- method) [68, 77, 85], quasi-linearization for multi-exponential decay curves fitting , etc.
The number of the exponential terms in Equation (39) depends on the accuracy of data fitting. Commonly fourth-order Prony series fit adequately the stress-relaxation data in cases of non-linear viscoelastic behaviors [62, 87–93], while long-term relaxation tests need 10-15 terms [76, 85, 94–100]. In general, the problem corresponds to identifications of the kernel of an integral Fredholm equation which, actually, is ill-posed problem [68, 101–103] requiring Tikhonov regularization . Comments on the required terms in the Prony decomposition are available in .
5.3. Examples of Interconversions With Exponential Terms
5.3.1. Example 3. Interconversion of a Single Exponential Model 
For a single exponential model (the Maxwell model) the relaxation and creep functions are
The condition (Equation 27) is automatically satisfied, while for t = 0 in (Equation 43) (and taking into account ;Equation 27) the results is
The Laplace transforms of (Equations 43, 44) lead to an equation that should be solved
Since (Equation 46) should be satisfied for all p it follows that τ1 = j1λ1 and this condition with (Equation 46) leads to Anderssen et al. 
184.108.40.206. Interconversion from relaxation to creep
220.127.116.11. Interconversion from creep to relaxation
More details related to the sensitivities and relative errors of these two interconversions are available in Anderssen et al. 
5.3.2. Example 4. Interconversion of a Double Exponential Model 
In this case the relaxation and creep functions are
The values of the relaxation times τi are zeros of
5.3.3. Example 5. Interconversion of Multi-Term Exponential Model 
By Prony decomposition the relaxation and creep functions can be generally presented through discrete spectra as Anderssen et al. 
The inverse times, that is the relaxation times τi and the retardation times λi satisfy
The advantage of this (Prony decomposition) approach is that the monotonicity is automatically satisfied 
5.3.4. Example 6. Interconversion of Prony Series of Relaxation and Creep
With Prony decomposition of the relaxation spectrum, the next step is to determine the creep modulus (compliance) J(t). The interconversion equation (Equation 56) simple gives [22, 70, 71, 105, 111, 112].
As a rule in the rheological studies, the discrete exponential models for G(t) and J(t), corresponding to associated relaxation and creep spectra H(τ) and L(τ)  (as sums of delta functions ) are used that assures the monotonic behavior automatically and we have
The constraint imposed [70, 105] by (Equations 26, 27) (see also Equation 56) allows H(τ) to be defined as function of L(τ) and vice versa. We will skip this general problem solved in  and to some extent demonstrated by Example 3 and 4, and will focus the attention on the discrete approximation of the relaxation (and compliance) spectra by Prony series.
The Prony series approximation of the stress relaxation (in a normalized form) is
The corresponding form of J(t) is
The inverse times, that is the relaxation times τi and the retardation times λi satisfy the inequalities (Equation 55). The attempts to solve (Equation 56) address many approaches for either J(t) or G(t) depending on the type of experimental data available [30, 106–109].
6. Caputo-Fabrizio Operator in the Constitutive Viscoelastic Equations
6.1. Relaxation Function in Terms of Caputo-Fabrizio Operator
Since we operate with a Prony series, that is with a finite sum of exponentials, then the inversion of the summation and the integral yields 
Thus, the memory effect from the convolution integral can be easy incorporated in each term of the Prony series ,namely
The finite sum of N+1 terms σi(t) leads directly to the generalized Maxwell model, that is at any time t we have
σi(t) is a product of the spring modulus Ei and its current strain ki(t), which to some extent could be considered as a hidden material variable. The clear physical meaning of this result is that the strain ki(t) at a given time t is expressed as convolution integral with an exponential kernel . Therefore, the model parameters that should be defined, following the approximation of the right-hand side of (Equation 62), are [14, 113] : the single separate spring stiffness E0 and the spring stiffness Ei as well as the relaxations time τi of each i−th Maxwell element.
6.2. Relationships of the Relaxation Time Spectrum and Fractional Order Spectrum
As it was mentioned at the beginning, the fractional parameter α is related to the dimensionless relation time as α = 1/(1 − τ/t0). Since we have a spectrum of relaxation times τi, then the spectrum of the fractional orders (parameters) is
Therefore, ki(t) can be expressed in a form related to the construction of the Caputo-Fabrizio operator, namely
6.2.1. Stress Relaxation in Terms of Caputo-Fabrizio Operator
Thus, the constitutive equation of the stress relaxation can be presented as Hristov 
Now, we turn on the determination of the spectrum of fractional orders αi = τi/t0 . From experimental data fittings there are a limited number of numerical values of relaxation times τi. Now, the principle problem at issue is the determination of the characteristic time t0. Since the experiments last limited times then we may assume that t0 equal the elapsed time te of the experiments. The literature data concerning Prony decompositions reveal (see the analysis in ) that the relaxation times form two groups (see the comments at the begging when the fractional order determination of the Caputo-Fabrizio operators was commented):: (1) relaxation times less then the elapsed time τi < te ⇒ τi/te ≤ 1 and (2) relaxation time greater than the elapsed time τi > te ⇒ τi/te ≥ 1 .
Consequently, addressing the fractional orders αi we have: i) fast relaxations for αi ∈ [0.5−1) corresponding to τi/te ≤ 1 and ii) slow relaxations for αi ∈ (0−0.5] when τi/te ≥ 1. Numerical examples supporting this estimates are reported in Hristov 
6.2.2. Creep Compliance in Terms of Caputo-Fabrizio Operator
By analogy of the stress relaxation expression in terms of Caputo-Fabrizio derivative, we may transform the Prony series decomposed compliance (Equation 54) as
where are the scaled (dimensionless) retardation times.
Hence, with the construction of the convolution integral describing the strain history we have
Thus, the strain history is
Then, the strain can be expressed as
Therefore,(Equation 73) defines a Caputo-Fabrizio operators with respect to σ(t) and the fractional order βi is related to the scaled retardation time ( from the spectrum) as
6.3. Description of Maxwell-Type Viscoelastic Media With Material Responses Modeled by Bessel Functions in Terms of Caputo-Fabrizio Operators
In the last two years (2016-2017) an interesting approach was developed by the group around Professor Mainardi [114–119] which could be considered as attempts to generalize the relaxation functions in the linear viscoelastic models, an approach also investigated in Colombaro et al.  and Guisti and Colombaro . The main idea comes from the possibility to represent the relaxation function in a viscoelastic Maxwell-type body by infinite discrete spectrum with times related to the zeros of Bessel functions of the first kind [118, 119]. Here we will present some key points of these studies,in the context of the ideas developed in this article, that is to show the incorporation of the Caputo-Fabrizio operators in the relaxation functions expressed as infinite Dirichlet series (or infinite Prony series as the authors of these studies defined them).
The analysis developed in Colombaro et al. , Colombaro et al. , Colombaro and Guisti , Colombaro et al. , Guisti and Mainardi , and Guisti and Mainardi  is based on the power series representation of the modified Bessel function of the fist kind as Abramowitz and Stegun 
with an asymptotic representation
Now, if the function Fν(t) has a Laplace transform 
It can be presented in the time domain as Guisti and Mainardi 
which is an expression as Dirichlet series, locally integrable, positive and increasing function for t > 0 , that is
because for s → ∞ from Equation (77) we have
If the is considered as a part of the Laplace transformation of the relaxation memory function Φν(t) ( here we use the original notations as that used in Guisti and Mainardi ), that is , we have
In the time domain the relaxation function Φν(t) is [recall equation (Equation 78)]
Now, let us turn on to the interconversion problem. From as it is suggested  that
In the time domain the relaxation modulus G(t) is
which follows from the interconversion relationships.
For further deep reading in this elegant mathematical studies we refer to Colombaro et al. , Colombaro and Guisti , Colombaro et al. , Guisti and Mainardi , and Guisti and Mainardi  where it is was demonstrated that the asymptotic behaviors of the viscoelastic responses in the so-called Bessel medium are (precisely in Colombaro et al. )
which actually follows from Equation (79).
The main idea to use a relaxation function expressed by modified Bessel functions of first kind comes from the possibility to calculate the sum of the infinite series of reciprocal positive zeros of the Bessel functions Jν , namely
With (Equation 89) the relaxation modulus becomes 
Hence, the relaxation modulus is represented as an infinite Dirichlet series which are absolutely converging, which actually represent a discrete spectrum where the stiffness of the nth element (spring-dashpot) is
(the prefix B means Bessel) and the relaxation times are . That is, in the context of the analysis done in this study, we have
Therefore, we get an expression in terms of an infinite sum of exponential functions, as in the Prony approximations, but now the relaxation times are defined by the zeros of the Bessel function. If a time scale t0 of modeled relaxation process exists, then nth dimensionless relaxation time is 0 < τn/t0 ≤ 1, as it was demonstrated several times with the Prony series. Now, let us consider that τn/t0 = (1 − βn)/βn ,where 0 < β < 1 and then we may re-write (Equation 91) as
Then, the construction the convolution integral of the stress relaxation with G(t) yields
Interchanging the orders of the summation and the integral in equation (Equation 93) , we get
Further, the step toward incorporation of the Caputo-Fabrizio operator is straightforward, namely
or in a compact form as
with the condition
Here denotes a Caputo-Fabrizio fractional operator with a fractional order βn based on the positive zeros of Jv.
Now, the naturally question coming to mind is: How this result can be applied to real data related to the stress relaxation and creep compliance of real viscoelastic materials? Unfortunately, no answers to this question exist in all works dealing with so-called Bessel media [114–119]. The first problem immediately appearing is: how the relaxation times can be related to the real data? As mentioned in several points of this article, the relaxation times corresponds to measurements taken at equidistantly distributed (in normal scale or logarithmic scale) points along the time axis. How, this real approach could be related to the zeros of the Bessel function is a question still remaining open. The second problem comes from the impossibility to work with infinite sum in (Equations 90 , 95 , and 97) . Actually, all computer simulations should use finite number of terms that immediately leads to truncated Dirichlet series obeying the condition (Equation 97). This immediately, transforms the approximation of the relaxation function to approximation through Prony series, but the first problem formulated above, still remains unanswered.
Actually, the approximation by Prony series or infinite Dirichlet series, generally speaking, is an approach considering approximations of Non-Debye responses (relaxations) by superpositions of sub-processes (as Debye relaxations) with different relaxation times (an approach widely applied in the relaxation processes in glass transitions , for instance).
7. Final Comments
The ideas and results developed in this work focused on the new fractional operator conceived in [1, 10] (2015) naturally appearing when we stayed on the shoulders of two classical results: the Prony series (1795)  and the Boltzmann superposition principle (1864) .
It was demonstrated that in many cases there are viscoelastic materials which experimental behaviors exhibit strong departures from the power-law. In such cases it is natural to rise the questions about the adequate modeling of the dynamic processes in such media and to ask for new fractional operators. This is, actually, the same question raised by Bagley and Torvik [49, 124] who in case of power-law media suggested the constitutive relationship
In (Equation 99) the derivatives have power-law kernels (Riemann-Liouville derivatives), which for N = 1 (one-term series of fractional derivatives) results in the simple expression
containing two fractional derivatives with different orders.
This article does not focus on development of different viscoelastic models based on the Caputo-Fabrizio fractional operator since this is out of its scope and draws new problems to be resolved. Despite this, if the kernels in the convolution integrals departure from the power-law we may consider a similar (formally) constitutive equation, namely
where, following the main idea of the Prony decompositions of the relaxation and the compliance, the series of fractional order of both sides of the equation have equal numbers of terms. Following [20, 107] the retardation times λi are satisfying the conditions
Therefore, from αi = 1/(1+τi/t0) and βi = 1/(1+λi/t0) it follows from Equation (102) that the fractional orders should satisfy the inequalities
The example of Renardy  (see also  and ), related to polymer rheology, reveals that a discrete relaxation spectrum with accumulation point at zero behaves as a power-law for short times, namely
Thus, in the asymptotic case (for t → 0 ) we may expect that the model (Equation 101) would converge to the model (Equation 100).
If we suggest only, for example, that N = 1, which actually is a departure of the main idea of Prony series, we have
This equation contains two fractional derivatives of different orders as in Equation (100). For instance, fractional viscoelastic models containing fractional derivatives and operators of two different order are thoroughly analyzed in Rossikhin and Shitikova  and Rossikhin and Shitikova  in the light of power-law materials. All these formal similarities and dissimilarities provoke new ideas that should be resolved.
To recapitulate, actually, we demonstrated that the Caputo-Fabrizio fractional operator naturally appears in the constitutive viscoelastic equations based on hereditary integrals when the material responses do not match the power-law. Moreover, it plays the same role as the fractional operators (derivatives) based on power-law memory kernels when they are inapplicable. At this moment the task of this study is accomplished and a lot of future works would be developed starting from results obtained here.
The author confirms being the sole contributor of this work and has approved it for publication.
Conflict of Interest Statement
The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
3. Hristov J. Derivatives with non-singular kernels: from the Caputo-Fabrizio definition and beyond: Appraising analysis with emphasis on diffusion models. In: Frontiers in Fractional Calculus, 2017 , Bhalekar S, editor. Sharjah: Bentham Science Publishers. (2017). pp. 269–342
4. Hristov J. Transient heat diffusion with a non-singular fading memory: from the cattaneo constitutive equation with Jeffrey's kernel to the caputo-fabrizio time-fractional derivative. Thermal Sci. (2016) 20:765–70. doi: 10.2298/TSCI160112019H
13. Hristov J. Derivation of fractional Dodson's equation and beyond: Transient mass diffusion with a non-singular memory and exponentially fading-out diffusivity. Progr Fract Differ Appl. (2017) 3: 255–70. doi: 10.18576/pfda/030402
24. Ciambella J, Paolone A, Vidoli S. A comparison of nonlinear integral-based viscoelastic models through compression tests on filled rubber. Mech Mater. (2010) 42:932–44. doi: 10.1016/j.mechmat.2010.07.007
25. Ciambella J, Paolone A, Vidoli S. Identification of the viscoelastic properties of soft materials at low frequencies: Performance, ill-conditioning and extrapolation capabilities of fractional and exponential models. J Mech Behav Biomed Mater. (2014) 37:286–98. doi: 10.1016/j.jmbbm.2014.05.031
37. Jozwiak B, Orczykowska M, Dziubinski M. Fractionalgeneralization of maxell and Kelvin-Voight models for biopolymer characterization. PLoS ONE (2015) 10:e0143090. doi: 10.1371./journal.pone.0143090
42. Schiessel H, Metzler R, Blumen A, Nonnenmacher TF. Generalized viscoelastic models: their fractional equations with solutions. J Phys A Math Gen. (1995) 28:6567–84. doi: 10.1088/0305-4470/28/23/012
43. Papagiannnopoulos A, Sotiropoulos K, Pispas S. Particle tracking microrheology of the power-law viscoelasticity of xantan solution. Food Hydrocolloids (2016) 61:201–10. doi: 10.1016/j.foodhyd.2016.05.020
56. Carillo S, Giorgi C. Non-classical memory kernels in linear viscoelasticity. In: Viscoelastic and Viscoplastic Materials, Mohamed Fathy El-Amin editor. Split: IntechOpen (2016). p. 295–331. doi: 10.5772/64251
57. Bradshaw RD, Brinson LC. A sing control method for fitting and interconverting material functions for linear viscoelastic solids. Mech Time Depend Mater. (1997) 1:85–108. doi: 10.1023/A:1009772018066
62. Phillips A, Pankaj P, May F, Taylor K, Howie C, Usmani A. Constitutive models for impacted morsellised cortico-cancelous bone. Biomateials (2006) 27:2162–70. doi: 10.1016/j.biomaterials.2005.10.034
63. Ravikumar N, Noble C, Cramphorn E, Taylor ZA. A short constitutive model for ballistic gelatine at surgical strain rates. J Mech Behav Biomed Mater. (2015) 47:87–94. doi: 10.1016/j.jmbbm.2015.03.011
64. Londono JG, Berger-Vergioat L, Waisman H. A Prony series type viscoelastic solid coupled with continuous damage law for polar ice modeling. Mech Mater. (2016) 98:81–97. doi: 10.1016/j.mechmat.2016.04.002
65. Babaei B, Davarian A, Pryse KM, Elson EL, Genin GM. Efficient and optimized identification of generalized viscoelastic relaxation spectra. J Mech Behav Biomed Mater. (2016) 55:32–41. doi: 10.1016/j.jmbbm.2015.10.008
68. Jaloha D, Constantinescu A, Neviere R. Revisiting the identification of the generalized Maxwell models from experimental results. Int J Solids Struct. (2015) 67–68: 169–81. doi: 10.1016/j.ijsolstr.2015.04.018
69. Sun Y, Chen J, Huang B. Characterization of asphalt concrete linear viscoelastic behaviour utilizing the Havriliak-Negami complex modulus model. Constr Build Mater. (2015) 99:226–34. doi: 10.1016/j.conbuildmat.2015.09.016
75. Prieto-Munoz P, Yin H, Testa R. Mechaniscs of an adhesive anchor system subjected to a pullout load, II: viscoelastic analysis. J Struc Eng. (2013) 140:04013053. doi: 10.1061/(ASCE)ST.1943-541X.0000822
81. Tayeb A, Arfaoui M, Zine A, Hamdi A, Benabdallah J, Ichchou M. On the non-linear viscoelastic behaviour of rubber-like materials: constitutive description and identification Int J Mech Sci. (2017) 130:437–47. doi: 10.1016/j.ijmecsci.2017.06.032
82. Vandenberghe E, Choucharina S, Luca S, De Ketelaere B, De Baerdemaeker J, Claes J. Spatio-temporal gradients of dry matter content and fundamental material parameters of Gouda cheese. J Food Eng. (2014) 142:31–48. doi: 10.1016/j.jfoodeng.2014.05.019
84. Mitra S, Mitra A. A genetic algorithms based techniques for computing the nonlinear least squares estimates of parameters of sum of exponential model. Expert Syst Appl. (2012) 39: 6370–9. doi: 10.1016/j.eswa.2011.12.033
85. Huang AW, Lu CH, Wu SC, Vinci RP, Brown WL, Lin MT. Viscoelastic mechanical properties measurement of thin AL and Al-Mg films using bulge testing. Thin Solid Films (2016) 618:2–7. doi: 10.1016/j.tsf.2016.03.064
87. DeHoff PH, Barrett AA, Lee RB, Anusavice KJ. Thermal compatibility of dental ceramic system using cylindrical an spherical geometries. Dental Mater. (2008) 24:744–52. doi: 10.1016/j.dental.2007.08.008
89. Andrews SHJ, Rattner JB, Shrive NG, Ronsky JL. Swelling significantly affects the material properties of the menisci in compression. J Biomech. (2015) 48:1485–9. doi: 10.1016/j.jbiomech.2015.02.001
93. Sliker LJ, Ciuti G, Rentschler M, Menxciasi A. Frictional resistance model for tissue-capsule endoscope sliding contact in the gastrointenstinal tract. Tribol Int. (2016) 2012:472–84. doi: 10.1016/j.triboint.2016.06.003
94. Cui T, Lin CW, Chien CH, Chao YJ, Van Zee JW. Service life estimation of liquid silicone rubber seals in polymer electrolyte membrane fuel cell environments. J Power Sources (2011) 196:1216–21. doi: 10.1016/j.jpowsour.2010.08.075
97. Pelayo F, Lamela-Rey M, Minuz-Calvente M, Lopez-Aenlle M, Alvarez-Vasquez A. Study of the time-temperature behaviour of PVB: Application of laminated glass elements. Thin Walled Struct.(2017) 119:324–31. doi: 10.1016/j.tws.2017.06.030
98. Schneider J, Hilcken J, Aronen A, Karvinen R, Olesen JF, Nielsen J. Stress relaxation in tempered glass caused by hear soak testing. Eng Struct. (2016) 122:42–9. doi: 10.1016/j.engstruct.2016.04.024
107. Luk-Cyr J, Crochon T, li C, Levesque M. Interconversion of linearly viscoelastic material functions expressed as Prony series: a closure. Mech Time Depend Mater (2013) 17:53–82. doi: 10.1007/s11043-012-9176-y
110. Hristov J. Steady-state heat conduction in a medium with spatial non-singular fading memory: derivation of caputo-fabrizio space-fractional derivative with Jeffrey's kernel and analytical solutions. Thermal Sci. (2017) 21:827–39. doi: 10.2298/TSCI160229115H
113. Canestrati F, Stimilli A, Bahia HA, Virgili A. Pseudo-variables method to calculate HMA relaxation modulus through low-temperature induced stress and strain. Mater Design (2015) 76:141–9. doi: 10.1016/j.matdes.2015.03.063
126. Rossikhin YuA, Shitikova M. Analysis of viscoelastic rod dynamics via models involving fractional derivatives or operators of different orders. Shock Vibration Digest (2004) 36:3–26. doi: 10.1177/0583102404039131
Keywords: linear viscoelasticity, response functions, interconversion, power-law response, non power-law responses, prony series, caputo-fabrizio operator
Citation: Hristov JY (2018) Linear Viscoelastic Responses: The Prony Decomposition Naturally Leads Into the Caputo-Fabrizio Fractional Operator. Front. Phys. 6:135. doi: 10.3389/fphy.2018.00135
Received: 27 September 2018; Accepted: 07 November 2018;
Published: 05 December 2018.
Edited by:Dumitru Baleanu, University of Craiova, Romania
Reviewed by:Yilun Shang, Northumbria University, United Kingdom
Aliyu Isa Aliyu, Federal University Dutse, Nigeria
Copyright © 2018 Hristov. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jordan Yankov Hristov, firstname.lastname@example.org