Two Integral Representations for the Relaxation Modulus of the Generalized Fractional Zener Model

: A class of generalized fractional Zener-type viscoelastic models with general fractional derivatives is considered. Two integral representations are derived for the corresponding relaxation modulus. The ﬁrst representation is established by applying the Laplace transform to the constitutive equation and using the Bernstein functions technique to justify the change of integration contour in the complex Laplace inversion formula. The second integral representation for the relaxation modulus is obtained by applying the subordination principle for the relaxation equation with generalized fractional derivatives. Two particular examples of the considered class of models are discussed in more detail: a model with fractional derivatives of uniformly distributed order and a model with general fractional derivatives, the kernel of which is a multinomial Mittag-Lefﬂer-type function. To illustrate the analytical results, some numerical examples are presented.


Introduction
Analytical methods play an essential role in the study of linear viscoelastic models [1,2]: this is, in general, due to the fact that in such models the rheological properties of a viscoelastic medium are described by a linear constitutive relation between stress σ and strain ε.Such a stress-strain relation can be studied analytically by employing the Laplace transform.
Let us consider the uniaxial case, in which σ = σ(x, t) and ε = ε(x, t), where x ∈ R and t ≥ 0 are spatial and temporal variables, respectively.For systems that are at rest before some starting time, t = 0, the relaxation modulus G(t) is defined via the convolutional relation [2] σ(x, t) = t 0 G(t − τ) ∂ε ∂τ (x, τ) dτ, (1) which allows us to study the relaxation modulus by the use of the Laplace transform with respect to time.
Fractional calculus has been extensively employed in linear viscoelasticity, due to its ability to model phenomena with memory [1][2][3].For instance, incorporating non-integer time derivatives in viscoelastic equations provides an appropriate framework for the description of properties of polymer solution and melts [4].One of the most used fractional viscoelastic models of solid-like behavior is the fractional Zener model, which is defined by the following stress-strain constitutive relation [1,2]: (1 + aD α t )σ(x, t) = (1 + bD α t )ε(x, t), 0 < a < b, 0 < α < 1. (2) Here, D α t is the fractional Riemann-Liouville derivative of order α ∈ (0, 1), defined by (see, e.g., [5]) In the limiting case α = 1, the stress-strain relation (2) represents the classical Zener model, which is also referred to as the Standard Linear Solid model [2].
The relaxation modulus G(t), in the case of the fractional Zener model (2), admits the explicit representation where E α (•) denotes the Mittag-Leffler function Different integral representations for the relaxation modulus (4) can be derived based on such representations for the Mittag-Leffler function.First, the integral representation formula for the Mittag-Leffler function in [6], Equation (A30), yields where the spectral function K(r) is defined as follows: Since K(r) ≥ 0 under the restrictions on the parameters in (2), representation (6) implies that the function G(t) is completely monotone.
Another type of integral representation follows from the formula (see, e.g., [2]) where M α (•) is the following function of Wright type: This function is usually referred to as the Mainardi function.Inserting ( 7) into (4) yields the integral representation where Various types of generalizations of the classical fractional Zener constitutive model (2) have been proposed and studied in the literature.In [7][8][9][10][11][12], different Zener-type models of distributed order are studied.In [7], the following stress-strain relation of distributed order was proposed: where the weight functions p σ (α) and p ε (α) are non-negative.By setting where 0 ≤ α 0 < α 1 < ... < α N < 1 and a n , b n > 0 for all n = 0, 1, ..., N, and δ(.) denotes the Dirac delta function, the following constitutive equation is obtained: The multi-term stress-strain relation (11) was studied in [8], and restrictions on the parameters a n , b n were found, which guarantee thermodynamic compatibility.When the weight functions in (10) have the form then the model is the so-called power-type distributed-order model.It was proved in [7] that the rheological model (10) with power weight functions ( 12) is thermodynamically compatible, provided a < b.The constitutive Equation ( 10) and the related mechanical models are studied in [9,10].Integral representations of the relaxation moduli, which generalize representation (6), were deduced in [11] for the model (11) and the model with continuous power-type distribution ( 10)-( 12) under the above-mentioned thermodynamic restrictions on the parameters.The obtained representations imply that the corresponding relaxation moduli are completely monotone functions.In [12], restrictions on the parameters guaranteeing thermodynamic compatibility were formulated for the distributed-order model (10), with weight functions of a more general form: where A(α) and B(α) are non-negative integrable functions defined for α ∈ [0, 1].In recent works [13,14], constitutive models of the Zener type, containing both fractional integrals and derivatives, have been formulated and their thermodynamical consistency has been studied.Relaxation moduli were calculated, and restrictions on model parameters narrowing thermodynamical constraints were posed, in order to ensure complete monotonicity.
A Zener-type model with general fractional derivatives was proposed in [15], and restrictions on the coefficients were derived that followed from the dissipation inequality.General fractional differential operators are a subject of increasing interest-see, e.g., [16][17][18], to mention only a few of many recent publications.
Motivated by [15], where a generalization of model ( 2) was proposed, with fractional derivative D α t replaced by a general convolutional derivative D (k) t , we consider the generalized fractional Zener model with the constitutive equation Here, D t is the convolutional derivative in the Riemann-Liouville sense, defined by where k(t) is a locally integrable memory kernel.Further assumptions on the kernel k(t) are specified later in Section 3.
Based on the properties of completely monotone, Stieltjes, and complete Bernstein functions, in [19] the complete monotonicity of the relaxation modulus for the generalized Zener model (13) was proved and the subordination principle for the corresponding Zenertype wave equation was established.
Integral representations for the relaxation moduli of viscoelastic models have proven useful for the analytical and numerical study of the corresponding models.Therefore, the aim of the present work is to establish integral representations for the relaxation modulus G(t) of the Zener model (13), which generalize representations (6) and (9).
First, we establish a representation of the form (6), based on deforming the contour of integration in the complex inversion formula for the Laplace transform.Such a technique has been extensively used in the literature.In [20], sufficient conditions were formulated under which the change of integration contour to the so-called Hankel path was justified.In [1,10], the method was applied to different fractional evolution equations.Applications to distributed-order relaxation equations can be found in [21,22].
To obtain an alternative integral representation for G(t), which generalises (9), we use the relation of this function and the solution of a relaxation equation with a general fractional Caputo derivative with kernel k(t), and we apply the subordination principle for this equation [16,23,24].The concept of subordination, originally introduced by Bochner in the theory of stochastic processes, has developed recently into a powerful tool in the study of anomalous relaxation and diffusion processes and the physics of complex systems (see, e.g., [25][26][27] and the recent review paper [28]).
The present work is organized as follows.Section 2 contains preliminaries on completely monotone, Stieltjes and complete Bernstein functions.In Section 3, the generalized fractional Zener model is formulated and some basic properties are discussed.In Section 4, the first integral representation for the relaxation modulus is derived.Section 5 contains two particular examples: the Zener model with fractional derivatives of uniformly distributed order and the one with general fractional derivatives, the kernel of which is a multinomial Mittag-Leffler-type function.Numerical results, based on the first integral representation, are given.Applying the subordination principle for the generalized fractional relaxation equation, a second integral representation for the relaxation modulus is obtained in Section 6. Concluding remarks are given in Section 7.

Preliminaries
The sets of real and complex numbers are denoted by R and C, respectively, and The Laplace transform of a function f (t) is denoted by f or L{ f }: The class of completely monotone functions, which we denote by CMF , consists of all real-valued functions f (t), defined on the half-line R + , such that According to Bernstein's theorem, a function is completely monotone if and only if it can be represented as the Laplace transform of a measure on R + .
The class SF of Stieltjes functions consists of all non-negative functions ϕ defined on R + that admit the representation (see [16]) where A ≥ 0, B ≥ 0, f (t) ∈ CMF and the Laplace transform of f exists for any s > 0.
The constants A and B in ( 16) are determined by the identities Any Stieltjes function is completely monotone.The class CBF of complete Bernstein functions consists of all non-negative functions ϕ defined on R + , which admit the representation [29] where A ≥ 0, B ≥ 0, and m(t) is a completely monotone function, such that The constants A and B in ( 18) are defined as follows: Next, some properties of the above special classes of functions are summarized for further use in this work.

Theorem 1 ([29]
). Suppose that ϕ is a non-negative function on R + and s > 0. The following statements hold true: (P1) The sets CMF , SF , and CBF are convex cones, closed under pointwise limits.
which is defined by the expression where A, B ≥ 0 are the constants defined in (19) and ν is a Borel measure on R + satisfying ) Any function ϕ from the classes CBF or SF admits an analytic extension to the complex plane cut along the negative real axis C\(−∞, 0], such that where * stands for complex conjugate.Moreover, if z ∈ C\(−∞, 0], then and (P6) Let ϕ be the analytic extension of a complete Bernstein function to C\(−∞, 0].Then, the limit exists and is real.
For details on the proofs, we refer to [29].
The following result is usually referred to as the Karamata-Feller theorem.
Here, and later in this work, the notation f (t) ∼ g(t) as t → t * stands for lim

Model Formulation
We consider the generalized fractional Zener model defined by constitutive Equation ( 13) with the convolutional derivative in the Riemann-Liouville sense (14).We assume that the Laplace transform of the kernel k(t) exists for all s > 0 and where SF denotes the class of Stieltjes functions (see (16)).Moreover, we suppose lim s→+∞ s k(s) = +∞ (27) in the Laplace domain (1), yielding Here, σ(x, s) and ε(x, s) denote the Laplace transforms of the functions σ(x, t) and ε(x, t) with respect to t, while x is considered as a parameter.Applying the Laplace transform to constitutive Equation ( 13), and using the identity we deduce the relation In this way, we obtain Let us mention some basic characteristics of the model.Based on (29) and assumption (27) in [19] we found that Furthermore, assumption ( 26) is equivalent to (see (P2)) which, according to property (P6), implies where s → 0+ denotes (0, ∞) s → 0. Therefore, applying the final value theorem for the Laplace transform, we obtain Since G(+∞) > 0, i.e., there is no full relaxation, constitutive Equation ( 13) is a model for solid-like viscoelastic behavior as the classical fractional one (2).
It was proven in [15] that condition 0 < a < b is sufficient for the thermodynamic compatibility of model (13).In [19], it was established that this condition implies a stronger property of complete monotonicity of the relaxation modulus.Moreover, it appears that the constraint a < b is also a necessary condition for the physical acceptability of model (13).Indeed, in a physically meaningful model, the function G(t) should be monotonically decreasing, in particular G(0+) > G(+∞): this, by taking into account (30) and (33), implies a < b.
In the rest of this work, we assume k 0 = 0, where k 0 is defined in (32).This is a normal assumption for the kernels of general fractional derivatives (see [16]).In the present work, this restriction is posed only for technical convenience.

First Integral Representation
To derive an integral representation for G(t), we use the derived Laplace transform G(s), given in (29), which we rewrite in the form where with g(s) = s k(s), as defined in (31).Since g(s) ∈ CBF , g(s) admits an analytic extension to C\(−∞, 0]-see (P5)-and g(s) is a non-negative and monotonically non-decreasing function for s > 0. Therefore, g(s) > 0 for s > 0, and (35) can be written in the form Applying properties (P1) and (P3), it follows that 1 + (ag(s)) −1 ∈ SF .Then, by (P2), Applying again (P3), we obtain F(s) ∈ SF .This, by (P5), in particular, implies that the function F(s) admits an analytic extension to the complex plane cut along the negative real axis C\(−∞, 0].Therefore, to derive a representation for G(t), we can apply the inverse Laplace integral formula to (34) and change the integration contour in the complex plane.The result is given in the next theorem.As usual, the notation φ → π− stands for φ → π, φ < π.Theorem 3. Assume 0 < a < b and let the kernel k ∈ L 1 loc (R + ) be such that k(s) ∈ SF and Suppose, moreover, that there exists ε > 0, such that where p(r) is a function, which does not depend on φ and p(r)e −ωr ∈ L 1 (R + ) for any ω > 0.
Assume the limit lim φ→π− k(re iφ ) exists for almost all r > 0.Then, the relaxation modulus G(t) for constitutive model ( 13) is a continuous function for t ≥ 0, infinitely differentiable and completely monotone for t > 0, and admits the integral representation Here, the function K(r) ≥ 0 is defined by the identity Proof.Applying the inversion formula for the Laplace transform to (34), and taking into account that L{1} = 1/s, we obtain where function F(s) is defined in (35).The function F(s) has the Stieltjes function representation ( 16) with A = B = 0, calculated by applying (17) and (37).Therefore, F(s) is the Laplace transform of a completely monotone function, denoted by f (t)-that is, F(s) = f (s).Then, G(t) = 1 + f (t) ∈ CMF and, in particular, it is infinitely differentiable for t > 0. It follows from (37) that Applying Theorem 2, (35) implies In this way, we prove that G(t) is continuous for t ≥ 0 and lim t→0+ G(t) = b/a.Next, our aim is to find a representation for f (t).

Examples
In this section, we consider two examples of kernels k in the definition (14) of the convolutional derivative D (k) t , and we show that for both examples the conditions of Theorem 3 are satisfied.Asymptotic expansions of the corresponding relaxation moduli are also discussed.

Fractional Derivative of Uniformly Distributed Order
Consider first the following particular case of the Zener-type model (13).
Example 1.Let ( 14) be the distributed-order fractional derivative with uniform distribution in the interval (0, 1).In this case, In the case of kernel (48), the functions k 1 (r) and k 2 (r) in (40) are given by From ( 35) and (48), we obtain which implies Therefore, assumption (38) is satisfied.
Let us establish the asymptotic expansion of the relaxation modulus G(t) for t → +∞ in the case of Example 1.We use the expansion which follows from (49).Applying the Karamata-Feller theorem (Theorem 2) with α = 1 and L(x) = (a + ln x) −1 , and taking into account relation (34), we derive Therefore, the function G(t) exhibits a slow logarithmic decay for large times and Theorem 3 is applied for numerical computation of G(t) for the model in Example 1.The results obtained for several values of the parameters a and b are plotted in Figure 1.
In this subsection, we consider the following special case of the generalized Zener model.

Example 2.
Assume the kernel in definition ( 14) is given by the completely monotone function of multinomial Mittag-Leffler type: where In this case, k(s) is given by (55), limits (37) hold, and k(s) ∈ SF as the Laplace transform of a completely monotone function.Inserting (55) into (29) yields A comparison to the fractional Zener model with multiple fractional derivatives (11) studied in [8] shows that the model of Example 2 is a special case of it and corresponds to the constitutive equation Furthermore, inserting g(s) = s k(s) with k(s) from (55) into (36) yields where c = b/a − 1.By the use of identity (55), we deduce that (58) is a Laplace transform of a multinomial Mittag-Leffler-type function, which leads to the following explicit representation Therefore, the relaxation modulus has the explicit representation G(t) = 1 + f (t) with f (t) given in (59), which is a generalization of Formula (4) for the classical fractional Zener model.
Let us check the rest of the conditions of Theorem 3. Since (58) implies This can be deduced from the property g(s) 1 1−β ∈ CBF (c.f., Proposition 3.1 in [32]) for the function This property implies by the use of ( 22) the inequality which leads to estimate (43) (in the same way as in the proof of Theorem 3), which, in turn, implies (42).This estimate guarantees the integrability of |F(s)|e −ω|s| for all s except those near the origin.Furthermore, in the case of kernel (57), the functions k 1 (r) and k 2 (r) in (40) are given by , where Let us establish the asymptotic behavior of the function G(t) as t → 0+ and t → +∞.For t → 0+ the series expansion (53) implies For t → ∞ we use the asymptotic expansion of the multinomial Mittag-Leffler function, given in [32], for the function f (t) in the representation (59), which yields Therefore, the relaxation modulus in the case of Example 2 exhibits an algebraic decay for large times to the final value G(+∞) = 1.
Although in the case of Example 2 there exists an explicit representation for the function G(t), in terms of the multinomial Mittag-Leffler-type function (59), the integral representation of Theorem 3 is still useful.For instance, it is convenient for numerical computation of the relaxation modulus G(t) and, respectively, the multinomial Mittag-Leffler-type function (59) (the infinite series (53) is not suitable for this purpose).
Theorem 3 is applied for the numerical computation of G(t) for the model with kernel (57) with m = 2, λ 1 = λ 2 = 1, α 1 = 0.5, α 2 = 0.3, β = 0.8. Figure 2 presents the numerical results obtained for several values of the parameters a and b.In order to enable comparison, we have chosen the same values of a and b as in Figure 1.Compared to the slow (logarithmic) decay in Figure 1, in Figure 2 we observe a faster (algebraic) decay for large times: this is in agreement with the asymptotic expansions (60) and (52).For all plots, the numerical initial value is in agreement with formula (30).

Second Integral Representation
Another integral representation for the function G(t) can be derived by applying the subordination principle for the generalized fractional relaxation equation [16,23,24]: where λ > 0 and C D t is the convolutional derivative in Caputo sense [16]: Problem (61) is studied in [16], see also [23,24].Detailed study in the case of distributed order derivatives can be found in [21,22].
The subordination relation for the solution to problem (61) is given next, adapted for use here (see [23,24]).
Theorem 4 ( [23,24]).Assume k ∈ L 1 loc (R + ) is such that k(s) ∈ SF and satisfies (37).Then, the solution u(t) to problem (61) admits the integral representation where the function ϕ(t, τ) is a unilateral probability density, i.e., which is defined by the identity From the relation between the convolutional derivatives in Riemann-Liouville and Caputo sense in (62) and identity (28), we obtain for the solution of problem (61) the following expression in the Laplace domain: Applying the uniqueness property of the Laplace transform, the obtained identity compared to (35) implies that the function f (t) in the representation G(t) = 1 + f (t) is the solution of the relaxation Equation (61), with convolutional derivative in Caputo sense and with the same kernel k as in (14), and Applying the subordination relation (63), we deduce the following: Theorem 5. Assume 0 < a < b, the kernel k ∈ L 1 loc (R + ) is such that k(s) ∈ SF and satisfies conditions (37).Then, the relaxation modulus G(t) to model ( 13) admits the integral representation where the function ϕ(t, τ) is defined in (65) and is a unilateral probability density function, i.e., it satisfies (64).
We close this section with an estimate, which follows from the established relation of G(t) and the solution of the generalized fractional relaxation Equation (61), by applying Theorem 5.8 in the dissertation [24].
Corollary 1. Assume the conditions of Theorem 5 are satisfied.Then, the relaxation modulus G(t) to model (13) obeys the estimate where the function l(t) is related to the kernel k(t) via the identity t 0 l(t − τ)k(τ) dτ = t.
For instance, in the case of Example 2, by applying the Laplace transform we obtain where the functions ω α (t) are defined in (3).

Concluding Remarks
In the present work, two integral representations are derived for the relaxation modulus of the generalized fractional Zener model (13).The results are formulated in Theorems 3 and 5.The main difficulty in the proof of Theorem 3 is to justify the change of the integration contour for the considered class of kernels.This is done based on the relevant properties of Stieltjes functions and related classes of functions.On the other hand, Theorem 5 is a straightforward implication of the subordination principle for the generalized fractional relaxation equation.
The first representation, (39), demonstrates that G(t) is a completely monotone function under the assumptions of Theorem 3: this is implied by the non-negativity of the spectral function K(r).For some models, deriving such a representation with a nonnegative spectral function is the most convenient way to establish complete monotonicity of the relaxation modulus (see e.g., [11]).In addition, representation (39) is appropriate for numerical computation of the relaxation modulus.
The second integral representation, (66), splits the relaxation modulus into two parts: a probability density function ϕ(t, τ), which depends only on the kernel k and is independent of the model parameters a and b, and the function (b/a − 1)e −τ/a .Such a decomposition makes representation (66) suitable for obtaining estimates for the relaxation modulus G(t), as well as for the study of the effect of the parameters a and b on its behavior, which will be a subject of future work.
The technique presented in this paper could be applied to other viscoelastic constitutive equations, to obtain integral representations for the corresponding relaxation moduli.Such representations would be useful for the analytical and numerical study of the models.

Figure 1 .
Figure 1.Relaxation modulus G(t) for the model of Example 1: (a) a = 0.3 and different values of b; (b) b = 3 and different values of a.
38) is satisfied.Let us note that in the case of kernel (57) estimate (42) is satisfied for all s ∈ C\(−∞, 0].