Determining damping terms in fractional wave equations

We consider the second order in time damped wave equation

$$\begin{equation}{u}_{tt}+{c}^{2}\mathcal{A}u+\sum\limits _{k=1}^{N}{b}_{k}{\mathrm{\partial }}_{t}^{{\alpha }_{k}}{\mathcal{A}}^{{\beta }_{k}}u=r,\end{equation} \tag{ 3 }$$

where $\mathcal{A}=-\mathbb{L}$ on ${\Omega}\subseteq {\mathbb{R}}^{d}$, equipped with homogeneous Dirichlet/Neumann/impedance boundary conditions, for an elliptic differential operator $\mathbb{L}$ and ${\alpha }_{k}\in \left(0,1\right]$, ${\beta }_{k}\in \left(\frac{1}{2},1\right]$.

Throughout this paper, we denote by ${\partial }_{t}^{\alpha }$ the (partial) Caputo–Djrbashian fractional time derivative of order α ∈ (n − 1, n) with $n\in \mathbb{N}$ by

$$\begin{equation*}{\partial }_{t}^{\alpha }u=\frac{1}{{\Gamma}(n-\alpha )}{\int }_{0}^{t}\frac{{\partial }_{t}^{n}u(s)}{{(t-s)}^{\alpha +1-n}}\,\mathrm{d}s,\end{equation*}$$

where ${\partial }_{t}^{n}$ denotes the nth integer order partial time derivative and for γ ∈ (0, 1), and ${I}_{t}^{\gamma }$ is the Abel integral operator defined by ${I}_{t}^{\gamma }[v](t)=\frac{1}{{\Gamma}(\gamma )}{\int }_{0}^{t}\frac{v(s)}{{(t-s)}^{1-\gamma }}\,\mathrm{d}s.$ For details on fractional differentiation and subdiffusion equations, we refer to, e.g., [8, 9, 27, 31, 32]. See also the tutorial paper on inverse problems for anomalous diffusion processes [19].

We assume that the differentiation orders with respect to time are distinct, that is

$$\begin{equation}0< {\alpha }_{1}< {\alpha }_{2}< \cdots < {\alpha }_{N}\leqslant 1,\qquad 0< {\gamma }_{1}< {\gamma }_{2}< \cdots < {\gamma }_{J}\leqslant 1\end{equation} \tag{ 4 }$$

a property that is crucial for distinguishing the different asymptotic terms in the solution u and its Laplace transform.

More generally than (3), in most of this paper we also take higher than second order time derivatives into account

$$\begin{equation}{u}_{tt}+{c}^{2}\mathcal{A}u+\sum\limits _{j=1}^{J}{d}_{j}{\partial }_{t}^{2+{\gamma }_{j}}u+\sum\limits _{k=1}^{N}{b}_{k}{\partial }_{t}^{{\alpha }_{k}}{\mathcal{A}}^{{\beta }_{k}}u=r.\end{equation} \tag{ 5 }$$

Note that the usual models from [21] are contained in (5), which is actually an extension of the general model from, e.g., the book [1] in the sense that it also allows for fractional powers of the operator $\mathcal{A}$.

Typically, in acoustics we just have $\mathbb{L}={\triangle}$, and then the operator $\mathcal{A}$ is known. To take into account a (possibly unknown) spatially varying speed of sound c(x), one can instead consider

$$\begin{equation}{u}_{tt}+{\mathcal{A}}_{c}u+\sum\limits _{j=1}^{J}{d}_{j}{\partial }_{t}^{2+{\gamma }_{j}}u+\sum\limits _{k=1}^{N}{b}_{k}{\partial }_{t}^{{\alpha }_{k}}{\mathcal{A}}_{c}^{{\beta }_{k}}u=r\end{equation} \tag{ 6 }$$

with ${\mathcal{A}}_{c}=-c{(x)}^{2}{\triangle}$; in order to get a selfadjoint operator ${\mathcal{A}}_{c}$ we then use the weighted L² inner product with weight function $\frac{1}{{c}^{2}}$.

Note that in order to assume (4) without loss of generality, so allowing for all possible combinations of α_k , β_k , we would have to consider

$$\begin{equation}{u}_{tt}+{c}^{2}\mathcal{A}u+\sum\limits _{k=1}^{N}{\partial }_{t}^{{\alpha }_{k}}\sum\limits _{\ell =1}^{{M}_{k}}{c}_{k\ell }{\mathcal{A}}^{{\beta }_{k\ell }}u=r\end{equation} \tag{ 7 }$$

or, including higher than second time derivatives such as in (5)

$$\begin{equation}{u}_{tt}+{c}^{2}\mathcal{A}u+\sum\limits _{j=1}^{J}{d}_{j}{\partial }_{t}^{2+{\gamma }_{j}}u+\sum\limits _{k=1}^{N}{\partial }_{t}^{{\alpha }_{k}}\sum\limits _{\ell =1}^{{M}_{k}}{c}_{k\ell }{\mathcal{A}}^{{\beta }_{k\ell }}u=r.\end{equation} \tag{ 8 }$$

However, (7) and (8) might be viewed as over-parameterised models and they actually lead to difficulties in proving uniqueness of solutions.

Important special cases of (6) are the Caputo–Wismer–Kelvin–Chen–Holm

$$\begin{equation}{u}_{tt}+{\mathcal{A}}_{c}u+b{\partial }_{t}^{\alpha }{\mathcal{A}}_{c}^{\beta }u=r\end{equation} \tag{ 9 }$$

and the fractional Zener fz model

$$\begin{equation}{u}_{tt}+{\mathcal{A}}_{c}u+d{\partial }_{t}^{2+\gamma }u+b{\partial }_{t}^{\alpha }{\mathcal{A}}_{c}u=r.\end{equation} \tag{ 10 }$$

These pde models will be considered on a time interval t ∈ (0, T) and driven by initial conditions and/or a separable source term

$$\begin{equation}\begin{aligned}\hfill & r(x,t)=\sigma (t)f(x),\quad x\in {\Omega},\quad t\in (0,T)\hfill \\ \hfill & u(x,0)={u}_{0}(x),\quad {u}_{t}(x,0)={u}_{1}(x),\quad ({u}_{tt}(x,0)={u}_{2}(x)\quad \text{if}\ {\gamma }_{J} > 0)\quad x\in {\Omega}.\hfill \end{aligned}\end{equation} \tag{ 11 }$$

Note that boundary conditions are already incorporated into the operator $\mathcal{A}$ or ${\mathcal{A}}_{c}$, respectively.

As in [18], we assume that data has been obtained from a quite general experimental setup in which the pde (5) with its unknown coefficients remains the same, while both excitation (11) and observation location may vary between the experiments. We number the experiments with an index i to denote the corresponding excitation by (u_0,i , u_1,i , f_i ), and the corresponding observation operator by B_i . Given I experiments, we thus assume to have the following I time trace observations for different driving by initial and/or source data:

$$\begin{equation}\begin{aligned}\hfill & {h}_{i}(t)=({B}_{i}{u}_{i})(t),\ t\in (0,T)\quad \text{where}\ {u}_{i}\,\text{solves}\;(5)\enspace \mathrm{a}\mathrm{n}\mathrm{d}\enspace (11)\hfill \\ \hfill & \qquad \text{with}\;({u}_{0},{u}_{1},f)=({u}_{0,i},{u}_{1,i},{f}_{i}),\enspace i=1,\dots ,I.\hfill \end{aligned}\end{equation} \tag{ 12 }$$

Examples of observation operators are

$$\begin{equation*}(\mathrm{a})\ {B}_{i}v=v({x}_{i})\quad \;\text{or}\;\quad (\mathrm{b})\ {B}_{i}v={\int }_{{{\Sigma}}_{i}}{\eta }_{i}(x)v(x)\mathrm{d}x\end{equation*}$$

for some points ${x}_{i}\in \bar{{\Omega}}$ and some weight functions η_i ∈ L^∞(ω), Σ_i ⊆ ∂Ω, provided that these evaluations are well-defined, that is, $u(t)\in C(\bar{{\Omega}})$ in case (a), or u(t) ∈ L¹(ω) in case (b), which according to theorem 2.1 is the case if ${\Omega}\subseteq {\mathbb{R}}^{d}$ with d = 1 for case (a) or d ∈ {1, 2, 3} for case (b).

The coefficients in (5) that we aim to in recover are

$$\begin{equation*}N,J,c,{b}_{k},{d}_{j},{\alpha }_{k},{\beta }_{k},{\gamma }_{j},\quad j=1,\dots ,J,\ k=1,\dots ,N.\end{equation*}$$

Moreover, we consider the problem of identifying, besides the differentiation orders, some space dependent quantities. In practice the most relevant ones are either the speed of sound c = c(x) in ultrasound tomography, or the initial data u₀ = u₀(x) while u₁ = 0 (equivalently the source term f = f(x), while u₀ = 0, u₁ = 0) in photoacoustic tomography pat. For the latter purpose, we will consider observations not only at single points (a) or patches (b), but over a surface Σ.

In this section we prove well-posedness of the forward problem (6). Just as in the inverse problem which is the main topic of this paper, we focus on the case of constant coefficients b_k , d_j , while the sound speed c = c(x) contained in the operator $\tilde{\mathcal{A}}$ may vary in space. We refer to chapter 7 of the author' upcoming book [22] for an analysis in the case when all coefficients are space dependent.

Our analysis will be based on energy estimates obtained by multiplying the pde (or actually its Galerkin semidiscretization with respect to space) with appropriate multipliers. For this purpose, we will make use of certain mapping properties of fractional derivative operators that can be quantified via upper and lower bounds.

We will use the following results that hold for $\theta \in \left[0,1\right)$.

• From [10, lemma 2.3]; see also [35, theorem 1]: for any w ∈ H^−θ/2(0, t),
  $$\begin{equation}{\int }_{0}^{t}{{\langle }_{0}{I}_{t}}^{\theta }w(s),w(s)\rangle \mathrm{d}s\geqslant \mathrm{cos}(\pi \theta /2){\Vert}w{{\Vert}}_{{H}^{-\theta /2}(0,t)}^{2}.\end{equation} \tag{ 13 }$$
• From [23, theorems 2.1, 2.2, 2.4 and proposition 2.1]: there exist constants $0< \underline{C}(\theta )< \bar{C}(\theta )$ such that for any w ∈ H^θ (0, t) (with w(0) = 0 if $\theta \geqslant \frac{1}{2}$), the equivalence estimates
  $$\begin{equation}\underline{C}(\theta ){\Vert}w{{\Vert}}_{{H}^{\theta }(0,t)}\leqslant {\Vert}{\partial }_{t}^{\theta }w{{\Vert}}_{{L}^{2}(0,t)}\leqslant \bar{C}(\theta ){\Vert}w{{\Vert}}_{{H}^{\theta }(0,t)}\end{equation} \tag{ 14 }$$
  hold. Likewise, for any w ∈ L²(0, t) we have ${\partial }_{t}^{\theta }w\in {H}^{-\theta }(0,t)$ and
  $$\begin{equation}\underline{C}(\theta ){\Vert}w{{\Vert}}_{{L}^{2}(0,T)}\leqslant {\Vert}{\partial }_{t}^{\theta }w{{\Vert}}_{{H}^{-\theta }(0,t)}\leqslant \bar{C}(\theta ){\Vert}w{{\Vert}}_{{L}^{2}(0,t)}.\end{equation} \tag{ 15 }$$
  Moreover, (by a duality argument) there exists a constant C(2θ) > 0 such that for any w ∈ L²(0, t)
  $$\begin{equation}{\int }_{0}^{t}{{\left(\right.}_{0}{I}_{t}}^{2\theta }w(s)\left.\right)w(s)\mathrm{d}s\leqslant C(2\theta ){\Vert}w{{\Vert}}_{{H}^{-\theta }(0,t)}^{2}.\end{equation} \tag{ 16 }$$
  Note that the latter two facts (15) and (16) hold true for θ = 1 as well.

We are now prepared to state and prove our main theorem on well-posedness of an initial value problem

$$\begin{equation}\left\{\begin{aligned}\hfill \sum\limits _{j=0}^{J}{d}_{j}& {\partial }_{t}^{2+{\gamma }_{j}}u+\sum\limits _{k=0}^{K}{b}_{k}{\partial }_{t}^{{\alpha }_{k}}\tilde{\mathcal{A}}u=r\quad \;\text{in}\;{\Omega}\times (0,T),\hfill \\ \hfill & u(x,0)={u}_{0}(x),\ {u}_{t}(x,0)={u}_{1}(x)\quad x\in {\Omega},\hfill \\ \hfill & (\text{if}\;{\gamma }_{J} > 0:\ {u}_{tt}(x,0)={u}_{2}(x)\quad x\in {\Omega}),\hfill \end{aligned}\right.\end{equation} \tag{ 17 }$$

which clearly covers the models from section 1.1 in case of β₁ = ⋯ = β_K . Here $\tilde{\mathcal{A}}$ is a selfadjoint (with respect to some possibly weighted space ${L}_{w}^{2}({\Omega})$ with inner product ⟨⋅, ⋅⟩) operator with eigenvalues and eigenfunctions ${({\lambda }_{i},{\varphi }_{i})}_{i\in \mathbb{N}}$, and ${\dot{H}}^{1}({\Omega})=\left\{v\in {L}_{w}^{2}({\Omega}):{\Vert}{\tilde{\mathcal{A}}}^{1/2}v{{\Vert}}_{{L}_{w}^{2}({\Omega})}< \infty \right\}$. This includes the cases $\tilde{\mathcal{A}}=-{\triangle}$, $\tilde{\mathcal{A}}=-c{(x)}^{2}{\triangle}$, $\tilde{\mathcal{A}}=-c{(x)}^{2}\nabla \cdot (\frac{1}{\rho (x)}\nabla )$ relevant for acoustics, where in the latter two cases we use the inner product $\langle {v}_{1},{v}_{2}\rangle ={\int }_{{\Omega}}\frac{1}{{c}^{2}(x)}{v}_{1}(x){v}_{2}(x)\mathrm{d}x$ in which $\;\tilde{\phantom{\rule{0ex}{0ex}}\mathcal{A}}$ is selfadjoint.

The coefficients and orders will be assumed to satisfy

$$\begin{equation}0\leqslant {\gamma }_{0}< {\gamma }_{1}< \cdots < {\gamma }_{J}\leqslant 1,\qquad 0\leqslant {\alpha }_{0}< {\alpha }_{1}< \cdots < {\alpha }_{K}\leqslant 1,\quad {\gamma }_{J}\leqslant {\alpha }_{K},\end{equation} \tag{ 18 }$$

and

$$\begin{equation}{d}_{j}\geqslant 0,\quad j\in \left\{1,\dots ,J\right\},\qquad {b}_{k}\geqslant 0,\quad k\in \left\{{k}_{\ast },\dots ,K\right\},\quad {d}_{J} > 0,\quad {b}_{K} > 0\end{equation} \tag{ 19 }$$

with k_* to be specified below. The condition γ_J ⩽ α_K has been shown to be necessary for thermodynamic consistency, see [1, lemma 3.1].

Note that the proof of theorem 2.1 also goes through without the restrictions γ_j ⩽ 1, α_k ⩽ 1, but orders larger than one would require introducing further initial conditions, which would make the exposition less transparent. Moreover, the focus in this paper is actually on wave type equations where this restriction of the orders is physically relevant.

We first of all consider the following time integrated and therefore weaker form with $\tilde{r}(x,t)={\int }_{0}^{t}r(x,s)\mathrm{d}s$

$$\begin{equation}\left\{\begin{aligned}\hfill & \sum\limits _{j=0}^{J}{d}_{j}\left({\partial }_{t}^{1+{\gamma }_{j}}u-\frac{{t}^{1-{\gamma }_{j}}}{{\Gamma}(2-{\gamma }_{j})}{u}_{2}(x)\right)\hfill \\ \hfill & +\sum\limits _{k=0}^{K}{b}_{k}{{\left(\right.}_{0}{I}_{t}}^{1-{\alpha }_{k}}\tilde{\mathcal{A}}u-\frac{{t}^{1-{\alpha }_{k}}}{{\Gamma}(2-{\alpha }_{k})}\tilde{\mathcal{A}}{u}_{0}\left.\right)=\tilde{r}\quad \;\text{in}\;{\Omega}\times (0,T),\hfill \\ \hfill & \qquad u(x,0)={u}_{0}(x),\ {u}_{t}(x,0)={u}_{1}(x)\quad x\in {\Omega}.\hfill \end{aligned}\right.\end{equation} \tag{ 20 }$$

The results are obtained by testing the pde in (17) with ${\partial }_{t}^{\tau }u$ for some τ > 0 that needs to be well chosen in order to exploit the estimates (13)–(16). Generalizing the common choice τ = 1 in case of the classical second order wave equation γ₀ = ⋯ = γ_J = α₀ = ⋯ = α_K = 0 and in view of the stability condition γ_J ⩽ α_K we choose τ ∈ [1 + γ_J , 1 + α_K ], actually, $\tau =1+{\alpha }_{{k}_{\ast }}$ with k_* according to

$$\begin{equation}{k}_{\ast }\in \left\{1,\dots ,K\right\}\quad \text{such}\;\text{that}\;{\alpha }_{{k}_{\ast }-1}< {\gamma }_{J}\leqslant {\alpha }_{{k}_{\ast }},\end{equation} \tag{ 21 }$$

if α₀ < γ_J and k_* = 0 otherwise. More precisely, we will assume that

$$\begin{equation}(\mathrm{i})\ {\alpha }_{0}< {\gamma }_{J},\ {k}_{\ast }\;\text{as}\;\text{in}\;(\text{21}),\ \;\text{or}\;\ (\mathrm{i}\mathrm{i})\ {\alpha }_{0}\geqslant {\gamma }_{J},\ {k}_{\ast }{:=}0\end{equation} \tag{ 22 }$$

holds. Additionally, we assume

$$\begin{equation}{\alpha }_{K}-1\leqslant {\alpha }_{{k}_{\ast }}\leqslant 1+\mathrm{min}\left\{{\gamma }_{0},{\alpha }_{0}\right\}\ \;\text{and}\;\text{either}\;\ {\gamma }_{J}< {\alpha }_{K}\;\text{or}\;\ {\gamma }_{J}={\alpha }_{K}={\alpha }_{0}.\end{equation} \tag{ 23 }$$

The upper and lower bounds on ${\alpha }_{{k}_{\ast }}$ in (23) are needed to guarantee $\theta \in \left[0,1\right)$ in several instances of θ used in (13)–(16) in the proof below.

Note that we are considering the case γ_J = α_K only in the specific setting of a single $\tilde{\mathcal{A}}$ term. Indeed, when admitting several $\tilde{\mathcal{A}}$ terms, well-posedness depends on specific coefficient constellations. As an example of this, let us mention the fractional Moore–Gibson–Thompson equations arising in acoustics

$$\begin{equation}{u}_{tt}+\mathfrak{T}{\partial }_{t}^{2+\alpha }u-{\mathfrak{c}}^{2}{\Delta}u-\mathfrak{b}{\partial }_{t}^{\alpha }{\Delta}u=f\end{equation} \tag{ 24 }$$

with positive parameters $\mathfrak{T},\mathfrak{b},\mathfrak{c},\alpha$. It is well-posed only if $\mathfrak{b}\geqslant {\mathfrak{c}}^{2}\mathfrak{T}$, according to [20, proposition 7.1]. As the analysis of this example also shows, even lower order time derivatives of $\tilde{\mathcal{A}}u$ (such as the $-{\mathfrak{c}}^{2}{\Delta}u$ term in (24)) may prohibit global in time boundedness of solutions. It will turn out in the proof that we therefore have to impose a bound on the magnitude of the lower order term coefficients in terms of the highest order one

$$\begin{equation}\begin{aligned}\hfill & \sum\limits _{k=0}^{{k}_{\ast }-1}C(1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k})\bar{C}((1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k})/2)\vert {b}_{k}\vert {{\Vert}}_{0}{I}_{t}^{({\alpha }_{K}-{\alpha }_{k})/2}{{\Vert}}_{{L}^{2}(0,T)\to {L}^{2}(0,T)}\hfill \\ \hfill & \qquad \leqslant \underline{C}((1+{\alpha }_{{k}_{\ast }}-{\alpha }_{K})/2){b}_{K}\,\mathrm{cos}(\pi (1+{\alpha }_{{k}_{\ast }}-{\alpha }_{K})/2)\hfill \end{aligned}\end{equation} \tag{ 25 }$$

that due to T dependence of ${{\Vert}}_{0}{I}_{t}^{({\alpha }_{K}-{\alpha }_{k})/2}{{\Vert}}_{{L}^{2}(0,T)\to {L}^{2}(0,T)}$ clearly restricts the maximal time T unless all b_k with index lower than k_* vanish, cf corollary 2.1 below.

By applying the test functions as mentioned above, we will derive estimates for the energy functionals (which are not supposed to represent physical energies but are just used for mathematical purposes)

$$\begin{align}{\mathcal{E}}_{\gamma }[u](t)& =\sum\limits _{j=0}^{J-1}\mathrm{c}\mathrm{o}\mathrm{s}(\pi (1+{\gamma }_{j}-{\alpha }_{{k}_{\ast }})/2){\Vert} \sqrt{{d}_{j}}{\mathrm{\partial }}_{t}^{(3+{\gamma }_{j}+{\alpha }_{{k}_{\ast }})/2}u{{\Vert} }_{{L}^{2}(0,t;{L}_{w}^{2}({\Omega}))}^{2}+{\underline{\mathcal{E}}}_{\gamma }[u](t),\end{align} \tag{ 26 }$$

where

$$\begin{equation}{\underline{\mathcal{E}}}_{\gamma }[u](t)=\begin{cases}\tilde{d}{\Vert}{\partial }_{t}^{(3+{\gamma }_{J}+{\alpha }_{{k}_{\ast }})/2}u{{\Vert}}_{{L}^{2}(0,t;{L}_{w}^{2}({\Omega}))}^{2}\quad \hfill & \;\text{if}\;{\gamma }_{J}< {\alpha }_{{k}_{\ast }}\hfill \\ \frac{{d}_{J}}{4}{\Vert}{\partial }_{t}^{1+{\gamma }_{J}}u{{\Vert}}_{{L}^{\infty }(0,t;{L}_{w}^{2}({\Omega}))}^{2}\quad \hfill & \;\text{if}\;{\gamma }_{J}={\alpha }_{{k}_{\ast }}\hfill \end{cases}\end{equation} \tag{ 27 }$$

where $\tilde{d}={d}_{J}\,\mathrm{cos}(\pi (1+{\gamma }_{J}-{\alpha }_{{k}_{\ast }})/2)\underline{C}{((1+{\gamma }_{J}-{\alpha }_{{k}_{\ast }})/2)}^{2}$ with $\underline{C}$ as in (15), $\theta =(1+{\gamma }_{J}-{\alpha }_{{k}_{\ast }})/2\in \left[0,1\right)$ and

$$\begin{equation}\begin{aligned}\hfill {\mathcal{E}}_{\alpha }[u](t)& =\frac{1}{4}{\Vert}\sqrt{{b}_{{k}_{\ast }}}{\partial }_{t}^{{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}u{{\Vert}}_{{L}^{\infty }(0,t);{L}_{w}^{2}({\Omega})}^{2}+\sum\limits _{k={k}_{\ast }+1}^{K}\mathrm{cos}(\pi (1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k})/2)\hfill \\ \hfill & \quad \times {\Vert}\sqrt{{b}_{k}}{\partial }_{t}^{(1+{\alpha }_{{k}_{\ast }}+{\alpha }_{k})/2}{\tilde{\mathcal{A}}}^{1/2}u{{\Vert}}_{{L}^{2}(0,t;{L}_{w}^{2}({\Omega}))}^{2},\hfill \end{aligned}\end{equation} \tag{ 28 }$$

where the sum in (28) is void in case k_* = K.

The solution space induced by these energies is U_γ ∩ U_α where

$$\begin{equation}\begin{aligned}\hfill {U}_{\gamma }& =\begin{cases}{H}^{(3+{\gamma }_{J}+{\alpha }_{{k}_{\ast }})/2}(0,T;{L}_{w}^{2}({\Omega}))\quad \hfill & \;\text{if}\;{\gamma }_{J}< {\alpha }_{{k}_{\ast }}\hfill \\ \left\{\right.v\in {L}^{2}(0,T;{L}_{w}^{2}({\Omega})):{\partial }^{1+{\gamma }_{J}}v\in {L}^{\infty }(0,T;{L}_{w}^{2}({\Omega}))\quad \hfill & \;\text{if}\;{\gamma }_{J}={\alpha }_{{k}_{\ast }}\hfill \end{cases}\hfill \\ \hfill {U}_{\alpha }& =\begin{cases}{H}^{(1+{\alpha }_{K}+{\alpha }_{{k}_{\ast }})/2}(0,T;\dot{H}({\Omega}))\quad \hfill & \;\text{if}\;K > {k}_{\ast }\hfill \\ \left\{\right.v\in {H}^{(1-\varepsilon )/2+{\alpha }_{K}}(0,T;\dot{H}({\Omega})):{\partial }^{{\alpha }_{K}}v\in {L}^{\infty }(0,T;\dot{H}({\Omega}))\quad \hfill & \;\text{if}\;K={k}_{\ast }\hfill \end{cases}\hfill \end{aligned}\end{equation} \tag{ 29 }$$

for any ɛ ∈ (0, α_K − max{α_K−1, γ_J }).

Theorem 2.1. Assume that the coefficients γ_j , α_k , d_j , b_k satisfy (18), (19), (22), (23), and that T and/or the coefficients b_k for k ⩽ k_* are small enough so that (25) holds.

Then for any ${u}_{0},{u}_{1}\in {\dot{H}}^{1}({\Omega})$, (${u}_{2}\in {L}_{w}^{2}({\Omega})$ if γ_J > 0), $r\in {H}^{{\alpha }_{{k}_{\ast }}-{\gamma }_{J}}(0,T;{L}_{w}^{2}({\Omega}))$, the time integrated initial boundary value problem (20) (to be understood in a weak ${\dot{H}}^{-1}({\Omega})$ sense with respect to space and an L²(0, T) sense with respect to time) has a unique solution u ∈ U_γ ∩ U_α cf (29). This solution satisfies the energy estimate

$$\begin{align}\hfill {\mathcal{E}}_{\gamma }[u](t)+{\mathcal{E}}_{\alpha }[u](t)& \leqslant C(t)\left({\Vert}{u}_{0}{{\Vert}}_{{\dot{H}}^{1}({\Omega})}^{2}+{C}_{0}+{\Vert}r{{\Vert}}_{{H}^{{\alpha }_{{k}_{\ast }}-{\gamma }_{J}}(0,t,{L}_{w}^{2}({\Omega}))}^{2}\right),\quad t\in (0,T)\hfill \end{align} \tag{ 30 }$$

for some constant C(t) depending on time, which is bounded as t → 0 but in general grows exponentially as t → ∞, ${\mathcal{E}}_{\gamma }$, ${\mathcal{E}}_{\alpha }$ defined as in (26), (28), and

$$\begin{equation}{C}_{0}=\begin{cases}{\Vert}{u}_{1}{{\Vert}}_{{\dot{H}}^{1}({\Omega})}^{2}+{\Vert}{u}_{2}{{\Vert}}_{{L}_{w}^{2}({\Omega})}^{2}\quad \hfill & \;\text{if}\;{\alpha }_{{k}_{\ast }} > 0\hfill \\ {\Vert}{u}_{1}{{\Vert}}_{{L}_{w}^{2}({\Omega})}^{2}\quad \hfill & \;\text{if}\;{\alpha }_{{k}_{\ast }}=0.\hfill \end{cases}\end{equation} \tag{ 31 }$$

Proof. We exclude the case α_K = 0 since this via (18) implies that also γ_J = 0 and we would therefore deal with the conventional second order wave equation whose analysis can be found in textbooks, cf e.g., [11]. Thus for the remainder of this proof we assume α_K > 0 to hold.

Step 1. Galerkin discretisation.

In order to prove existence and uniqueness of solutions to (17), we apply the usual Faedo–Galerkin approach of discretisation in space with eigenfunctions φ_i corresponding to eigenvalues λ_i of $\tilde{\mathcal{A}}$, $u(x,t)\approx {u}^{L}(x,t)={\sum }_{i=1}^{L}\,{u}_{i}^{L}(t){\varphi }_{i}(x)$ and testing with φ_j , that is,

$$\begin{equation}\begin{aligned}\hfill & \begin{cases}\langle \sum\limits _{j=0}^{J}{d}_{j}{\partial }_{t}^{2+{\gamma }_{j}}{u}^{L}(t)+\sum\limits _{k=0}^{K}{b}_{k}{\partial }_{t}^{{\alpha }_{k}}\tilde{\mathcal{A}}{u}^{L}(t)-r(t),{v\rangle }_{{L}_{w}^{2}({\Omega})}=0\quad t\in (0,T)\quad \hfill \\ \langle {u}^{L}(0)-{u}_{0},v\rangle =\langle {u}_{t}^{L}(0)-{u}_{1},v\rangle =0,\quad (\text{if}\;{\gamma }_{J} > 0:\langle {u}_{tt}-{u}_{2},v\rangle =0)\quad \hfill \end{cases}\hfill \\ \hfill & \;\text{for}\;\text{all}\;v\in \text{span}({\varphi }_{1},\dots ,{\varphi }_{L}).\hfill \end{aligned}\end{equation} \tag{ 32 }$$

To prove existence of a solution ${\underline{u}}^{L}\in {H}^{2+{\gamma }_{J}}(0,T;{\mathbb{R}}^{L})$ to (32), we rewrite it as a system of Volterra integral equations for

$$\begin{equation*}{\underline{\xi }}^{L}{:=}{\partial }_{t}^{2+{\gamma }_{J}}{\underline{u}}^{L}=\begin{cases}{}_{0}I_{t}^{1-{\gamma }_{J}}{\underline{u}}_{ttt}^{L}\quad \hfill & \;\text{if}\;{\gamma }_{J} > 0\hfill \\ {\underline{u}}_{tt}^{L}\quad \hfill & \;\text{if}\;{\gamma }_{J}=0.\hfill \end{cases}\end{equation*}$$

Unique solvability of the resulting system in L²(0, T^L ) follows from [16, theorem 4.2, p 241 in section 9]. Then from

$$\begin{equation*}\left\{{\partial }_{t}^{{\gamma }_{J}}{\underline{u}}_{tt}^{L}=\underline{\xi }\in {L}^{2}(0,{T}^{L}),\quad {\underline{u}}_{tt}^{L}(0)={\underline{u}}_{2}^{L}\right.\end{equation*}$$

in case γ_J > 0 we have a unique ${\underline{u}}_{tt}^{L}\in {H}^{{\gamma }_{J}}(0,{T}^{L})$ (cf [23, section 3.3]); the same trivially holds true in case γ_J = 0.

From the energy estimates below, this solution actually exists for all times $t\in \left[0,T\right)$, so the maximal time horizon of existence is actually T^L = T, independently of the discretisation level L.

Step 2. Energy estimates.

With k_* defined by (21) in case (i) and set to k_* = 0 in case (ii), we use $v={\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{u}^{L}(s)$ as a test function in (32) (with t replaced by s) and integrate for s from 0 to t. The resulting terms can be estimated as follows.

Inequality (13) with $\theta =1+{\gamma }_{j}-{\alpha }_{{k}_{\ast }}\in \left[0,1\right)$ yields

$$\begin{equation*}\begin{aligned}\hfill & {\int }_{0}^{t}\langle {d}_{j}{\partial }_{t}^{2+{\gamma }_{j}}{u}^{L}(s),{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{u}^{L}(s)\rangle \mathrm{d}s={\int }_{0}^{t}\langle {d}_{j}{\partial }_{t}^{2+{\gamma }_{j}}{u}^{L}{(s),}_{0}{{I}_{t}}^{1+{\gamma }_{j}-{\alpha }_{{k}_{\ast }}}{\partial }_{t}^{2+{\gamma }_{j}}{u}^{L}(s)\rangle \mathrm{d}s\hfill \\ \hfill & \qquad \geqslant \mathrm{cos}(\pi (1+{\gamma }_{j}-{\alpha }_{{k}_{\ast }})/2){\Vert}\sqrt{{d}_{j}}{\partial }_{t}^{2+{\gamma }_{j}}{u}^{L}{{\Vert}}_{{H}^{-(1+{\gamma }_{j}-{\alpha }_{{k}_{\ast }})/2}(0,t;{L}_{w}^{2}({\Omega}))}^{2}\quad \text{for}\,{\gamma }_{j}< {\alpha }_{{k}_{\ast }}\hfill \end{aligned}\end{equation*}$$

and with $\theta =1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k}\in \left[0,1\right)$

$$\begin{equation*}\begin{aligned}\hfill & {\int }_{0}^{t}\langle {b}_{k}{\partial }_{t}^{{\alpha }_{k}}\tilde{\mathcal{A}}{u}^{L}(s),{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{u}^{L}(s)\rangle \mathrm{d}s={\int }_{0}^{t}\langle {b}_{k}{\partial }_{t}^{{\alpha }_{k}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}(s),{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}(s)\rangle \mathrm{d}s\hfill \\ \hfill & \qquad ={\int }_{0}^{t}\langle {b}_{k0}{{I}_{t}}^{1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k}}{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}(s),{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}(s)\rangle \mathrm{d}s\hfill \\ \hfill & \qquad \geqslant \mathrm{cos}(\pi (1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k})/2){\Vert}\sqrt{{b}_{k}}{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}{{\Vert}}_{{H}^{-(1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k})/2}(0,t;{L}_{w}^{2}({\Omega}))}^{2}\quad \text{for}\,k > {k}_{\ast }.\hfill \end{aligned}\end{equation*}$$

In particular for k = K > k_* (which implies $1+{\alpha }_{{k}_{\ast }}-{\alpha }_{K}< 1$; the case k_* = K will be considered below), by (15) with $\theta =(1+{\alpha }_{{k}_{\ast }}-{\alpha }_{K})/2\in \left[0,1\right)$, noting that ${\partial }_{t}^{(1+{\alpha }_{{k}_{\ast }}+{\alpha }_{K})/2}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}$ vanishes at t = 0 and that ${\partial }_{t}^{(1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k})/2}{\partial }_{t}^{(1+{\alpha }_{{k}_{\ast }}+{\alpha }_{K})/2}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}={\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}$,

$$\begin{align}\hfill & {\int }_{0}^{t}\langle {b}_{K}{\partial }_{t}^{{\alpha }_{K}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}(s),{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}(s)\rangle \mathrm{d}s\geqslant {\underline{C}}_{K,{\alpha }_{{k}_{\ast }}}{\Vert}{\partial }_{t}^{(1+{\alpha }_{{k}_{\ast }}+{\alpha }_{K})/2}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}{{\Vert}}_{{L}^{2}(0,t;{L}_{w}^{2}({\Omega}))}^{2}\hfill \end{align} \tag{ 33 }$$

with

$$\begin{equation*}{\underline{C}}_{K,{\alpha }_{{k}_{\ast }}}=\underline{C}((1+{\alpha }_{{k}_{\ast }}-{\alpha }_{K})/2){b}_{K}\,\mathrm{cos}(\pi (1+{\alpha }_{{k}_{\ast }}-{\alpha }_{K})/2).\end{equation*}$$

From the identity ${\int }_{0}^{t}\langle {\partial }_{t}v(s),v(s)\rangle \mathrm{d}s=\frac{1}{2}{\Vert}v(t){{\Vert}}^{2}-\frac{1}{2}{\Vert}v(0){{\Vert}}^{2}$ and Young's inequality in the case ${\gamma }_{J}={\alpha }_{{k}_{\ast }}$ we get

$$\begin{equation*}\begin{aligned}\hfill & {\int }_{0}^{t}\langle {d}_{J}{\partial }_{t}^{2+{\gamma }_{J}}{u}^{L}(s),{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{u}^{L}(s)\rangle \mathrm{d}s={\int }_{0}^{t}\langle {d}_{J}\left({\partial }_{t}{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{u}^{L}(s)-\frac{{s}^{-{\alpha }_{{k}_{\ast }}}}{{\Gamma}(1-{\alpha }_{{k}_{\ast }})}{u}_{2}^{L}\right)\!,{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{u}^{L}(s)\rangle \mathrm{d}s\hfill \\ \hfill & \geqslant \frac{1}{2}{\Vert}\sqrt{{d}_{J}}{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{u}^{L}(t){{\Vert}}_{{L}_{w}^{2}({\Omega})}^{2}-\frac{{t}^{2(1-{\alpha }_{{k}_{\ast }})}{\Vert}\sqrt{{d}_{J}}{u}_{2}^{L}{{\Vert}}_{{L}_{w}^{2}({\Omega})}^{2}}{{\Gamma}{(2-{\alpha }_{{k}_{\ast }})}^{2}}-\frac{1}{4}{\Vert}\sqrt{{d}_{J}}{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{u}^{L}{{\Vert}}_{{L}^{\infty }(0,t);{L}_{w}^{2}({\Omega})}^{2}\hfill \\ \hfill & \;\text{if}\;{\gamma }_{J}={\alpha }_{{k}_{\ast }}\in (0,1)\hfill \end{aligned}\end{equation*}$$

$$\begin{equation*}\begin{aligned}\hfill {\int }_{0}^{t}\langle {d}_{J}{\partial }_{t}^{2+{\gamma }_{J}}{u}^{L}(s),{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{u}^{L}(s)\rangle \mathrm{d}s& =\frac{1}{2}{\Vert}\sqrt{{d}_{J}}{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{u}^{L}\left.\right)(t){{\Vert}}_{{L}_{w}^{2}({\Omega})}^{2}-\frac{1}{2}{\Vert}\sqrt{{d}_{J}}{u}_{1+{\alpha }_{{k}_{\ast }}}^{L}{{\Vert}}_{{L}_{w}^{2}({\Omega})}^{2}\hfill \\ \hfill & \quad \text{if}\;{\gamma }_{J}={\alpha }_{{k}_{\ast }}\in \left\{0,1\right\}\hfill \end{aligned}\end{equation*}$$

and

$$\begin{equation*}\begin{aligned}\hfill & {\int }_{0}^{t}\langle {b}_{{k}_{\ast }}{\partial }_{t}^{{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}(s),{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}(s)\rangle \mathrm{d}s\hfill \\ \hfill & ={\int }_{0}^{t}\langle {b}_{{k}_{\ast }}{\partial }_{t}{\partial }_{t}^{{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}(s)-\frac{{s}^{-{\alpha }_{{k}_{\ast }}}}{{\Gamma}(1-{\alpha }_{{k}_{\ast }})}{\tilde{\mathcal{A}}}^{1/2}{u}_{1}^{L},{\partial }_{t}^{{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}(s)\rangle \mathrm{d}s\hfill \\ \hfill & \geqslant \frac{1}{2}{\Vert}\sqrt{{b}_{{k}_{\ast }}}{\partial }_{t}^{{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}\left.\right)(t){{\Vert}}_{{L}_{w}^{2}({\Omega})}^{2}-\frac{{t}^{2(1-{\alpha }_{{k}_{\ast }})}{\Vert}\sqrt{{b}_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}_{1}^{L}{{\Vert}}_{{L}_{w}^{2}({\Omega})}^{2}}{{\Gamma}{(2-{\alpha }_{{k}_{\ast }})}^{2}}\hfill \\ \hfill & \qquad -\frac{1}{4}{\Vert}\sqrt{{b}_{{k}_{\ast }}}{\partial }_{t}^{{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}{{\Vert}}_{{L}^{\infty }(0,t);{L}_{w}^{2}({\Omega})}^{2}\ \;\text{if}\;{\alpha }_{{k}_{\ast }}\in (0,1),\hfill \end{aligned}\end{equation*}$$

$$\begin{equation*}\begin{aligned}\hfill & {\int }_{0}^{t}\langle {b}_{{k}_{\ast }}{\partial }_{t}^{{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}(s),{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}(s)\rangle \mathrm{d}s\hfill \\ \hfill & =\frac{1}{2}{\Vert}\sqrt{{b}_{{k}_{\ast }}}{\partial }_{t}^{{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}\left.\right)(t){{\Vert}}_{{L}_{w}^{2}({\Omega})}^{2}-\frac{1}{2}{\Vert}\sqrt{{b}_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}_{{\alpha }_{{k}_{\ast }}}^{L}{{\Vert}}_{{L}_{w}^{2}({\Omega})}^{2}\hfill \\ \hfill & \qquad \qquad \text{if}\;{\alpha }_{{k}_{\ast }}\in \left\{0,1\right\}.\hfill \end{aligned}\end{equation*}$$

Hence in the case k_* = K, where due to our assumption that α_K > 0 at the beginning of the proof, ${\partial }_{t}^{{\alpha }_{K}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}(0)=0$,

$$\begin{equation*}\begin{aligned}\hfill & \underset{{t}^{\prime }\in (0,t)}{\mathrm{sup}}{\int }_{0}^{{t}^{\prime }}\langle {b}_{K}{\partial }_{t}^{{\alpha }_{K}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}(s),{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}(s)\rangle \mathrm{d}s\hfill \\ \hfill & \geqslant \frac{{b}_{K}}{4}{\Vert}{\partial }_{t}^{{\alpha }_{K}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}{{\Vert}}_{{L}^{\infty }(0,t);{L}_{w}^{2}({\Omega})}^{2}-\frac{{t}^{2(1-{\alpha }_{K})}{\Vert}\sqrt{{b}_{K}}{\tilde{\mathcal{A}}}^{1/2}{u}_{1}^{L}{{\Vert}}_{{L}_{w}^{2}({\Omega})}^{2}}{{\Gamma}{(2-{\alpha }_{K})}^{2}}.\hfill \end{aligned}\end{equation*}$$

Since the difference $1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k}$ is larger than one for k < k_*, we have to employ (16) in that case (for this purpose, note that still, according to our assumptions, $\theta =(1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k})/2\in [0,1]$). From this and (15), again with $\theta =(1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k})/2\in [0,1]$, we obtain

$$\begin{equation}\begin{aligned}\hfill & {\int }_{0}^{t}\langle {b}_{k}{\partial }_{t}^{{\alpha }_{k}}{u}^{L}(s),{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{u}^{L}(s)\rangle \mathrm{d}s\hfill \\ \hfill & ={\int }_{0}^{t}\langle {b}_{k0}{{I}_{t}}^{1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k}}{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}(s),{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}(s)\rangle \mathrm{d}s\hfill \\ \hfill & \geqslant -C(1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k}){\Vert}\sqrt{{b}_{k}}{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}{{\Vert}}_{{H}^{-(1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k})/2}(0,t;{L}_{w}^{2}({\Omega}))}^{2}\hfill \\ \hfill & \geqslant -C(1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k})\bar{C}((1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k})/2){\Vert}\sqrt{{b}_{k}}{\partial }_{t}^{(1+{\alpha }_{{k}_{\ast }}+{\alpha }_{k})/2}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}{{\Vert}}_{{L}^{2}(0,t;{L}_{w}^{2}({\Omega}))}^{2}\hfill \\ \hfill & \;\text{for}\;k< {k}_{\ast }.\hfill \end{aligned}\end{equation} \tag{ 34 }$$

Note that in case k_* < K, the potentially negative contributions arising from the k < k_* terms (34) are dominated by means of the leading energy term from (28), that is, ${\Vert}{\partial }_{t}^{(1+{\alpha }_{{k}_{\ast }}+{\alpha }_{K})/2}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}{{\Vert}}_{{L}^{2}(0,t;{L}_{w}^{2}({\Omega}))}^{2}$ from (33), due to the fact that $(1+{\alpha }_{{k}_{\ast }}+{\alpha }_{k})/2< (1+{\alpha }_{{k}_{\ast }}+{\alpha }_{K})/2$ for k < k_*.

To achieve this also in case k_* = K with α_K−1 < γ_J < α_K , we additionally test with $v={\partial }_{t}^{1+{\alpha }_{K}-\varepsilon }{u}^{L}(s)$ in (32) where 0 < ɛ < α_K − max{α_K−1, γ_J }. Repeating the above estimates with ${\alpha }_{{k}_{\ast }}$ replaced by α_K − ɛ one sees that the leading energy term then is ${\Vert}{\partial }_{t}^{{\alpha }_{K}+(1-\varepsilon )/2}{\tilde{\mathcal{A}}}^{1/2}{u}^{L}{{\Vert}}_{{L}^{2}(0,t;{L}_{w}^{2}({\Omega}))}^{2}$ from (33). This still dominates the potentially negative contributions arising from the k < K terms according to (34) (with α_K − ɛ in place of ${\alpha }_{{k}_{\ast }}$), since (1 + α_K − ɛ + α_k )/2 < α_K + (1 − ɛ)/2 for k < K.

Note that in the case γ_J = α_K = α₀, no potentially negative terms (34) arise.

Finally, the right-hand side term can be estimated by means of (15) with $\theta ={\alpha }_{{k}_{\ast }}-{\gamma }_{J}\in \left[0,1\right)$ and Young's inequality

$$\begin{equation*}\begin{aligned}\hfill & {\int }_{0}^{t}\langle r(s),{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{u}^{L}(s)\rangle \mathrm{d}s\leqslant {\Vert}r{{\Vert}}_{{H}^{{\alpha }_{{k}_{\ast }}-{\gamma }_{J}}(0,T,{L}_{w}^{2}({\Omega}))}{\Vert}{\partial }_{t}^{1+{\alpha }_{{k}_{\ast }}}{u}^{L}(s){{\Vert}}_{{H}^{-({\alpha }_{{k}_{\ast }}-{\gamma }_{J})}(0,t,{L}_{w}^{2}({\Omega}))}\hfill \\ \hfill & \leqslant \bar{d}{\Vert}{\partial }_{t}^{1+{\gamma }_{J}}{u}^{L}(s){{\Vert}}_{{L}^{2}(0,t,{L}_{w}^{2}({\Omega}))}^{2}+\frac{1}{4\bar{d}\bar{C}{({\alpha }_{{k}_{\ast }}-{\gamma }_{J})}^{2}}{\Vert}r{{\Vert}}_{{H}^{{\alpha }_{{k}_{\ast }}-{\gamma }_{J}}(0,T,{L}_{w}^{2}({\Omega}))}^{2}\hfill \\ \hfill & \leqslant {\int }_{0}^{t}{\underline{\mathcal{E}}}_{\gamma }[{u}^{L}](s)\mathrm{d}s+C{\Vert}r{{\Vert}}_{{H}^{{\alpha }_{{k}_{\ast }}-{\gamma }_{J}}(0,T,{L}_{w}^{2}({\Omega}))}^{2}\hfill \end{aligned}\end{equation*}$$

with some constants $\bar{d},C > 0$.

Altogether we arrive at the following estimate

$$\begin{equation}{\mathcal{E}}_{\gamma }[{u}^{L}](t)+{\mathcal{E}}_{\alpha }[{u}^{L}](t)\leqslant {\int }_{0}^{t}{\underline{\mathcal{E}}}_{\gamma }[{u}^{L}](s)\mathrm{d}s+\text{rhs}[{u}^{L}](t)+{\text{rhs}}_{0r}^{L}(t),\end{equation} \tag{ 35 }$$

where ${\mathcal{E}}_{\gamma }$, ${\underline{\mathcal{E}}}_{\gamma }$, ${\mathcal{E}}_{\alpha }$ are defined as in (26)–(28),

$$\begin{equation*}\begin{aligned}\hfill & \text{rhs}[u](t)=\sum\limits _{k=0}^{{k}_{\ast }-1}C(1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k})\bar{C}((1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k})/2){\Vert}\sqrt{{b}_{k}}{\partial }_{t}^{(1+{\alpha }_{{k}_{\ast }}+{\alpha }_{k})/2}{\tilde{\mathcal{A}}}^{1/2}u{{\Vert}}_{{L}^{2}(0,t;{L}_{w}^{2}({\Omega}))}^{2}\hfill \\ \hfill & {\text{rhs}}_{0,r}^{L}(t)={1}_{{\gamma }_{J}={\alpha }_{{k}_{\ast }}\in (0,1)}\frac{{t}^{2(1-{\alpha }_{{k}_{\ast }})}{\Vert}\sqrt{{d}_{J}}{u}_{2}^{L}{{\Vert}}_{{L}_{w}^{2}({\Omega})}^{2}}{{\Gamma}{(2-{\alpha }_{{k}_{\ast }})}^{2}}+{1}_{{\alpha }_{{k}_{\ast }}\in (0,1)}\frac{{t}^{2(1-{\alpha }_{{k}_{\ast }})}{\Vert}\sqrt{{b}_{{k}_{\ast }}}{\tilde{\mathcal{A}}}^{1/2}{u}_{1}^{L}{{\Vert}}_{{L}_{w}^{2}({\Omega})}^{2}}{{\Gamma}{(2-{\alpha }_{{k}_{\ast }})}^{2}}\hfill \\ \hfill & \qquad +\frac{1}{2\tilde{d}\bar{C}{({\alpha }_{{k}_{\ast }}-{\gamma }_{J})}^{2}}{\Vert}r{{\Vert}}_{{H}^{{\alpha }_{{k}_{\ast }}-{\gamma }_{J}}(0,T,{L}_{w}^{2}({\Omega}))}^{2},\hfill \end{aligned}\end{equation*}$$

where the Boolean variable 1_B is equal to one if B is true and vanishes otherwise.

We now substitute ${\mathcal{E}}_{\gamma }[u](t)$, ${\mathcal{E}}_{\alpha }[u](t)$ by lower bounds containing only an estimate from below of their leading order terms ${\underline{\mathcal{E}}}_{\gamma }[u](t)$ as defined in (27) and

$$\begin{equation}{\underline{\mathcal{E}}}_{\alpha }[u](t){:=}\begin{cases}{\underline{C}}_{K,{\alpha }_{{k}_{\ast }}}{\Vert}{\partial }_{t}^{(1+{\alpha }_{{k}_{\ast }}+{\alpha }_{K})/2}{\tilde{\mathcal{A}}}^{1/2}u{{\Vert}}_{{L}^{2}(0,t;{L}_{w}^{2}({\Omega}))}^{2}\quad \hfill & \;\text{if}\;{k}_{\ast }< K\hfill \\ {\underline{C}}_{K,{\alpha }_{K}-\varepsilon }{\Vert}{\partial }_{t}^{{\alpha }_{K}+(1-\varepsilon )/2}{\tilde{\mathcal{A}}}^{1/2}u{{\Vert}}_{{L}^{2}(0,t;{L}_{w}^{2}({\Omega}))}^{2}\quad \hfill & \;\text{if}\;{k}_{\ast }=K.\hfill \end{cases}\end{equation} \tag{ 36 }$$

We then obtain from (35)

$$\begin{equation}{\underline{\mathcal{E}}}_{\gamma }[{u}^{L}](t)+{\underline{\mathcal{E}}}_{\alpha }[{u}^{L}](t)\leqslant {\int }_{0}^{t}{\underline{\mathcal{E}}}_{\gamma }[{u}^{L}](s)\mathrm{d}s+\text{rhs}[{u}^{L}](t)+{\text{rhs}}_{0,r}^{L}(t)+{\tilde{\text{rhs}}}_{0}^{L}(t),\end{equation} \tag{ 37 }$$

where

$$\begin{equation*}{\tilde{\text{rhs}}}_{0}^{L}(t)={1}_{{k}_{\ast }=K,{\alpha }_{K}\in (0,1)}\frac{{t}^{2(1-{\alpha }_{K})}{\Vert}\sqrt{{b}_{K}}{\tilde{\mathcal{A}}}^{1/2}{u}_{1}^{L}{{\Vert}}_{{L}_{w}^{2}({\Omega})}^{2}}{{\Gamma}{(2-{\alpha }_{K})}^{2}}.\end{equation*}$$

It is readily checked that $\text{rhs}[u](t)\leqslant {c}_{b,\alpha }(t){\underline{\mathcal{E}}}_{\alpha }[u](t)$ for

$$\begin{align*}\hfill {c}_{b,\alpha }(t)& =\frac{1}{{\underline{C}}_{K,{\alpha }_{{k}_{\ast }}}}\sum\limits _{k=0}^{{k}_{\ast }-1}C(1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k})\bar{C}((1+{\alpha }_{{k}_{\ast }}-{\alpha }_{k})/2)\vert {b}_{k}\vert {{\Vert}}_{0}{I}_{t}^{({\alpha }_{K}-{\alpha }_{k})/2}{{\Vert}}_{{L}^{2}(0,t)\to {L}^{2}(0,t)}.\hfill \end{align*}$$

Therefore, under the smallness assumption (25), that is, c_b,α (T) < 1 estimate (37) yields

$$\begin{equation*}\underline{{\mathcal{E}}_{\gamma }}[{u}^{L}](t)+(1-{c}_{b,\alpha }(T)){\underline{\mathcal{E}}}_{\alpha }[{u}^{L}](t)\leqslant \tilde{C}\left({\int }_{0}^{t}{\underline{\mathcal{E}}}_{\gamma }[{u}^{L}](s)\mathrm{d}s+{\Vert}{u}_{0}{{\Vert}}_{{\dot{H}}^{1}({\Omega})}^{2}+{C}_{0}+{\Vert}r{{\Vert}}_{{H}^{{\alpha }_{{k}_{\ast }}-{\gamma }_{J}}(0,T,{L}_{w}^{2}({\Omega}))}^{2}\right)\end{equation*}$$

for some $\tilde{C}$ independent of t, L, with C₀ as in (31). Applying Gronwall's inequality and re-inserting into (35) we end up with

$$\begin{equation}{\mathcal{E}}_{\gamma }[{u}^{L}](t)+{\mathcal{E}}_{\alpha }[{u}^{L}](t)\leqslant C(T)\left({\Vert}{u}_{0}{{\Vert}}_{{\dot{H}}^{1}({\Omega})}^{2}+{C}_{0}+{\Vert}r{{\Vert}}_{{H}^{{\alpha }_{{k}_{\ast }}-{\gamma }_{J}}(0,T,{L}_{w}^{2}({\Omega}))}^{2}\right).\end{equation} \tag{ 38 }$$

Step 3. Weak limits.

As a consequence of (38), the Galerkin solutions u^L are uniformly bounded in U_γ ∩ U_α . Therefore the sequence ${({u}^{L})}_{L\in \mathbb{N}}$ has a weakly(^*) convergent subsequence ${({u}^{{L}_{k}})}_{k\in \mathbb{N}}$ with limit u ∈ U_γ ∩ U_α . To see that u satisfies (20), we revisit (32), integrate it with respect to time and take the L²(0, T) product with arbitrary smooth test functions ψ to conclude that

$$\begin{equation*}\begin{aligned}\hfill & {\int }_{0}^{T}\psi (t)\left\{\langle \sum\limits _{j=0}^{J}{d}_{j}\left({\partial }_{t}^{1+{\gamma }_{j}}u(t)-\frac{{t}^{1-{\gamma }_{j}}}{{\Gamma}(2-{\gamma }_{j})}{u}_{2}\right),{v\rangle }_{{L}_{w}^{2}({\Omega})}\right.\hfill \\ \hfill & \qquad \qquad +\langle \sum\limits _{k=0}^{K}{b}_{k}{{\left(\right.}_{0}{I}_{t}}^{1-{\alpha }_{k}}\tilde{\mathcal{A}}u(t)-\frac{{t}^{1-{\alpha }_{k}}}{{\Gamma}(2-{\alpha }_{k})}\tilde{\mathcal{A}}{u}_{0}\left.\right)-r(t),{v\rangle }_{\dot{H}{({\Omega})}^{\ast },\dot{H}({\Omega})}\left.\right\}\mathrm{d}t\hfill \\ \hfill & \qquad \quad ={\int }_{0}^{T}\psi (t)\left\{\langle \sum\limits _{j=0}^{J}{d}_{j}{\partial }_{t}^{1+{\gamma }_{j}}(u-{u}^{{L}_{k}})(t),{v\rangle }_{{L}_{w}^{2}({\Omega})}\right.\hfill \\ \hfill & \left.\qquad \qquad +\langle \sum\limits _{k=0}^{K}{b}_{k0}{{I}_{t}}^{1-{\alpha }_{k}}\tilde{\mathcal{A}}(u-{u}^{{L}_{k}})(t),{v\rangle }_{\dot{H}{({\Omega})}^{\ast },\dot{H}({\Omega})}\right\}\mathrm{d}t\hfill \end{aligned}\end{equation*}$$

for all v ∈ span(φ₁, ..., φ_K ), K ⩽ L_k and any $\psi \in {C}_{c}^{\infty }(0,T)$. Taking the limit k → ∞ by using the previously mentioned weak(^*) limit of ${u}^{{L}_{k}}$ we conclude

$$\begin{equation*}\begin{aligned}\hfill & {\int }_{0}^{T}\psi (t)\left\{\langle \sum\limits _{j=0}^{J}{d}_{j}\left({\partial }_{t}^{1+{\gamma }_{j}}u(t)-\frac{{t}^{1-{\gamma }_{j}}}{{\Gamma}(2-{\gamma }_{j})}{u}_{2}\right),{v\rangle }_{{L}_{w}^{2}({\Omega})}\right.\hfill \\ \hfill & \qquad \qquad +\langle \sum\limits _{k=0}^{K}{b}_{k}{{\left(\right.}_{0}{I}_{t}}^{1-{\alpha }_{k}}\tilde{\mathcal{A}}u(t)-\frac{{t}^{1-{\alpha }_{k}}}{{\Gamma}(2-{\alpha }_{k})}\tilde{\mathcal{A}}{u}_{0}\left.\right)-r(t),{v\rangle }_{\dot{H}{({\Omega})}^{\ast },\dot{H}({\Omega})}\left.\right\}\mathrm{d}t=0\hfill \end{aligned}\end{equation*}$$

for all v ∈ span(φ₁, ..., φ_K ) with arbitrary $K\in \mathbb{N}$, and any $\psi \in {C}_{c}^{\infty }(0,T)$. Therefore the weak limit u indeed satisfies the time integrated pde and we have proven the existence part of the theorem.

To verify the initial conditions, we use the fact that due to our assumption α_K > 0, U_α continuously embeds into $C([0,T];\dot{H}({\Omega}))$ and therefore u(0) = u₀ is attained in an $\dot{H}({\Omega})$ sense. Also, U_γ continuously embeds into ${C}^{1}([0,T];{L}_{w}^{2}({\Omega}))$ in case ${\gamma }_{J}< {\alpha }_{{k}_{\ast }}$ or γ_J > 0. In the remaining case ${\gamma }_{J}={\alpha }_{{k}_{\ast }}=0$, attainment of u_t (0) = u₁ in an ${L}_{w}^{2}({\Omega})$ sense can be shown analogously to the conventional second order wave equation, cf, e.g., [11, theorem 3 in section 7.2].

Moreover, taking weak limits in the energy estimate (38) together with weak lower semicontinuity of the norms contained in the definitions of ${\mathcal{E}}_{\gamma }$, ${\mathcal{E}}_{\alpha }$, implies (30).

Step 4. Uniqueness.

The manipulations carried out in step 2 of the proof are also feasible with the Galerkin approximation u^L replaced by a solution u of the pde itself, and lead to the energy estimate (30) independently of the Galerkin approximation procedure. From this we conclude that the pde with vanishing right-hand side and initial data only has the zero solution, which due to linearity of the pde implies uniqueness.□

Corollary 2.1. Under the assumptions of theorem 2.1 with b_k = 0, k = 0, ..., k_* − 1, the assertion extends to hold globally in time, that is, for all t ∈ (0, ∞). If additionally, r = 0, then the energy is globally bounded

$$\begin{equation*}{\mathcal{E}}_{\gamma }[u](t)+{\mathcal{E}}_{\alpha }[u](t)\leqslant C\left({\Vert}{u}_{0}{{\Vert}}_{{\dot{H}}^{1}({\Omega})}^{2}+{C}_{0}\right),\quad t\in (0,\infty )\end{equation*}$$

with some constant C > 0 independent of time.

Remark 2.1. Using the so-called multinomial (or multivariate) Mittag–Leffler functions and separation of variables in principle enables a solution representation by separation of variables, cf [26, theorem 4.1]. However, proving convergence of the infinite sums in this representation would require extensive estimates of these Mittag–Leffler functions. Thus, in order to keep the exposition transparent without having to introduce too much additional machinery, we chose to remain with the energy argument in the proof above.

In order to establish sufficient regularity of the highest order term ${\partial }_{t}^{2+{\gamma }_{J}}{u}^{L}$ so that also the original equation (17) holds, we make use of the pde itself.

Corollary 2.2. Under the conditions of theorem 2.1 the solution u of (20) satisfies ${\partial }_{t}^{2+{\gamma }_{J}}u\in {L}^{2}(0,T;\dot{H}{({\Omega})}^{\ast })$ and is therefore also a solution to the original pde (17) in an ${L}^{2}(0,T;\dot{H}{({\Omega})}^{\ast })$ sense.

Proof. From the Galerkin approximation (32) of (17), due to invariance of the space span(φ₁, ..., φ_L ) we can substitute v by ${\tilde{\mathcal{A}}}^{-1/2}v$ there and move ${\tilde{\mathcal{A}}}^{-1/2}$ to the left-hand side of the inner product by taking the adjoint to arrive at

$$\begin{equation}\begin{aligned}\hfill & \langle \sum\limits _{j=0}^{J}{\partial }_{t}^{2+{\gamma }_{j}}{\tilde{\mathcal{A}}}^{-1/2}[{d}_{j}{u}^{L}(t)]+\sum\limits _{k=0}^{K}{\partial }_{t}^{{\alpha }_{k}}{\tilde{\mathcal{A}}}^{-1/2}[{b}_{k}\tilde{\mathcal{A}}{u}^{L}(t)]-{\tilde{\mathcal{A}}}^{-1/2}r(t),v{\rangle }_{{L}_{w}^{2}({\Omega})}=0\hfill \\ \hfill & \qquad \quad t\in (0,T)\quad \text{for}\,\text{all}\;v\in \text{span}({\varphi }_{1},\dots ,{\varphi }_{L}).\hfill \end{aligned}\end{equation} \tag{ 39 }$$

From theorem 2.1 we know that $\tilde{r}=-{\sum }_{k=0}^{K}\,{\partial }_{t}^{{\alpha }_{k}}{\tilde{\mathcal{A}}}^{-1/2}[{b}_{k}\tilde{\mathcal{A}}u(t)]+{\tilde{\mathcal{A}}}^{-1/2}r(t)$ is contained in ${L}^{2}(0,T;{L}_{w}^{2}({\Omega}))$. Thus (39) can be viewed as a Galerkin discretisation of a Volterra integral equation for $\zeta ={\partial }_{t}^{2+{\gamma }_{J}}{\tilde{\mathcal{A}}}^{-1/2}u(t)$ with $\tilde{r}$ as right-hand side. The corresponding coefficient vector ${\underline{\zeta }}^{L}$ of the solution ζ^L is therefore bounded by

$$\begin{equation*}{\Vert}{\underline{\zeta }}^{L}{{\Vert}}_{{\ell }^{2}({\mathbb{R}}^{L})}\leqslant C({\Vert}\tilde{r}{{\Vert}}_{{L}^{2}(0,T;{L}_{w}^{2}({\Omega}))}+{\Vert}\zeta (0){{\Vert}}_{{L}_{w}^{2}({\Omega})})\end{equation*}$$

with C independent of L. From this we deduce

$$\begin{equation*}{\Vert}{\zeta }^{L}{{\Vert}}_{{L}^{2}(0,T;{L}_{w}^{2}({\Omega}))}\leqslant C({\Vert}\tilde{r}{{\Vert}}_{{L}^{2}(0,T;{L}_{w}^{2}({\Omega}))}+{\Vert}\zeta (0){{\Vert}}_{{L}_{w}^{2}({\Omega})})\end{equation*}$$

and thus, by taking the limit L → ∞ the same bound for ${\Vert}{\zeta }^{L}{{\Vert}}_{{L}^{2}(0,T;{L}_{w}^{2}({\Omega}))}$ itself.□

There are several methods available for solving multi-term fractional ordinary differential equations. For our situation we require a solver that is accurate over the entire time interval $\left[0,\infty \right)$ since one problem we consider is the recovery of the components of the operator from large time values only. Our method of choice here is to convert the system with multiple fractional derivatives to a single ode of the form ${D}_{t}^{\gamma }u=Au$ where A is a square matrix. The seminal starting point here is the paper by Diethelm and Ford [7], which allows one to consider the equation ${D}^{{\alpha }_{0}}y(t)=f\left(t,y(t),{D}^{{\alpha }_{1}},\dots ,{D}^{{\alpha }_{n}}\right)$ subject to ${y}^{(k)}(0)={y}_{0}^{k}$ for k = 0, 1, ... ⌈α₀⌉ − 1, and its conversion to the form ${D}^{\gamma }Y(t)=g\left(t,Y(t)\right)$, with Y(0) = Y₀. Existence and uniqueness results are obtained and a numerical method for the solution formulated. This was later refined by Garrappa and Popolizio [15], in the case of linear systems into a solver whose accuracy approaches machine precision. This is based on using the Mittag–Leffler function with a matrix argument A that uses different approaches depending of the location of the eigenvalues of A in the complex plane and combined into an algorithm that works for matrices with any eigenvalue locations. See also [28]. We describe this above approach briefly since it comes with a certain restriction that should be clarified.

The following result is proven in [28]:

Theorem 2.2. Consider the initial value problem

$$\begin{align}\hfill {D}^{{\alpha }_{0}}& =f\left(t,y(t),{D}^{{\alpha }_{1}}y(t),\dots ,{D}^{{\alpha }_{n-1}}y(t)\right),\hfill \\ \hfill {y}^{(j)}(0)& ={y}_{0}(j),\quad j=0,\dots ,[{\alpha }_{0}]-1\hfill \end{align} \tag{ 40 }$$

where α₀ > α₁ > ⋯ α_k−1 and 0 < α_k−1 < 1. Assume that each ${\alpha }_{j}\in \mathbb{Q}$ and let M be the least common denominator of α₀, ..., α_k−1. Set γ = 1/M and N = Mα₀. Then the initial value problem (41) is equivalent to the system of equations

$$\begin{align}{D}^{\gamma }{y}_{0}(t)& ={y}_{1}(t),\qquad {D}^{\gamma }{y}_{1}(t)={y}_{2}(t),\\ {D}^{\gamma }{y}_{N-1}(t)& =f(t,{y}_{0}(t),{y}_{{\alpha }_{1}/\gamma }(t),\dots ,{y}_{{\alpha }_{k-1}/\gamma }(t))\end{align} \tag{ 41 }$$

together with the initial conditions

$$\begin{equation}{y}_{j}(0)={y}_{0}^{(j/M)}\quad \text{if}\;j/M\in {\mathbb{Z}}_{0},\enspace \text{else}\ 0.\end{equation} \tag{ 42 }$$

As a result of the above theorem the linear system can be reformulated to the form

$$\begin{equation}{D}^{\gamma }Y(t)=AY(t)+{e}_{N}f(t),\qquad Y(0)={Y}_{0},\quad {e}_{N}={(0,0,\dots ,1)}^{\mathrm{T}}\in {\mathbb{R}}^{N},\end{equation} \tag{ 43 }$$

where the coefficient matrix $A\in {\mathbb{R}}^{N\times N}$ is of companion type.

Then it is possible to express the solution to (43) in terms of Mittag–Leffler functions with matrix arguments as

$$\begin{equation}Y(t)={E}_{\alpha ,1}({t}^{\gamma }A){Y}_{0}+{\int }_{0}^{t}{(t-s)}^{\gamma -1}{E}_{\gamma ,\gamma }\left({(t-s)}^{\gamma }A\right){e}_{N}f(s)\mathrm{d}s.\end{equation} \tag{ 44 }$$

For general matrix functions the Schur–Parlett algorithm is the method of choice and to this end the final component needed is an accurate evaluation of the matrix Mittag–Leffler function and its derivatives. This is both nontrivial and quite technical and we simply refer to the above references [15, 28], and also to [14] for the details. Code to implement what is needed for this step is available from Roberto Garrappa's homepage [13].

Note that in our case we have α₀ = 2 and also the requirement that each of the exponents α_j = q_j /p_j must be rational numbers. The above restrictions can lead to a quite large system of size lcm(p₁, ..., p_k−1) and hence a small value of γ. However, we use this solver only for the purpose of providing simulated data and our inversion techniques are independent of its use.

Since we are using asymptotic formulas it is essential to use a reliable and accurate solver that is capable to cover a wide range of time values. The code from [14, 15, 28] provably fulfils these requirements and our numerical experiments confirmed this.

Separation of variables based on the eigensystem of $\mathcal{A}$ yields a solution representation

$$\begin{align*}\hfill u(x,t)& =\sum\limits _{n=1}^{\infty }\left(\sigma (t){w}_{fn}(t)\langle f,{\varphi }_{n}\rangle +{w}_{2n}(t)\langle {u}_{2},{\varphi }_{n}\rangle +{w}_{1n}(t)\langle {u}_{1},{\varphi }_{n}\rangle \right.+\left.{w}_{0n}(t)\langle {u}_{0},{\varphi }_{n}\rangle \right){\varphi }_{n}(x)\hfill \end{align*}$$

and convergence of the sums is guaranteed by the energy estimates from theorem 2.1. Here w_fn , w_2n , w_1n , w_0n denote the solutions of the resolvent equation

$$\begin{equation*}\begin{aligned}\hfill & {w}_{n,tt}+{c}^{2}{\lambda }_{n}{w}_{n}+\sum\limits _{j=1}^{J}{d}_{j}{D}_{t}^{2+{\gamma }_{j}}{w}_{n}+\sum\limits _{k=1}^{N}{b}_{k}{\lambda }_{n}^{{\beta }_{k}}{D}_{t}^{{\alpha }_{k}}{w}_{n}=\begin{cases}1\ \;\text{for}\;{w}_{fn}\quad \hfill \\ 0\ \;\text{else}\quad \hfill \end{cases}\hfill \\ \hfill & {w}_{n,tt}(0)=\begin{cases}1\ \;\text{for}\;{w}_{2n}\quad \hfill \\ 0\ \;\text{else}\quad \hfill \end{cases},\quad {w}_{n,t}(0)=\begin{cases}1\ \;\text{for}\;{w}_{1n}\quad \hfill \\ 0\ \;\text{else}\quad \hfill \end{cases},\quad {w}_{n}(0)=\begin{cases}1\ \;\text{for}\;{w}_{0n}\quad \hfill \\ 0\ \;\text{else}\quad \hfill \end{cases}.\hfill \end{aligned}\end{equation*}$$

From the Laplace transformed resolvent equations we obtain the resolvent solutions

$$\begin{equation*}\begin{aligned}\hfill & {\hat{w}}_{fn}(s)=\frac{1}{\omega (s;{\lambda }_{n})},\hfill \\ \hfill & {\hat{w}}_{0n}(s)=\frac{s+{\sum }_{j=1}^{J}{\,d}_{j}\,{s}^{{\gamma }_{j}+1}+{\sum }_{k=1}^{N}\,{s}^{{\alpha }_{k}-1}{b}_{k}{\lambda }_{n}^{{\beta }_{k}}}{\omega (s;{\lambda }_{n})}=\frac{\omega (s,{\lambda }_{n})-{c}^{2}{\lambda }_{n}}{s\omega (s;{\lambda }_{n})},\hfill \\ \hfill & {\hat{w}}_{1n}(s)=\frac{1+{\sum }_{j=1}^{J}{\,d}_{j}{s}^{{\gamma }_{j}}}{\omega (s;{\lambda }_{n})},\qquad {\hat{w}}_{2n}(s)=\frac{{\sum }_{j=1}^{J}{\,d}_{j}{s}^{{\gamma }_{j}-1}}{\omega (s;{\lambda }_{n})}\hfill \end{aligned}\end{equation*}$$

with

$$\begin{equation*}\omega (s;{\lambda }_{n})={s}^{2}+{c}^{2}{\lambda }_{n}+\sum\limits _{j=1}^{J}{d}_{j}{s}^{{\gamma }_{j}+2}+\sum\limits _{k=1}^{N}{s}^{{\alpha }_{k}}{b}_{k}{\lambda }_{n}^{{\beta }_{k}},\end{equation*}$$

where $\hat{w}$ denotes the Laplace transform of w.

Formally we can write

$$\begin{align*}\hfill \hat{u}(x,s)& =\sum\limits _{n=1}^{\infty }\left(\hat{\sigma }(s){\hat{w}}_{fn}(s)\langle f,{\varphi }_{n}\rangle +{\hat{w}}_{2n}(s)\langle {u}_{2},{\varphi }_{n}\rangle +{\hat{w}}_{1n}(s)\langle {u}_{1},{\varphi }_{n}\rangle \right.+\left.{\hat{w}}_{0n}(s)\langle {u}_{0},{\varphi }_{n}\rangle \right){\varphi }_{n}(x)\hfill \end{align*}$$

and therefore

$$\begin{align*}\hfill \hat{{B}_{i}u}(s)& =\sum\limits _{n=1}^{\infty }\left(\hat{\sigma }(s){\hat{w}}_{fn}(s)\langle f,{\varphi }_{n}\rangle +{\hat{w}}_{2n}(s)\langle {u}_{2},{\varphi }_{n}\rangle +{\hat{w}}_{1n}(s)\langle {u}_{1},{\varphi }_{n}\rangle \right.+\left.{\hat{w}}_{0n}(s)\langle {u}_{0},{\varphi }_{n}\rangle \right){B}_{i}{\varphi }_{n}.\hfill \end{align*}$$

We make use of the Weyl estimate λ_n ≈ C_d n^2/d in ${\mathbb{R}}^{d}$ as well as additional regularity of the excitation to approximate the sum over n (46) of fractions by a single fraction of the form

$$\begin{align*}\hfill \sum\limits _{n=1}^{\infty }\frac{\langle {f}_{i},{\varphi }_{n}\rangle {B}_{i}{\varphi }_{n}}{\omega (s,{\lambda }_{n})}& \approx \frac{{\sum }_{n=1}^{\infty }\langle {f}_{i},{\varphi }_{n}\rangle {B}_{i}{\varphi }_{n}}{\omega (s,{\lambda }_{1})}\text{with}\;\omega (s,\lambda )={s}^{2}+{c}^{2}\lambda +\sum\limits _{j=1}^{J}{d}_{j}{s}^{{\gamma }_{j}+2}+\sum\limits _{k=1}^{N}{s}^{{\alpha }_{k}}{b}_{k}{\lambda }^{{\beta }_{k}}\hfill \end{align*}$$

for large real positive s. We focus exposition on the source excitation case u_ℓ = 0, ℓ = 0, 1, 2, f ≠ 0 from (46), the other cases follow similar. Also, since we only need one excitation here I = 1, we just skip the subscript i. Indeed, for $s\in \mathbb{R}$, s ⩾ 1 we can estimate the difference

$$\begin{equation*}\begin{aligned}\hfill D{:=}& \left\vert \frac{{\sum }_{n=1}^{\infty }\langle f,{\varphi }_{n}\rangle B{\varphi }_{n}}{\omega (s,{\lambda }_{1})}-\sum\limits _{n=1}^{\infty }\frac{\langle f,{\varphi }_{n}\rangle B{\varphi }_{n}}{\omega (s,{\lambda }_{n})}\right\vert \hfill \\ \hfill & =\left\vert \sum\limits _{n=1}^{\infty }\frac{\langle f,{\varphi }_{n}\rangle B{\varphi }_{n}\left({c}^{2}({\lambda }_{n}-{\lambda }_{1})+{\sum }_{k=1}^{N}\,{s}^{{\alpha }_{k}}{b}_{k}({({\lambda }_{n})}^{{\beta }_{k}}-{({\lambda }_{1})}^{{\beta }_{k}})\right)}{\omega (s,{\lambda }_{n})\omega (s,{\lambda }_{1})}\right\vert \hfill \\ \hfill & \leqslant \underset{n\in \mathbb{N}}{\mathrm{sup}}\left\{{\lambda }_{n}^{-\nu }\vert B{\varphi }_{n}\vert \right\}{\Vert}f{{\Vert}}_{{\dot{H}}^{2\mu +2\nu }({\Omega})}{\left(\sum\limits _{n=1}^{\infty }{\left(\frac{{c}^{2}{\lambda }_{n}+{\sum }_{k=1}^{N}\,{s}^{{\alpha }_{k}}{b}_{k}{({\lambda }_{n})}^{{\beta }_{k}}}{{\lambda }_{n}^{\mu }{s}^{2}({s}^{2}+{c}^{2}{\lambda }_{n})}\right)}^{2}\right)}^{1/2}\hfill \\ \hfill & \leqslant \underset{n\in \mathbb{N}}{\mathrm{sup}}\left\{{\lambda }_{n}^{-\nu }\vert B{\varphi }_{n}\vert \right\}{\Vert}f{{\Vert}}_{{\dot{H}}^{2\mu +2\nu }({\Omega})}{s}^{-(2-{\alpha }_{N}+2/p)}\frac{{c}^{2}+{\sum }_{k=1}^{N}\,{b}_{k}{({\lambda }_{1})}^{{\beta }_{k}-1}}{{p}^{1/p}{(pp-1{c}^{2})}^{\frac{p-1}{p}}}{\left(\sum\limits _{n=1}^{\infty }{\lambda }_{n}^{-2(\mu -1/p)}\right)}^{1/2},\hfill \end{aligned}\end{equation*}$$

where we have used Young's inequality to estimate

$$\begin{equation*}{s}^{2}+{c}^{2}{\lambda }_{n}\geqslant {s}^{2/p}{p}^{1/p}{(\frac{p}{p-1}{c}^{2}{\lambda }_{n})}^{\frac{p-1}{p}}.\end{equation*}$$

Since λ_n ∼ n^2/d , the sum converges for $2(\mu -\frac{1}{p}) > \frac{d}{2}$. On the other hand, to get an $O({s}^{-(2+{\gamma }_{J}+\varepsilon )})$ estimate for D with ɛ > 0, we need 2 − α_N + 2/p > 2 + γ_J and therefore choose

$$\begin{equation}\mu > \frac{d}{4}+\frac{\mathrm{max}\left\{{\alpha }_{N}+{\gamma }_{J},{\tilde{\alpha }}_{\tilde{N}}+{\tilde{\gamma }}_{\tilde{J}}\right\}}{2},\qquad p{:=}\frac{2}{\mathrm{max}\left\{{\alpha }_{N}+{\gamma }_{J},{\tilde{\alpha }}_{\tilde{N}}+{\tilde{\gamma }}_{\tilde{J}}\right\}}\end{equation} \tag{ 47 }$$

(taking into account the tilde version of the above estimate as well). This yields

$$\begin{equation*}\begin{aligned}\hfill & \frac{{\sum }_{n=1}^{\infty }\langle {f}_{i}{\varphi }_{n}\rangle {B}_{i}{\varphi }_{n}}{{s}^{2}+{c}^{2}{\lambda }_{1}+{\sum }_{j=1}^{J}{\,d}_{j}{s}^{{\gamma }_{j}+2}+{\sum }_{k=1}^{N}\,{s}^{{\alpha }_{k}}{b}_{k}{({\lambda }_{1})}^{{\beta }_{k}}}\hfill \\ \hfill & \qquad =\frac{{\sum }_{n=1}^{\infty }\langle {\tilde{f}}_{i}{\tilde{\varphi }}_{n}\rangle {\tilde{B}}_{i}{\tilde{\varphi }}_{n}}{{s}^{2}+{\tilde{c}}^{2}{\tilde{\lambda }}_{1}+{\sum }_{j=1}^{\tilde{J}}{\tilde{\,d}}_{j}\,{s}^{{\tilde{\gamma }}_{j}+2}+{\sum }_{k=1}^{\tilde{N}}\,{s}^{{\tilde{\alpha }}_{k}}{\tilde{b}}_{k}{({\tilde{\lambda }}_{1})}^{{\tilde{\beta }}_{k}}}+O({s}^{-(2+\mathrm{max}\left\{{\gamma }_{J},{\tilde{\gamma }}_{\tilde{J}}\right\}+\varepsilon )}),\quad s\in [1,\infty )\hfill \end{aligned}\end{equation*}$$

with $\varepsilon =\mu -\frac{d}{4}-\frac{\mathrm{max}\left\{{\alpha }_{N}+{\gamma }_{J},{\tilde{\alpha }}_{\tilde{N}}+{\tilde{\gamma }}_{\tilde{J}}\right\}}{2} > 0$. Thus, our assumptions on μ, ν are

$$\begin{align}\hfill & \mu > \frac{d}{4}+\frac{\mathrm{max}\left\{{\alpha }_{N}+{\gamma }_{J},{\tilde{\alpha }}_{\tilde{N}}+{\tilde{\gamma }}_{\tilde{J}}\right\}}{2},\enspace \underset{n\in \mathbb{N}}{\mathrm{sup}}\left\{{\lambda }_{n}^{-\nu }\vert B{\varphi }_{n}\vert \right\},\enspace \underset{n\in \mathbb{N}}{\mathrm{sup}}\left\{{\tilde{\lambda }}_{n}^{-\nu }\vert B{\tilde{\varphi }}_{n}\vert \right\}< \infty .\hfill \end{align} \tag{ 48 }$$

In case of rational powers

$$\begin{equation}\begin{aligned}\hfill {\alpha }_{k}& =\frac{{p}_{k}}{{q}_{k}},\qquad {\tilde{\alpha }}_{k}=\frac{{\tilde{p}}_{k}}{{\tilde{q}}_{k}},\qquad {\gamma }_{j}=\frac{{n}_{j}}{{m}_{j}},\qquad {\tilde{\gamma }}_{j}=\frac{{\tilde{n}}_{j}}{{\tilde{m}}_{j}},\hfill \\ \hfill \bar{q}& =\mathrm{lcm}({q}_{1},\dots ,{q}_{N},{\tilde{q}}_{1},\dots ,{\tilde{q}}_{\tilde{N}},{m}_{1},\dots ,{m}_{J},{\tilde{m}}_{1},\dots ,{\tilde{m}}_{\tilde{J}}),\hfill \end{aligned}\end{equation} \tag{ 49 }$$

setting $z={s}^{1/\bar{q}}$ we read (51) as an equality between two rational functions of z, whose coefficients and powers therefore need to coincide (provided $\left\{{s}^{1/\bar{q}}:s\in {\mathbb{C}}_{\theta }\right\}\supseteq [\underline{z},\infty )$ for some $\underline{z}\geqslant 0$). In the general case of real powers one can use the induction proof from [18, theorem 1.1] to obtain the same uniqueness. This yields

$$\begin{equation*}\begin{aligned}\hfill & (\mathrm{I})\qquad \hfill & \hfill & N=\tilde{N},\qquad J=\tilde{J},\hfill \\ \hfill & (\mathrm{I}\mathrm{I})\hfill & \hfill & {c}^{2}{\lambda }_{1}={\tilde{c}}^{2}{\tilde{\lambda }}_{1},\hfill \\ \hfill & (\mathrm{I}\mathrm{I}\mathrm{I})\hfill & \hfill & {\alpha }_{k}={\tilde{\alpha }}_{k},\ {\gamma }_{j}={\tilde{\gamma }}_{j},\hfill & \hfill & k=1,\dots ,N,\ j=1,\dots ,J\hfill \\ \hfill & (\mathrm{I}\mathrm{V})\hfill & \hfill & {b}_{k}{({\lambda }_{1})}^{{\beta }_{k}}={\tilde{b}}_{k}{({\tilde{\lambda }}_{1})}^{{\tilde{\beta }}_{k}},\ {d}_{j}={\tilde{d}}_{j},\quad \hfill & \hfill & \ k=1,\dots ,N,\ j=1,\dots ,J.\hfill \end{aligned}\end{equation*}$$

From this it becomes clear that we get uniqueness of

$$\begin{equation*}N,\ J,\ c,\ {b}_{k},\ {\alpha }_{k},\ {d}_{j},\ {\gamma }_{j},\end{equation*}$$

provided all ${\beta }_{k}={\tilde{\beta }}_{k}$ are known and ${\lambda }_{1}={\tilde{\lambda }}_{1}$ (but still possibly unknown) since we can then divide by λ₁ and ${\lambda }_{1}^{{\beta }_{k}}$ in (II) and (IV), respectively. However, there seems to be no chance to simultaneously obtain b_k and β_k from (IV). Therefore, we will consider a setting with two different excitations in section 3.2 for this purpose.

The same derivation can be made for the variable c = c(x) model (6) in place of (5), with c, $\mathcal{A}$, L²(Ω), ${\Vert}\varphi {{\Vert}}_{{L}^{2}}=1$ replaced by 1, ${\mathcal{A}}_{c}$, ${L}_{1/{c}^{2}}^{2}({\Omega})$, ${\Vert}\varphi {{\Vert}}_{{L}_{1/{c}^{2}}^{2}}=1$. From (II) with c, $\tilde{c}$ replaced by unity, we get ${\lambda }_{1}={\tilde{\lambda }}_{1}$ so that we do not need to assume this.

Theorem 3.1. Let $f,\tilde{f}\in {\dot{H}}^{2\mu +2\nu }({\Omega})$ for μ, ν sufficiently large so that (48) holds.

Then

$$\begin{equation*}Bu(t)=\tilde{B}\tilde{u}(t)\quad t\in (0,T).\end{equation*}$$

• (a)  
  For the solutions u, $\tilde{u}$ of (5) with vanishing initial data and known equal ${\beta }_{k}={\tilde{\beta }}_{k}$, possibly unknown but equal ${\lambda }_{1}={\tilde{\lambda }}_{1}$, possibly unknown (rest of) $\mathcal{A}$, $\tilde{\mathcal{A}}$, and possibly unknown $f,\tilde{f}$ implies

  $$\begin{align*}\hfill & N=\tilde{N},\qquad J=\tilde{J},\qquad c=\tilde{c},\qquad {b}_{k}={\tilde{b}}_{k},\qquad {\alpha }_{k}={\tilde{\alpha }}_{k},\qquad {d}_{j}={\tilde{d}}_{j},\qquad {\gamma }_{j}={\tilde{\gamma }}_{j},\hfill \\ \hfill & \quad k\in \left\{1,\dots ,N\right\},\enspace j\in \left\{1,\dots ,J\right\}.\hfill \end{align*}$$
• (b)  
  For the solutions u, $\tilde{u}$ of (6) with vanishing initial data and known equal ${\beta }_{k}={\tilde{\beta }}_{k}$, and possibly unknown ${\mathcal{A}}_{c}$, ${\tilde{\mathcal{A}}}_{\tilde{c}}$ implies

  $$\begin{align*}\hfill & N=\tilde{N},\qquad J=\tilde{J},\qquad {b}_{k}={\tilde{b}}_{k},\qquad {\alpha }_{k}={\tilde{\alpha }}_{k},\qquad {d}_{j}={\tilde{d}}_{j},\qquad {\gamma }_{j}={\tilde{\gamma }}_{j},\hfill \\ \hfill & \quad k\in \left\{1,\dots ,N\right\},\enspace j\in \left\{1,\dots ,J\right\}.\hfill \end{align*}$$

Remark 3.1. This uniqueness approach also extends to different excitation combinations such as u₀ ≠ 0, u_j = 0, j = 1, 2, f = 0; u₁ ≠ 0, u_j = 0, j = 0, 2, f = 0; u₂ ≠ 0, u_j = 0, j = 0, 1, f = 0; (the latter only in case γ_J > 0).

Consider for example the case u₀ ≠ 0, u_j = 0, j = 1, 2, f = 0 that will be relevant for additionally recovering u₀, cf remark 3.6 and which is the one that differs most from the setting above in view of the fact that the numerator of the resolvent solution depends on s and λ_n . The crucial estimate on the error that is made by replacing λ_n by λ₁ then becomes, for $s\in \mathbb{R}$, s ⩾ 1, and using the decomposition ω(s, λ) = s² + ω_d (s) + ω_b (s, λ) with ${\omega }_{d}(s)={\sum }_{j=1}^{J}{\,d}_{j}\,{s}^{{\gamma }_{j}+2}$, ${\omega }_{b}(s,\lambda )={\sum }_{k=1}^{N}\,{s}^{{\alpha }_{k}}{b}_{k}{\lambda }^{{\beta }_{k}}$

$$\begin{equation*}\begin{aligned}\hfill D{:=}& \left\vert \sum\limits _{n=1}^{\infty }\frac{\langle {u}_{0},{\varphi }_{n}\rangle B{\varphi }_{n}(\omega (s,{\lambda }_{1})-{c}^{2}{\lambda }_{1})}{s\omega (s,{\lambda }_{1})}-\sum\limits _{n=1}^{\infty }\frac{\langle {u}_{0},{\varphi }_{n}\rangle B{\varphi }_{n}(\omega (s,{\lambda }_{n})-{c}^{2}{\lambda }_{n})}{s\omega (s,{\lambda }_{n})}\right\vert \hfill \\ \hfill & =\left\vert \sum\limits _{n=1}^{\infty }\frac{\langle {u}_{0},{\varphi }_{n}\rangle B{\varphi }_{n}{c}^{2}\left({\lambda }_{n}({s}^{2}+{\omega }_{d}(s)+{\omega }_{b}(s,{\lambda }_{1}))-{\lambda }_{1}({s}^{2}+{\omega }_{d}(s)+{\omega }_{b}(s,{\lambda }_{n}))\right)}{s\omega (s,{\lambda }_{n})\omega (s,{\lambda }_{1})}\right\vert \hfill \\ \hfill & \leqslant \sum\limits _{n=1}^{\infty }\frac{\vert \langle {u}_{0},{\varphi }_{n}\rangle B{\varphi }_{n}\vert {c}^{2}\left(({\lambda }_{n}-{\lambda }_{1})({s}^{2}+{\omega }_{d}(s))+{\lambda }_{n}{\omega }_{b}(s,{\lambda }_{1})-{\lambda }_{1}{\omega }_{b}(s,{\lambda }_{n})\right)}{s\omega (s,{\lambda }_{n})\omega (s,{\lambda }_{1})}\hfill \\ \hfill & \leqslant \underset{n\in \mathbb{N}}{\mathrm{sup}}\left\{{\lambda }_{n}^{-\nu }\vert B{\varphi }_{n}\vert \right\}{\Vert}{u}_{0}{{\Vert}}_{{\dot{H}}^{2\mu +2\nu }({\Omega})}{\left(\sum\limits _{n=1}^{\infty }{\left(\frac{({s}^{2}+\bar{d}{s}^{2+{\gamma }_{J}}+\bar{b}{s}^{{\alpha }_{N}}){c}^{2}{\lambda }_{n}}{{\lambda }_{n}^{\mu }{s}^{3}({s}^{2}+{c}^{2}{\lambda }_{n})}\right)}^{2}\right)}^{1/2}\hfill \end{aligned}\end{equation*}$$

for $\bar{d}={\sum }_{j=1}^{J}{\,d}_{j}$, $\bar{b}={\sum }_{k=1}^{N}\,{b}_{k}\,\mathrm{max}\left\{1,{\lambda }_{1}\right\}$, where we have used (λ_n − λ₁)(s² + ω_d (s)) ⩾ 0 and ${\lambda }_{n}{\omega }_{b}(s,{\lambda }_{1})-{\lambda }_{1}{\omega }_{b}(s,{\lambda }_{n})={\sum }_{k=1}^{N}\,{s}^{{\alpha }_{k}}{b}_{k}\left({\lambda }_{n}^{{\beta }_{k}}{\lambda }_{1}^{{\beta }_{k}}({\lambda }_{n}^{1-{\beta }_{k}}-{\lambda }_{1}^{1-{\beta }_{k}})\right)\geqslant 0$.

To achieve an $O({s}^{-(2+\mathrm{max}\left\{{\gamma }_{J},{\tilde{\gamma }}_{\tilde{J}}\right\}+\varepsilon )})$ estimate, instead of (48) we thus assume

$$\begin{align}\hfill & \mu > \frac{d}{4}+1,\qquad \mathrm{max}\left\{{\gamma }_{J},{\tilde{\gamma }}_{\tilde{J}}\right\}< \frac{1}{2}\qquad \underset{n\in \mathbb{N}}{\mathrm{sup}}\left\{{\lambda }_{n}^{-\nu }\vert B{\varphi }_{n}\vert \right\},\qquad \underset{n\in \mathbb{N}}{\mathrm{sup}}\left\{{\tilde{\lambda }}_{n}^{-\nu }\vert B{\tilde{\varphi }}_{n}\vert \right\}< \infty .\hfill \end{align} \tag{ 50 }$$

to recover theorem 3.1 with f, $\tilde{f}$ replaced by u₀, ${\tilde{u}}_{0}$.

Remark 3.2. Note that in view of the fact that $\mathcal{A}$, ${\mathcal{A}}_{c}$ do not need to be known (up to the fact that a single eigenvalue is supposed to coincide), we are able to identify all relevant information on the damping model in an unknown medium.

Assume that we know at least some of the eigenfunctions and can use them as excitations ${f}_{i}={\varphi }_{{n}_{i}}$, ${\tilde{f}}_{i}={\tilde{\varphi }}_{{\tilde{n}}_{i}}$, i = 1, ..., I, so that (46) becomes

$$\begin{equation}\begin{aligned}\hfill & \frac{\langle {f}_{i},{\varphi }_{{n}_{i}}\rangle {B}_{i}{\varphi }_{{n}_{i}}}{{s}^{2}+{c}^{2}{\lambda }_{{n}_{i}}+{\sum }_{j=1}^{J}{\,d}_{j}\hspace{0.5pt}{s}^{{\gamma }_{j}+2}+{\sum }_{k=1}^{N}\,{s}^{{\alpha }_{k}}{b}_{k}{({\lambda }_{{n}_{i}})}^{{\beta }_{k}}}\hfill \\ \hfill & \qquad =\frac{\langle {\tilde{f}}_{i},{\tilde{\varphi }}_{{\tilde{n}}_{i}}\rangle {\tilde{B}}_{i}{\tilde{\varphi }}_{{\tilde{n}}_{i}}}{{s}^{2}+{\tilde{c}}^{2}{\tilde{\lambda }}_{{\tilde{n}}_{i}}+{\sum }_{j=1}^{\tilde{J}}{\tilde{\,d}}_{j}\hspace{0.5pt}{s}^{{\tilde{\gamma }}_{j}+2}+{\sum }_{k=1}^{\tilde{N}}\,{s}^{{\tilde{\alpha }}_{k}}{\tilde{b}}_{k}{({\tilde{\lambda }}_{{\tilde{n}}_{i}})}^{{\tilde{\beta }}_{k}}},\quad i=1,\dots ,I,\enspace s\in {\mathbb{C}}_{\theta }.\hfill \end{aligned}\end{equation} \tag{ 51 }$$

We can use the rational function or induction argument as above to compare powers of s and conclude the following:

$$\begin{equation*}\begin{aligned}\hfill & (\mathrm{I})\qquad \hfill & \hfill & N=\tilde{N},\qquad J=\tilde{J}\hfill \\ \hfill & (\mathrm{I}\mathrm{I})\hfill & \hfill & {c}^{2}{\lambda }_{{n}_{i}}={\tilde{c}}^{2}{\tilde{\lambda }}_{{\tilde{n}}_{i}},\quad \hfill & \hfill & i=1,\dots ,I\hfill \\ \hfill & (\mathrm{I}\mathrm{I}\mathrm{I})\hfill & \hfill & {\alpha }_{k}={\tilde{\alpha }}_{k},\ {\gamma }_{j}={\tilde{\gamma }}_{j},\hfill & \hfill & k=1,\dots ,N,\ j=1,\dots ,J\hfill \\ \hfill & (\mathrm{I}\mathrm{V})\hfill & \hfill & {b}_{k}{({\lambda }_{{n}_{i}})}^{{\beta }_{k}}={\tilde{b}}_{k}{({\tilde{\lambda }}_{{\tilde{n}}_{i}})}^{{\tilde{\beta }}_{k}},\ {d}_{j}={\tilde{d}}_{j},\quad \hfill & \hfill & i=1,\dots ,I,\ k=1,\dots ,N,\ j=1,\dots ,J.\hfill \end{aligned}\end{equation*}$$

To extract b_k and β_k , which can be done separately for each k, (skipping the subscript k) we use I = 2 excitations, assume that the corresponding eigenvalues of $\mathcal{A}$ and $\tilde{\mathcal{A}}$ are known and equal and define

$$\begin{equation*}F(b,\beta )={(b{({\lambda }_{{n}_{1}})}^{\beta },b{({\lambda }_{{n}_{2}})}^{\beta })}^{\mathrm{T}}.\end{equation*}$$

One easily sees that the 2 × 2 matrix F'(b, β) is regular for ${\lambda }_{{n}_{1}}\ne {\lambda }_{{n}_{2}}$ and b ≠ 0 and thus by the inverse function theorem F is injective. Thus from (I)–(IV) we obtain all the constants in (45).

Theorem 3.2. Let ${f}_{1}={\varphi }_{{n}_{1}}$, ${f}_{2}={\varphi }_{{n}_{2}}$, ${\tilde{f}}_{1}={\tilde{\varphi }}_{{\tilde{n}}_{1}}$, ${\tilde{f}}_{2}={\tilde{\varphi }}_{{\tilde{n}}_{2}}$ for some n₁ ≠ n₂, ${\tilde{n}}_{1}\ne {\tilde{n}}_{2}\in \mathbb{N}$. Then

$$\begin{equation*}{B}_{i}{u}_{i}(t)={\tilde{B}}_{i}{\tilde{u}}_{i}(t),\quad t\in (0,T),\enspace i=1,2\end{equation*}$$

for the solutions u_i , ${\tilde{u}}_{i}$ of (5) with vanishing initial data, known and equal ${\lambda }_{i}={\tilde{\lambda }}_{i}$, i = 1, 2 and possibly unknown (rest of) $\mathcal{A}$, $\tilde{\mathcal{A}}$ implies

$$\begin{align*}\hfill & N=\tilde{N},\quad J=\tilde{J},\quad c=\tilde{c},\quad {b}_{k}={\tilde{b}}_{k},\quad {\alpha }_{k}={\tilde{\alpha }}_{k},\quad {\beta }_{k}={\tilde{\beta }}_{k},\quad {d}_{j}={\tilde{d}}_{j},\quad {\gamma }_{j}={\tilde{\gamma }}_{j},\hfill \\ \hfill & \qquad \qquad k\in \left\{1,\dots ,N\right\},\enspace j\in \left\{1,\dots ,J\right\}.\hfill \end{align*}$$

Remark 3.3. Scaling with unity in the excitation has been chosen here just for simplicity of exposition. Clearly also multiples of eigenfunctions can be used.

Again, uniqueness can be shown analogously with eigenfunction excitation by initial data rather than by sources.

We briefly consider the additional recovery of a spatially varying wave speed c and for this purpose look at the ch model (9) as an example. With ${\mathcal{A}}_{c}=-c{(x)}^{2}{\triangle}$, (I)–(IV) (with N = 1, c = 1, β = 1, b_k = b) yields

$$\begin{equation*}\alpha =\tilde{\alpha },\qquad {\lambda }_{{n}_{i}}={\tilde{\lambda }}_{{\tilde{n}}_{i}},\qquad b{\lambda }_{{n}_{i}}=\tilde{b}{\tilde{\lambda }}_{{\tilde{n}}_{i}},\quad i=1,\dots ,I.\end{equation*}$$

Using inverse Sturm–Liouville theory, one can recover the spatially varying coefficient c(x) in one space dimension from knowledge of all eigenvalues (actually just two sets of eigenvalues corresponding to different boundary conditions), see, [6, 29]. However, in order to do so via this single mode excitation approach, we would need infinitely many measurements $\left\{{n}_{i}:i=1,\dots ,I\right\}=\mathbb{N}$.

The goal of obtaining knowledge on the eigenvalues can be much better achieved by going back to (46) and considering equality of poles (and residues), which according to [21, section 4] gives us λ_n and $b{\lambda }_{n}^{\beta }$ (note that we have set c = 1 since c(x) is contained in ${\mathcal{A}}_{c}$) from a single excitation u₁ with nonvanishing Fourier coefficients. We will now consider this approach for the more general setting (3) in section 3.3.

In this section we focus on the second order in time case (3) and achieve separation of summands in the eigenfunction expansion by using poles (and residues) of $\hat{h}(s)$. For this purpose we have to prove that each summand λ_n gives rise to at least one pole (in case of ch we know there are two complex conjugates) and that these poles differ for different λ_n . This is also known for ch and fz, see [21, section 4]. Indeed the results from [21, lemma 4.1, remark 4.1] on ch carry over to the model (3) in the setting

$$\begin{equation}{\beta }_{1}=\cdots ={\beta }_{k}=\beta ,\qquad {b}_{1},\dots ,{b}_{k}\geqslant 0\end{equation} \tag{ 52 }$$

that is, (3) with a single β and dissipative behaviour (as is natural to demand from a damping model). Existence of at least one pole for each λ_n can even be shown for the most general model (7).

Lemma 3.1. For each λ_n there exists at least one root of ${\omega }_{n}^{\text{mg}}(s)={s}^{2}+{c}^{2}{\lambda }_{n}+{\sum }_{k=1}^{N}{s}^{{\alpha }_{k}}{\sum }_{\ell =1}^{{M}_{k}}{c}_{k\ell }{\lambda }_{n}^{{\beta }_{k\ell }}$, that is, a pole of the Laplace transformed relaxation solution ${\hat{w}}_{fn}^{\text{mg}}$, ${\hat{w}}_{0n}^{\text{mg}}$, or ${\hat{w}}_{1n}^{\text{mg}}$. Moreover, the poles are single.

In case of (52), the poles of ${\hat{w}}_{fn}^{\text{g}}$, ${\hat{w}}_{0n}^{\text{g}}$, ${\hat{w}}_{1n}^{\text{g}}$, which are the roots of ${\omega }_{n}^{\text{g}}(s)={s}^{2}+{c}^{2}{\lambda }_{n}+{\sum }_{k=1}^{N}{b}_{k}{s}^{{\alpha }_{k}}{\lambda }_{n}^{\beta }$, lie in the left half plane and differ for different λ_n .

Proof. Let f_n (z) = z² + c² λ_n , ${g}_{n}(z)={\sum }_{k=1}^{N}{z}^{{\alpha }_{k}}{\sum }_{\ell =1}^{{M}_{k}}{c}_{k\ell }{\lambda }_{n}^{{\beta }_{k\ell }}$. Let C_R be the circle radius R, centre at the origin. Then, due to α_k ⩽ 1, for a sufficiently large R > c² λ_n , the estimate |g_n (z)| < |f_n (z)| holds on C_R and so Rouché's theorem shows that f_n (z) and ${\omega }_{n}^{\text{mg}}(z)=({f}_{n}+{g}_{n})(z)$ have the same number of roots, counted with multiplicity, within C_R . For f_n these are only at $z=\pm i\sqrt{{\lambda }_{n}}c$ and so ${\omega }_{n}^{\text{mg}}$ has precisely one single root in the third and in the fourth quadrant, respectively.

The fact that the poles of ${\hat{w}}_{fn}^{\text{g}}$, ${\hat{w}}_{0n}^{\text{g}}$, ${\hat{w}}_{1n}^{\text{g}}$ lie in the left half plane in case of (52) follows by a standard energy argument. Suppose now that ω^g has a root at p = r e^iθ , where π/2 ⩽ θ < π, for both λ_n and λ_m . Then for p = r e^iθ , subtracting the equations ${\omega }_{n}^{\text{g}}(p)=0$ and ${\omega }_{m}^{\text{g}}(p)=0$ we obtain

$$\begin{equation}\frac{{c}^{2}({\lambda }_{n}-{\lambda }_{m})}{{\lambda }_{n}^{\beta }-{\lambda }_{m}^{\beta }}=-\sum\limits _{k=1}^{N}{b}_{k}{p}^{{\alpha }_{k}}=\frac{{p}^{2}+{c}^{2}{\lambda }_{n}}{{\lambda }_{n}^{\beta }}\end{equation} \tag{ 53 }$$

thus

$$\begin{equation*}\frac{{p}^{2}}{{c}^{2}}=\frac{{\lambda }_{n}^{\beta }({\lambda }_{n}-{\lambda }_{m})}{{\lambda }_{n}^{\beta }-{\lambda }_{m}^{\beta }}-{\lambda }_{n}=\frac{{\lambda }_{n}^{\beta }{\lambda }_{m}^{\beta }({\lambda }_{n}^{1-\beta }-{\lambda }_{m}^{1-\beta })}{{\lambda }_{n}^{\beta }-{\lambda }_{m}^{\beta }}.\end{equation*}$$

Now if λ_n ≠ λ_m then the right-hand side is positive and real and so 2θ = π. This means that ${\sum }_{k=1}^{N}{b}_{k}{s}^{{\alpha }_{k}}={\sum }_{k=1}^{N}{b}_{k}{r}^{{\alpha }_{k}\pi /2}$ (where all b_k have the same sign) has nonvanishing imaginary part, a contradiction to the fact that the left-hand side of (53) is real.□

Since we will separate the summands by taking residues on both sides of the identity (46), the following identity is also crucial.

Lemma 3.2. For each λ_n , $n\in \mathbb{N}$, the residues of ${\hat{w}}_{fn}^{\text{g}}$, ${\hat{w}}_{0n}^{\text{g}}$, ${\hat{w}}_{1n}^{\text{g}}$ do not vanish.

Proof. By l'Hospital's rule and due to the fact that the poles p_n are single, we get for $\omega (s)={\omega }_{fn}^{\text{mg}}(s)={\omega }_{1n}^{\text{mg}}(s)$

$$\begin{equation*}\text{Res}({\hat{w}}_{n}^{\text{mg}};{p}_{n})=\underset{s\to {p}_{n}}{\mathrm{lim}}\,\frac{(s-{p}_{n})}{\omega (s)}=\underset{s\to {p}_{n}}{\mathrm{lim}}\,\frac{1}{{\omega }^{\prime }(s)}\ne 0.\end{equation*}$$

Moreover, with $\omega (s)={\omega }_{0n}^{\text{mg}}(s)$, $d(s)=\frac{\omega (s)-c{\lambda }^{2}}{s}$,

$$\begin{equation*}\begin{aligned}\hfill \text{Res}({\hat{w}}_{n}^{\text{mg}};{p}_{n})& =\underset{s\to {p}_{n}}{\mathrm{lim}}\,\frac{(s-{p}_{n})d(s)}{\omega (s)}=\underset{s\to {p}_{n}}{\mathrm{lim}}\,\frac{d(s)+(s-{p}_{n}){d}^{\prime }(s)}{{\omega }^{\prime }(s)}=\frac{d({p}_{n})}{{\omega }^{\prime }({p}_{n})}\hfill \\ \hfill & =\frac{\omega ({p}_{n})-c{\lambda }_{n}^{2}}{{p}_{n}{\omega }^{\prime }({p}_{n})}=-\frac{c{\lambda }_{n}^{2}}{{p}_{n}{\omega }^{\prime }({p}_{n})}\ne 0.\hfill \end{aligned}\end{equation*}$$

□

As a consequence of lemma 3.1 and 3.2, by taking residues at the poles p_n of the Laplace transformed time trace data $\hat{h}(s)$, we obtain from (46)

$$\begin{equation}\begin{aligned}\hfill \text{Res}(\hat{h};{p}_{n})& =\frac{1}{2{p}_{n}+{\sum }_{k=1}^{N}{\alpha }_{k}{b}_{k}{p}_{n}^{{\alpha }_{k}-1}\,{({\lambda }_{n})}^{{\beta }_{k}}}\sum\limits _{k\in {K}^{{\lambda }_{n}}}\langle f,{\varphi }_{n,k}\rangle B{\varphi }_{n,k}\hfill \\ \hfill & =\frac{1}{2{p}_{n}+{\sum }_{k=1}^{\tilde{N}}{\tilde{\alpha }}_{k}{\tilde{b}}_{k}{p}_{n}^{{\tilde{\alpha }}_{k}-1}{({\tilde{\lambda }}_{n})}^{{\tilde{\beta }}_{k}}}\sum\limits _{k\in {K}^{{\tilde{\lambda }}_{n}}}\langle \tilde{f},{\tilde{\varphi }}_{n,k}\rangle \tilde{B}{\tilde{\varphi }}_{n,k},\hfill \\ \hfill & \hspace{-3pt}i=1,\dots ,I,\enspace n\in \mathbb{N}\hfill \end{aligned}\end{equation} \tag{ 54 }$$

where the cardinality of the index set ${K}^{{\lambda }_{n}}$ is equal to the multiplicity of the eigenvalue λ_n .

In order to achieve nonzero numerators in (46), (51) and (54), we assume

$$\begin{equation}B{\varphi }_{n}\ne 0,\qquad \tilde{B}{\tilde{\varphi }}_{n}\ne 0\quad \;\text{for}\;\text{all}\;n\in \mathbb{N},\end{equation} \tag{ 55 }$$

and

$$\begin{equation}\langle f,{\varphi }_{n}\rangle \ne 0,\qquad \langle \tilde{f},{\tilde{\varphi }}_{n}\rangle \ne 0\quad \;\text{for}\;\text{all}\;n\in \mathbb{N},\end{equation} \tag{ 56 }$$

see remark 3.4 below.

Assume that,

$$\begin{equation}\sum\limits _{k\in {K}^{{\lambda }_{n}}}\langle f,{\varphi }_{n,k}\rangle B{\varphi }_{n,k}=\sum\limits _{k\in {K}^{{\tilde{\lambda }}_{n}}}\langle \tilde{f},{\tilde{\varphi }}_{n,k}\rangle \tilde{B}{\tilde{\varphi }}_{n,k}\ne 0,\quad \text{for}\,\text{all}\;n\in \mathbb{N}\end{equation} \tag{ 57 }$$

(e.g. by assuming $f=\tilde{f}$, $B=\tilde{B}$, $\mathcal{A}=\tilde{\mathcal{A}}$). Then from the reciprocals of (54) we get, after division by ${\sum }_{k\in {K}^{{\lambda }_{n}}}\langle f,{\varphi }_{n,k}\rangle B{\varphi }_{n,k}$ and subtracting 2p_n on both sides

$$\begin{equation}\sum\limits _{k=1}^{N}{\alpha }_{k}{b}_{k}{p}_{n}^{{\alpha }_{k}-1}{({\lambda }_{n})}^{{\beta }_{k}}=\sum\limits _{k=1}^{\tilde{N}}{\tilde{\alpha }}_{k}{\tilde{b}}_{k}{p}_{n}^{{\tilde{\alpha }}_{k}-1}{(\tilde{{\lambda }_{n}})}^{{\tilde{\beta }}_{k}}.\end{equation} \tag{ 58 }$$

In case $\beta ={\beta }_{k}={\tilde{\beta }}_{k}=\tilde{\beta }$, ${\lambda }_{n}={\tilde{\lambda }}_{n}$, this can be divided by ${({\lambda }_{n})}^{\beta }$ and simplifies to

$$\begin{equation*}\sum\limits _{k=1}^{N}{\alpha }_{k}{b}_{k}{p}_{n}^{{\alpha }_{k}-1}=\sum\limits _{k=1}^{\tilde{N}}{\tilde{\alpha }}_{k}{\tilde{b}}_{k}{p}_{n}^{{\tilde{\alpha }}_{k}-1}\end{equation*}$$

that is, in case of rational powers (49) an equality between two polynomials of $z={s}^{1/\bar{q}}$ at infinitely many different (according to the second part of lemma 3.1) points ${p}_{n}^{1/\bar{q}}$. This yields

$$\begin{equation*}N=\tilde{N},\qquad {\alpha }_{k}={\tilde{\alpha }}_{k},\qquad {b}_{k}={\tilde{b}}_{k},\quad k=1,\dots ,N.\end{equation*}$$

Thus we recover the results of theorem 3.1 under essentially more restrictive assumptions—however we can drop the smoothness assumption on the excitation f, so that f = δ is admissible. Thus the use of poles and residues has no clear advantage over the smooth excitation approach as long as we only aim at recovering the damping term constants. However, it allows to gain additional information on spatially varying coefficients c = c(x) or f = f(x) or u₀ = u₀(x).

3.3.1. Additional recovery of c(x)

First of all, just using knowledge of the pole locations of the data, lemma 3.1 allows to obtain the following corollary of theorem 3.1(b) on uniqueness of the numbers N, J, coefficients b_k , orders α_k , as well as the function c(x) in one space dimension. Note that the assumption f = δ is not compatible with the smoothness assumptions from theorem 3.1, thus we cannot apply the strategy from [21] of obtaining, besides the eigenvalues, also the values φ(x₀) and therewith the norming constants of the eigenfunctions. We thus apply a different result from Sturm–Liouville theory, according to which a pair of eigenfunction sequences corresponding to two different boundary impedances suffices to recover the potential q of the differential operator −△ + q, see [6, 30]. Using the transformation of variables

$$\begin{align*}\hfill \xi (x)& ={\int }_{0}^{x}c{(s)}^{-1}\mathrm{d}s,\qquad \eta (\xi )={c}^{-1/2}(x)y(x),\quad q(\xi )=\frac{1}{4}({c}^{\prime }{(x)}^{2}-2c(x){c}^{\prime \prime }(x)),\hfill \end{align*}$$

so that

$$\begin{align*}\hfill -{c}^{2}(x){y}_{xx}(x)& =\lambda y(x),\quad x\in (0,L){\Leftrightarrow}-{\eta }_{\xi \xi }(\xi )+q(\xi )\eta (\xi )=\lambda \eta (\xi ),\quad \xi \in \left(0,{\int }_{0}^{L}c{(s)}^{-1}\mathrm{d}s\right)\hfill \end{align*}$$

this transfers to the recovery of the sound speed c(x).

Since these results hold in one space dimension only, in this subsection we consider

$$\begin{equation}\begin{aligned}\hfill & {u}_{tt}-c(x){\triangle}u\sum\limits _{k=1}^{N}{b}_{k}{\partial }_{t}^{{\alpha }_{k}}{(-c(x){\triangle}u)}^{{\beta }_{k}}u=\sigma (t)f\quad \hfill & \hfill & x\in (0,L),\ t\in (0,T)\hfill \\ \hfill & {u}^{\prime }(0,t)-{H}_{0}u(0,t)=0,\quad {u}^{\prime }(L,t)+{H}_{L}u(L,t)=0\quad \hfill & \hfill & t\in (0,T)\hfill \\ \hfill & u(x,0)={u}_{0}(x),\quad {u}_{t}(x,0)={u}_{1}(x)\quad {u}_{tt}(x,0)={u}_{2}(x)\quad \hfill & \hfill & x\in (0,L)\hfill \end{aligned}\end{equation} \tag{ 59 }$$

for some L > 0 and some nonnegative impedance constants H₀, H_L . For two different right hand boundary impedance values H_L,i , i = 1, 2 we denote by uⁱ the solution of (59) with H_L = H_L,i and by ${({\varphi }_{n,i},{\lambda }_{n,i})}_{n\in \mathbb{N}}$ the eigenpairs of −c(x)△ with boundary conditions ${\varphi }_{n,i}^{\prime }(0)-{H}_{0}{\varphi }_{n,i}(0)=0$, ${\varphi }_{n,i}^{\prime }(L)+{H}_{L,i}{\varphi }_{n,i}(L)=0$, and weighted L² normalization ${\int }_{0}^{L}\frac{1}{c{(x)}^{2}}{\varphi }_{n,i}{(x)}^{2}\mathrm{d}x=1$, i = 1, 2.

For constructing sufficiently smooth excitations with nonvanishing coefficients with respect to the eigenfunction basis, we point to remark 3.4 below.

Moreover, we assume that the observations satisfy

$$\begin{equation}{B}_{i}{\varphi }_{n,i}\ne 0,\qquad {\tilde{B}}_{i}{\tilde{\varphi }}_{n,i}\ne 0\quad \text{for}\,\text{all}\;n\in \mathbb{N},\enspace i=1,2,\end{equation} \tag{ 60 }$$

which is, e.g., the case with the choice B_i v := v(x₀), x₀ ∈ {0, L}, see remark 3.4 below.

Analogous notation is used for (59) with c, b_k , α_k , β_k , f, H₀, H_L replaced by their corresponding tilde counterparts.

From theorem 3.1 and uniqueness of the eigenvalues from the poles we obtain the following result on combined identification of the coefficients in the damping model and the space-dependent sound speed.

Corollary 3.1. Let $f,\tilde{f}\in {\dot{H}}^{2\mu +2\nu }(0,L)$ for μ, ν satisfying (48) and let (56) and (60) hold.

Then

$$\begin{equation*}{B}_{i}{u}_{i}(t)={\tilde{B}}_{i}{\tilde{u}}_{i}(t)\quad t\in (0,T),\enspace i=1,2\end{equation*}$$

for the solutions u_i , ${\tilde{u}}_{i}$ of (59) with vanishing initial data, H_L = H_L,i , ${\tilde{H}}_{L}={\tilde{H}}_{L,i}$, known equal ${\beta }_{k}={\tilde{\beta }}_{k}$, and unknown c, $\tilde{c}\in {C}^{2}(0,L)$, implies

$$\begin{equation*}\begin{aligned}\hfill & N=\tilde{N},\qquad {b}_{k}={\tilde{b}}_{k},\qquad {\alpha }_{k}={\tilde{\alpha }}_{k},\quad k\in \left\{1,\dots ,N\right\},\hfill \\ \hfill & c(x)=\tilde{c}(x),\quad x\in (0,L).\hfill \end{aligned}\end{equation*}$$

The same result can be obtained without smoothness assumption on f provided (57), and $\beta ={\beta }_{k}={\tilde{\beta }}_{k}=\tilde{\beta }$ by a poles/residual argument, cf (58).

Remark 3.4. Some comments on the verification of conditions (56), (60) and (57) are in order, in particular in view of the fact that c(x), $\tilde{c}(x)$ and therefore also the eigenfunctions are unknown. First of all, note that choosing $f,\tilde{f}={\delta }_{{x}_{0}}$ formally satisfies this condition provided we can make sure that none of the eigenfunctions vanish at the point x₀. This is, e.g., the case if x₀ ∈ {0, L} is a boundary point for then in case x₀ = L the impedance conditions would imply ${\varphi }_{n}^{\prime }(0)-{H}_{0}{\varphi }_{n}(0)=0$, ${\varphi }_{n}^{\prime }(L)=0$, φ_n (L) = 0, which by uniqueness of solutions to the second order ODE $-c(x){\varphi }_{n}^{\prime \prime }(x)-{\lambda }_{n}{\varphi }_{n}(x)=0$ with these boundary conditions would imply φ_n ≡ 0, a contradiction to φ_n being an eigenfunction. Likewise for x₀ = 0.

Note that setting B_i v := v(x₀), x₀ ∈ {0, L}, the same argument yields (60) as well as (57).

However, clearly $f={\delta }_{{x}_{0}}$ does not satisfy the regularity requirements from the first part of the corollary, thus we also indicate a possibility of constructing a smooth f that satisfies (56). For any s > 0, ɛ > 0, the function $f(x)=c{(x)}^{2}{\sum }_{n=1}^{\infty }{\lambda }_{n}^{-s}{n}^{-(\frac{1}{2}+\varepsilon )}{\varphi }_{n}(x)$ lies in ${\dot{H}}^{s}$ and satisfies ${\langle f,{\varphi }_{k}\rangle }_{{L}_{1/{c}^{2}}^{2}}\geqslant {\lambda }_{n}^{-s}{n}^{-(\frac{1}{2}+\varepsilon )}$ for all $k\in \mathbb{N}$. Unfortunately this is only a result on existence of f, since the reconstruction requires knowledge of c.

3.3.2. Additional recovery of f(x) (or u₀(x) or u₁(x))

In order to reconstruct the possibly unknown excitation f (analogously for u₁, u₀, the latter relevant in PAT) we use measurements not only at finitely many points but on a surface Σ and make the linear independence assumption

$$\begin{equation}\left(\sum\limits _{k\in {K}^{{\lambda }_{n}}}{m}_{k}{\varphi }_{k}(x)=0\enspace \text{for}\;\text{all}\;x\in {\Sigma}\right)\Longrightarrow\left({m}_{k}=0\enspace \text{for}\;\text{all}\;k\in {K}^{{\lambda }_{n}}\right)\quad n\in \mathbb{N}.\end{equation} \tag{ 61 }$$

Moreover, for the purpose of recovering f (or u₁, or u₀) we assume c(x) to be known and therefore ${\varphi }_{n}={\tilde{\varphi }}_{n}$, ${\lambda }_{n}={\tilde{\lambda }}_{n}$. Thus in place of (54) we have

$$\begin{equation}\begin{aligned}\hfill & \frac{1}{2{p}_{n}+{\sum }_{k=1}^{N}{\alpha }_{k}{b}_{k}{p}_{n}^{{\alpha }_{k}-1}{({\lambda }_{n})}^{\beta }}\sum\limits _{k\in {K}^{{\lambda }_{n}}}\langle f,{\varphi }_{n,k}\rangle {\varphi }_{n,k}(x)\hfill \\ \hfill & \qquad =\frac{1}{2{p}_{n}+{\sum }_{k=1}^{\tilde{N}}{\tilde{\alpha }}_{k}{\tilde{b}}_{k}{p}_{n}^{{\tilde{\alpha }}_{k}-1}{({\lambda }_{n})}^{\tilde{\beta }}}\sum\limits _{k\in {K}^{{\tilde{\lambda }}_{n}}}\langle \tilde{f},{\varphi }_{n,k}\rangle {\varphi }_{n,k}(x),\quad x\in {\Sigma},\enspace n\in \mathbb{N}.\hfill \end{aligned}\end{equation} \tag{ 62 }$$

Again, we can obtain a combined identification result from theorem 3.1, uniqueness of the eigenvalues from the poles, and (62).

Corollary 3.2. Let $f,\tilde{f}\in {\dot{H}}^{2\mu +2\nu }({\Omega})$ for μ, ν satisfying (48) and let (62) hold.

Then

$$\begin{equation*}u(x,t)=\tilde{u}(x,t)\quad x\in {\Sigma},\ t\in (0,T),\end{equation*}$$

for the solutions u, $\tilde{u}$ of (6) with vanishing initial data known equal ${\beta }_{k}={\tilde{\beta }}_{k}$, and unknown f, $\tilde{f}\in {L}^{2}({\Omega})$, implies

$$\begin{equation*}\begin{aligned}\hfill N& =\tilde{N},\qquad {b}_{k}={\tilde{b}}_{k},\qquad {\alpha }_{k}={\tilde{\alpha }}_{k},\quad k\in \left\{1,\dots ,N\right\},\hfill \\ \hfill f(x)& =\tilde{f}(x),\quad x\in {\Omega}.\hfill \end{aligned}\end{equation*}$$

Remark 3.5. The same result can be obtained without smoothness assumption on f provided (57), and $\beta ={\beta }_{k}={\tilde{\beta }}_{k}=\tilde{\beta }$ by a poles/residual argument, cf (58).

Remark 3.6. The result carries over almost without change to the setting f = 0, u₀ = 0, where u₁ = u₁(x) ≠ 0 is the unknown to be recovered.

In the setting f = 0, u₁ = 0, of recovering u₀ = u₀(x) ≠ 0 (as relevant in PAT), we can rely on remark 3.1 to conclude that Corollary 3.2 remains valid with f, $\tilde{f}$, (48) replaced by u₀, ${\tilde{u}}_{0}$, (50).

Starting from the equation

$$\begin{align*}\hfill \hat{h}(s)& =\frac{\langle f,\varphi \rangle B\varphi }{{s}^{2}+{c}^{2}\lambda +{\sum }_{k=1}^{N}{s}^{{\alpha }_{k}}{b}_{k}{\lambda }^{{\beta }_{k}}}\begin{cases}1\,\hfill & \quad \text{in}\;\text{case}\;\,{u}_{0}=0,\ \,{u}_{1}=\varphi ,\ \,f=0\hfill \\ \hat{\sigma }(s)\,\hfill & \quad \text{in}\;\text{case}\;\,{u}_{0}=0,\ \,{u}_{1}=0,\ \,f=\varphi \hfill \\ s+\sum\limits _{k=1}^{N}{s}^{{\alpha }_{k}-1}{b}_{k}{\lambda }^{{\beta }_{k}}\,\hfill & \quad \text{in}\;\text{case}\;\,{u}_{0}=\varphi ,\ \,{u}_{1}=0,\ \,f=0\hfill \end{cases}\hfill \end{align*}$$

(cf (51) for the slightly simpler setting of (5)) and noting that taking the reciprocal and rearranging makes the problem somewhat less nonlinear and the derivative easier to compute, we consider

$$\begin{equation*}F({b}_{1},\dots ,{b}_{N},{\beta }_{1},\dots ,{\beta }_{N},{\alpha }_{1},\dots ,{\alpha }_{N};s)=G(s),\quad s\in \left\{{s}_{1},\dots ,{s}_{M}\right\}\end{equation*}$$

with

$$\begin{equation*}F({b}_{1},\dots ,{b}_{N},{\beta }_{1},\dots ,{\beta }_{N},{\alpha }_{1},\dots ,{\alpha }_{N};s)=\sum\limits _{k=1}^{N}{b}_{k}{\lambda }^{{\beta }_{k}}{s}^{{\alpha }_{k}}\end{equation*}$$

and

$$\begin{equation*}G(s)=\begin{cases}\frac{B\varphi }{\hat{h}(s)}-{s}^{2}-{c}^{2}\lambda \quad \hfill & \;\text{in}\;\text{case}\;\,{u}_{0}=0,\ \,{u}_{1}=\varphi ,\ \,f=0\hfill \\ \frac{B\varphi \hat{\sigma }(s)}{\hat{h}(s)}-{s}^{2}-{c}^{2}\lambda \quad \hfill & \;\text{in}\;\text{case}\;\,{u}_{0}=0,\ \,{u}_{1}=0,\ \,f=\varphi \hfill \\ \frac{{c}^{2}\lambda }{1-\frac{B\varphi }{s\hat{h}(s)}}-{s}^{2}-{c}^{2}\lambda \quad \hfill & \;\text{in}\;\text{case}\;\,{u}_{0}=\varphi ,\ \,{u}_{1}=0,\ \,f=0.\hfill \end{cases}\end{equation*}$$

For applying Newton's method, the Jacobian is given by

$$\begin{align}& \frac{\mathrm{\partial }F}{\mathrm{\partial }{b}_{i}}({b}_{1},\dots ,{b}_{N},{\beta }_{1},\dots ,{\beta }_{N},{\alpha }_{1},\dots ,{\alpha }_{N};s):={\lambda }^{{\beta }_{i}}{s}^{{\alpha }_{i}}\\\ & \frac{\mathrm{\partial }F}{\mathrm{\partial }{\beta }_{i}}({b}_{1},\dots ,{b}_{N},{\beta }_{1},\dots ,{\beta }_{N},{\alpha }_{1},\dots ,{\alpha }_{N};s):={b}_{i}\;\mathrm{l}\mathrm{n}(\lambda )\;{\lambda }^{{\beta }_{i}}{s}^{{\alpha }_{i}}\\ & \frac{\mathrm{\partial }F}{\mathrm{\partial }{\alpha }_{i}}({b}_{1},\dots ,{b}_{N},{\beta }_{1},\dots ,{\beta }_{N},{\alpha }_{1},\dots ,{\alpha }_{N};s):={b}_{i}{\lambda }^{{\beta }_{i}}\;\mathrm{l}\mathrm{n}(s)\;{s}^{{\alpha }_{i}}.\end{align}$$

Taking the logarithm makes the equation less nonlinear with respect to α and β by considering (in case N = 1)

$$\begin{equation*}\begin{aligned}\hfill {F}_{1}(b,\beta ,\alpha )& {:=}b{\lambda }^{\beta }{s}_{1}^{\alpha }=\frac{B\varphi \hat{\sigma }({s}_{1})}{\hat{h}({s}_{1})}-{s}_{1}^{2}-{c}^{2}\lambda \hfill \\ \hfill {F}_{m}(b,\beta ,\alpha )& {:=}\mathrm{log}\,b+\beta \,\mathrm{log}\,\lambda +\alpha \,\mathrm{log}\,{s}_{m}=\mathrm{log}\left(\frac{B\varphi \hat{\sigma }({s}_{m})}{\hat{h}({s}_{m})}-{s}_{m}^{2}-{c}^{2}\lambda \right),\enspace m=2,\dots ,M.\hfill \end{aligned}\end{equation*}$$

Then F₁ is linear with respect to b and the F_m are linear with respect to α and β. This helps to enlarge the convergence radius of Newton's method, as we experienced in our numerical tests.

In case u₀ = φ, u₁ = 0, f = 0 we have

$$\begin{equation}\hat{h}(s)=\frac{s+{\sum }_{k=1}^{N}{b}_{k}{\lambda }^{{\beta }_{k}}{s}^{{\alpha }_{k}-1}}{{s}^{2}+{\sum }_{k=1}^{N}{b}_{k}{\lambda }^{{\beta }_{k}}{s}^{{\alpha }_{k}}+{c}^{2}\lambda }B\varphi \sim \frac{B\varphi }{{c}^{2}}\sum\limits _{k=1}^{N}{b}_{k}{\lambda }^{{\beta }_{k}-1}{s}^{{\alpha }_{k}-1}\quad \;\text{as}\;s\to 0.\end{equation} \tag{ 63 }$$

We will heavily make use of a Tauberian theorem (for s → 0, t → ∞, see [12, theorem 2 in chapter XIII.5]) which we here quote for the convenience of the reader.

Theorem 4.1. Let $w:\left[0,\infty \right)\to \mathbb{R}$ be a monotone function and assume its Laplace–Stieltjes transform $\hat{w}(s)={\int }_{0}^{\infty }{\mathrm{e}}^{-st}\,\mathrm{d}w(t)$ exists for all s in the right hand complex plane. Then for ρ ⩾ 0

$$\begin{equation*}\hat{w}(s)\sim \frac{C}{{s}^{\rho }}\quad \text{as}\ s\to 0\end{equation*}$$

if and only if

$$\begin{equation*}w(t)\sim \frac{C}{{\Gamma}(\rho +1)}{t}^{\rho }\quad \text{as}\ t\to \infty .\end{equation*}$$

Note that here dw is actually a measure, which in case of an absolutely continuous function w can be written as dw(t) = w'(t)dt.

Applying theorem 4.1 to (63) we get

$$\begin{equation}h(t)=\frac{B\varphi }{{c}^{2}}\sum\limits _{k=1}^{N}\frac{{b}_{k}{\lambda }^{{\beta }_{k}-1}}{{\Gamma}(1-{\alpha }_{k})}{t}^{-{\alpha }_{k}}+O({t}^{-2{\alpha }_{1}})\quad \text{as}\,t\to \infty .\end{equation} \tag{ 64 }$$

Therefore we get an asymptotic formula (like the one in [17] for the subdiffusion equation) for the smallest order

$$\begin{equation}{\alpha }_{1}=-\underset{t\to \infty }{\mathrm{lim}}\,\frac{\mathrm{log}\,h(t)}{\mathrm{log}\,t}.\end{equation} \tag{ 65 }$$

By l'Hospital's rule and actually thereby removing the constant we can instead also use

$$\begin{equation*}{\alpha }_{1}=-\underset{t\to \infty }{\mathrm{lim}}\,\frac{\frac{\mathrm{d}}{\mathrm{d}t}(\mathrm{log}\,h(t))}{\frac{\mathrm{d}}{\mathrm{d}t}(\mathrm{log}\,t)}=-\underset{t\to \infty }{\mathrm{lim}}\,\frac{t{h}^{\prime }(t)}{h(t)}.\end{equation*}$$

After having determined α₁ this way, we can also compute b₁ as a limit

$$\begin{equation}{b}_{1}=\underset{t\to \infty }{\mathrm{lim}}(h(t){t}^{{\alpha }_{1}}){c}^{2}{\lambda }^{1-{\beta }_{1}}{\Gamma}(1-{\alpha }_{1}).\end{equation} \tag{ 66 }$$

So we can successively (by the above procedure and subtracting the recovered terms one after another) reconstruct those terms k for which α_k < 2α₁, see the remainder term in (64). However, if there are terms with α_k ⩾ 2α₁, they seem to get masked by the $O({t}^{-2{\alpha }_{1}})$ remainder.

The same holds true if we do the excitation by u₀ = 0, u₁ = φ, f = 0, or u₀ = 0, u₁ = 0, f = φ, where (63) changes to $\hat{h}(s)=\frac{1}{{s}^{2}+{\sum }_{k=1}^{N}{b}_{k}{\lambda }^{{\beta }_{k}}{s}^{{\alpha }_{k}}+{c}^{2}\lambda }B\varphi$.

In order to avoid the restriction α_k < 2α₁, we thus refine the expansion of $\hat{h}(s)$ in terms of powers of s and retain the singular ones, that is, those with negative powers, since we are looking at the limiting case s → 0. Using the geometric series formula and the multinomial theorem, with the abbreviations ${\Sigma}={\sum }_{k=1}^{N}{b}_{k}{\lambda }^{{\beta }_{k}}{s}^{{\alpha }_{k}}$, $q=\frac{{s}^{2}+{\Sigma}}{{c}^{2}\lambda }$, yields

$$\begin{equation*}\begin{aligned}\hfill \hat{h}(s)& =\frac{s+{\Sigma}/s}{{s}^{2}+{\Sigma}+{c}^{2}\lambda }=\frac{1}{s}\frac{q}{q+1}=-\frac{1}{s}\sum\limits _{m=1}^{\infty }{(-q)}^{m}=-\frac{1}{s}\sum\limits _{m=1}^{\infty }{(-\frac{1}{{c}^{2}\lambda })}^{m}\sum\limits _{\ell =0}^{m}{s}^{2(m-\ell )}{{\Sigma}}^{\ell }\hfill \\ \hfill & =-\sum\limits _{m=1}^{\infty }{(-\frac{1}{{c}^{2}\lambda })}^{m}\sum\limits _{\ell =0}^{m}{s}^{2(m-\ell )-1}\sum\limits _{{i}_{1}+\cdots +{i}_{N}=\ell }\left(\genfrac{}{}{0pt}{}{\ell }{{i}_{1},\dots ,{i}_{N}}\right)\prod\limits _{j=1}^{N}{b}_{j}^{{i}_{j}}{\lambda }^{\sum\limits _{j=1}^{N}{\beta }_{j}{i}_{j}}{s}^{\sum\limits _{j=1}^{N}{\alpha }_{j}{i}_{j}}.\hfill \end{aligned}\end{equation*}$$

The terms corresponding to ℓ < m are obviously O(s), so to extract the singularities it suffices to consider ℓ = m. The set of indices leading to singular terms is then

$$\begin{align*}\hfill {I}_{m}={I}_{m}({\alpha }_{1},\dots ,{\alpha }_{N})& =\left\{\vec{i}=({i}_{1},\dots ,{i}_{N}):{i}_{1}+\cdots +{i}_{N}=m,\sum\limits _{j=1}^{N}{\alpha }_{j}{i}_{j}< 1\right\},\enspace m\leqslant {m}_{\mathrm{max}}=[\frac{1}{{\alpha }_{1}}]\hfill \end{align*}$$

and the singular part of $\hat{h}$ can be written as

$$\begin{align}\hfill {\hat{h}}_{\text{sing}}(s)& =-\sum\limits _{m=1}^{{m}_{\mathrm{max}}}{(-\frac{1}{{c}^{2}\lambda })}^{m}\sum\limits _{\vec{i}\in I}{\tilde{b}}_{m,\vec{i}}{s}^{\sum\limits _{j=1}^{N}{\alpha }_{j}{i}_{j}-1}\enspace \text{with}\,{\tilde{b}}_{\vec{i}}\,=\left(\genfrac{}{}{0pt}{}{m}{{i}_{1},\dots ,{i}_{N}}\right)\prod\limits _{j=1}^{N}{b}_{j}^{{i}_{j}}{\lambda }^{\sum\limits _{j=1}^{N}{\beta }_{j}{i}_{j}}.\hfill \end{align} \tag{ 67 }$$

In case N = 2 with i = i₁, i₂ = m − i, this reads as

$$\begin{align*}\hfill {\hat{h}}_{\text{sing}}(s)& =-\sum\limits _{m=1}^{{m}_{\mathrm{max}}}{(-\frac{1}{{c}^{2}\lambda })}^{m}\sum\limits _{i={i}_{m,\mathrm{min}}}^{m}{\tilde{b}}_{m,i}{s}^{{\alpha }_{1}i+{\alpha }_{2}(m-i)-1}\text{with}\;{\tilde{b}}_{m,i}=\left(\genfrac{}{}{0pt}{}{m}{i}\right){b}_{1}^{i}{b}_{2}^{m-i}{\lambda }^{{\beta }_{1}i+{\beta }_{2}(m-i)}\hfill \end{align*}$$

since

$$\begin{align*}\hfill {I}_{m}& =\left\{i\in \left\{0,\dots ,m\right\}:{\alpha }_{1}i+{\alpha }_{2}(m-i)< 1\right\}=\left\{{i}_{m,\mathrm{min}},\dots ,m\right\}\;\text{with}\;{i}_{m,\mathrm{min}}=\left[\frac{m{\alpha }_{2}-1}{{\alpha }_{2}-{\alpha }_{1}}\right].\hfill \end{align*}$$

For (67), by theorem 4.1

$$\begin{equation*}{\hat{h}}_{\text{sing}}(t)=-\sum\limits _{m=1}^{{m}_{\mathrm{max}}}{(-\frac{1}{{c}^{2}\lambda })}^{m}\sum\limits _{\vec{i}\in I}\frac{{\tilde{b}}_{m,\vec{i}}}{{\Gamma}\left(1-{\sum }_{j=1}^{N}{\alpha }_{j}{i}_{j}\right)}{t}^{-\sum\limits _{j=1}^{N}{\alpha }_{j}{i}_{j}}\end{equation*}$$

which confirms (65) and (66).

For the first damping term, we also get an asymptotic formula for the smallest order in Laplace domain—the small s counterpart to (65): due to

$$\begin{equation}\mathrm{log}\,\hat{h}(s)\approx \mathrm{log}\,{b}_{1}+({\alpha }_{1}-1)\mathrm{log}\,s-\mathrm{log}({c}^{2}{\lambda }^{1-{\beta }_{1}})\quad \text{as}\,s\to 0\end{equation} \tag{ 68 }$$

we have

$$\begin{equation}1-{\alpha }_{1}=-\underset{s\to 0}{\mathrm{lim}}\,\frac{\mathrm{log}\,\hat{h}(s)}{\mathrm{log}\,s}.\end{equation} \tag{ 69 }$$

One might realise this limit by extrapolation: fitting the values $\frac{\mathrm{log}\,\hat{h}({s}_{m})}{\mathrm{log}\,{s}_{m}}$ at M sample points s₁, ..., s_M to a regression line or low order polynomial r(s) and setting α₁ = 1 + r(0). Next, recover b₁ from

$$\begin{equation*}\mathrm{log}\,{b}_{1}\approx \mathrm{log}\,\hat{h}(s)+\mathrm{log}({c}^{2}{\lambda }^{1-{\beta }_{1}})-({\alpha }_{1}-1)\mathrm{log}\,s.\end{equation*}$$

The formula (68) also shows why recovering b₁ and α₁ simultaneously appears to be so hard: the factor multiplied with log b₁ is unity while the factor multiplied with 1 − α₁ tends to infinity.

In case u₀ = φ, u₁ = 0, f = 0 where the Laplace transformed observations are given by

$$\begin{equation}\hat{h}(s){:=}\frac{s+{\sum }_{k=1}^{N}{b}_{k}{\lambda }^{{\beta }_{k}}{s}^{{\alpha }_{k}-1}}{{s}^{2}+{\sum }_{k=1}^{N}{b}_{k}{\lambda }^{{\beta }_{k}}{s}^{{\alpha }_{k}}+{c}^{2}\lambda }\end{equation} \tag{ 70 }$$

large s/small t asymptotics cannot be exploited for the following reason. As s → ∞ then $\hat{h}(s)\to 0$ and in fact $s\hat{h}(s)\to 1$. More precisely, for

$$\begin{align*}\hfill d(s)& {:=}s\hat{h}(s)-1=\frac{{s}^{2}+{\sum }_{k=1}^{N}{b}_{k}{\lambda }^{{\beta }_{k}}{s}^{{\alpha }_{k}}}{{s}^{2}+{\sum }_{k=1}^{N}{b}_{k}{\lambda }^{{\beta }_{k}}{s}^{{\alpha }_{k}}+{c}^{2}\lambda }-1=-\frac{{c}^{2}\lambda }{{s}^{2}+{\sum }_{k=1}^{N}{b}_{k}{\lambda }^{{\beta }_{k}}{s}^{{\alpha }_{k}}+{c}^{2}\lambda }.\hfill \end{align*}$$

Thus clearly ${s}^{\gamma }[s\hat{h}(s)-1]\to 0$ for any $\gamma \in \left[0,2\right)$.

Not even large λ values (or some combined asymptotics as s → 0, λ_j → ∞) seems to work because then the c² λ term takes over and again masks the terms containing α_k .

However, in case u₀ = 0, u₁ = φ, f = 0 the situation is different since then

$$\begin{equation}\hat{h}(s){:=}\frac{1}{{s}^{2}+{\sum }_{k=1}^{N}{b}_{k}{\lambda }^{{\beta }_{k}}{s}^{{\alpha }_{k}}+{c}^{2}\lambda }\end{equation} \tag{ 71 }$$

and therefore

$$\begin{equation}1-({s}^{2}+{c}^{2}\lambda )\hat{h}(s){:=}\underset{\to 1}{\underbrace{\frac{{s}^{2}}{{s}^{2}+{\sum }_{k=1}^{N}{b}_{k}{\lambda }^{{\beta }_{k}}{s}^{{\alpha }_{k}}+{c}^{2}\lambda }}}\left(\sum\limits _{k=1}^{N}{b}_{k}{\lambda }^{{\beta }_{k}}{s}^{{\alpha }_{k}-2}\right)\end{equation} \tag{ 72 }$$

which via theorem 4.1 (where due to [12, theorem 3 in chapter XIII.5] we can replace s → 0, t → ∞ by s → ∞, t → 0) yields

$$\begin{equation*}-{c}^{2}\lambda h(t)-{h}^{\prime \prime }(t)\sim \sum\limits _{k=1}^{N}\frac{{b}_{k}{\lambda }^{{\beta }_{k}}}{{\Gamma}(2-{\alpha }_{k})}{t}^{1-{\alpha }_{k}}\quad \text{as}\,t\to 0,\end{equation*}$$

where we have used the fact that h(0) = 0, h'(0) = Bφ because of u₀ = 0, u₁ = φ.

This yields the asymptotic formulas

$$\begin{equation*}1-{\alpha }_{N}=\underset{t\to 0}{\mathrm{l}\mathrm{i}\mathrm{m}}\;\frac{\mathrm{l}\mathrm{n}(-{c}^{2}\lambda h(t)-{h}^{\mathrm{\prime }\mathrm{\prime }}(t))}{\mathrm{l}\mathrm{n}(t)}=\underset{t\to 0}{\mathrm{l}\mathrm{i}\mathrm{m}}\vert \frac{t({c}^{2}\lambda h(t)+{h}^{\mathrm{\prime }\mathrm{\prime }}(t))}{{c}^{2}\lambda h(t)+{h}^{\mathrm{\prime }\mathrm{\prime }}(t)}\vert ,\end{equation*}$$

where we have used l'Hospital's rule and

$$\begin{equation*}{b}_{N}=-\underset{t\to 0}{\mathrm{lim}}({c}^{2}\lambda h(t)+{h}^{\prime \prime }(t)){t}^{{\alpha }_{N}-1}\left.\right){\lambda }^{-{\beta }_{N}}{\Gamma}(2-{\alpha }_{N}).\end{equation*}$$

Again, in order to also recover α_N−1, ..., α₁, one may also take into account mixed terms in the expansion of (72) analogously to the large t asymptotics case. In case of two damping terms, this results in

$$\begin{equation}\begin{aligned}\hfill & -{c}^{2}\lambda h(t)-{h}^{\prime \prime }(t)\hfill \\ \hfill & =\frac{{b}_{1}{\lambda }^{{\beta }_{1}}}{{\Gamma}(2-{\alpha }_{1})}{t}^{1-{\alpha }_{1}}+\frac{{b}_{2}{\lambda }^{{\beta }_{2}}}{{\Gamma}(2-{\alpha }_{2})}{t}^{1-{\alpha }_{2}}+\frac{{b}_{1}{\lambda }^{{\beta }_{1}}{c}^{2}\lambda }{{\Gamma}(4-{\alpha }_{1})}{t}^{3-{\alpha }_{1}}+\frac{{b}_{2}{\lambda }^{{\beta }_{2}}{c}^{2}\lambda }{{\Gamma}(4-{\alpha }_{2})}{t}^{3-{\alpha }_{2}}\hfill \\ \hfill & \quad +\frac{{b}_{1}^{2}{\lambda }^{2{\beta }_{1}}}{{\Gamma}(4-2{\alpha }_{1})}{t}^{3-2{\alpha }_{1}}+\frac{{b}_{1}{b}_{2}{\lambda }^{{\beta }_{1}+{\beta }_{2}}}{{\Gamma}(4-{\alpha }_{1}-{\alpha }_{2})}{t}^{3-{\alpha }_{1}-{\alpha }_{2}}+\frac{{b}_{2}^{2}{\lambda }^{2{\beta }_{2}}}{{\Gamma}(4-2{\alpha }_{2})}{t}^{3-2{\alpha }_{2}}\hfill \\ \hfill & \quad +O({t}^{5-3{\alpha }_{2}})\quad \text{as}\ t\to 0.\hfill \end{aligned}\end{equation} \tag{ 73 }$$

The leftmost graph in figure 1 shows the actual time trace of h(t) over the range 0 ⩽ t ⩽ 40 for a damping model with three terms ${b}_{i}{\partial }_{t}^{{\alpha }_{i}}$, i = 1, 2, 3. Note the damped oscillatory behaviour evident and the range shown is the one we are labelling as 'full-time' despite the fact that we will later look at the important case of the very long term behaviour of the solution where the time values are substantially outside of this range. This time range terminology is highly dependent on the values of critical constants in the model. These include the size of the domain (which determines the scaling of the eigenvalues {λ_n }), the wave speed c and of course the strength of the damping given by the coefficients {b_i }. In our reconstructions we will typically take these constants to be of approximate order unity but note the dependence on these quantities and the effect that a substantial change would make.

Zoom In Zoom Out Reset image size

To obtain this h(t) our damping constants are of order $\sim 0.1-0.2$ and the combined term Λ = c² λ is taken to be unity. Changes to the former would modify the time-decrease in h and changes to the latter would alter the frequency of the oscillations.

The rightmost graph in figure 1 shows the logarithm of the Laplace transform $H(s)=\hat{h}(s)$ together with the values of ${\hat{h}}_{i}(s)$ after iteration i of the scheme.

This demonstrates two important facts. First, the very fast convergence of the scheme in the sense of the convergence of the target Laplace transform function as the parameters {α_i , b_i } are resolved. The actual $\hat{h}(s)$ is shown by the solid black line, the initial approximation by the dotted red curve and the first iteration by the dashed green curve which, at this scale is already almost indistinguishable from the actual $\hat{h}$. Second, given this it should be quite feasible to obtain reasonably accurate values the parameters {α_i , b_i }. In addition, assuming we excite the system with an initial condition equal to an eigenfunction corresponding to an actual eigenvalue λ, we will be able to reconstruct the value of c² from the composite Λ = c² λ.

In this situation, to reconstruct the parameters {α_i , b_i , c²} we work directly with the representation $\hat{h}(s)$ or more exactly with its logarithm which is a more convenient form for computing the Jacobian

$$\begin{equation}\mathrm{log}\left(\hat{h}(s)\right)=\mathrm{log}\left(s+\sum\limits _{1}^{n}{b}_{i}{s}^{{\alpha }_{i}-1}\right)-\mathrm{log}\left({s}^{2}+\sum\limits _{1}^{n}{b}_{i}{s}^{{\alpha }_{i}}+{c}^{2}\lambda \right).\end{equation} \tag{ 74 }$$

Of course, to obtain $\hat{h}(s)$ we must approximate the integral $\hat{h}(s)={\int }_{0}^{\infty }{\mathrm{e}}^{-st}\,h(t)\mathrm{d}t$ and this indeed does require full time measurements. However, we do not actually need the values of the analytic function $\hat{h}(s)$ for all s for the Newton scheme and so if we only use s values with s > s₀ then this allows us to ignore very large time measurements and indeed we truncated these to t ⩽ 40.

The reconstruction of the components in the function $\hat{h}(s)$ is shown graphically in figure 2 and in tabular form in table 1. The exact values were α = {0.25, 0.5, 0.75}, b = {0.2, 0.25, 0.1}, Λ = 4. Note that we have included the value of the residual here and keeping track of this value allows the iteration scheme to terminate when saturation has occurred.

Zoom In Zoom Out Reset image size

Here we are trying to simulate the asymptotic values of the constituent powers of t occurring in the data function h(t). This was achieved by using a sample of points between t_min and t_max.

A Newton based approach

We apply Newton's method to recover the constants {p_i , c_k,ℓ } in

$$\begin{equation*}\begin{aligned}\hfill h(t)& ={c}_{1,1}{t}^{-{p}_{1}}+{c}_{1,2}{t}^{-{p}_{2}}+{c}_{1,3}{t}^{-{p}_{3}}+{c}_{2,1}{t}^{-2{p}_{1}}+{c}_{2,2}{t}^{-2{p}_{2}}+{c}_{2,3}{t}^{-2{p}_{3}}\hfill \\ \hfill & \quad +{c}_{2,4}{t}^{-{p}_{1}-{p}_{2}}+{c}_{2,5}{t}^{-{p}_{2}-{p}_{3}}+{c}_{2,6}{t}^{-{p}_{3}-{p}_{1}}+{c}_{3,1}{t}^{-3{p}_{1}}+{c}_{3,2}{t}^{-3{p}_{2}}\hfill \\ \hfill & \quad +{c}_{3,3}{t}^{-2{p}_{1}-{p}_{2}}+{c}_{3,4}{t}^{-{p}_{1}-2{p}_{2}}+{c}_{3,5}{t}^{-2{p}_{1}-{p}_{3}}+{c}_{3,6}{t}^{-{p}_{1}-2{p}_{3}}\hfill \\ \hfill & \quad +{c}_{3,7}{t}^{-2{p}_{2}-{p}_{3}}+{c}_{3,8}{t}^{-{p}_{2}-2{p}_{3}}+O({t}^{-3{p}_{3}})\hfill \end{aligned}\end{equation*}$$

which results from the expansion (67) in case of three terms, and then recover α_i = p_i and b_i = c_1,i Γ(1 − p_i ). In the above we have neglected the term ${t}^{-3{\alpha }_{3}}$ as this would not arise from the Tauberian theorem in the case that the largest power ${\alpha }_{3}\geqslant \frac{1}{3}$. We may also have to exclude other terms such as ${t}^{-2{\alpha }_{1}-{\alpha }_{2}}$ if 2α₁ + α₂ ⩾ 1. In practice, during the iteration process, terms should be included or excluded in the code depending on this criterion: we did so by checking if the argument passed to the Γ function would be negative in which case the term is deleted from use for that iteration step.

Values for the table shown below were t_min = 5 × 10⁴ and t_max = 2 × 10⁵. As a general rule, terms with small α values can be resolved with a smaller value of t_max, but for, say the recovery of a pair of damping terms with α_i > 0.8, a larger value of t_max with commensurate accuracy will be needed.

The case of three damping terms ${\left\{{b}_{i}{\partial }_{t}^{{\alpha }_{i}}\right\}}_{i=1}^{3}$ with $\alpha =\left\{\frac{1}{4},\frac{1}{3},\frac{2}{3}\right\}$ and b_i = 0.1 is shown in table 2 below. The initial starting values were taken to be between ten and thirty percent of the actual. These are shown in the line corresponding to iteration 0.

There are features here that are typical of such reconstructions. The reconstruction method resolves the lowest fractional power α₁ and its coefficient b₁ quickly as this term is the most persistent one for large times: essential numerical convergence for {α₁, b₁} is obtained by the third iteration. The next lowest power and coefficient lags behind; here {α₂, b₂} is already at the stated accuracy by the fifth iteration. In each case the power is resolved faster and more accurately than its coefficient. The third term also illustrates this; the power is essentially resolved by iteration 8, but in fact its coefficient b₃ is not resolved to the third decimal place until iteration number 30. This is seen quite clearly in the singular values of the Jacobian: the largest singular values correspond to the lowest α-values and the smallest to the coefficients of the largest α-powers.

As might be expected, resolving terms whose powers are quite close is in general more difficult. This is relatively insignificant for low α values. For example, with α₁ = 0.2 and 0.22 < α₂ < 0.25 say, correct resolution will be obtained although the coefficients will take longer to resolve than indicated in table 2. On the other hand if, say ${\alpha }_{1}=\frac{1}{4}$ and α₂ = 0.85, α₃ = 0.9, then with the indicated range of time values used the code will fail to recover this last pair. If this is sensed and now only a single second power is requested this will give a good estimate for α₂ but its coefficient will be overestimated. Also, as indicated previously, α values close to one require an extended time measurement range to stay closer to the asymptotic regime of u(t).

Note that the starting values were at least ten percent away from the exact values and one cannot expect the convergence radius of Newton's methods to be much larger for a problem exhibiting as high nonlinearity as the one at hand. Still, let us point to the fact that in the single term case we do not need any starting guesses but obtain very good results directly from the asymptotic formulas (65) and (66). This might give the idea of using the asymptotics to construct starting guesses in the multi term case and then apply Newton. However, it is unclear how the asymptotics could yield starting values for terms other than the first one. A (theoretical, as it turns out) possibility for exploiting asymptotics in the multi term case is described below.

Successive sequential use of asymptotic formulas

A few words are in order about an approach that from the above discussion might seem a good or even a better alternative.

Since each damping term ${b}_{i}{\partial }_{t}^{{\alpha }_{i}}$ contributes a time trace term with large time behaviour ${c}_{i}{t}^{-{\alpha }_{i}}$, it is feasible to take T sufficiently large so that ${c}_{i}{t}^{-{\alpha }_{i}}\ll {c}_{1}{t}^{-{\alpha }_{1}}$ for t > T, that is, all but the smallest damping power is negligible, and this can then be recovered. In successive steps then we subtract this from the data h to get ${h}_{1}(t)=h(t)-{c}_{1}{t}^{-{\alpha }_{1}}$ and now seek to recover the next lowest α power from the large time values of h₁(t) in a range (δT, T) for 0 < δ < 1. Then these steps can be repeated until there is no discernible signal remaining in the sample interval t_min, t_max.

This indeed works well under the right circumstances for recovering two α values but the coefficients {b_i } are less well resolved. It also requires a delicate splitting of the time interval and gives a much poorer resolution of the two terms in the case where, say α₁ = 0.2 and α₂ = 0.25 than that recovered from the Newton scheme. For the recovery of three damping terms this was in general quite ineffective.

Every time an α_i has been recovered, the remaining signal is significantly smaller than the previous, leading to an equally significant drop in effective accuracy. Also, even if we just make a small error in the coefficient b_i , the relative error that is caused by this becomes completely dominant for large times.

In short, this is an elegant and seemingly constructive approach to showing uniqueness for a finite number of damping terms. However, it has limited value from a numerical recovery perspective when used under a wide range of parameter values.

In this case we are simulating measurements taken over a very limited initial time range: in fact we take the measurement interval to be $t\in \left[0,0.1\right)$. The line of attack is to use the known form of $\hat{h}(s)$ for large values of s and convert the powers of s appearing into powers of t for small times using the Tauberian theorem. In case of two damping terms this gives

$$\begin{equation*}\begin{aligned}\hfill -{c}^{2}\lambda {h}_{\text{small}}(t)-{h}_{\text{small}}^{\prime \prime }(t)& ={c}_{1,1}{t}^{1-{\alpha }_{1}}+{c}_{1,2}{t}^{1-{\alpha }_{2}}+{c}_{2,1}{t}^{3-{\alpha }_{1}}+{c}_{2,2}{t}^{3-{\alpha }_{2}}\hfill \\ \hfill & \quad +{c}_{2,3}{t}^{3-2{\alpha }_{1}}+{c}_{2,4}{t}^{3-{\alpha }_{1}-{\alpha }_{2}}+{c}_{2,5}{t}^{3-2{\alpha }_{2}},\hfill \end{aligned}\end{equation*}$$

where each term c_k,ℓ is computed in terms of {α_i }, {b_i } and λ, cf (73).

The values of {α_i , c_k,ℓ } were then computed from the data by a Newton scheme then finally converted back to the derived values of {α_i } and {b_i }.

The exact values chosen were α = {0.25, 0.2}, b = {0.1, 0.1} and the initial starting guesses were α = {0.3, 0.16}, b = {0.08, 0.12}. We show the progression of the iteration in table 3.

While theory predicts reconstructibility of an arbitrary number of terms in both cases, there is a clear difference in ability to reconstruct terms between the small time and the large time asymptotics. First of all, the method we described effectively only recovers two terms with small time measurements, as compared to three in the large time. The b_i coefficients, which are always harder to obtain than the α_i exponents, are much worse than in the small time than in the large time regime. This is partly explained by the higher degree of ill-posedness due to the necessity of differentiating the data twice.