Lack of superstable trajectories in linear viscoelasticity: a numerical approach

Given a positive operator A on some Hilbert space, and a nonnegative decreasing summable function μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document}, we consider the abstract equation with memory u¨(t)+Au(t)-∫0tμ(s)Au(t-s)ds=0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \ddot{u}(t)+ A u(t)- \int _0^t \mu (s)Au(t-s) ds=0 \end{aligned}$$\end{document}modeling the dynamics of linearly viscoelastic solids. The purpose of this work is to provide numerical evidence of the fact that the energy E(t)=(1-∫0tμ(s)ds)‖u(t)‖12+‖u˙(t)‖2+∫0tμ(s)‖u(t)-u(t-s)‖12ds\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {{\textsf{E}}}(t)=\Big (1-\int _0^t\mu (s)ds\Big )\Vert u(t)\Vert ^2_1+\Vert \dot{u}(t)\Vert ^2 +\int _0^t\mu (s)\Vert u(t)-u(t-s)\Vert ^2_1ds \end{aligned}$$\end{document}of any nontrivial solution cannot decay faster than exponential, no matter how fast might be the decay of the memory kernel μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document}. This will be accomplished by simulating the integro-differential equation for different choices of the memory kernel μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} and of the initial data.

μ(s)Au(t − s)ds = 0 modeling the dynamics of linearly viscoelastic solids. The purpose of this work is to provide numerical evidence of the fact that the energy E(t) = 1 − t 0 μ(s)ds u(t) 2 1 + u(t) 2 + t 0 μ(s) u(t) − u(t − s) 2 1 ds of any nontrivial solution cannot decay faster than exponential, no matter how fast might be the decay of the memory kernel μ. This will be accomplished by simulating the integro-differential equation for different choices of the memory kernel μ and of the initial data.

The physical model
Within the theory of uniaxial deformations in isothermal viscoelasticity, the dynamics of a homogeneous isotropic linearly viscoelastic solid, occupying a bounded volume ⊂ R 3 at rest, is described by the integrodifferential evolution equation (see, e.g., [3,10,18,23,25,29,37]) The unknown u = u(x, t) : × [0, ∞) → R describes the axial displacement field relative to the reference configuration . Throughout the paper, the space variable x will be always omitted. The function κ appearing in the convolution integral has the form where κ ∞ > 0, and the memory kernel μ is a (nonnull) nonnegative, decreasing, piecewise smooth, and summable function defined on R + = (0, ∞). In particular, we allow μ to have a finite number of discontinuities (downward jumps). These assumptions imply the convexity of κ. The values κ(0) > κ(∞) > 0 represent the instantaneous elastic modulus, and the relaxation modulus of the material, respectively. In more generality, in this work we will consider an abstract version of (1.1), obtained by replacing the operator − with Dirichlet boundary conditions, acting on the Hilbert space L 2 ( ), with a strictly positive selfadjoint linear operator A, acting on some real Hilbert space H , with dense domain D(A) ⊂ H .

Notation
We agree to denote by ·, · and · the inner product and norm on H , respectively. We also introduce the Hilbert space V = D(A 1/2 ), endowed with the inner product and norm u, v 1 = A 1/2 u, A 1/2 v and u 1 = A 1/2 u .
Observe that for the particular case of (1.1)-(1.2), the space V is nothing but the Sobolev space H 1 0 ( ). Setting for simplicity κ(0) = 1, which translates into the constraint we end up with the abstract Volterra evolution equation of which the boundary value problem (1.1)-(1.2) is just a particular instance. Nonetheless, the abstract equation (1.4) covers also different models: for example, the one employed in the description of the vibrations of thin viscoelastic rods, namely, or of viscoelastic plates, namely, with ν > 0 small (see [29,31]). Both equations can be given the form (1.4) by setting in the first case and in the second one Incidentally, the first occurrence of A is (more precisely, extends to) a bounded linear operator on L 2 ( ). The abstract equation (1.4) is complemented with the initial conditions given at the initial time t = 0 where u 0 and v 0 are assigned data. It is well known that for all initial data (u 0 , v 0 ) ∈ V × H , problem (1.4)-(1.5) has a unique weak solution Such a solution is actually a strong one, namely, whenever the initial data (u 0 , v 0 ) belong to the more regular space D(A) × V . The natural energy of the system is defined by and it can be interpreted as the sum of the mechanical energy 6) and the memory energy which is ultimately responsible for the dissipation. Note that, for large t The dissipativity of (1.4) is witnessed by the fact that the E(t) is a decreasing function of t (see, e.g., [8,10,18]; but see also the next formula (1.11)). Equation (1.4) is sometimes referred to as a model of (linear) viscoelasticity of Volterra type. A different and perhaps more accurate model is the one of viscoelasticity with infinite memory, which amounts to considering the convolution integral on the whole R + , to wit,ü (1.8) In this case, however, in order to compute the convolution integral one has to know the values of u for all times before the actual time t. Accordingly, being the initial time conventionally set at t = 0, the function u for negative times is assumed to be a known initial datum, which need not solve the equation. Indeed, (1.8) is assumed to hold only for t > 0. As pointed out in the seminal work [12] (see also [22]), one way to handle (1.8) is to introduce the auxiliary variable which keeps track of the past history of u, and is assigned as an initial datum at t = 0. Hence, assuming the initial condition where η 0 = η 0 (s) is an assigned datum, we immediately see that (1.9) With this choice, and recalling (1.3), equation (1.8) transforms intö the Lebesgue space of square summable V -valued functions on R + with respect to the measure μ(s)ds, we introduce the Hilbert space 2 1 ds.
Then, system (1.9)-(1.10) is known to generate a strongly continuous linear semigroup S(t) of solutions 1 acting on H, in the sense of [15,35]. This means that, for any initial datum u 0 = (u 0 , v 0 , η 0 ) ∈ H, there is a unique solution

Remark 1.1 The infinitesimal generator of S(t) is the linear operator A on H acting as
with domain Here, T is the infinitesimal generator of the right-translation semigroup on M, defined as where the prime denotes the distributional derivative. Accordingly, system (1.9)-(1.10) can be equivalently viewed as the differential equation on H d dt The energy at time t corresponding to the initial datum u 0 reads and is a decreasing function of t, for it satisfies (for any regular initial datum) the differential inequality where s n are the discontinuity points (if any) of μ (see [12,33]). In the semigroup language, this means that S(t) is a contraction semigroup, i.e., the operator norm of the semigroup is less than or equal to 1 for all t ≥ 0. The analysis of S(t) and its asymptotic features has been carried out by several authors since the Seventies of the last century. To mention just a few of them, we quote the works [2, 12, 19-21, 30, 32-34] (stability and exponential stability issues), [4,13] (spectral properties), [6,7] (models with time dependent kernels), [14,16,17] (modeling aspects), [26][27][28] (semigroup approach). The link between the semigroup S(t) and the original Volterra equation (1.4) is expressed by the following proposition.

Indeed, defining the function
it is readily seen that system (1.9)-(1.10) can be equivalently written as (1.12) This is nothing but (1.4), in presence of a time-dependent forcing term. Notice that the function G is identically zero whenever u 0 = (u 0 , v 0 , u 0 ). In this case, we recover viscoelasticity of Volterra type. Besides, observing that the third (or memory) component we obtain (as it should be) the equality E(t) = F(t).
2 The Decay rate of the energy

Exponential stability
Being S(t) a contraction semigroup, either its operator norm is always equal to 1, or it goes to zero as t → ∞, and the decay is necessarily (at least) of exponential type (see [15,35]). In terms of the energy F, this means that there exist constants M ≥ 1 and ω > 0 such that, for all initial data in H, When (2.1) holds, the semigroup is said to be exponentially stable. The issue of the exponential stability of S(t) has attracted the attention of several authors since the appearance of [12], and it has been first proved within the sufficient condition that the memory kernel μ, besides our general hypotheses, satisfies the differential inequality for some δ > 0 (see [17,26,27,30]). Later, in [2], it has been shown that a necessary condition in order for S(t) to be exponentially stable is that for some C ≥ 1 and δ > 0, which turns out to be equivalent to (2.2) when C = 1. The necessary and sufficient condition has been established only in the more recent work [34], where the following theorem is proved.

Remark 2.2
If the operator A is bounded, the result is slightly different: exponential stability occurs if and only if (i) holds, with the only exception of a certain class of flat kernels, called resonant, where trajectories with conserved energy arise (see [2,34]).
Summarizing, with regard to the exponential stability of the semigroup the picture is nowadays completely clear. It is worth noting that (2.3) implies an exponential decay of the kernel μ of rate δ (at least). The situation is a little bit more complicated if we restrict our attention to the decay of the energy E of the Volterra equation (1.4). Here, lacking a semigroup structure, decay types other than exponential are possible. Without claiming to give the result in full detail, roughly speaking what happens is that the energy E and the kernel μ share the same decay type, but this occurs up to the decay of exponential type (see [5,8]). For faster decays, no results are available, as remarked in [5].

Superstability
We now come to the core of the present work, that is, the existence of nontrivial trajectories which converge to zero faster than any exponential. To this end, we define the (exponential) decay rate of the energy 2 of the semigroup S(t) to be the nonnegative number Thus, S(t) is exponentially stable if and only if ω > 0, whereas if ω = 0 then it has unitary operator norm for all t. [4], it is actually possible to give a quantitative version of Theorem 2.1, which says that if (2.1) holds for some ω > 0, then condition (2.3) is necessarily satisfied for some δ ≥ ω.

Remark 2.4 As shown in
This motivates the following definition.
Accordingly, in view of the remark above, if S(t) is superstable, then the memory kernel μ is superexponential. Paradigmatic examples of superexponential kernels are μ(s) = e −s 2 and any compactly supported kernel (complying with the general assumptions of the preceding section).
The main question now is: Are there superexponential μ which render the corresponding S(t) superstable? The answer to this question is negative, and this is perhaps the main result of the very recent paper [4]. From a physical viewpoint, this fact somehow tells something about the persistence of the elastic force versus viscous effects. Notwithstanding, it might as well be possible that there do exist individual nontrivial trajectories that go to zero in a superexponential fashion. In light of (1.12), the trajectories which are more likely to decay fast are the ones corresponding to G ≡ 0, i.e., the solutions to the Volterra equation (1.4). Our claim is that, for any nonzero initial data, the corresponding E cannot decay arbitrarily fast. More precisely, we have the following conjecture.

Conjecture 2.6
For any given memory kernel μ, there exists a structural constant κ > 0 with the following property: for any (u 0 , v 0 ) = (0, 0), the corresponding E fulfills the relation Unfortunately, an analytic proof of such a conjecture seems to be out of reach. Indeed, although we have refined techniques to study the global behavior of the semigroup, we do not have at disposal significant tools that enable us to discuss in detail the behavior of a single trajectory. This is the point of the story where the Numerics steps in.

The equation
We assume for simplicity that the domain D(A) of A is compactly embedded into H . This assumption, which is satisfied in the concrete case of the Laplace-Dirichlet operator occurring in viscoelasticity, guarantees that the spectrum of A is made by eigenvalues only. Besides, denoting the eigenvalues of A (each counted with its own multiplicity) by 0 < λ 1 ≤ λ 2 ≤ · · · ≤ λ j → ∞, the (normalized) eigenvectors w j of λ j form a complete orthonormal basis of H (see, e.g., [38]). This allows us to decompose the solution of (1.4) into the infinite sum where u j (t) is the solution to the ordinary differential equation obtained by projecting (1.4) along with the initial data (1.5) onto the eigenspace relative to the eigenvalue λ j , that is (up to a straightforward change of variables) with initial values By the classical Plancherel Theorem, the energy E(t) of (1.4) splits into the infinite sum is the energy corresponding to the jth-mode. Accordingly, This means that the decay rate of the energy of a particular trajectory is certainly smaller than or equal to the decay rate of the fastest (nonzero) mode. Reason why we can shift our focus to the study of the equation (3.1), ruling the evolution of the single mode associated to λ j . Theorem 2.1 implies that, if the memory kernel complies with assumption (2.3), then E j (t) decays at least exponentially for every j ∈ N. The validity of our Conjecture 2.6 translates into the fact that every E j (t), for different choices of the initial data, should not decay at an exponential rate higher than a certain threshold κ > 0. To this end, we will solve numerically equation (3.1) for different values of λ j over a time lapse [0, T ], with T sufficiently large for the energy to stabilize on a specific decay pattern. Clearly, we are interested in choosing a kernel μ(s) which is superexponential, in the sense of Definition 2.5. Here a problem arises, namely, the fact that after a certain period of time (possibly much smaller than T ) a superexponential kernel becomes so small that some approximation issues may occur in the numerical scheme. To overcome this difficulty, we introduce the truncated kernel where ε > 0 will be chosen suitably small. Then, we consider the problem 3) The following result, borrowed from [9], ensures that the solutions u ε j of converge to u j as ε → 0.

Proposition 3.1 Let u j and u ε
j be solutions to respectively (3.1) and (3.3) with the same initial data u 0, j , v 0, j . Then there exists a constant C > 0 such that where To summarize, having fixed a superexponential kernel μ(s), for suitably chosen T and ε we will perform a numerical simulation of the ordinary integrodifferential equation (where we write just u in place of u ε j and λ in place of λ j ) and study the behavior of the related energy in dependence of λ and of the initial data u 0 , v 0 .

The numerical scheme
In this section we briefly describe the approach to numerically approximate a secondorder integrodifferential equation of the form for given λ, and for integration limits a(t) and b(t) such that 0 ≤ a(t) ≤ b(t) ≤ T . Notice that (3.4) can be written in the form (3.5) by choosing Over the years many numerical techniques have been developed to tackle problems of this kind. We recall, among the others, perturbation methods, as well as spectral and Chebyshev collocation methods (see, e.g., [1,11,40] and references therein), but the list is far from being exhaustive. Here we employ a numerical approach that suitably combines an ODE solver with an iterative scheme, as proposed in [24]. For a given integer N h , let t n , with n = 0, . . . , N h , be a set of time snapshots within the interval [0, T ], with the convention that t 0 = 0. We define y i to be the approximate solution of u(t) at t = t i , i.e., To set up the algorithm, we consider an initial guess y (0) n . This is obtained by solving numerically (3.5) and neglecting the integral term on the right hand-side, by using for example, the variable-step, variable-order Adams-Bashforth-Moulton (ABM) predictor-corrector scheme (see, e.g., [36,39]). Next, we compute y (k+1) n ← y (k) n as described in Algorithm 1.
Notice that in the first step of Algorithm 1, the integral appearing on the right hand side is solved numerically, employing a Lobatto quadrature algorithm. Furthermore, α is a suitable smoothing parameter, which is is equivalent to a relaxation parameter in the context of iterative methods for linear systems. The iterative algorithm is terminated whenever two consecutive iterates differ up to a prescribed (user-defined) tolerance TOL, i.e., or whenever a maximum number of iterations N it is reached.

Numerical verification
For a fixed value of the tolerance TOL = 10 −8 , we are interested in validating the numerical scheme, showing that the approximate solution y n converges to the exact solution u(t n ) as N h → ∞ (or, equivalently, as h → 0). To this end, we consider the following Cauchy problem where t ∈ [0, 1]. This is the test case considered in [24]. For such an equation, the exact solution is known and is equal to

Fig. 1 Computed err h at the discretization points for different choices of N h
We compute the numerical solution y n of (3.6) using an increasingly fine equispaced approximation of the time span [0, 1] made by N h = 50·2 j intervals, for j = 6, 8, 10.
Then, for each of the N h , we evaluate the pointwise error at the discretization nodes In Fig. 1 we report err h at the discretization points for our choices of N h . As we can see, the numerical performance of the scheme improves as N h increases.

Numerical results
We investigate numerically (3.5) for three different superexponential kernels, namely, where χ denotes the characteristic function. Without claiming to be complete, we believe that these functions provide a somewhat exhaustive overview of the most important features which a superexponential kernel can enjoy. Indeed, μ 1 (s) is the prototype of superexponential kernel, μ 2 (s) is a function which decays extremely fast (faster than any function of the form e −s p , with p > 0), and μ 3 (s) is superexponential simply because it is compactly supported.

Remark 4.1
We also tested other types of memory kernels, obtaining comparable results. Among others, we mention piecewise continuous kernels exhibiting a finite number of (downward) jumps, as well as kernels possessing a weak singularity at zero, such as with p ∈ (0, 1) and k p sufficiently small in order to guarantee that the total mass is less than 1, in compliance with (1.3).
For all the choices of superexponential kernels μ i (s), i = 1, 2, 3 we consider two test cases.

Test Case 1
In the first set of numerical experiments, we fix the initial data (u 0 , v 0 ) and compute numerically the energy E j for the choice λ j = j, with j = 1, . . . , 20. Then, we plot log(E j ) against the time t. If Conjecture 2.6 were true, we would expect to observe the following facts: 1. The single energy profile associated to a certain λ j is not log-concave. Accordingly, since by Theorem 2.1 we know that the energy decays (at least) exponentially, its profile should eventually become a straight line in the logarithmic plot. 2. The exponential decay rates, which in this context correspond to the slopes of the lines in the logarithmic plot, are bounded from above. This is witnessed by the fact that either such rates reach a maximum (recall that the rate is positive by definition), or they stabilize on a certain value as λ j increases.

Test Case 2
In the second set of numerical experiments, we investigate what happens by changing the initial data. For every choice of (u 0 , v 0 ) we compute the exponential decay rate of E j for increasing choices of λ j , and seek for differences in the behavior depending on the initial data. Once again, if Conjecture 2.6 were true, barring minor distinctions, we should observe that for every choice of the initial data the decay rates either reach a maximum or stabilize for large λ j .
All the numerical computations have been performed in MATLAB™, version R2021a. Throughout the section the parameter α appearing in Algorithm 1 is chosen to be α = 1/2.

Remark 4.2
In our Test Cases, we assumed the first eigenvalue λ 1 of A to be equal to 1. No difference occurs in the results if one starts from a different (but of course strictly positive) λ 1 .
Before getting into the details of the results, we highlight a crucial issue. We would like to perform a sufficiently accurate simulation of (3.5), up to a final time T large enough to observe how the energy E j behaves asymptotically. However, since we are dealing with exponentially and superexponentially decaying objects, in our computations we need to pay particular attention not to reach machine epsilon precision, as well as to avoid possible round-off errors. In particular, whenever T becomes too large both the energy and the superexponential kernel might become too small, making the simulation inaccurate. This is precisely the reason why we introduced the truncated kernel μ ε (s). Now let ε > 0 be small enough that where L is a user-defined tolerance. In practice we choose L equal to 10 −16 . Now, let μ ε (s) be defined as in (3.2). Then, recalling Proposition 3.1, we have In the case of μ 1 (s), we have where we have used (1.3) in the last inequality. Besides, since e s > s 2 for every s ∈ R + , it is clear that the latter inequality holds also for μ 2 (s). In light of Proposition 3.1, this translates into an error of the order of √ L. As we will see, this numerical error is of a lower order than the computed energy, and therefore it does not affect the final result of the numerical simulations. We do not need to define an ε for μ 3 (s), as the function is already compactly supported (and the support is reasonably small).

The kernel 1
For the first kernel, we choose ε = 1/6 and T = 12. In particular, we truncate the function at t = 6. In the first numerical test we work with fixed initial data and we compute the energy E j (t) for λ j = 1, . . . , 20. The tolerance in Algorithm 1 is set at TOL = 10 −8 and N h = 2000. In Fig. 2 we report the computed E j as a function of the time (semi-logarithmic scale). For visualization purposes, we only selected the profiles corresponding to the eigenvalues λ j = 2, 5, 9, 14, 20. We refer to Table 1 for all the (least-square) computed exponential decay rates. We observe that the energy decays exponentially fast for every λ j , but not superexponentially, as highlighted by the fact that none of the decay profiles is log-concave. Furthermore, we see that, as the λ j increases, the rate settles on a value of about 1.05, whereas the fastest decay, about 1.1, is obtained in correspondence of λ 11 . Also observe that, in accordance with our theoretical expectations, the energy computed numerically is indeed decreasing.
In the second numerical investigation (Test Case 2), we perform the same test as before, except that now we try different initial data. The only difference in the numerical parameters is the number of intervals, which is now set at N h = 1600. In  Computed exponential decay rates for initial data defined as in (4.1) view of the linearity of (3.5), we restrict ourselves to initial data of Euclidean norm equal to 1. Specifically, we choose where For every initial datum, we verify that the energy profile is not log-concave (hence it is log-linear). Then, we compute the decay rate (i.e., the slope) of the energy E j for each λ j . In Fig. 3 we can see the plot of the results. To wit, on the x-axis we have the different initial data, while on the y-axis we have λ j . For small λ j , the decay rate is also small, for every choice of the initial data. Then, as we increase λ j , the rate starts to increase and stabilizes on a certain profile. This limit profile attains its minimum for the choice (u 0 , v 0 ) = (0, 1).

The kernel 2
Here we take ε = 5/18 and T = 12. Again, we choose TOL = 10 −8 , and N h = 2000 for Test Case 1, whereas N h = 1600 for Test Case 2. We perform the same set of numerical experiments as before, obtaining similar results. Indeed, in Fig. 4, we see that the energy relative to the initial data in (4.1) has an exponential decay. Moreover, from the results reported in Table 2 there is the evidence that the decay rate settles around 0.38 for large λ j , while Fig. 5 shows that the situation does not change for different initial data. In fact, for this particular kernel, the experimental results exhibit a less accentuated dependence on the initial data.

The kernel 3
The last set of numerical experiments features a compactly supported kernel. As already mentioned, in this case there is no need use Proposition 3.1 and, therefore, to define . For these experiments we set T = 6, which is already enough for the energy to stabilize. The other numerical parameters are exactly the same as before.
Concerning Test Case 1, we see that every single energy profile is, again, log-linear. The main difference with respect to the preceding kernels is that now the exponential decay rate does not stabilize for large λ j . Indeed, it reaches its maximum at λ = 13, and then starts to decrease for larger values (see Fig. 6 and Table 3). Such a behavior is also observed in Test Case 2. Here, with very minor differences depending on the initial data, the decay rate decreases as λ j increases. To confirm  Computed decay rate for variable initial data defined as in (4.2) this pattern, we actually performed other experiments with the values of λ j up to 40, without noting significant differences (see Fig. 7).

Completely flat kernels
We conclude this work by exploiting our numerical approach to examine an issue that, although not pertinent to the verification of Conjecture 2.6, has an independent interest in the context of the theory of viscoelasticity. As shown in [2], for a completely flat kernel μ (i.e., with μ = 0 almost everywhere) the contraction semigroup S(t) generated by (1.10) is not exponentially stable. This is the case, for instance, if we take the step kernel μ(s) = χ [0,1/2] (s), We aim to demonstrate, by employing our arguments, that the same feature reflects on the corresponding energy E(t) of the Volterra equation. To this end, we perform a final numerical simulation by taking μ as in (5.1) with ρ = 1, ρ = 1 8 , ρ = 1 32 , ρ = 0.
We expect the energy to exponentially decay for ρ = 0, but not for ρ = 0. We simulate equation (3.5), with initial data (4.1), for the chosen kernels for λ = 2, N h = 2000 and T = 12. The plot of the logarithm of the computed energies can be seen in Fig. 8. Notably, the profile is a decreasing line whose angular coefficient tends to zero as ρ → 0.

Conclusions
For several choices of initial data, we have shown that the (numerical) energy of the solution is never log-concave. Furthermore, its decay rate remains bounded as λ j increases. Such numerical results provide a strong evidence in support of Conjecture 2.6. Another interesting question concerns with the role played by each of the components E 0 (t) and E 1 (t) of the energy, defined in (1.6)-(1.7), in the decay rate. One might argue that the lack of superstable trajectories suggested by our numerical simulations is, in a sense, artificial, and only due to the presence of the memory energy E 1 (t), while the "effective" mechanical energy E 0 (t) might still exhibit a superexponential decay.
In fact, we believe that a stronger statement than Conjecture 2.6 holds, namely, that E 0 (t) itself cannot decay faster than an exponential. In this case, the numerical analysis is slightly more complicated, since the mechanical energy is not a decreasing function of time any longer, which introduces some technical difficulties in the interpretation of the results.
In Fig. 9 we can see the logarithmic plot of E 0, j . Although the behavior of the (logarithm of the) mechanical energy is not as regular as the one of E j (t), there is a clear stabilization of the profile for every λ j . Although a single example cannot be taken as an ultimate evidence, this suggests that Conjecture 2.6 holds in the stronger form.