Oscillations and damping in the fractional Maxwell materials

This paper examines the oscillatory behaviour of complex viscoelastic systems with power law-like relaxation behaviour. Specifically, we use the fractional Maxwell model, consisting of a spring and fractional dashpot in series, which produces a power-law creep behaviour and a relaxation law following the Mittag-Leffler function. The fractional dashpot is characterised by a parameter beta, continuously moving from the pure viscous behaviour when beta=1 to the purely elastic response when beta=0. In this work, we study the general response function and focus on the oscillatory behaviour of a fractional Maxwell system in four regimes: stress impulse, strain impulse, step stress, and driven oscillations. The solutions are presented in a format analogous to the classical oscillator, showing how the fractional nature of relaxation changes the long-time equilibrium behaviour and the short-time transient solutions. We specifically test the critical damping conditions in the fractional regime, since these have a particular relevance in biomechanics.


I. INTRODUCTION
The damping of free oscillations is a subject with a long and rich history. However, there is less work available looking at how the damping of oscillations depends specifically on the viscoelasticity of a system, and in particular, the types of viscoelasticity seen in biological tissues and complex composite materials. The relaxation function in many such systems is distinctly nonexponential, which often referred to as a distribution or spectrum of internal relaxation times. An alternative way to describe such complex viscoelasticity is via a viscoelastic element characterized by a fractional time derivative. In this work, we examine the oscillatory behavior of a viscoelastic element characterized by a fractional Maxwell model [1][2][3]. The fractional Maxwell model is a generalization of the classical model (based on the elastic and damping elements in series) and successfully describes the relaxation function of polymeric materials [4][5][6][7], biological tissues [8][9][10], cells [11,12], and foods [13][14][15] over a large range of time scales. It is readily generalized to a fractional Zener model [16,17] for the description of viscoelastic solids (since the Maxwell model ultimately produces a plastic flow in response to a constant force). It is useful to be able to relate the well-understood characteristics of viscoelastic behavior to how a system dissipates energy in free oscillations. This in particular can give insight into the role biological tissues have in dissipating energy from their relaxation function.
This work draws from the rich field of fractional oscillators, with earlier rigorous work of Rossikhin and Shitikova looking at its applications in viscoelasticity [18,19]. Fractional dynamics involves taking the integer order differential operators and replacing them with operators of a fractional order [20,21]; however, in many papers which use fractional dynamics there is no direct physical justification for this. This is particularly true for harmonic oscillators where the inertial term is taken as fractional (with arguably unclear physical foundations). Here, we retain a standard inertia effects, but include the fractional-Maxwell relaxation function that models realistic viscoelasticity, and use it as the origin of elasticity and damping to derive oscillatory response as a result to several representative initial conditions.
In this paper, we first review the linear viscoelasticity and the corresponding relaxation function to establish a reference point, and then justify the use of the fractional Maxwell model and where it might apply. We then derive the general response function of a viscoelastic element characterized by a given relaxation function, G(t). Sections III A-C discuss the response of the fractional Maxwell model to a stress impulse, strain impulse, step stress, and driven oscillations, while drawing an analogy to the classical case. Sections IV and V discuss the frequency and decay of free oscillations for the fractional Maxwell model. The interest here is to give insight into the damping of oscillations in biological materials, which typically do have a more complicated viscoelasticity of a fractional type. One example would be a tendon, an essential mechanical element of large organisms, which should have low damping and respond in a similar way to all frequencies of oscillation to fulfill its biological function. On the other hand, for other soft connective tissues it would be preferable to be positioned close to critical damping conditions to prevent injury by efficiently dissipating any excess mechanical energy from impacts that could not be fully absorbed by muscle (such as in running or landing from a jump).

II. VISCOELASTIC MODELS
In this section, we briefly review basic viscoelasticity and examine relaxation functions. We also demonstrate the a) crossover between some of the phenomenological models and how often it can be difficult to distinguish among them. In linear viscoelasticity, the stress r arising in response to some arbitrary strain, eðtÞ, is determined by the Boltzmann superposition integral [22] rðtÞ ¼ ð t

À1
Gðt À t 0 Þ_ eðt 0 Þdt 0 ; (1) where G(t) is the characteristic linear response function of the material, often called the relaxation function. We can see that G(t) is measured by the system response to a step strain, where _ eðtÞ ¼ De Á dðt 0 Þ and Eq. (1) becomes rðtÞ ¼ GðtÞ Á De. The Boltzmann superposition integral itself is the limit sum over infinitesimal step strains determined by de ¼ ðde=dsÞds. Figure 1(a) shows a schematic of typical relaxation curve, which includes the key features of: (a) The instantaneous response to step strain is often called the glass modulus; (b) the modulus relaxes over some characteristic time, s, to an equilibrium modulus G e ; (c) the amount of modulus relaxation is labeled as G r .
The relaxation time s is only clearly defined in the case of simple exponential relaxation GðtÞ ¼ G e þ G r e Àt=s . However, it is very rare that a practical material can be described with a single exponential relaxation, except in the cases of transient networks with a single rate of crosslink breaking [23]. The attempts to rationalize a more complex relaxation response often apply a "brute force" fit with multiple relaxation times, the so called generalized Maxwell model. Another common relaxation function that fits well with materials with a fractal or hierarchical structure is the power law. This is often called scale-free rheology as there is no end to relaxation that follows At Àb at all probed time scales, with only the magnitude A changing. An issue with power law relaxation is that there is no finite modulus at t ¼ 0, an issue usually avoided by using the asymptotic power law 1=ð1 þ At b Þ.
Viscoelasticity is often described through phenomenological models built of springs and dashpot elements of the basic Maxwell model; the simplest alternatives being the Kelvin-Voigt model where the spring and dashpot are connected in parallel and the Zener model where there is one additional elastic element in series with the Voigt frame. However, to fit realistic relaxation functions a great many dashpots and springs are often required, and with the number of parameters increasing the basic physical sense rapidly dissipates. The ideas of fractional viscoelasticity have been developed to address the power-law relaxation observations while retaining a fixed small number of fitting parameters. In the words of [2], it achieves a "mimicry of memory." In this case, the dashpots are replaced with fractional dashpots (sometimes called Scott-Blair elements), which on their own follow power-law relaxation [24]. The use of fractional dashpot changes the constitutive equation from r ¼ gde=dt to where the fractional derivative is a formal generalization of a derivative to noninteger order (which in practice is only meaningful in the reciprocal representation of Fourier or Laplace transformations). Sometimes, the fractional Maxwell model is presented as two fractional dashpots in series, which gives an extra fitting parameter; however, for simplicity we focus on the case with just one fractional dashpot. For a good introduction to fractional derivatives in viscoelasticity, see [3,25], while the general matters of fractional calculus can be reviewed in [26]. Table I shows several model viscoelastic relaxation functions and their Laplace transforms. The fractional Maxwell model is useful to model a variety of relaxation processes by adding only one additional parameter, b. This is opposed to more complicated models, such as multirelaxation models which may need a great many parameters before a good fit is found. It also differs from other relaxation functions, such as the stretched-exponential and nonasymptotic power law, which do not have phenomenological models to underpin the empirical expressions. The step strain response of an isolated fractional dashpot element is simply a power law of the form At Àb =C½1 À b, where C is gamma function. When b ¼ 0 the fractional dashpot acts as a perfectly elastic solid; as b increases, the element becomes more viscous, until the point of b ¼ 1 is reached, when ðAt À1 =C½0Þ ¼ AdðtÞ becomes the relaxation of a pure Newtonian fluid. The fractional Maxwell model has its relaxation described by the Mittag-Leffler function, E b ½Àðt=sÞ b (see Appendix A for details), which behaves very similarly to the nonasymptotic power law, 1=ð1 þ ðt=sÞ b Þ, being equal to one at t ¼ 0 and converging on the power law at long times. For small b, the Mittag-Leffler function is fitted exceptionally well by the asymptotic power law, which is also roughly equivalent to the stretched exponential with 1=2 exp½1 À ðt=sÞ b=2 . In many situations, these functions are equivalent, or practically indistinguishable over typical experimental time scales [27], as shown in Fig. 2.

Model
Response, G(t) (Pa)GðsÞ (Pa s) It is clear that in a great many situations the Mittag-Leffler function is as good a fit, if not better, than the powerlaw, asymptotic power law, or a stretched exponential. It has the key benefits of a finite value at zero time, having few parameters, and having a compact Laplace transform which is useful for any analytical calculations. The Mittag-Leffler function is frequently used to fit the relaxation in polymeric materials [28]. In particular, in biological tissues [29][30][31] and cells [32][33][34], the stress relaxation often exhibits low power laws, which lend themselves being fit to the Mittag-Leffler function. For example, the Mittag-Leffler function has been used to fit human breast tissue cells relaxation function with b % 0:5 [35].
The following sections examine the oscillatory response of the fractional Maxwell model in response to an impulse, impulses strain, step strain, and driven oscillations. In each case, we first solve the classical Maxwell model to highlight the differences and similarities.

III. GENERAL SOLUTION
Let us examine the simple model of a mass attached to a generic viscoelastic element, with the other end fixed, and derive the general solution for its motion, see Fig. 1(b). The position of the mass is governed by its inertia, an external force it is subjected to, F(t), and the relaxation function of the viscoelastic material, G(t). Thus, the equation of motion takes the form where x ¼ x 0 eðtÞ is the position of the mass, A is the crosssection area, and e the strain in the viscoelastic element. This is a balance between inertial forces, the stress in the viscoelastic element given through the Boltzmann superposition integral, and the external force. As dx=dt ¼ x 0 de=dt, Eq. (2) can be written in terms of stress and strain where l ¼ mx 0 =A is the reduced mass and rðtÞ ¼ FðtÞ=A is the external stress, respectively. The Laplace transformation of Eq. (3) gives a general relation We assume that until t ¼ 0 the mass was stationary, such that eðtÞ ¼ 0 and _ eðtÞ ¼ 0 (no infinite acceleration is allowed). The nontrivial initial condition that initiates motion is that at t ¼ 0 a stress rðtÞ is applied. Under these conditions, we can write the general solution in Laplace spacẽ e s ð Þ ¼r s ð Þ ls 2 þ sG s ð Þ : To find the strain response as a function of time, we must first substitute into Eq. (5) the Laplace transform of the external stress, which reflects the specific mode of dynamic experiment, and also the specific viscoelastic relaxation function, which reflects the physical nature of the material. After this, the inverse Laplace transformation has to be performed. Table II shows the four characteristic experiments, which produce initial conditions on the externally applied stress that we look at in this paper, along with their Laplace transforms. The characteristic of an impact is a "stress impulse" (this solution is also called the transfer function). In this case, the external condition is a step-function change in momentum (or velocity) centered at t ¼ 0, which manifests itself as a delta-function stress. A step stress is equivalent to having a weight suddenly attached, allowing the system to TABLE II. Four external stress conditions and their corresponding Laplace transforms. The first three cases arise from a discontinuity that occurs between the ambient state before t ¼ 0 and the solutions for t > 0.

Initial condition
rðtÞ ðPaÞrðsÞ ðPa sÞ Impulse stress R 0 dðtÞ R 0 Step stress r 0 r 0 =s Impulse strain lDeðddðtÞ=dtÞ lDe Driven oscillation r 0 sinðxtÞ r 0 x=ðx 2 þ s 2 Þ find a new equilibrium position. The strain impulse is the act of "plucking" the mass, that is, moving it instantaneously from zero strain to some finite strain and then releasing: The stress function used in this initial condition may seem strange, but the same general solution in Laplace space can be found by having no initial stress and setting eð0Þ to De. In this case, an instantaneous change in strain is reflected by a step-function change in position at t ¼ 0, which corresponds to the derivative of the delta function in force. The final initial condition is the response to driven oscillations, which is another focus of this paper.
In each case, before solving for the fractional Maxwell model for each case of the initial condition, we first solve for the classical Voigt and Maxwell models. The Voigt model has the same solutions as the linear damped oscillator, which we use to remind the reader of the familiar classical results. The solution of the classical Maxwell model (b ¼ 1) allows to draw a useful analogy and highlight the origin of the different terms in the fractional Maxwell solutions. For clarity, Table III lists the general solutions in Laplace space, for the initial conditions and models looked at in this paper. The expressions in this table use shorthand notations, which are important to list here to avoid ambiguity, as they will be frequently used below. In the Voigt model, In the Maxwell models, the commonly used parameters take the form:

A. Response to a stress impulse
In this case, the momentum is instantaneously changed at t ¼ 0 from 0 to some finite value, i.e., from stationary the mass is given an initial velocity via impact. The initial stress condition can be written as rðtÞ ¼ R 0 dðtÞ where R 0 has units of kg m À1 s À1 . The equation of motion is Laplace space is thenẽ The classical case of damped oscillations is equivalent to using the Voigt model, where an elastic restoring force and a viscous damping force act in parallel (viscous dissipation could either be due to internal or external friction). Taking the inverse Laplace transform of the solution in Table III gives the three modes of response depending on whether x is zero, real, or imaginary. These modes are shown in Table IV. These three solutions correspond to the classical cases of under-damped, over-damped, and critically damped oscillation, respectively. The response in this case of a Voigt model, as expected, settles to an equilibrium at zero strain, and has no long-term memory of it being subjected to an impulse. Now looking at the Maxwell model, and taking the inverse Laplace transformation of the solution listed in Table  III, we obtain see Appendix B for derivation. At t ¼ 0, the terms in the bracket cancel exactly, as we would expect. At long times, the oscillatory part decays and the strain settles at the equilibrium position e 1 ¼ R 0 =g [after simplifying the prefactor in Eq. (7)]. This is the unrecoverable deformation that happened through the dashpot creep in a Maxwell model. As with the Voigt model, we can identify three distinct regimes of underdamped oscillations for x 1 > k, slow overdamped relaxation for x 1 < k, and the fastest approach to equilibrium for x 1 ¼ k (critical damping). Now turning to the factional Maxwell model, taking the inverse Laplace transformation of the corresponding solution in Table III is more difficult due to noninteger powers; the full derivation is discussed in Appendix C. Once solved, the response of this viscoelastic system to a stress impulse takes the form where the values of shorthand parameters A, /, and the function qðxÞ are listed in Appendix C. The values of x and k are determined by the roots of s 2 þ s À1 s 2Àb þ x 2 1 ¼ 0, because for the inverse Laplace transformation of a function with simple poles, the only time dependent part is the exponential, i.e., the locations of the poles directly determine the frequency and decay rate of the oscillations. When b ¼ 1, we recover the classical Maxwell model. The integral term vanishes quickly and is mainly responsible for stabilizing the behavior at short times, in particular the asymptotic behavior of the power law as t ! 0. Equation (8) actually proves a conjecture made in [36], referring to a criterion for the crossover between liquidand solidlike behavior of a viscoelastic system.
There is a small technical point to be made here. In general, the integral with qðxÞ in Eq. (8) cannot be evaluated ModelẽðsÞ: Impulse stressẽðsÞ: Impulse strainẽðsÞ: Step stress analytically. In most cases, it converges well and can easily be found numerically, and the dependence on key parameters extracted explicitly. However, in some cases (usually when xs ( 1) the convergence is poor and the result could be less clear. To address this issue, here and in similar expressions below, we deliberately isolate the integral part to play as small a role as possible, bringing out the explicit long-time contributions which can be solved exactly. Figure 3 shows how the oscillatory behavior changes with the index b, reminding that at small fraction b the damping element is close to an elastic unit in its response (hence the significant oscillations are retained). In this plot, the parameters are chosen such that x 1 ¼ 1=2s, corresponding to the critical damping regime. Hence when b is close to one, the relaxation to the equilibrium strain is fast. The power-law determines the long time behavior and the equilibrium strain, with the decaying oscillations happening about this value. At short times, all solutions for a given b converge, as all start with the same velocity and have the same spring constant, G r . The system differs from the classical case in two key ways, first it has a long time memory of the impact owing to the power-law creep, whereas in the classical case the system is at rest as soon as the oscillations decay. The second difference is the system will always oscillate and there is no critical damping or overdamped situation in the strict sense, a consequence of there being no purely real solutions to s 2 þ s À1 s 2Àb þ x 2 1 ¼ 0 when 0 < b < 1. However, there are still values of x 2 1 and s which give the fastest decay of oscillations; this is discussed in more detail in Sec. V.

B. Response to strain impulse
This regime practically is analogous to plucking a string or a suspended mass, and is perhaps less common in standard rheological testing. The inverse Laplace transformation of the Voigt model solution illustrated in Table III is given by the standard decaying oscillation where / ¼ cos À1 ðx=x eq Þ, and as before: k ¼ 2g=l; Here, we will once again have the underdamped, overdamped, and critical damped regimes when x is real, imaginary, and zero, respectively. The strain rate can be simplified to a compact expression There is a subtle point here that the rate is not exactly out of phase with the strain, which is a direct consequence of the damping; as k ! 0, then cos À1 / ! 0 and they become properly 90 out of phase.
For the classical Maxwell model, the strain response is again Here, / ¼ cos À1 ðx=x 1 Þ, and k ¼ 1=2s; x 2 1 ¼ G r =l, and The expressions have the same form as in the Voigt model, although the k term has a different origin. It is clear why the Voigt model settles to zero strain in equilibrium, as there will always be a restoring force with the spring and dashpot in parallel. In this specific case of strain impulse, the dashpot of the Maxwell model does not extend in response to the initial instant deformation, which is taken by the spring; the system will then oscillate about zero strain, dissipating energy evenly above and below zero strain, eventually settling at zero strain.
The solution of the fractional Maxwell model is where the detailed derivation, solutions for x and k, and the expressions for parameters A, /, and qðxÞ are given in Appendix D. As with the classical Maxwell model, the solution returns to zero strain at equilibrium for the same fundamental reason of specific response to the strain impulse. In this case, the integral term has to decay slowly enough that it forces the strain remain below De at early times, countering the larger than one amplitude of the oscillations. However, as in the earlier case of stress impulse, this term must decay more quickly than oscillations themselves, leaving the decaying oscillations the only relevant response at longer times. Figure 4 shows that at b close to one, that is, for an almost Newtonian dashpot, the response is close to the  x real R 0 lx sinðxtÞe Àkt critical damping. Whereas at small b the fractional element is closer to an elastic unit by its nature, and the response is only a weakly decaying oscillation.

C. Response to step stress
This represents one of the common rheological experiments or situations, often called the creep compliance test. Taking the straightforward inverse Laplace transformation of the Voigt model solution for this initial condition (Table III), we obtain the decaying oscillation with respect to the equilibrium elastic deformation where the parameters are the same as in the earlier cases of the Voigt model. The solution to the classical Maxwell model reflects the equilibrium creep under constant stress e t ð Þ ¼ _ e r t þ e r 1 À 1 Again, the parameters are the same as in the earlier examples of Maxwell model. The characteristic relaxation time after the onset of stress is s ¼ g=G r ; e r ¼ r 0 =G r is the extension of the spring element in the Maxwell model, and _ e r ¼ r 0 =g is the extension rate of the dashpot element in this creep experiment. As we might expect, the transient regime has damped oscillations on top of the standard creep behavior of the Maxwell model. Now examining the fractional case of the step stress response of a material, whose intrinsic relaxation is described by the Mittag-Leffler function, the solution for the time-dependent deformation is As usual, the details of derivation, expressions for x and k, as well as A, /, and qðxÞ are given in the corresponding Appendix E. As in Eq. (8), we evaluate the integral containing qðxÞ numerically. To address a possibility of poor convergence, we isolate this integral part to play as small a role as possible, bringing out the explicit long-time contributions which can be solved exactly.
Using the analogy with the classical case, we denote e r ¼ r 0 =G r to be the constant contribution to strain from the extension of the spring under constant stress. Figure 5 illustrates the role of the fractional dashpot. Rather than linear (Newtonian) creep, we now find the power-law growth of deformation under the constant stress, following the t b law at long times. The second power-law term in Eq. (15), proportional to t 2ðbÀ1Þ , represents a transient decay (since b < 1). This term, together with the integral term of qðxÞ, controls the initial regime and account for the time it takes for the mass to accelerate until the baseline creep reaches the equilibrium regime. As with the classical case, this is a more significant effect when the mass is large and the viscosity is low, as it takes a long time to accelerate to the equilibrium velocity; however, this is only relevant when b is close to one. The response of a material with low b is closer to a purely elastic element, with oscillations taking a very long time to decay. This type of oscillatory behavior is commonly seen in creep experiments where the instrument inertia cannot be avoided, and has been observed in materials with fractional Maxwell behavior [13,15], where they fit the oscillations with a numerical model.

IV. RESPONSE TO DRIVEN OSCILLATIONS
In this section, we look at the response of our model viscoelastic materials to driven oscillations. As is typical in this case, we will focus on the equilibrium response, rather than transient solutions. At equilibrium, the solution for the amplitude of oscillations driven by an external stress r 0 sinðxtÞ takes the form As in other sections of this paper, we use notations: . Figure 6 shows the normalized strain oscillation amplitude, A, as a function of scaled frequency, x=x 1 for different values of fraction b. Plots (a) and (b) show the result for c ¼ 0:01 and c ¼ 100, respectively. One could think of low c meaning low viscosity, but high mass and/or relaxation modulus; high c represents low inertia spring and high viscosity situation. However, such interpretation is only strictly valid when b ¼ 1.
For small c, the dashpot is relatively weak and extends very easily, and as the spring and fractional dashpot are in series this means the spring plays a minor role in the oscillations. When b is close to one, with lower frequency the amplitude gets much larger as the dashpot can flow much further for a given stress in with a longer period of oscillation. This is essentially the power law response of a fluid. Consequently, as frequency tends to zero the amplitude will tend to infinity, just as we would expect for constant stress on the classical Maxwell model [see Eq. (14)]. This is different in the Voigt model, or any model with a spring in parallel to the dashpot, where the amplitude would tend to the value of the spring with no damping elements. As b gets lower, a resonance peak starts to rise, as both the value c b gets larger and the dashpot becomes a weak solid contributing to the elasticity, which also allows the spring to play a greater role as the magnitudes become comparable. As we have a quasispring and an ordinary spring in series, their elastic contributions will also add in series, so the effective spring constant will be lower than the contribution of the weakest element, which in this case is the fractional dashpot due to the small c. As b decreases and c b increases, the dashpot becomes more solid and the resonant frequency will increase, converging when b ¼ 0. At b ¼ 0, the fractional dashpot becomes a spring with modulus G r , and the resonant frequency becomes x 1 = ffiffi ffi 2 p and the amplitude is double that of a single spring on its own.
In the case when c is large and b is close to one, the dashpot will play relatively less role when compared to the spring and we see a typical spring resonance curve centered at x 1 . However, for very low frequencies, the amplitude still becomes large as even though the dashpot has a high viscosity-it is still free to extend given enough time. As b decreases, so will c b , and the fractional dashpot will become more comparable to the spring element, that is, contributing elastically. So we see both a lowering of the resonance peak due to greater movement of dashpot and higher damping, as well as a lowering of the resonant frequency. Once again, as b tends to zero the fractional dashpot will become a spring of modulus G r and the resonant frequency is x 1 = ffiffi ffi 2 p . The shift of resonant frequency with b is shown in greater detail in Fig. 7. Here, we plot the resonance frequency (defined as a position of the peak amplitude) as a continuous function of b for several values of c. In all cases, as b tends to zero, the resonant frequency becomes x 1 = ffiffi ffi 2 p , here c b ¼ 1 and we left with two springs of the same modulus in series, irrespective of the value of viscosity g. For c close to, or below the critical damping condition (c ¼ 1=2), there is no  resonant peak for large b, that is, when the viscous element in series is sufficiently small (perhaps counter-intuitively, in the Maxwell model higher viscosity means less damping, due to the series assembly). The value of b needed to produce a resonant peak becomes smaller with c, as the fractional dashpot needs to become more elastic before the spring in series can start storing energy, rather than energy simply being dissipating through the fractional dashpot.

V. CRITICAL DAMPING
In contrast to the case of steady driven oscillation above, in this section we discuss the frequency and decay rate of oscillations during the transient regime. We specifically focus on how much damping occurs in different viscoelastic models in comparable conditions. All initial conditions for the imposed stress give essentially the same transient response in terms of oscillation frequency and damping rate. What differs between each case is the amplitude and phase difference of the oscillations, as well as the equilibrium strain. The oscillation frequency and the damping rate are determined by the location of the poles on the complex plane of the Laplace transformation (i.e., the solutions of s 2 þ s Àb s 2Àb þ x 2 1 ¼ 0): There are always two solutions, which are complex conjugates of each other. The real part, Àk, determines the decay rate of oscillations, while the imaginary part of these solutions, x, determines the frequency of the transient oscillations (see Appendix C for details). Figure 8 shows the plots of the dimensionless decay rate, jkj=x 1 , and the dimensionless frequency, x=x 1 , as a function of the dimensionless variable c ¼ x 1 s. We have seen in Sec. IV that c essentially measures the relative "strength" of the dissipative element compared to the inertial spring element of the mechanical model.
The two plots in Fig. 8 show how, for a given natural frequency x 1 ¼ ffiffiffiffiffiffiffiffiffiffi ffi G r =l p , the decay rate and oscillation frequency change with changing the relaxation time s for a given b. Remember that in terms of the fractional Maxwell model, s ¼ ðb=G r Þ 1=b (see Appendix A for details). For the classical Maxwell model, when b ¼ 1, the peak decay rate is reached exactly at c ¼ 1=2 (i.e., the viscosity 2g ¼ ffiffiffiffiffiffiffiffi G r l p ). The value of the peak decay rate is k ¼ x 1 , meaning that transient oscillations decay over the half-period of the natural frequency. For c < 1=2, the oscillation frequency is zero, and there are no oscillations with an increasingly long steady decay rate (the classically overdamped regime). Conversely, for c > 1=2, the oscillation frequency rapidly converges to x 1 with an increasingly longer decay rate, corresponding to the classically underdamped regime.
However, these intuitive benchmarks become more obscure in the fractional relaxation case. On decreasing b, the decay rate has a much less defined peak, which becomes lower but also more spread out. This means that even though the maximum decay rate is lower, the significant damping occurs over a broad range of parameters (e.g., for systems with a different effective mass). As with the classical Maxwell case, the oscillation frequency increases with increasing c; however, with b 6 ¼ 1 there is no classical overdamped region-there always is a residual oscillation in a fractional viscoelastic system. Figure 9 shows the critical damping factor c Ã ¼ s Ã x 1 , at which there is a peak in the largest decay rate jkj in Fig. 8(a). In the classical case, this point also corresponds to the sharp transition between the overdamped and underdamped regimes. As b decreases, the damping peak moves to lower c. In the limiting case below b % 0:1 the peak is at very low  9. The value of c at which the "critical damping" occurs: c Ã ¼ s Ã x 1 , i.e., the position of peaks in Fig. 8(a), plotted as a function of continuously changing fraction b. Note that at low b the damping peaks become so diffuse that the whole notion of critical damping becomes ill defined. values of c, with a simple scaling relationship of c Ã ' ð1=2Þ 1=b . In general, for low b the peak value of damping becomes irrelevant, as the system acts primarily elastically with light damping that changes little over a broad range of xs.
The interesting behavior occurs for the intermediate range of fractionality, 0:1 < b < 0:9, where we still have significant damping but over a broader range of parameters. A viscoelastic material following the classical Maxwell relaxation behavior would be quite limited in terms of damping oscillations, as for a given s and G r , the reduced mass, l, would have be tuned very carefully to give significant damping (the peak around the critical damping point is very sharp). In the fractional Maxwell model, for a given s, G r , and b, the requirements on reduced mass are much less stringent for significant damping of oscillations to be present. This argument could be tuned to any other parameter from the relevant set, e.g., for a range of G r values for a fixed mass. This behavior is useful where the reduced mass (i.e., mass, area, and initial length of the construction) can change in the process, but damping is still required at a similar level.
A good example would be in the case of biological connective tissues, such as ligaments, tendons, or fascia. Each of these tissues has their own specific role and associated viscoelastic properties. For instance, tendon exhibits a lowexponent power-law relaxation behavior [37], and as such fits well with fractional Maxwell with a low b < 0:05. This agrees with the perceived biological function of a tendon, which is to passively store elastic energy allowing the muscle and soft tissues to dissipate energy at their own preferred rate without over-stretching [38,39]. It would be detrimental for a tendon function to have b much greater than 0.1: It needs to act as an almost purely elastic element for a whole range of parameters (reduced mass, modulus, viscosity, rate of deformation), and that would be difficult to achieve in a construct made of inherently dissipative soft-matter materials. But by making b very low, it can remain in this universally nondamping mode.
In contrast, soft tissues (such as fascia) have relaxation functions with higher values of b $ 0:2 À 0:4, meaning that they can damp significantly over a relatively broad range of reduced mass, without the need to directly alter their mechanical properties. This is useful in running or in response to impact, where oscillations after a stress impulse must be rapidly dissipated in soft tissue (e.g., before the next stride) to avoid resonance and reinforcement of oscillation amplitude [40]. Similarly, natural rubber typically has a power-law decay of relaxation function with b $ 0:6, and not surprisingly it is widely used in broad impact or vibration damping situations (from earthquake protection of buildings and resonance-proofing of bridges to vibration-insulation of cooling fans on circuit boards).

VI. CONCLUSIONS
Fractional viscoelasticity is a relatively established field, driven by the practical need to describe materials with complex rheological response in a universal manner, using as few fitting parameters as possible [3,25]. We were especially motivated by the work of Rossikhin and Shitikova [18,19], who pioneered many steps that we had to follow in this work. Specifically, they have correctly solved the problem of stress impulse [obtaining the transfer function in Eq. (8) and Appendix C]. We reproduced their solution here with a particular focus on underlying physical processes, separating oscillating and nonoscillating parts, and highlighting relevant regimes. Several other types of rheological experiment (initial conditions) studied here, and the analysis of driven oscillation and critical damping is new in this field. Each of these problems has relevance in a different experimental setting, and these were discussed in the respective sections above.
In this paper, we focused on one specific aspect of complex viscoelasticity, which involves inertial effects and damping; in particular, it was having qualitative and significant effects in damping of oscillations. There are several other areas where the fractional calculus (and specifically-Maxwell model) plays a big role, for instance in the problems of anomalous diffusion in such media [41,42] and in associated approach to ageing [43].
One may ask why we have paid relatively less attention to the versions of Kelvin-Voigt model, which many would associate with a basic viscoelastic solid (as opposed to the classical Maxwell model that inherently shows plastic flow at long times). First of all, such intuitive understanding of Maxwell model representing viscoelastic creep is no longer valid once b < 1: We have seen how the fractional viscoelastic element effectively represents the internal elasticity of the system. There is also a subtle and little known issue with Voigt model in the fractional case, as both Naber [20] and Rossikhin [18] have noticed: Once a fractional damping element is inserted, there is a very abrupt transition in behavior (between b ¼ 1 and b ¼ 0:99). In the nominally overdamped regime, one expects [and indeed finds in the Maxwell model, cf. Fig. 8(a)] the rate of decay to be slow. However, as soon as the small fractionality is added to the Voigt model, the rate of decay becomes increasingly faster as one increases the magnitude of g relative to ffiffiffiffiffiffiffiffi G r l p . At the same time, the fractional Voigt model is entirely valid and well-behaved (and thus corresponds to relevant experiments) in the underdamped regime. Our main point is that in fractional-relaxation systems, the Maxwell model is adequate for viscoelastic solids and predicts physically correct effects in oscillation and energy damping.

ACKNOWLEDGMENTS
The authors are grateful to Fanlong Meng and Samuel Bell for many useful discussions. This research was supported by the EPSRC Critical Mass Grant for Cambridge Theoretical Condensed Matter (EP/J017639).

APPENDIX A: RELAXATION FUNCTION OF THE FRACTIONAL MAXWELL MODEL
The relationship between stress and strain in the fractional Maxwell model arises from the series combination of the elastic element governed by r 1 ¼ G r e 1 , and the fractional dashpot element governed by r 2 ¼ bd b e 2 =dt b . As the elements are in series e ¼ e 1 þ e 2 and r ¼ r 1 ¼ r 2 in equilibrium, giving the stress strain relationship of where s b ¼ b=G r . The dimensionality of the parameter b must be a function of index b for the units to balance. The relaxation function is the stress response to a step strain, and it can be shown either through Laplace transform or directly through the analysis of fractional derivatives [25] that this is given by where s ¼ ðb=G r Þ 1=b is the characteristic relaxation time scale and E b ½z is the Mittag-Leffler function named in honor of the Swedish mathematician G€ osta M. Mittag-Leffler. The function is defined as the series expansion, which is convergent for the whole of the complex plane This is a generalization of the basic exponential function, since when b ¼ 1 one recovers Cðn þ 1Þ ¼ n! For a more detailed discussion, see Appendix E in the textbook [25].

APPENDIX B: IMPACT IN A MAXWELL MODEL: INVERSE LAPLACE SOLUTION
"Impact" is the condition that we referred to as stress impulse above. Starting from the equation of motion in Laplace space we first rearrange it to isolate singularities in denominator. The first step is The second term can be split further 1 s s þ k ð Þ 2 þ x 2 1 Finally, we obtain the equation in a form where we can easily look up standard Laplace transformsẽ The three terms above correspond to the following parts in the equation of motion in time: which can be assembled in compact form by joining the oscillating functions where / ¼ cos À1 ð2kx=x 2 1 Þ.

APPENDIX C: IMPACT IN A FRACTIONAL MAXWELL MODEL: INVERSE LAPLACE SOLUTION
Here, we consider the stress impulse situation again. The equation of motion in Laplace space can be manipulated into a form with isolated dimensionless singularities in denominator where 0 b 1; z ¼ ss and c ¼ x 1 s. The solution can then be written as the sum of two functions each representing the inverse Laplace transform of a term in Eq. (C1). For both functions, the inverse Laplace transform is part of the contour integral shown in Fig. 10, which due to the Cauchy theorem must be equivalent to the sum of enclosed residues. The paths along the negative real axis are branch cuts which account for the fractional power in the functions. The transformation involves calculating the following contour integrals: It can be shown that the functionf 1 ðsÞ has no residues, so only the branch cuts contribute. The branch cuts are the two path integrals from À1 to 0 and 0 to À1, where we substitute in z ! Àz ¼ z expðipÞ for the first integral and z ! Àz ¼ z expðÀipÞ for the second integral where we used the substitution x ¼ zt=s, and the identities C½a ¼ Ð 1 0 x aÀ1 e Àx dx and p=sin½ap ¼ C½1 À aC½a. It can be quickly checked that a Laplace transform of this result gives 1=z b .
The second part of the inverse Laplace transform off 2 ðzÞ has complicated residues and branch cuts and requires careful analysis. As before, the branch cuts are given by substituting z ¼ Àz ¼ z expðipÞ to the integral from À1 to 0 and z ¼ Àz ¼ z expðÀipÞ to the integral from 0 to À1 such that X Branchcuts¼ We are left with a purely real integral over x, with a function in the integrand given by Þþ 2x 2Àb c 2 þ x 2 ð Þcos pb ½ : (C6) Changing variables gives X Branch cuts ¼ À sin pb ½ t bÀ1 ps b ð 1 0 q ys t y Àb e Ày dy ; (C7) multiplying the numerator and denominator of q ys t À Á by t 4 . Note that as t ! 0, the function takes a limit q ! 1 and the branch cuts simplify to Àt bÀ1 =ðs b C½bÞ canceling the identical term in f 1 ðtÞ. This is a delicate but important point that prevents an unphysical divergence of inverse power-law solutions at t ! 0. As t increases, the value of q will decrease for any given y, and due to the exp½Ày, this results in the integrand getting small rapidly. When x > c or y > x 1 s 2 =t, q will once again tend to 1 owing to the y 4 term dominating over c. However, by this time exp½Ày will be very small. In summary, the role of the branch cuts here is to counter the large value of f 1 ðtÞ at early times, but then it quickly vanishes and has no influence on the equilibrium position or relevant regime of damped oscillations.
To calculate the residues, we must first find the location of the poles by solving z 2 þ z 2Àb þ c 2 ¼ 0. This is most easily done by substituting in polar form z ¼ re ih and separating into real and imaginary parts r 2 cos½2h þ r 2Àb cos½ð2 À bÞh þ c 2 ¼ 0; r 2 sin½2h þ r 2Àb sin½ð2 À bÞh ¼ 0: (C8) These equations can only be simultaneously be satisfied if the solution is to the left of the imaginary axis, and there are clearly two solutions which are complex conjugates of each other (i.e., at h and Àh). The pair of simultaneous equations, for ðr; hÞ, are more usefully represented as and r ¼ À In the region of p=2 < h < p; sin½bh is always positive and sin½2h is always negative, hence Eq. (C9) is only satisfied if p=2 < h < p=ð2 À bÞ where sin½ð2 À bÞh will be positive. The poles can be found to be first solving Eq. (C9) for the argument, h, for given values of b, s, and c, this can then be substituted in Eq. (C10) to find the absolute value. With the two poles (z 1 and z 2 ) known, the residues can now be calculated, as these are simple poles. The result is where k ¼ 1=2s; x 2 ¼ x 2 1 À k 2 ; x 2 1 ¼ G r =l, and s ¼ g=G r . e r ¼ r 0 =G r and is the extension of the spring element in the Maxwell model in a creep experiment, and _ e r ¼ r 0 =g and is the extension rate of the dashpot element in a creep experiment. The inverse is given by e t ð Þ ¼ _ e r t þ e r 1 À 1 x 2 1 s 2 1 À cos xt ½ e Àt=2s À Á À 1 xs 3 À 1 x 2 1 s 2 sin xt ½ e Àt=2s ¼ _ e r t þ e r 1 À 1 x 2 1 s 2 À e r x 1 x cos xt þ / ½ e Àt=2s ; where e t 0 B B B @ 1 C C C A ; where x ¼ x=x 1 and c ¼ x 1 s ¼ g= ffiffiffiffiffiffiffiffi G r l p .