Analysis of models for viscoelastic wave propagation

We consider the problem of waves propagating in a viscoelastic solid. For the material properties of the solid we consider both classical and fractional differentiation in time versions of the Zener, Maxwell, and Voigt models, where the coupling of different models within the same solid are covered as well. Stability of each model is investigated in the Laplace domain, and these are then translated to time-domain estimates. With the use of semigroup theory, some time-domain results are also given which avoid using the Laplace transform and give sharper estimates. We take the time to develop and explain the theory necessary to understand the relation between the equations we solve in the Laplace domain and those in the time-domain which are written using the language of causal tempered distributions. Finally we offer some numerical experiments that highlight some of the differences between the models and how different parameters effect the results.


Introduction
This paper offers a thorough introduction to mathematical tools to describe wave propagation in solids modeled with a wide collection of viscoelastic laws. Before we even attempt a general description of the models we will be addressing, let us emphasize what our goals are and what has not been tackled in the present paper. We aim for a unified mathematical description of a wide collection of known viscoelastic models, including basic well-posedness results. The models will include all classical viscoelastic wave models, fractional versions thereof, and couplings of different models in different subregions. The techniques that we will employ in the first part of the article (Sections 3 through 7) are those of Laplace transforms, understanding the influence of volume forces, normal stresses, and given displacements, as well as the strain-stress relation as transfer functions describing linear distributional processes in the time domain. These techniques are classical, although we will use them in a language that is borrowed from the recent literature of time domain integral equations. In a second part of the paper (Sections 8 to 10), we will introduce and use tools from the theory of strongly continuous semigroups to analyze the three classical models and some transitional situations where, for instance, classical Zener-style viscoelasticity coexists with pure linear elasticity in different subdomains, with smooth or abrupt transition regions. Our goals for the current piece of work are not in the realm of the modeling: we will analyze but not discuss known models, and we will not deal with physical justifications thereof. A particular issue where we will be very restrictive is the fact that we will only deal with solids moving from equilibrium (no displacement, strain, or stress) at time zero. There are practical reasons for this choice (since stress remembers past strain, it is not entirely justifiable to start the clock with known displacement and stress), but we are also restricted because of our analysis goals. In both parts (transfer function analysis and semigroup analysis) we will give a hint at how to deal with initial conditions.
Barring initial and boundary conditions that are needed to fully describe the model, our goal is the study of a linear elasticity equation in a bounded domain of d-dimensional space. Here u is the displacement field, upper dots denote time differentiation, σ is the stress tensor and f represents the volumetric forces. Linear strain ε = 1 2 (∇u + (∇u) ) determines stress through a generic convolutional law (we only display the time variable in the following formulas) where D is a time-dependent, possibly distributional, tensor-valued kernel. Following the careful description given by Francesco Mainardi in his monograph [19], the four classical models of viscoelasticity are named after Zener, Voigt, Maxwell, and Newton. The strainto-stress relationship is given by a differential equation as follows: (Linear elasticity) Here a is a non-negative function and C 0 and C 1 are four-index tensors satisfying hypotheses similar to the tensor that is used to describe linear elasticity (we will be covering heterogeneous anisotropic solids), with some additional conditions that will be introduced when we explain the models in detail. Formally speaking, the last three models can be considered as particular examples of Zener's model. However, they have very different properties. Newton's model is equivalent to a parabolic equation, as can be seen by substituting the formula for σ = C 1ε in the equation of conservation of momentum and integrating once in the time variable. This model therefore does not produce waves and will be ignored in the sequel. Voigt's model also gives an explicit strain-to-stress relation: if we substitute σ = C 0 ε + C 1ε in the dynamical equation we can see that we end up with a PDE of order three (there are terms with two derivatives in space and one in time). The models of Zener, Maxwell, and Voigt dissipate energy.
Let us now give a quick literature review, including some relevant work from the modeling and mathematical communities. Some of the monographs [18,19,1] contain a large collection of references that can be used for a deeper introduction to this fascinating area. For a very early introduction to linear viscoelasticity, see [10]. A generalized model that encompasses all classical models describe above four is examined in [1,Chapter 3]. An overview of the physics of viscoelasticity and its relation to rheology (fluid and soft solid flow) can be found in [29] and [23], where the latter also develops numerical methods for solving problems associated to viscoelasticity, while an overview of the mathematical theories and techniques, including the problem of waves propagating in viscoelastic media can be found in [8]. Viscoelastic models have also been formulated in the language of integral equations (see [32,8,12]), and electrostatic models like the Cole-Cole model in [16].
We wish to consider waves propagating in viscoelastic solids using both the classical models and those using fractional derivatives. The relation between fractional derivatives and viscoelastic models (including introductions to the Mittag-Leffler functions) can be found in [19,22], and [20] offers a short survey of the history of development of this theory. Mainardi, who it can easily be seen is at the center of much of the development and communication of waves in viscoelasticity, shows in [21] that the fractional relaxation process with constant coefficients is equivalent to a similar process governed by variable coefficient ODE. This relationship has also been explored in [16], where fractional derivatives are bypassed by using a method of Yuan and Agrawal to convert the equations to a system of first order ODE. The Zener model has been extensively studied in the context of waves with both classical [2,5,14] and fractional [26] time derivatives, as well as in the quasistatic case [30,28]. Other authors have explored the Maxwell and fractional Maxwell models [7], different variations of Maxwell models [6,9], and the Voigt model [5,12,15].
For some of our stability results, we make rigorous use of Laplace transforms for vector-valued distrivutions. Laplace transforms also appear in [2,8,9,10,26,28], used in the context of showing existence and uniqueness of solutions, justifications of models, exploration of the constitutive relations, or to simplify numerical implementation. In the context of stability analysis, a Green's function representation of the solution to the three dimensional wave problem is used in [8]. While we obtain some estimates in the time domain by making use of an inversion theorem for the Laplace transformation (in the spirit of the Payley-Wiener theorems), in some cases can make use of semigroup theory to obtain estimates directly. Similar analysis has been performed for waves in unbounded domains [5] and in bounded domains [14], where semigroup analysis is used to show the existence and uniqueness of solutions for waves in a Zener model.
While not explored in detail in the present work, we are also interested in the analysis and implementation of numerical schemes for the simulation of waves in viscoelastic materials. As mentioned earlier, [23] contains an overview of numerical methods for problems in viscoelasticity including finite element, boundary element, and finite volume formulations. Numerical implementation with finite elements has also been explored in [12] and [35] for the simulation and comparison to real world data, specifically blood flow in [35]. Also in the context of blood flow, [28] uses discontinuous Galerkin methods to simulate a quasistatic nonlinear 1D fractional Zener model. A DG method for a general linear quasistatic viscoelastic model is proposed in [30] and a priori error estimates are derived. Convergence of finite element methods for viscoelasticity is explored [13] and [15,14], where the first reference focuses on convergence in time while the latter two are concerned with optimal order of convergence in space. Coupling of elastic and viscoelastic subdomains is examined in [24] where boundary elements for the viscoelastic subdomain are coupled with finite elements for the elastic components, whereas in [34] a scheme involving only finite elements are used for the same problem and the two schemes are compared. Our paper is structured as follows. After introducing the general model (Section 2), we give a general framework for the viscoelastic material law as a transfer function (Section 3) and then move on to prove that the main classical models (Zener, Maxwell, Voigt), fractional versions of them, and combinations of different models in different subdomains, fit in our general framework (Section 4). Sections 5-7 contain the Laplace domain analysis of the model carried out as follows: first we do a transfer function analysis, then we give the general theory of how to understand transfer functions as Laplace transforms of distributional convolutions in the time variable, and finally we give estimates for the case of smooth (in time) data. In Sections 8-10, we start with a semigroup analysis of the classical models. Because we are striving for generality, we make an effort to include models where the classical viscoelastic models can degenerate into classical elasticity, which motivates a careful discussion of the closure process of a normed space with respect to a certain seminorm and how this affects the action of some operators. Section 9 gives a detailed treatment of Zener's model, using the tools of the previous section and wellknown results of the theory of strongly continuous semigroups in Hilbert spaces. Section 10 sketches the main changes that need to be made to the preceding analysis to study Voigt's and Maxwell's models. Finally, and just for the sake of illustration, we show some simulations for one and three dimensional models.
Before we proceed with the work at hand, let us give here some quick notational pointers. While the transient models we will be describing and analysing in this paper take values on spaces of real valued functions, the transfer function analysis will require the introduction of complex variables and complex-valued functions. To be on the safe side, all brackets and forms considered in this paper will be linear or bilinear, never conjugate linear or sesquilinear. We will write with no conjugation involved. The upperscript will be used for transposition of matrices, without conjugation. Note that A : B = 0 for all A ∈ C d×d sym := {A ∈ C d×d : A = A}, and B ∈ C d×d skw := {B ∈ C d×d : B = −B}. We will write M 2 = M : M. Given two Banach spaces X and Y , we will consider the space B(X, Y ) of bounded linear maps from X to Y with the operator norm. We will shorten B(X) := B(X, X).

An introduction to the model problem
The wave propagation problem will be given in an open bounded domain Ω ⊂ R d , whose boundary is denoted by Γ. In order to have a well defined trace operator in classical Sobolev spaces, we will assume that Ω is locally a Lipschitz hypograph, although this hypothesis can be relaxed as long as we have a trace operator. We will assume that Γ is decomposed into Dirichlet and Neumann parts, Γ D and Γ N , satisfying The inner products in the Lebesgue spaces (note that the latter is a space of symmetric-matrix-valued functions) will be respectively denoted by a well-defined Gelfand triple. The reciprocal duality product of the two fractional spaces on Γ N will be denoted with the angled bracket ·, · Γ N . We will consider the space for symmetric-tensor-valued functions In the above definition, the divergence operator is applied to the rows of the matrix valued function S. Following well-known results on Sobolev spaces, we can define a bounded linear and surjective operator γ N : H(div, Ω) → H −1/2 (Γ N ) so that the following weak formulation of Betti's formula holds.
Pending a precise introduction of the material law, which we will give in the Laplace domain in Section 3, we are now ready to give a functional form for the viscoelastic wave propagation problem. We look for u : [0, ∞) → H 1 (Ω) and σ : [0, ∞) → H(div, Ω) satisfying for all t ≥ 0 (2.1c) Here ρ ∈ L ∞ (Ω), with ρ ≥ ρ 0 almost everywhere for some positive constant ρ 0 , models the mass density in the solid, which at the initial time t = 0 is at rest on the reference configuration Ω u(0) = 0,u(0) = 0.
Upper dots are used to denote time derivatives. The data are functions What characterizes the viscoelastic model (for small deformations where the strain at time t can be described by ε(u(t))) is the existence of a strain-stress relation of the form The viscoelastic law (2.1e) is formally written in terms of a convolutional kernel D. In a first approximation, this kernel can be considered as a fourth order tensor (with some symmetric properties) depending on the time variable. As we will see later on, the most interesting examples arise when the causal convolution (2.1e) is a distributional one and D is described as a causal tensor-valued distribution of the real variable.

Viscoelastic laws in the Laplace domain
In this section we are going to give a precise meaning to the general viscoelastic law (2.1e). Instead of writing the convolutional law (2.1e) in the time domain, we are going to introduce a Laplace transformed model which we will analyze in detail. This model will then be used to justify a family of distributional models in the time domain. At this point, we consider the formal Laplace transform of the convolutional process (2.1e) and introduce C(s) := s L{D}(s).
(Note that multiplication by s takes care of time differentiation.) Laplace transforms will be defined in the complex half-plane The viscoelastic material tensor can be described as a transfer function through a holomorphic map Given M ∈ C d×d we thus have C + s → C(s)M ∈ L ∞ (Ω; C d×d ).  For each s ∈ C + , we take C(s) to be the smallest real number so that, almost everywhere in Ω, C(s)M ≤ C(s) M ∀M ∈ C d×d sym .

Hypothesis 3 (Boundedness).
There exists an integer r ≥ 0 and non-increasing and almost everywhere in Ω We can equivalently introduce this material tensor with a collection of holomorphic functions (the material coefficients) In this case, we just define C ijkl (s)m kl i, j = 1, . . . , d, and notice that (3.1) is equivalent to (3.4). When we write the material tensor in terms of coefficients, we can take as an upper bound of C(s) almost everywhere. With this point of view C(s) can be considered as an element of L ∞ (Ω; C d×d×d×d ). The general viscoelastic material law in the Laplace domain is σ(u) = C(s)ε(u). (3.5) Here and in the sequel, we will identify C(s) with the associated 'multiplication' operator C(s) ∈ B(L 2 (Ω; C d×d sym )), noticing that The associated bilinear form is a(u, w; s) := (C(s)ε(u), ε(w)) Ω .

Examples
Before detailing the main examples covered with our theory, let us introduce a definition that will make our exposition simpler. Let where c > 0 is a constant. For simplicity, in the future, we will write C ≥ c to refer to the last inequality and we will say that C is a steady Hookean material model, when conditions (4.1) are satisfied. The constant c > 0 will be called a lower bound for the model. Note that we can apply the model to complex-valued matrices C(M re + ıM im ) := CM re + ıCM im .
When hypotheses (4.1) are satisfied with c = 0, we will call C a non-negative Hookean model. For a steady Hookean model C, we will write  (d) If a ∈ L ∞ (Ω) satisfies a ≥ a 0 > 0 almost everywhere, then a C is a steady Hookean model.
Proof. Note first that CM ∈ C d×d sym almost everywhere for all M ∈ C d×d . Therefore, if N ∈ C d×d skw , then 0 = CM : N = M : CN ∀M ∈ C d×d , which implies CN = 0. Property (b) follows from (a). Properties (c) and (d) are straightforward.

Elastic models
In the case where we take a Hookean model C 0 and we consider the constant function C(s) ≡ C 0 , it is simple to see that the

Zener's classical viscoelastic model
We now consider the material law where a ∈ L ∞ (Ω) is strictly positive, C 0 and C 1 are steady Hookean material models, with C diff := C 1 − a C 0 ≥ 0, that is, almost everywhere This hypothesis makes C diff a non-strict Hookean model. As we will see in the direct time-domain analysis of Section 9, C diff is the diffusive part of the elastic model, while C 0 acts as a base or ground elastic model. We will make use of the formula The Laplace domain stress-strain relation can be written in implicit form σ + a s σ = C 0 ε(u) + sC 1 ε(u), corresponding to the differential relation in the time domain σ(t) + aσ(t) = C 0 ε(u(t)) + C 1 ε(u(t)), σ(0) = 0.
where a 0 = 1/ a −1 L ∞ (Ω) is a lower bound for a.
Proof. To prove (a), note first that and therefore where all the bracketed quantities in the right hand side are real. Taking real parts and noticing that Re (1 + as) −1 = 1 |1 + as| 2 (1 + aRe s) ≥ 0, the result follows.
Let then g 0 = C 0 ijkl and g 1 = C 1 ijkl . An easy computation shows that we easily estimate which proves the result.
The above exposition of Zener's model allows for full anisotropy. The isotropic viscoelastic model can be easily described with two variable coefficients. To do that we let λ, µ : C + → L ∞ (Ω) be holomorphic functions with the following properties being satisfied almost everywhere in Ω and for all s ∈ C + : We then define Examples of functions satisfying (4.4) can be found using a variant of Zener's model for viscoelasticity: let a, m µ , b µ , m λ , b λ ∈ L ∞ (Ω) be strictly positive (bounded below by a positive number) and such that Then (4.4). To prove the lower bound (4.4b), note that Re (sµ(s)) = Re

Fractional Zener models
In this section we explore models of the form C(s ν ) where C(s) is a Zener model and ν ∈ (0, 1). For fractional powers in the complex plane we will always take the principal determination of the argument, i.e., the one with a branch cut at the negative real axis. A fractional Zener model has the form where a, C 0 , and C 1 satisfy the same hypotheses as in Section 4.2. In the time domain, this corresponds to is a Caputo fractional derivative of order ν. Note that this fractional derivative coincides with the Riemann-Liouville fractional derivative of the same order, if we are assuming homogeneous initial conditions for all variables. This fractional derivative can also be defined as a distributional fractional derivative in the entire real line. Proof. To prove (a), using the same idea as in Proposition 4.2, we write We thus only need to show To see that, first observe which proves the result. The proof of (b) is a direct consequence of Proposition 4.2(b) and Lemma 4.4.

Maxwell's model
Maxwell's model is given by where C 1 is a steady Hookean model (with lower bound c 1 > 0) and a ∈ L ∞ (Ω) satisfies a ≥ a 0 > 0 almost everywhere, for some constant a 0 . In the time domain this gives again an implicit strain-to-stress relation Proof. Hypothesis 1 is easy to verify. Hypothesis 3 can be verified using the proof of Proposition 4.2 taking C 0 = 0. To prove Hypothesis 2 note that almost everywhere Note now that since and therefore, almost everywhere and for all s ∈ C + which proves the result. Proof. Hypothesis 1 is straightforward and Hypothesis 3 follows from the fact that as follows from Lemma 4.4. To verify Hypothesis 2 we first estimate almost everywhere. Using (4.6) and Lemma 4.4, we can easily bound At the same time (we have used Lemma 4.4 again) and the proof is finished.

Voigt's model
as a viscoelastic parameter model, where C 0 is a steady Hookean material model and C 1 is a non-negative Hookean model. In those parts of the domain where C 1 = 0, Voigt's model reduces to classical linear elasticity. In the time domain, this model gives an explicit differential expression for the strain-to-stress relationship Note that this can be plugged into the momentum equation yielding which shows that this model is a third order differential equation, although we admit the possibility that the third order terms vanish in some regions.
where c 0 is the lower bound of C 0 .
Proof. It is straightforward proof using the type of inequalities of the proof of Proposition 4.2.
Let Ω 1 , . . . , Ω J be non-overlapping subdomains of Ω such that Ω = ∪ J j=1 Ω j . Assume that C j is a viscoelastic model in the domain Ω j , satisfying the Hypotheses 1-3. Then defines a viscoelastic model in the full domain Ω.
This means, in particular, that we can combine all the above models (Zener, Maxwell, and Voigt) in their differential or fractional versions, with different fractional orders in different subdomains subdomains. Because all our definitions are distributional, whenever there is an interface (smooth or not) we are implicitly imposing continuity of the displacement (this is done by assuming u too take values in H 1 (Ω)) and of the normal stress (since σ takes values in H(div, Ω)).
Newton's model is a fourth choice among classical viscoelastic models, given by the law where C 1 is a steady Hookean model. Since the viscoelastic law in the time domain becomes σ(t) = C 1 ε(u(t)), the model can be simplified to (here g is an antiderivative of f and takes care of non-vanishing initial conditions if needed). Therefore, Newton's model becomes a parabolic equation for linear elasticity and does not produce waves. Since the interest of this paper is wave models, we will not investigate this simple model any further.

Transfer function analysis
We consider the energy norm, tagged in a parameter c > 0: Recall that the mass density function is ρ ∈ L ∞ (Ω) is strictly positive. By (4.3), it follows that Proposition 5.1. If ψ, φ and r are the functions and integer appearing in Hypotheses 2 and 3 and we define and (5.2b) follows from the fact that 1 ≤ |s|/Re s for all s ∈ C + . The asymptotic bounds (5.3) can be proved easily.
be a bounded extension operator and consider the solution of the elliptic boundary value problem These are d uncoupled scalar problems for each of the components of u. Using the Bamberger-HaDuong lifting lemma (originally stated in [3], see [31, Proposition 2.5.1] for a rephrasing in the current language), it follows that However, the solution of (5.4) minimizes ||| u||| c among all u satisfying γ D u = α. Therefore This finishes the proof.
The main theorem of this section studies the operator associated to the Laplace transform of problem (2.1). In fact, problem (5.6) below is the Laplace transform of (2.1) for data that as functions of time are Dirac masses at time equal to zero.
for a certain constant C depending on ρ and the geometry.
Proof. The solution of (5.6) for data (f , β, α) can be decomposed as the sum of the solutions for (f , 0, 0), (0, β, 0) and (0, 0, α). For the first one, we note that u ∈ H 1 D (Ω; C) satisfies b(u, w; s) = (f , w) Ω ∀w ∈ H 1 D (Ω; C), and we can use w = su as test function. Applying (5.2a), it follows that For the second one, we use the same argument to bound where we have used Korn's inequality and (5.1). For the final one, we write u = u + u 0 , where u = u(α, |s|) is the solution of (5.4) with c = |s| and u 0 ∈ H 1 D (Ω; C). Then Taking w = su 0 above, and using (5.2b), we can bound Therefore, Using (4.3) and the result follows.

Distributional propagation of viscoelastic waves
In this section we show how the transfer function studied in Section 5 (specifically in Theorem 5.3) is the Laplace domain transform of the solution operator for a distributional version of the viscoelastic wave propagation problem (6.5) (cf. Proposition 6.4 below). We start with some language about vector-valued distributions, borrowed from [31].

Background on operator valued distributions
Given a real Hilbert space X, its complexification X C := X + ıX is a complex Hilbert space that is isometric to X × X with the product by complex scalars defined in the natural way. The complexification X C has a naturally defined conjugation, which is a conjugate linear isometric involution in X C . The Lebesgue and Sobolev spaces of complex-valued functions that we have used in Section 5 are complexifications of the corresponding real spaces. If X and Y are real Hilbert spaces, the space of bounded linear operators B(X C , Y C ) can be understood as the subspace of B(X 2 , Y 2 ) formed by matrices of operators of the form This is easily seen to be isomorphic (with an equivalent but not equal norm) to the complexification of the Banach space B(X, Y ). (For the problem of the many possible equivalent complexifications of real Banach spaces, see [25]. where in the right hand side we use the natural conjugations of X C and Y C . With this definition Ax = Ax. An X-valued tempered distribution is a continuous linear map from the Schwartz class S(R) to X. We say that the X-valued tempered distribution h is causal, when the action of h on any element of S(R) supported in (−∞, 0) is zero. Causal tempered distributions have a well defined Laplace transform, which is a holomorphic function where the conjugation on the left-hand side is the one in X C . We will write h ∈ TD(X) whenever h is an X-valued causal tempered distribution whose Laplace transform satisfies where µ ∈ R and ψ : (0, ∞) → (0, ∞) is non-increasing and at worst rational at the origin, i.e., sup 0<x<1 x k ψ(x) < ∞ for some k ≥ 0. When, instead of a Hilbert space X and its complexification X C , we are dealing with bounded linear operators B(X C , Y C ), the conjugation in (6.1) is the one for operators between complexified spaces. Some pertinent observations and results: (a) If h ∈ TD(X) and T ∈ B(X, Y ) is a 'steady-state' operator, then T h ∈ TD(Y ). We can reverse the roles of T and h and show that if an operator valued distribution T ∈ TD(B(X, Y )) acts on a constant h ∈ X, it defines T h ∈ TD(Y ). (c) If h ∈ TD(R) and a ∈ X, then the tensor product a ⊗ h defines a distribution in TD(X).
If X and Y are Banach spaces and h ∈ TD(B(X, Y )), then the convolution product h * λ is well defined for any X-valued causal distribution λ, independently on whether it is tempered or not [33]. In the simpler case where λ ∈ TD(X), we have the convolution theorem L{h * λ}(s) = L{h}(s)L{λ}(s) ∀s ∈ C + , which can be used as an equivalent (and simple) definition of the convolution of h * λ, and we also have h * λ ∈ TD(Y ), as can be easily proved from the definition. is well defined.
Proof. Let us first recall that for s ∈ C + , we have defined the operator and that we have (see (3.6)) This means that C : C + → B(L 2 (Ω; C d×d sym )) satisfies the conditions (6.1) and (6.2) and therefore there exists C ∈ TD(B(L 2 (Ω))) such that L{C} = C. If we fix M ∈ R d×d sym , we can easily show that the Laplace transform of CM ∈ TD(L 2 (Ω)) is C(s)M. Finally, if u ∈ TD(H 1 (Ω)), then ε(u) ∈ TD(L 2 (Ω)) and the convolution C * ε(u) is well defined.
This general result about the weak (distributional) version of the viscoelastic wave propagation problem includes the possibility of adding non-homogeneous initial conditions for the displacement and the velocity, since they are just included in the volume forcing function f . These non-zero conditions will make the distributional solution of (6.5) nonsmooth at time t = 0 and therefore the estimates of Section 7 will not be valid.

Estimates in time
In this section we translate the estimates for the transfer function given in Theorem 5.3 into time-domain estimates. The following theorem rephrases [31, Proposition 3.2.2], which is an inversion theorem for the Laplace transform of causal convolution operators.
If λ ∈ C m+1 (R; X) is a causal function such that λ (m+2) is integrable, then the unique Y -valued causal function u such that L{u} = FL{λ} is continuous and satisfies The integrability of the (m + 2)-th derivative of λ can be relaxed to local integrability, but in that case we need to add a hypothesis concerning the growth of λ (m+2) (t) X , since we need to be able to take Laplace transforms in C + . In any case, since all the processes that we are dealing with are causal (the solution depends on the past values of the data, and never on the future ones), the result holds in finite intervals of time assuming only local integrability of the last derivative. Note also that if X and Y are complexifications of real spaces, the process in the time domain is real-valued for real-valued data when we have F(s) = F(s) for all s.
We now consider Sobolev spaces where X is any Hilbert space.
for some constants C 1 , C 2 , C 3 . If additionally then σ ∈ C([0, ∞); L 2 (Ω)) and we have the bound Proof. The result is a more or less direct consequence of Corollary 7.2, using Theorem 5.3 for the Laplace domain estimates, and going through the language of Propositions 6.2 and 6.3. Since the problem is linear, it can be decomposed as the sum of problems with data (f , 0, 0), (0, α, 0), and (0, 0, β). We follow the notation of Proposition 6.2 and define nine instances of spaces and operators to apply Corollary 7.2: we list the spaces X and Y (before complexification), as well as the value of m and µ and the function ϕ in the estimate (7.1). We separate the operator S(s) in Proposition 6.2 as a sum of three operators S(s) = S f (s) + S α (s) + S β (s).
We will also use the embedding of H 1 (Ω) into L 2 (Ω). The following table lists all nine operators: This proves the continuity properties for u and σ as well as the estimates (7.3a) and (7.3c). Note that zero initial values hold whenever the output function is continuous (see Corollary 7.2). To obtain the bounds for u(t) Ω we can use a simple shifting argument, since if (f , α, β) →u (i.e., we use the operators multiplied by s and with values in L 2 (Ω)), then (f (−1) , α (−1) , β (−1) ) → u.
In general (see Proposition 5.1)

Technical work towards a time domain analysis
In Section 9, we give a different analysis, based on the theory of C 0 -semigroups of operators, of the viscoelastic wave propagation for Zener's model, including situations where part of the domain is described with a purely elastic material. This requires a certain amount of preparatory work, which we will present in full detail. To avoid keeping track of constants, in this section and in the next we will use the symbol to absorb constants in inequalities that we do not want to display. Let C diff ∈ L ∞ (Ω; R (d×d)×(d×d) ) be a non-negative Hookean model, which will play the role of the diffusive part of the viscoelastic model. We will identify the tensor with the operator which is bounded, selfadjoint, and positive semidefinite. We also consider a strictly positive function a ∈ L ∞ (Ω), and the associated multiplication operator Note that m a C diff = C diff m a . We next consider the following seminorm in L 2 (Ω) Proposition 8.1. The following properties hold: (a) |S| c S Ω for all S ∈ L 2 (Ω).
(c) |S| c = 0 if and only if S ∈ ker C diff .
Proof. Since m a and C diff are bounded linear operators in L 2 (Ω), the proof of (a) is straightforward |S| 2 c = (aC diff S, S) Ω ≤ m a C diff S Ω S Ω S 2 Ω . To prove (b) note first that since m 1/a is bounded, then and therefore, using (a) This finishes the proof.
The dynamics of the viscoelastic model will allow for subdomains where energy is conserved and subdomains where the diffusive part of the model is active. The fact that we allow for transition areas between purely elastic and strictly diffusive models forces us to take the completion of the space L 2 (Ω) with respect to the seminorm | · | c . This standard process requires first to eliminate the kernel of the seminorm and move to its orthogonal complement. We thus consider the subspace which is clearly closed in L 2 (Ω). If T ∈ M and |T| c = 0, then T ∈ (ker C diff ) ⊥ ∩ ker C diff = {0} (Proposition 8.1(c)) and therefore | · | c is a norm in M . Note also that if S ∈ L 2 (Ω), we can decompose An easy argument about extensions and the fact that M = (ker C diff ) ⊥ prove that We now consider b ∈ L ∞ (Ω) and the associated multiplication operator m b : L 2 (Ω) → L 2 (Ω) that we recall was given by (c) (S, RE) c = ( C diff m a S, E) Ω for all S ∈ M and E ∈ L 2 (Ω).
Proof. It is simple to see that Πm b = m b | M Π and therefore which proves (a). If S ∈ M and note that We can always write T = χ Ω 1 T + χ Ω 2 T and note that by (8.3b) If C diff T = 0, then χ Ω 1 T : T = 0 by (8.3a) and therefore T = 0 almost everywhere in Ω 1 .
However, now, due to (8.3b), we have and therefore M = M . In this case R : L 2 (Ω) → M is just the restriction to Ω 1 of matrixvalued functions defined on Ω, C diff is the restriction of the action of C diff to functions defined on L 2 (Ω 1 ; R d×d sym ), and the same happens to multiplication operators m b . This simple situation arises when the domain is subdivided into two parts: in one part we will deal with a purely elastic medium (C diff = 0), while in the other part we will handle a viscoelastic medium that is 'strictly' diffusive, as expressed in (8.3a).

Semigroup analysis of a general Zener model
In the coming two sections we will use some basic results on C 0 -semigroups in Hilbert spaces, namely the Lumer-Philips theorem (characterizing the generators of contractive semigroups in Hilbert spaces as maximal dissipative operators) and existence theorems for strong solutions of equations of the forṁ where A : D(A) ⊂ H → H is maximal dissipative and F : [0, ∞) → H is a sufficiently smooth function. These results are well known and the reader is referred to any classical book dealing with semigroups (for instance, Pazy's popular monograph [27]) or to the chapters on semigroups in many textbooks on functional analysis.
Consider now the space Take now (f , F, G) ∈ H, solve the coercive variational problem and define (see Proposition 8.2(a)) E :=ε(u) + F, so that S + m a S = Rε(u) + m a G and therefore At the same time, by (8.2) and Proposition 8.2(a), we have which is equivalent to Summing up, we have (u, E, S) ∈ D(A) and (u, E, S) = A(u, E, S) + (f , F, G).
Corollary 9.3. If α, β, and f satisfy the conditions of Theorem 9.2, (u, E, S) is the solution of (9.2), and we define the pair (u, σ) satisfies the bounds for all t ≥ 0 As a consequence σ(0) = 0.
Proof. If we solve problem (9.2) with (ḟ ,α,β) as data, and we integrate from 0 to t, we obtain the solution of (9.2) which is therefore an element of C 2 ([0, ∞); H). This is enough to prove everything else.
The estimates of Corollary 9.3 greatly improve those of Section 7 (see Theorem 7.3) in two aspects: less regularity required for the data, and constants independent of the time variable. There are three particular cases included in the analysis of this section that are worth paying special attention to.
there is no need to use the completion process since M = M = L 2 (Ω) and then C diff = C diff , m a = m a , and R is the identity operator. The space is now endowed with the norm which is equivalent to the usual norm. This makes the analysis of this strictly diffusive viscoelastic problem much simpler.
(c) When C diff = 0, we have M = M = {0} and the third equation and unknown do not play any role. In this case the operators ±A are maximal dissipative and therefore, A is the generator of a group of isometries in H, i.e., this model is conservative. This should not be a surprise, since in this case we recover a first order formulation of the classical linear elastic wave equation.
The introduction of non-zero initial conditions for the most general version of this model is not trivial. When the model is strictly diffusive (case (b) in the above discussion, i.e., when (9.14) holds), we are allowed to impose initial conditions which would be the natural ones for the formulation (9.11). This is done by modifying the system (9.2), using initial conditions u(0) = u 0 , E(0) = 0, S(0) = 0, and adding ρ v 0 to the right-hand side of (9.2a) and G 0 , with to the right-hand-side of (9.2c). The general case is much more complicated, given that the initial conditions for σ(0) must match C 0 ε(u 0 ) in purely elastic subregions.

More semigroup analysis
We are now going to take advantage of the preparatory work of Section 8 to give a quick view of the formulations and estimates that can be obtained for the Maxwell and Voigt models.

Maxwell's model
Maxwell's model can be understood as the particular case of Zener's model when C 0 = 0 and C diff is strictly positive. However, its analysis is not included in the treatment given in Section 9, due to the fact that C 0 was used to define the norm of the space H. We thus start again, with a new space which is equivalent to the usual norm. The domain of the operator The operator A is maximal dissipative. The proof of surjectivity of U → U − AU starts with the solution of the coercive problem for given (f , G) ∈ H. (Note that the strict positivity of C diff is key for this argument to hold.) This is followed by the definition of S = (1 + a) −1 (ε(u) + aG).
Using the operator A in an equation of the form (9.6), we can prove that the hypotheses of Theorem 9.2 are sufficient to provide a solution (u, S) ∈ C 1 ([0, ∞); H) of the problem with vanishing initial conditions. Introducing the stress tensor we have a solution of a.e. − t, with vanishing initial conditions. The estimates of Corollary 9.3 hold for this model as well.
A combination of Zener's and Maxwell's models is also available. It requires an even more general framework so that C 0 and C diff can vanish on separate parts of the domain as long as a certain combination stays strictly positive. (See the variational problem (9.1) that is solved as a starting step to prove maximal dissipativity. As long as C 0 + m −1 1+a C diff is a Hookean model, everything else will work.) The analysis would require a completion process with respect to the seminorm (C 0 ·, ·) 1/2 Ω and the corresponding restriction operator. This is a simple (while a little cumbersome) extension that the reader can do to prove their handle of the techniques developed above.

Voigt's model
The analysis of Voigt's viscoelastic model (Zener's model with a = 0, C 0 strictly positive and C diff ≥ 0), including areas transitioning to classical linear elasticity, follows from a simple modification of the ideas of Section 9. In Voigt's model C diff = C 1 plays the role of a dissipative term. The semigroup analysis of this model is slightly different in involving a second order differential operator. Like in Maxwell's model, there is no need to involve a completion process to handle transitions to classical linear elasticity. We now consider the following ingredients: A simple argument shows that If we take (f , F) ∈ H, solve the coercive problem and define E := ε(u) + F, it is easy to prove that (u, E) ∈ D(A) and (u, E) − A(u, E) = (f , F). Therefore, A is maximal dissipative. Going carefully over the proof of Theorem 9.2, it is easy to see that with the same hypotheses on the data f , α, and β, we have a unique (u, E) ∈ C 1 ([0, ∞); H), vanishing at zero, and solving The associated stress tensor is defined by σ(t) := C 0Ė (t) + C diff ε(u(t)) = C 0 ε(u(t)) + C diff ε(u(t)) and the resulting pair (u, σ) is a solution to (9.11) (with a = 0) satisfying also the estimates of Corollary 9.3.

Some experiments
We now present some numerical experiments of the various viscoelastic models that we described in Section 4. We use finite elements for space discretization and a trapezoidal rule-based convolution quadrature (TRCQ) for time discretization [4,11,17]. The numerical experiments will be divided into three groups. First, we investigate 1D uniaxial wave propagation. Through these experiments, we observe how the viscoelastic behavior is dependent upon the choosing of parameters in the constitutive equations. In the second group of experiments we compare the behaviors of 1D uniaxial wave propagation in elastic, classical viscoelastic, fractional viscoelastic, and heterogeneous models by plotting their 2D space-time contour graphs. In the heterogeneous model, we decompose the region into two subdomains with different viscoelastic models, where the reflection and refraction of waves can be observed at the transition interface. Finally, we present the 3D simulation of a viscoelastic rod. Similar to the heterogeneous domain in the previous experiments, the rod is decomposed into two different subdomains. The snapshots we present show how the simulation accurately captures the memory and relaxation effects of the rod under a sudden change in displacement. The numerical analysis of the discretization schemes employed in this section is the goal of future research. Tests have been performed in sufficiently refined space-and-time meshes to obtain some sort of eye-ball convergence to a solution.

1D experiments
For simplicity, in the one-dimensional examples we will use traditional PDE notation, as opposed to the notation of evolutionary equations used throughout the paper. We first present numerical experiments for different fractional models in one dimension where the constitutive relation that determines the model is defined through σ and ρ. We use the window function displayed in the left of Figure 1 as Dirichlet boundary data at x = 0, while we take a homogeneous Neumann boundary condition at x = 1. The strain-to-stress relation is given by a general formula σ + a ∂ ν t σ = C 0 u x + C 1 ∂ ν t u x , (11.1e) for parameters C 0 , C 1 , a and ν to be determined. For the discretization in space we use P 4 finite elements on a mesh with 513 subintervals of equal size. Discretization in time is carried out using a TRCQ with 10,240 time-steps of equal size in the interval [0, 40]. For the first two sets of experiments, we implement the Dirichlet boundary condition g(t) in Figure 1 at x = 0 and observe the evolution of u(1, t). We use: (a) the values given in Table 1, varying C 1 , to create the results of Figure 2, and (b) the values in Table  2, varying ν, for the results of  The parameters used to create the plots given in Figure 2. The first choice of C 1 for the Zener and Voigt models reduce the model to linear elasticity.  Table 1.
In Table 1 the values for C 1 are chosen so that the parameter C diff , which controls the diffusion, takes same values for the three different models, except for Maxwell, where C diff = 0 results in a model that is identically zero. To avoid this while still getting comparable results, we use C 1 = 0.05 for the first value in the Maxwell model. In Figure  2, as we increase C 1 , all three models show a faster energy dissipation. Compared to the Zener and Voigt models, the Maxwell model exhibits less oscillations as a response to the Dirichlet boundary condition.

Zener
Maxwell Voigt C 0 1.5 0 1.5 C  The parameters used to create the plots given in Figure 3.  Figure 3: Effect of changing ν using the parameters given in Table 2.
In Figure 3, we observe that the decreasing the fractional power ν leads to a slower rate of energy dissipation. We also notice that the change of ν does not have any obvious effect on the frequency of the oscillations.
We now change the boundary conditions (11.1b) and (11.1c) to where h is the function in the right of Figure 1. Once again, we vary ν to see the effects that the fractional order of the derivative has on the Zener, Maxwell, and Voigt models respectively. The parameters we choose for this experiment is given in Table 3 and we again plot u(1, t) in Figure 4.   Here, both Zener and Voigt models show exponential rate response to the suddenly applied traction h(t) and then converge to some equilibrium state. The Maxwell model, as expected, shows a linear rate of creep when the traction is constant in time. As the fractional power ν decreases, we observe that the amplitude of oscillations is larger for the Zener and Voigt models, on the other hand we see the creep rate is slower for the Maxwell model.

Spacetime plots
We now study space-time contour plots of the displacement solution of (11.1) with Zener, Maxwell and Voigt models. We show three simulations focusing on one of these models where in each experiment we compare it with its fractional version, and observe its behavior in a heterogenous domain coupled with an elastic model. When we work on a heterogenous domain we split the interval [0, 1] into [0, 1/2) and [1/2, 1], where the first half is elastic and the second half is one of the models we are comparing: Zener, Maxwell or Voigt. We show space-time plots of the displacement corresponding to different models side by side where elastic model is included in all cases for the sake of comparison. We implement two different signals as Dirichlet boundary condition: a single pulse and a train of pulses (see Figure 5). For each experiment we display eight plots where the first four are the results of a single pulse while the last four are the results of the same experiment but for a train of pulses. Our first test focuses on the Zener model where we utilize the parameters given in Table 4, and Figure 6 shows the outcome of this experiment. Here, the contrast between the energy conservative elastic model and the dissipative Zener model can easily be seen. We also notice that fractional Zener model displays slower dissipation than the classical model. In the heterogenous domain we observe a reflection and refraction of waves at the interface x = 1/2. Although the elastic model is conservative, in the case of a heterogenous domain we see the dissipation of the Zener model affecting the coupled system with a loss in wave amplitude. We also observe that in the results whose Dirichlet boundary condition is a train of pulses (the four panels on the right), the solution enters quickly into a time-harmonic regime.    Table 4 where first and last four subplots correspond to the signals shown on the left and right of Figure 5.
We perform a similar comparison for the Maxwell model using the parameters given in Table 5 and collecting the results in Figure 7. When a single pulse is used, we notice that the waves in the Maxwell model exhibit dissipation. Moreover, the fractional Maxwell reveals less dissipation with a little dispersion. In the heterogenous domain the reflections at x = 1/2 are less obvious for both of the input signals when comparing to the corresponding Zener simulations.     Table 5 where first and last four subplots correspond to the signals shown on the left and right of Figure 5.
Lastly, we demonstrate the space-time plots corresponding to the Voigt model in Figure  8 using the parameters from Table 6. Comparing to Zener and Maxwell, we observe that this model displays more dispersion. In particular, this dispersion is on dramatic display in the heterogenous domain where the reflections are clearly seen in the elastic part, while the dispersion occurs in the Voigt subdomain.     Table 5 where first and last four subplots correspond to the signals shown on the left and right of Figure 5.

3D numerical simulation
We now present a numerical simulation for viscoelastic waves propagating in the parallelepiped Ω = (0, 1) × (0, 10) × (0, 1) with a Dirichlet boundary on one of the small faces Therefore when y < 5, the model is purely elastic and when y ≥ 5 the model is a Zener viscoelastic model. For the discretization in space we use P 2 finite elements on a mesh of 30,720 tetrahedra obtained by partitioning a uniform quadrilateral mesh of 8 × 80 × 8 elements. Discretization in time is done by TRCQ using 500 time-steps over the interval [0, 50]. Figures 9 and 10 show snapshots of the displacement from the simulation, where the coloring in the first one exhibits the norm of the stress (averaged on each tetrahedron).
We remark that in both figures, while we choose the same snapshots, the snapshots are not uniform in time. This is so we can highlight some of the more interesting aspects of the simulation which occur earlier in the time interval.  In the top rows of these figures, we observe the elastic waves, generated by the sudden deformation at y = 0, propagating along the rod in the y direction. We also note that the elastic part of the rod (y < 5) responds quickly to the sudden deformation and adjusts to the new displacement (enforced by the Dirichlet boundary condition) in about 70 time-steps, while the viscoelastic part of the rod (y > 5) shows a much slower response to the change of the displacement taking around 400 time-steps to adjust.