1 Introduction

This work is stimulated by the physical models studied in [12, 13], where longitudinal elastic waves of a membrane are coupled to viscous fluid flow in the enclosing half space. The aims are to understand the damping of the elastic waves through the coupling to the viscous fluid, on the one hand, and to explain the appearance of the non-classical dispersion relation, on the other hand. Denoting by \(k \in \mathbb {R}^{d-1}\) the horizontal wave vector and \(\omega \in \mathbb {R}\) the angular frequency, the classical elastic wave satisfies a dispersion relation \(\omega ^2 \approx |k|^2\), while the longitudinal pressure waves, here referred to as Lucassen waves (cf. [15]), satisfy \(|\omega |^{3/2} \approx |k|^2\) such that the wave speed \(c(k)=\omega (k)/|k|\) depends on k.

The class of Lucassen waves attracted considerable attention over the last decade due to its biophysical relevance in living organisms, where the transmission of information over biologically relevant distances and timescales is fundamental. The standard model describes the propagation of signals on the vast network of nerve cells via a purely electrical mechanism, unable to explain a number of non-electric phenomena, like the effectiveness of anesthetics scaling with their solubility in lipid membranes [16, 18] or the lower heat dissipation of a nerve in contrast to an electrical cable [22]. As it is known that, alongside the electrical signal, a mechanical displacement travels along the nerve fiber [3, 10], there is a need for a more complete model incorporating these aspects. On the one hand, experimental scientists (see [5, 21]), using a lipid monolayer spread at the air–water interface as a minimal model, have shown that indeed interface-localized pressure waves can propagate in such systems. On the other hand, from a more theoretical viewpoint, all possible surface wave solutions for a viscoelastic membrane atop a half space of viscous fluid have been determined in [9], including the experimentally observed Lucassen waves and their dispersion relations of the type \(|\omega |^{3/2}\approx |k|^2\). In particular, a fractionally damped wave equation was derived for describing the Lucassen waves efficiently in [12, 13]. The biophysical relevance of Lucassen waves is demonstrated by the fact that the wave solutions depend directly on the lateral membrane compressibility \(\kappa \). For example, adsorption of lipophilic substances, like anesthetics, into the membrane presumably alters \(\kappa \) and, as a consequence, is expected to change the wave propagation properties. In addition to that, at large amplitudes, the pressure pulse locally modifies the compressibility \(\kappa \) and thereby significantly increases the propagation distance [21]. This nonlinear property suggests an all-or-none behavior, which indeed is observed in nerve pulse propagation.

Here we want to understand this phenomenon using the mathematically most simple model, which is given by the following coupled system:

$$\begin{aligned} {\rho _{\mathrm{memb}}\ddot{U} = \kappa \Delta _x U - \mu \partial _z v\quad {\text {for }} t>0, \ x\in \Sigma ,} \end{aligned}$$
(1.1a)
$$\begin{aligned} {\dot{U}(t,x)= v(t,x,0)\quad {\text {for }} t>0, \ x\in \Sigma ,} \end{aligned}$$
(1.1b)
$$\begin{aligned} {\rho _{\mathrm{bulk}}\dot{v} = \mu \Delta _{x,z} v \quad {\text {for }} t>0,\ (x,z)\in \Omega \,{:=}\, \Sigma {\times } {]{-}\infty ,0[}.} \end{aligned}$$
(1.1c)

In the physical setup of [12, 13], the domain \(\Sigma = \mathbb {R}^1\) denotes the membrane and \(U(t,x)\in \mathbb {R}\) denotes the horizontal displacement (longitudinal motion) of the membrane. The half space \({\Omega = \Sigma \times {]{-}\infty ,0[}\subset \mathbb {R}^d}\) is filled by a viscous fluid whose horizontal velocity component is \(v(t,x,z) \in \mathbb {R}\) (pure shear flow). Condition (1.1b) is a no-slip condition for the fluid along the membrane, while the induced stress of the sheared fluid is included in (1.1a) via \({{-} \mu \partial _z v(t,x,0)}\).

For mathematical purposes, we can allow \(\Sigma \subset \mathbb {R}^{d-1}\), but to avoid any complications with boundary conditions we assume that \(\Sigma \) is of the form

$$\begin{aligned} {\Sigma = \mathbb {R}^k \times \big (\mathbb {R}_{\!/(\ell \mathbb {Z})}\big )^n \quad {\text { with }} \ell >0 {\text { and }} k+n=d{-}1.} \end{aligned}$$
(1.2)

In particular, \(\Sigma \) is an additive group and (1.1) is translation invariant. In Sect. 2, we first show that the system has the natural energy

$$\begin{aligned} {\mathbb {E}}(U,\dot{U},v)= \int _\Sigma \Big \{\frac{\rho _{\mathrm{memb}}}{2} \dot{U}{}^2 + \frac{\kappa }{2} |\nabla U|^2\Big \} \;\!\mathrm {d}x + \int _\Omega \frac{\rho _{\mathrm{bulk}}}{2} v^2 \;\!\mathrm {d}z \;\!\mathrm {d}x \end{aligned}$$
(1.3)

as a Liapunov function. Thus, the function space \({{\mathbf{H}}:= {\mathrm H}^1(\Sigma )\times {\mathrm L}^2(\Sigma ) \times {\mathrm L}^2(\Omega )}\) is the natural state space. Note that this includes periodic boundary conditions for \({x \in \Sigma = \mathbb {R}^k \times \big (\mathbb {R}_{\!/(\ell \mathbb {Z})}\big )^n}\).

Moreover, we discuss suitable scalings of time t, the horizontal variable \(x\in \mathbb {R}^{d-1}\), and the vertical variable \(z\in {]{-}\infty ,0[}\). We can renormalize all constants such that the system of equations takes the form

$$\begin{aligned} \ddot{U}= \Delta _x U - \partial _zv|_{z=0} \quad {\text {for }} t>0, \ x\in \Sigma , \end{aligned}$$
(1.4a)
$$\begin{aligned} \dot{U}= v|_{z=0}\quad {\text {for }} t>0, \ x\in \Sigma , \end{aligned}$$
(1.4b)
$$\begin{aligned} \dot{v} = \varepsilon ^2 \Delta _x v + \partial _z^2 v \quad {\text {for }} t>0,\ (x,z)\in \Omega , \end{aligned}$$
(1.4c)

with the parameter \(\varepsilon = \mu /\sqrt{\rho _{\mathrm{memb}}k\,}\) . The essential point is here that the scaling of the horizontal variable \(x\in \mathbb {R}^{d-1}\) is different from the vertical variable \(z \in {]{-}\infty ,0[}\), thus breaking the isotropy of the diffusion \(\mu \Delta _{x,z} v\) in (1.1c).

Our interest lies in the case \(\varepsilon \rightarrow 0\). Indeed, in Sect. 2.4 we simply set \(\varepsilon =0\) in (1.4c) and show that this limit allows us to solve the scalar, one-dimensional heat equation \({\dot{v} = \partial _z^2 v}\) on \({]{-}\infty ,0[}\) for each \(x\in \Sigma \) independently. Assuming the initial condition \(v(0,x,z)=0\) and using the Dirichlet boundary condition \(v(t,x,0)=\Phi (t,x)\), the stress \({\partial _z v(t,x,0)}\) can be explicitly expressed via the heat kernel, namely

$$\begin{aligned} {\partial _z v(t,x,0)= \int _0^t \frac{1}{\sqrt{\pi (t{-}\tau )}} \, \dot{\Phi }(\tau ,x) \;\!\mathrm {d}\tau ,} \end{aligned}$$
(1.5)

see (4.3). It is this one-dimensional parabolic Dirichlet-to-Neumann map that introduces the fractional damping into the wave equation. In particular, denoting the fractional (Caputo) derivative of order \(\alpha \in {]0,1[}\) of the function g with \(g(0)=0\) by

$$\begin{aligned} {}^{\mathrm C}{\mathrm D}^\alpha g: \ t \ \mapsto \ \ \frac{1}{\Gamma (1{-}\alpha )} \int _0^1 \frac{1}{(t{-}\tau )^\alpha } \, \dot{g}(\tau )\;\!\mathrm {d}\tau , \end{aligned}$$
(1.6)

we see that the mapping in (1.5) takes the form \({\partial _z v(t,x,0)= (\,{}^{\mathrm C}{\mathrm D}^{1/2} \Phi )(t,x)}\), which is a Caputo derivative of order 1/2.

Indeed, if we solve (1.4) with \(\varepsilon =0\) and \(v(0,x,z)=0\), then we can eliminate v totally by exploiting (1.5) with \(\Phi =\dot{U}\), and U has to solve

$$\begin{aligned} \ddot{U}(t,x) + \int _0^t \frac{1}{\sqrt{\pi (t{-}\tau )}} \, \ddot{U}(\tau ,x) \;\!\mathrm {d}\tau = \Delta _x U(t,x) \quad {\text {for }} t>0, \ x\in \Sigma . \end{aligned}$$
(1.7)

This is a fractionally damped wave equation where the damping is generated by a fractional Caputo derivative of order 3/2, and this fractional derivative acts locally with respect to the space variable \(x \in \Sigma \).

In Sect. 2.5, we follow the approach in [9, 12, 12] and discuss the dispersion relations for our normalized system (1.4) and show that, for the limit case \(\varepsilon =0\), the dispersion relation reads \((\mathrm {i}\omega )^2 + (\mathrm {i}\omega )^{3/2} + |k|^2=0\), where \(\mathop {\mathrm {Im}}\omega \ge 0\) is enforced by the stability through the Liapunov function \({\mathbb {E}}\). Hence, for small |k| we obtain the new dispersion relation

$$\begin{aligned} \omega (k)=\Big (\pm \frac{\sqrt{3}}{2} + \frac{\mathrm {i}}{2} \Big ) \,|k|^{4/3} \ + \ O(|k|^2)_{k\rightarrow 0}. \end{aligned}$$

There is a rich mathematical literature on linear and nonlinear partial differential equations involving fractional time derivatives, see, e.g., [1, 11, 20, 23, 24, 26]. Moreover, the appearance of fractional derivatives via Dirichlet-to-Neumann maps is a well-known phenomenon, see, e.g., [6, Ch. 1.1, Eqn. (1.2)] for an early reference or [17] for applications in electrochemistry. Indeed, there is an extensive mathematical literature on the so-called harmonic extension method, where fractional operators are expressed via suitable Dirichlet-to-Neumann maps, see [2] and the recent generalization in [7] allowing for fractional time derivatives.

Our focus is different, because we want to show that (1.7) appears as a rigorous limit for \(\varepsilon \rightarrow 0^+\) in the coupled system (1.4). For this, in Sect. 3 we develop the linear semigroup theory by showing that the semigroups \(\mathrm {e}^{t A_\varepsilon }:{\mathbf{H}}\rightarrow {\mathbf{H}}\) exist for all \(\varepsilon \ge 0\) and are bounded in norm by \(C(1{+}t)\). In Theorem 3.3, we establish the strong convergence \(\mathrm {e}^{t A_\varepsilon } w_0 \rightarrow \mathrm {e}^{t A_0} w_0\) for \(\varepsilon \rightarrow 0^+\), which holds for all \(t>0\) and \(w_0 \in {\mathbf{H}}\). For more regular initial conditions \(w_0\), we obtain the quantitative estimate

$$\begin{aligned} \Vert \mathrm {e}^{t A_\varepsilon } w_0 - \mathrm {e}^{t A_0} w_0\Vert _{\mathbf{H}}\le \sqrt{\varepsilon \,t}\, (2.3{+}t)^2 \, \big ( \Vert w_0\Vert _{{\mathbf{H}}} + \Vert \nabla _x w_0\Vert _{{\mathbf{H}}}\big ). \end{aligned}$$

In Sect. 4, we return to the energetics and the dissipation for the damped wave equation. By starting from the natural energy and dissipation in the PDE system (1.4) with \(\varepsilon =0\) and the explicit solution for v(t, x, z) in terms of \(\dot{U}(\tau ,x)\), we obtain a natural energy functional \({\mathbf{E}}\) for the fractionally damped wave equation that is non-local in time:

$$\begin{aligned} {\mathbf{E}}(U(t), \big [\dot{U}(\cdot )\big ]_{[0,t]}) = \int _\Sigma \Bigg \{&\frac{1}{2}\dot{U}(t,x)^2 + \frac{1}{2}|\nabla U(t,x)|^2\nonumber \\&+\int _0^t\!\int _0^t \! \frac{1}{4\sqrt{\pi }(2t{-}r{-}s)^{3/2}} \, \dot{U}(r,x) \dot{U}(s,x) \;\!\mathrm {d}s \;\!\mathrm {d}r \Bigg \} \;\!\mathrm {d}x, \end{aligned}$$
(1.8a)

where \(\big [\dot{U}(\cdot )\big ]_{[0,t]}\) indicates the dependence on \(\dot{U}(s)\) for \(s\in [0,t]\). For solutions U of (1.7), we obtain an energy–dissipation balance with a non-local dissipation:

$$\begin{aligned} \frac{{\mathrm d}}{{\mathrm d}t} {\mathbf{E}}(U(t), \big [\dot{U}(\cdot )\big ]_{[0,t]}) = - \int _\Sigma \int _0^t\int _0^t \frac{1}{\sqrt{\pi }(2t{-}r{-}s)^{1/2}} \,\ddot{U}(r,x) \ddot{U}(s,x) \;\!\mathrm {d}s \;\!\mathrm {d}r \, \;\!\mathrm {d}x. \end{aligned}$$
(1.8b)

Related results are obtained in [8, 23, 25], but there typically only energy–dissipation inequalities are derived. It is surprising to see that the two non-local kernels in (1.8) that depend on \(t{-}r,\ t{-}s\in [0,T]\) are only depending on the sum \((t{-}r) + (t{-}s) = 2t{-}r{-}s\), which derives from very specific scaling properties of the heat kernel.

2 The formal modeling

In this section, we describe the formal modeling, including the energy functional, the scalings and the derivation of the fractionally damped wave equation as the scaling limit.

2.1 The energy functional and the state space

We return to the full system (1.1) and observe that it has the form of a damped Hamiltonian system with the total energy \({\mathbb {E}}(U,\dot{U}, v)\) given in (1.3). Indeed, taking the time derivative along solutions \(t\mapsto (U(t),v(t))\) of (1.1) we find

$$\begin{aligned}&\frac{{\mathrm d}}{{\mathrm d}t} {\mathbb {E}}(U(t),\dot{U}(t), v(t)) = \int _\Sigma \big (\rho _{\mathrm{memb}}\dot{U} \ddot{U} +\kappa \nabla _{\!x} \dot{U} \cdot \nabla _{\!x} U \big ) \;\!\mathrm {d}x + \int _\Omega \rho _{\mathrm{bulk}}v\,\dot{v} \;\!\mathrm {d}x \;\!\mathrm {d}z \\&\overset{\text {(1.1a),(1.1c)}}{=}\int _\Sigma U_t \mu v_z \;\!\mathrm {d}x + \int _\Omega \mu \,v\,\Delta _{x,z} v \;\!\mathrm {d}x \;\!\mathrm {d}z \overset{\text {(1.1b)}}{=} - \int _\Omega \mu ( |\nabla _{\!x} v|^2 {+}v_z^2) \;\!\mathrm {d}x \;\!\mathrm {d}z \le 0. \end{aligned}$$

Here we used that the integration by parts \(\int _\Sigma \nabla _{\!x} \dot{U} \cdot \nabla _{\!x} U \;\!\mathrm {d}x = - \int _\Sigma \dot{U} \Delta _x U\;\!\mathrm {d}x \) does not generate boundary terms because \(\Sigma \) has the form (1.2).

Thus, \({\mathbb {E}}\) acts as a Liapunov function and it is a bounded quadratic form on the Hilbert space \({{\mathbf{H}}= {\mathrm H}^1(\Sigma )\times {\mathrm L}^2(\Sigma )\times {\mathrm L}^2(\Omega )}\), which we consider as the basic state space for our problem. In Sect. 3 we will show that (1.1) has a unique solution for each initial value \(w_0=(U(0),\dot{U}(0),v(0)) \in {\mathbf{H}}\).

More precisely, the system (1.1) can be written as a damped Hamiltonian system for the states \(X=(U,P,p)\) where \(P=\rho _{\mathrm{memb}}\dot{U}\) and \(p=\rho _{\mathrm{bulk}}v\). With \({\mathbb {E}}(U,P,p)={\mathbb {E}}(U,\frac{1}{\rho _{\mathrm{memb}}}P,\frac{1}{\rho _{\mathrm{bulk}}} p)\) we have

$$\begin{aligned} {\left( \!\!\begin{array}{c}\dot{U}\\ \dot{P}\\ \dot{p} \end{array}\!\!\right) = \big ( {\mathbb {J}} {-} \mathbb K\big ) {\mathrm D}{\mathbb {E}}(U,P,p) \ {\text {with }}{\mathbb {J}}= \left( \begin{array}{ccc}0&{}I&{}0 \\ \!\!\!{-}I&{}0&{}0\\ 0&{}0&{}0 \end{array}\right) {\text { and }}{\mathbb {K}} = \left( \begin{array}{ccc}0&{}0&{}0 \\ 0&{}*&{}\mu \partial _z\Box |_{z=0}\!\!\\ 0&{}*&{}{-}\mu \Delta _{x,z}\! \end{array}\right) \!\!,} \end{aligned}$$

where \({{\mathbb {K}}: \mathrm {dom}({\mathbb {K}}) \subset {\mathbf{X}}\rightarrow {\mathbf{X}}:={\mathrm L}^2(\Sigma ) \times {\mathrm L}^2(\Sigma )\times {\mathrm L}^2(\Omega )}\) is defined as the self-adjoint, unbounded operator induced by the quadratic dissipation potential

$$\begin{aligned} {\mathcal R}^*(\Pi ,\Xi ,\xi )&:=\frac{\mu }{2} \int _\Omega |\nabla _{x,z}\xi |^2 \;\!\mathrm {d}x \;\!\mathrm {d}z + \chi _*(\Xi ,\xi ), \\&\qquad {\text { with }}\chi _*(\Xi ,\xi )= \left\{ \begin{array}{c@{\quad }l} 0 &{} {\text {if }}\Xi =\xi |_{z=0}, \\ \infty &{}\text {else}. \end{array}\right. \end{aligned}$$

2.2 A long-wave scaling

To obtain a first understanding of the different scaling of horizontal and vertical spatial variables, we study the long-wave scaling for the wave equation. This mean that we scale the horizontal space variable x and the time variable t with the same factor \(\delta >0\). For the moment, we assume that the membrane constants \(\rho _{\mathrm{memb}}\) and k are given and of order 1, while \(\rho _{\mathrm{bulk}}\) and \(\mu \) are much smaller. More general scalings are discussed in the following subsection.

Without loss of generality, we keep U fixed and obtain velocities \(\dot{U}\) of order \(\delta \). Hence, to keep the no-slip condition, we also need to rescale v by a factor \(\delta \). The main point is that we want z to be rescaled by a smaller factor, let us say \(\delta ^\alpha \) with \(\alpha \in {[0,1[}\). This implies that \({\partial _z v}\) scales like \(\delta ^{1+\alpha }\). Thus, to treat the coupling term \({\mu \partial _z v|_{z=0}}\) of the same order as \(\ddot{U}\) and \(\nabla _x U\), we need to assume that \(\mu \) also scales with \(\delta \), namely like \(\delta ^{1-\alpha }\). Finally, we also assume the appropriate scaling for \(\rho _{\mathrm{bulk}}\), namely

$$\begin{aligned} {\hat{x}} = \delta \,x, \ \ \ {\hat{t}} = \delta \, t, \ \ \ {\hat{z}} = \delta ^\alpha \, z, \ \ \ {\hat{U}} = U, \ \ \ v = \delta \, {\hat{v}}, \ \ \ \mu = \delta ^{1-\alpha } {\hat{\mu }}, \ \ \ \rho _{\mathrm{bulk}}=\delta ^\alpha \hat{\ }\rho _{\mathrm{bulk}}. \end{aligned}$$

Hence, this long-wave scaling with small \(\delta \) is indeed suitable, if the bulk quantities \(\rho _{\mathrm{bulk}}\) and \(\mu \) are much smaller than the membrane quantities \(\rho _{\mathrm{memb}}\) and stiffness k.

Inserting these scalings (and dropping the hats), we find the transformed system

$$\begin{aligned} {\begin{aligned}&\rho _{\mathrm{memb}}\ddot{U}= \kappa \Delta _x U - \mu \partial _z v|_{z=0} {\text { and }} \dot{U}= v|_{z=0}&{\text { on }}&\Sigma ,\\&\rho _{\mathrm{bulk}}\dot{v} = \mu (\delta ^{2{-}2\alpha } \Delta _x {+} \partial _z^2)v&{\text { in }}&\Omega . \end{aligned}} \end{aligned}$$
(2.1)

Here, the case of small \(\delta \) is relevant, and in the limit \(\delta \rightarrow 0^+\) we obtain the fractionally damped wave equation.

2.3 Non-dimensionalizing by a general scaling

We fully non-dimensionalize the system by considering general rescalings, where we scale x, z, and t independently:

$$\begin{aligned} {\hat{x}} = a \,x, \quad {\hat{t}} = b \,t, \quad {\text {and }} {\hat{z}} = c\, z, \end{aligned}$$

but do not assume any scaling on the material parameters \(\rho _{\mathrm{memb}},\ \kappa ,\ \rho _{\mathrm{bulk}}\), and \(\mu \). We keep \({\hat{U}}=U\) (which is always possible by linearity), but need to rescale \(v= b\,{\hat{v}}\) to transform the no-slip condition \(\dot{U} = v|_{z=0}\) into \({\partial _{{\hat{t}}} {\hat{U}} = \hat{v}|_{{\hat{z}}=0}}\). The transformed equations read (after dropping the hats) as follows

$$\begin{aligned} \rho _{\mathrm{memb}}b^2 \ddot{U}&= \kappa a^2 \Delta _{x}U - \mu bc\,\partial _z v|_{z=0} \ \ {\text { on }}\Sigma , \\ \rho _{\mathrm{bulk}}b^2 \dot{v}&= \mu bc^2 \partial _z^2 v + \mu a^2b \Delta _x v \qquad {\text { in }}\Omega . \end{aligned}$$

Dividing the equations by \(\rho _{\mathrm{memb}}b^2\) and \(\rho _{\mathrm{bulk}}b^2\) respectively, we can equate the first three of the four coefficients to 1, namely

$$\begin{aligned} \frac{\kappa a^2}{\rho _{\mathrm{memb}}b^2} = 1,\quad \frac{\mu c}{\rho _{\mathrm{memb}}b} = 1, \quad \frac{\mu c^2 }{\rho _{\mathrm{bulk}}b}=1. \end{aligned}$$

We obtain the solution

$$\begin{aligned} a^2 =\frac{\rho _{\mathrm{bulk}}^2\mu ^2}{\rho _{\mathrm{memb}}^3 \kappa }, \quad b = \frac{\rho _{\mathrm{bulk}}\mu }{\rho _{\mathrm{memb}}^2}, \quad {\text {and }} c = \frac{\rho _{\mathrm{bulk}}}{\rho _{\mathrm{memb}}}, \end{aligned}$$

and the remaining fourth coefficient reads

$$\begin{aligned} \varepsilon := \Big (\frac{\mu a^2}{\rho _{\mathrm{bulk}}b}\Big )^{1/2} = \frac{a}{c}= \frac{\mu }{\sqrt{\rho _{\mathrm{memb}}\kappa }}. \end{aligned}$$

The non-dimensionalized coupled system now reads

$$\begin{aligned} {\ddot{U} = \Delta _{x}U - \partial _z v|_{z=0} {\text { and }} \dot{U}=v|_{z=0} {\text { on }}\Sigma , \quad \dot{v} = \varepsilon ^2 \Delta _x v +\partial _z^2 v {\text { in }}\Omega ,} \end{aligned}$$

which is exactly the renormalized system (1.4), which is studied subsequently.

Hence, the system (1.1) has a unique non-dimensional parameter \(\varepsilon = \mu /\sqrt{\rho _{\mathrm{memb}}\, k\,}\) that describes the effective anisotropy of the diffusion in the bulk \(\Omega \). Subsequently, we are interested in the case of very small \(\varepsilon \) and indeed in the limit \(\varepsilon \rightarrow 0^+\).

We can interpret \(\varepsilon \) as the ratio of three different length scales. Choosing an arbitrary timescale \(t_0>0\), we have the diffusion length \(\ell _\mathrm {diff}\), the “equivalent membrane thickness” \(\ell _\mathrm {thick}\), and the membrane travel length \(\ell _\mathrm {trav}\) given by

$$\begin{aligned} \ell _\mathrm {diff}(t_0) = \Big ( \!\frac{\mu t_*}{\rho _{\mathrm{bulk}}}\!\Big )^{1/2}\!, \ \ \ell _\mathrm {thick}= \frac{\rho _{\mathrm{memb}}}{\rho _{\mathrm{bulk}}}, \ \ \ell _\mathrm {trav}(t_0)= t_0 c_\mathrm {memb} = t_0 \Big (\!\frac{\kappa }{\rho _{\mathrm{memb}}}\!\Big )^{1/2}, \end{aligned}$$

where \(c_\mathrm {memb}\) is the wave speed in the undamped membrane. Now our dimensionless parameter \(\varepsilon \) is given by

$$\begin{aligned} \varepsilon = \frac{\ell _\mathrm {diff}(t_0)^2}{\ell _\mathrm {thick}\; \ell _\mathrm {trav}(t_0)} \quad {\text {for all }}\; t_0>0. \end{aligned}$$

To make the definition even more intrinsic, we may choose \(t_0\) as a characteristic time \(t_*\) for the system. We ask that the time \(t_*\) is chosen such that the corresponding diffusion length scale equals the equivalent membrane thickness, viz. \(\ell _\mathrm {diff}(t_*)=\ell _\mathrm {thick}\). This yields

$$\begin{aligned} t_*= \frac{\rho _{\mathrm{memb}}^2}{ \mu \rho _{\mathrm{bulk}}} \quad \text {and} \quad \ell _\mathrm {trav}(t_*) = t_* c_\mathrm {memb} = \frac{ \rho _{\mathrm{memb}}^{3/2} \kappa ^{1/2}}{\mu \rho _{\mathrm{bulk}}}. \end{aligned}$$

The scalings of time and horizontal and vertical lengths are now given as

$$\begin{aligned} t = t_* \,{\hat{t}}, \quad x = \ell _\mathrm {trav}(t_*)\,\hat{x},\quad z = \ell _\mathrm {thick} \,{\hat{z}}. \end{aligned}$$

This leads to the final relation

$$\begin{aligned} \varepsilon = \frac{\ell _\mathrm {thick}}{ \ell _\mathrm {trav}( t_*) } =\frac{ \mu }{ \sqrt{\rho _{\mathrm{memb}}\kappa }}. \end{aligned}$$

Typical parameters for the experimental setup consisting of a lipid monolayer, such as DPPC at the water–air interface, are \(\rho _{\mathrm{memb}}= 10^{-6}~\mathrm{kg / m^2}\), \(\mu = 10^{-3}~\mathrm{Pa\,s} = 10^{-3}~\mathrm{kg/(s\,m)}\), and \(\kappa = 10^{-2}~\mathrm{N / m}\), where for \(\rho _{\mathrm{memb}}\) the surface excess mass density was used. These parameters yield \(\varepsilon = 10\). Although this value is not small, it does not contradict our argumentation. As shown in [9], different waves can coexist in such a system, the longitudinal capillary waves with \(|\omega |^{3/2} \approx |k|^2\) being only one of them. In particular, it is interesting that the dispersion of this wave, which has been known in the literature since [15], follows from our general calculation as a rigorous limit.

2.4 The limit model and the fractionally damped wave equation

We now study the limit equation by setting \(\varepsilon =0\) in the rescaled system (1.4). The justification of taking this limit is given in the following section.

After setting \(\varepsilon =0\), we obtain the system

$$\begin{aligned} \ddot{U} (t,x)= \Delta _x U(t,x) - \partial _z v(t,x,0)&{\text {for }} t>0, \ x\in \Sigma , \end{aligned}$$
(2.2a)
$$\begin{aligned} \dot{U}(t,x) = v(t,x,0)&{\text {for }} t>0, \ x\in \Sigma , \end{aligned}$$
(2.2b)
$$\begin{aligned} \dot{v}(t,x,z) = \partial _z^2 v(t,x,z)&{\text {for }} t>0,\ (x,z) \in \Omega . \end{aligned}$$
(2.2c)

The important point is now that the equation for (2.2c) can be solved explicitly by the use of the properly rescaled one-dimensional heat kernel \(H(t,y)=(4\pi t)^{-1/2} \mathrm {e}^{-y^2/(4t)}\). Note that \(x\in \Omega \) appears now as a parameter only, since the diffusion in x-direction is lost.

The solution of (2.2c) with the boundary condition (2.2b) and the initial condition \(v(0,x,z)=0\) takes the explicit form

$$\begin{aligned} {v(t,x,z)= \int _0^t 2 \partial _z H(t{-}\tau ,z) \, \dot{U}(\tau ,x) \;\!\mathrm {d}\tau ,} \end{aligned}$$

see Sect. 4.1 for a derivation. Taking the derivative with respect to z and using that the heat kernel H satisfies \({\partial _z^2 H= \partial _t H}\), we obtain

$$\begin{aligned} {\partial _z v(t,x,z)= \int _0^t 2\partial _t H(t{-}\tau ,z)\,\dot{U}(\tau ,x) \;\!\mathrm {d}\tau = \int _0^t 2 H(t{-}\tau ,z)\,\ddot{U}(\tau ,x) \;\!\mathrm {d}\tau ,} \end{aligned}$$

where for the integration by parts in the last identity we exploited \(\dot{U}(t,x,0)=v(t,x,0)=0\) and \(H(0,z)=0\) for \(z<0\). Thus, evaluation at \(z=0\) and using \(H(t,0)=(4\pi t)^{-1/2}\), the coupling term in (2.2a) reduces to

$$\begin{aligned} {\partial _z v(t,x,0)= \int _0^t 2H(t{-}\tau ,0) \ddot{U}(\tau ,x)\;\!\mathrm {d}\tau = \int _0^\tau \frac{1}{\sqrt{\pi (t{-}\tau )}} \, \ddot{U}(\tau ,x) \;\!\mathrm {d}\tau .} \end{aligned}$$
(2.3)

Through these formulas, we see how kinetic energy is moved from the membrane via the no-slip condition (2.2b) into the one-dimensional diffusion equation. Through the memory kernel in (2.3) the energy is restored partially in a delayed fashion, which leads to a fractional damping, here of order 3/2. Because in the bulk \(\Omega = \Sigma \times {]{-}\infty ,0[}\) there is no coupling between different points \(x \in \Sigma \), this damping is non-local in time but local with respect to \(x\in \Sigma \).

Joining everything, we see that the limiting system (2.2) contains the fractionally damped wave equation

$$\begin{aligned} \ddot{U}(t,x) + \int _0^t \frac{1}{\sqrt{\pi (t{-}\tau )}}\, \ddot{U}(\tau ,x)\;\!\mathrm {d}\tau = \Delta U(t,x) \quad \text {on } \Sigma . \end{aligned}$$
(2.4)

An analysis concerning existence of solutions and concerning the energetics is given in Sect. 4.

2.5 The dispersion relations

Following [12], we consider special solutions of (3.1) obtained by a Fourier ansatz. For the temporal growth factor \(\mu =\mathrm {i}\omega \in {\mathbb C}\) with \(\mathop {\mathrm {Re}}\mu \le 0\) and the wave vector \(k \in \mathbb {R}^{d-1}\) we set

$$\begin{aligned} (U(t,x),V(t,x), v(t,x,z)) = \mathrm {e}^{\mu t + \mathrm {i}k \cdot x} \big ( a,b, c\,\mathrm {e}^{\gamma z} \big ) \end{aligned}$$

with \(a,b,c,\gamma \in {\mathbb C}\) and \(\mathop {\mathrm {Re}}\gamma >0\). From \(V=U_t\), we obtain \(b=a\mu \), while \(v|_{z=0}=V\) implies \(c=b=a\mu \). Finally, we have to satisfy the membrane equation \(U_{tt} = \Delta _{x} U - v_z|_{z=0}\) and the diffusion equation \({\partial _t v= \varepsilon ^2 \Delta _x v+\partial _z^2 v}\), which leads to the algebraic relations (since only \(a\ne 0\) is interesting) \( \mu ^2 = -| k |^2 - \mu \gamma \) and \(\mu = -\varepsilon ^2 | k |^2 + \gamma ^2\). As in [12], we eliminate the variable \(\gamma \) and obtain the dispersion relation

$$\begin{aligned} 0=\Gamma (\mu , k )= \big (\mu ^2{+}|k|^2\big )^2 - \varepsilon ^2 \mu ^2|k|^2 - \mu ^3, \end{aligned}$$

where we still need to be careful to satisfy \(\mathop {\mathrm {Re}}\mu \le 0\) and \(\mathop {\mathrm {Re}}\gamma >0\) with \(\gamma ^2=\mu {+} \varepsilon ^2|k|^2\).

For short waves, i.e., \(|k|\gg 1\), we obtain the expansion

$$\begin{aligned} \mu = - \mathrm {i}|k| - \frac{|k|}{2} \,\Big ( \varepsilon ^2 - \frac{\mathrm {i}}{|k|}\Big )^{1/2} + \text {h.o.t.} \end{aligned}$$

In the case \(\varepsilon >0\), this means that short waves travel at speed 1, but are damped proportional to |k|. The limit \(\varepsilon =0\) leads to a significantly smaller damping, namely one of order \(|k|^{1/2}\).

As expected due to the scaling discussed in the previous subsections, the case of long waves, i.e., \(|k| \ll 1\), is not so sensitive with respect to \(\varepsilon \). For all \(\varepsilon \ge 0\), we find the expansion

$$\begin{aligned} \mu = - \Big ( \frac{1}{2} \pm \frac{\mathrm {i}\sqrt{3}}{2}\Big ) |k|^{4/3} + \text {h.o.t.} \end{aligned}$$

In particular, we find that the waves slow down for \(|k|\rightarrow 0\), because the wave speed takes the form \(c(k) =\mathop {\mathrm {Im}}(\mu (k)/|k|) = \pm |k|^{1/3}\sqrt{3}/2 + \)h.o.t. Moreover, the damping is very low, because it is proportional to \(|k|^{4/3} \).

3 Convergence result for the semigroup

From now on, it suffices to consider the rescaled system, where \(\varepsilon >0\) appears as the only small parameter:

$$\begin{aligned} \ddot{U}= \Delta _x U - \partial _zv|_{z=0}&{\text {for }} t>0, \ x\in \Sigma , \end{aligned}$$
(3.1a)
$$\begin{aligned} \dot{U}= v|_{z=0}&{\text {for }} t>0, \ x\in \Sigma , \end{aligned}$$
(3.1b)
$$\begin{aligned} \dot{v} = \varepsilon ^2 \Delta _x v + \partial _z^2 v&{\text {for }} t>0,\ (x,z)\in \Omega =\Sigma {\times } {]{-}\infty ,0[}. \end{aligned}$$
(3.1c)

In this section, we first prove existence of solutions for the initial boundary value problem and then show that in the limit \(\varepsilon \rightarrow 0\) the corresponding solutions \(t\mapsto w_\varepsilon (t) \in {\mathbf{H}}\) converge strongly to \(t\mapsto w_0(t)\) in the Hilbert space \({\mathbf{H}}\). For this, it is sufficient to employ the classical theory of Trotter and Kato, see, e.g., [19, Sec. 3.3], where convergence of the resolvent implies convergence of the semigroup.

3.1 Formulation via strongly continuous semigroups

By introducing the variable \(V=\dot{U}\) and setting \(w=(U,V,v)\), we rewrite system (3.1) in the form \(\dot{w} = A_\varepsilon w\) and will show that the solutions w can be obtained in the form \(w(t)= \mathrm {e}^{t A_\varepsilon } w_0\), i.e., we have to show that \(A_\varepsilon \) is the generator of a strongly continuous semigroup on the space \({\mathbf{H}}= {\mathrm H}^1(\Sigma )\times {\mathrm L}^2(\Sigma )\times {\mathrm L}^2(\Omega )\). We define the unbounded linear operators \(A_\varepsilon : D(A_\varepsilon ) \subset {\mathbf{H}}\rightarrow {\mathbf{H}}\) via

$$\begin{aligned}&D(A_\varepsilon ) = \big \{\, (U,V,v)\in {\mathrm H}^2(\Sigma ) \times {\mathrm H}^1(\Sigma ) \times \big (X_1^\varepsilon (\Omega ){\cup } Y^\varepsilon (\Omega )\big ) \, \big | \, v|_{z=0}=V \text { on } \Sigma \,\big \} , \\& A_\varepsilon \left( \begin{array}{c} U\\ V\\ v \end{array}\right) = \left( \begin{array}{c} V \\ \Delta _x U - (\partial _z v)|_{z=0} \\ \varepsilon ^2 \Delta _x v + \partial _z^2 v \end{array}\right) . \end{aligned}$$

Here the spaces \(X_\lambda ^\varepsilon (\Omega )\) with \(\lambda >0\) and \(Y^\varepsilon (\Omega )\) are defined via

$$\begin{aligned} X_\lambda ^\varepsilon (\Omega )&:= \big \{\, v \in {\mathrm L}^2 (\Omega ) \, \big | \, \varepsilon ^2\Delta _x v +\partial _z^2 v - \lambda v=0, \ v|_{z\in 0}\in {\mathrm H}^1(\Sigma ) \,\big \} \ \text { and} \\ Y^\varepsilon (\Omega )&:= \big \{\, v\in {\mathrm L}^2(\Omega ) \, \big | \, \varepsilon ^2 \Delta _x v + \partial _z^2 v \in {\mathrm L}^2(\Omega ),\ v|_{z=0}=0 \,\big \} \,. \end{aligned}$$

We emphasize that the domain for \(\varepsilon =0\) is different from the domains for \(\varepsilon >0\), because of the missing x-derivatives for v in the first case. Nevertheless, the trace of v at \(z=0\) is well defined in \({\mathrm L}^2(\Sigma )\) because \({\partial _z v}\) lies in \({\mathrm L}^2(\Sigma , {\mathrm H}^1({]{-}\infty ,0[}))\) and \({\mathrm H}^1( {]{-}\infty ,0[})\) embeds continuously into \({\mathrm C}^0( {]{-}\infty ,0[})\).

More precisely, for \(\varepsilon >0\) we may apply the classical elliptic regularity theory from [14] which shows that \(Y^\varepsilon (\Omega )\) is a closed subspace of \(H^2(\Omega )\) whereas \(X_\lambda ^\varepsilon (\Omega )\) is only contained in \(H^{3/2}(\Omega )\) but not in \({\mathrm H}^2(\Omega )\). In Step 3 of the proof below, we will show that for \(\varepsilon >0\) we have

$$\begin{aligned} X_\lambda ^\varepsilon (\Omega ) \cup Y^\varepsilon (\Omega ) = X_1^\varepsilon (\Omega ) \cup Y^\varepsilon (\Omega ) \text { for all }\lambda >0. \end{aligned}$$
(3.2)

For \(\varepsilon =0\), the spaces \(X_\lambda ^0(\Omega )\) and \(Y^0(\Omega )\) have lower regularity in \(x \in \Sigma \), namely

$$\begin{aligned} X_\lambda ^0(\Omega )&= \big \{\, v\in {\mathrm H}^1(\Omega ) \, \big | \, v(x,z)=\mathrm {e}^{\sqrt{\lambda } z}v(x,0),\ v(\cdot ,0)\in {\mathrm H}^1(\Sigma ) \,\big \} , \\ Y^0(\Omega )&:= \big \{\, v \in {\mathrm L}^2\big (\Sigma ;{\mathrm H}^2({]{-}\infty ,0[})\big ) \, \big | \, v(x,0)=0 \text { a.e.\ in } \Sigma \,\big \} . \end{aligned}$$

Since \(z\mapsto \mathrm {e}^z - \mathrm {e}^{\sqrt{\lambda }z}\) lies in \(H^2( {]{-}\infty ,0[}) \) and vanishes at \(z=0\), we easily see \(X_\lambda ^0(\Omega ) \cup Y^0(\Omega ) = X_1^0(\Omega ) \cup Y^0(\Omega )\) for all \(\lambda >0\).

Our first result in this section shows that for each \(\varepsilon \ge 0\) the operator \(A_\varepsilon \) generates a strongly continuous semigroup \((\mathrm {e}^{tA_\varepsilon })_{t\ge 0} \) on \({\mathbf{H}}\) with a uniform growth rate 1.

Theorem 3.1

(Generation of semigroups). For all \(\varepsilon \ge 0\), the operators \(A_\varepsilon \) defined above are closed. For real \( \lambda >0\), the resolvents \((A_\varepsilon {-}\lambda I)^{-1}: {\mathbf{H}}\rightarrow D(A_\varepsilon )\subset {\mathbf{H}}\) exist and satisfy the estimate

$$\begin{aligned} \big \Vert (A_\varepsilon {-}\lambda I)^{-1} \big \Vert _{{\mathbf{H}}\rightarrow {\mathbf{H}}} \le \frac{1}{\lambda -1} \text { for } \lambda >1. \end{aligned}$$
(3.3)

In particular, \(A_\varepsilon \) is the generator of the strongly continuous semigroup \(\mathrm {e}^{t A_\varepsilon }:{\mathbf{H}}\rightarrow {\mathbf{H}}\) satisfying \(\Vert \mathrm {e}^{t A_\varepsilon } \Vert _{{\mathbf{H}}\rightarrow {\mathbf{H}}} \le 1+t\,\mathrm {e}/2\) for \(t\ge 0\). Moreover, the functional energy

$$\begin{aligned} {\mathbb {E}}_0(U,\dot{U},v)= \int _\Sigma \Big \{\frac{1}{2} \dot{U}(x)^2 + \frac{1}{2} |\nabla U(x)|^2 + \int _{-\infty }^0 \frac{1}{2}\,v(x,z)^2 \;\!\mathrm {d}z \Big \} \;\!\mathrm {d}x \end{aligned}$$
(3.4)

is a Liapunov function, i.e., along solutions we have the estimate

$$\begin{aligned} {\mathbb {E}}_0(U(t),\dot{U}(t),v(t)) \le {\mathbb {E}}_0(U(s),\dot{U}(s),v(s)) \quad \text { for } t > s \ge 0. \end{aligned}$$

Proof

We first treat the case \(\varepsilon >0\) in Steps 1 to 3 and then discuss the differences for the case \(\varepsilon =0\) in Step 4.

Step 1: A priori estimate. For \(\alpha >0\) we use the norm \(|\cdot |_\alpha \) on \({\mathrm H}^1(\Sigma )\) defined via \(|U|_\alpha ^2 = \alpha ^2\Vert U\Vert _2^2 + \Vert \nabla _xU\Vert _2^2\), where \(\Vert \cdot \Vert _2\) is the standard \({\mathrm L}^2\) norm on \(\Sigma \). For \(\widetilde{w}=(\widetilde{U},\widetilde{V},\widetilde{v})\) and \(\widehat{w}=(\widehat{U},\widehat{V},\widehat{v})\) in \({\mathbf{H}}\), we define \(\big \langle \!\big \langle \widetilde{w},\widehat{w} \big \rangle \!\big \rangle _\alpha := \langle \widetilde{U},\widehat{U}\rangle _\alpha +\int _\Sigma \widetilde{V}\,\widehat{V}\;\!\mathrm {d}x + \int _\Omega \widetilde{v}\,\widehat{v}\;\!\mathrm {d}z\;\!\mathrm {d}x\).

For \(w=(U,V,v)\in D(A_\varepsilon )\) and \(\alpha \ge 0\), we obtain the estimate

$$\begin{aligned} {\begin{aligned} \big \langle \!\big \langle A_\varepsilon w,w \big \rangle \!\big \rangle _\alpha&:=\langle V,U\rangle _\alpha + \int _\Sigma \big (\Delta _xU{-}\partial _z v|_{z=0}\big ) V \;\!\mathrm {d}x + \int _\Omega \big ( \varepsilon ^2 \Delta _xv{+}\partial _z^2 v\big ) v \;\!\mathrm {d}z \;\!\mathrm {d}x\\&= \int _\Sigma \Big \{ \alpha ^2 UV+\nabla U\cdot \nabla V - \nabla U\cdot \nabla V - (\partial _z v\, v)|_{z=0}\\& - \int _{-\infty }^0\big (\varepsilon ^2 |\nabla _{x}v|^2 - (\partial _z v)^2\big ) \;\!\mathrm {d}z + (\partial _z v\, v)|_{z=0} \Big \} \;\!\mathrm {d}x \\&= \int _\Sigma \alpha ^2 UV \;\!\mathrm {d}x - \int _\Omega \big (\varepsilon ^2 |\nabla _{x}v|^2 - (\partial _z v)^2\big ) \;\!\mathrm {d}z \;\!\mathrm {d}x \\&\le \frac{\alpha }{2} \big ( \alpha ^2\Vert U\Vert _2^2 + \Vert V\Vert _2^2 \big ) \le \frac{\alpha }{2} \big \langle \!\big \langle w,w \big \rangle \!\big \rangle _\alpha = \frac{\alpha }{2} \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}w\varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_\alpha ^2, \end{aligned}} \end{aligned}$$
(3.5)

where we used the norm \(\varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}w\varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_\alpha = \big \langle \!\big \langle w,w\big \rangle \!\big \rangle ^{1/2}_\alpha \). For \(\lambda >0\) and \(F=(A_\varepsilon {-}\lambda I)w\) we obtain the estimate

$$\begin{aligned} \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}F \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_\alpha \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}w \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_\alpha&\ge -\big \langle \!\big \langle F,w \big \rangle \!\big \rangle _\alpha = - \big \langle \!\big \langle (A_\varepsilon {-}\lambda I) w,w \big \rangle \!\big \rangle _\alpha \\&\ge \big ( \lambda - \frac{\alpha }{2} \big ) \big \langle \!\big \langle w,w \big \rangle \!\big \rangle _\alpha =\big ( \lambda - \frac{\alpha }{2} \big ) \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}w\varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_\alpha ^2. \end{aligned}$$

Thus, for \(\alpha >0\) we have established the estimate

$$\begin{aligned} \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}(A_\varepsilon {-}\lambda I)^{-1} F\varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_\alpha \le \frac{1}{\lambda -\alpha /2} \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}F\varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_\alpha \quad \text {for all } \lambda >\alpha /2. \end{aligned}$$
(3.6)

Because of \(\varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}\cdot \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_{\alpha =1}=\Vert \cdot \Vert _{{\mathbf{H}}}\), we obtain

In particular, we have shown that the bounded linear operators \(A_\varepsilon {-}\lambda I: D(A_\varepsilon ) \rightarrow {\mathbf{H}}\) are injective. The following steps show that these operators are also surjective, i.e., the resolvent equations have a solution in \(D(A_\varepsilon )\).

Step 2: Reduction of resolvent equation. It remains to show that for all \(\lambda >0\) the resolvent equation \((A_\varepsilon {-} \lambda I)w=F \in {\mathbf{H}}\) has a solution in \(w=(U,V,v)\in D(A_\varepsilon )\). In this step, we reduce the problem to an equation for U alone.

Writing \(F=(G,H,f)\) the system reads

$$\begin{aligned} \text {on }\Sigma :\quad&V-\lambda U= G,\quad \Delta _x U-\partial _z v|_{z=0} - \lambda V=H, \quad V = v|_{z=0}, \end{aligned}$$
(3.7a)
$$\begin{aligned} \text {in } \Omega :\quad&\varepsilon ^2 \Delta _x v + \partial _z^2 v - \lambda v = f. \end{aligned}$$
(3.7b)

Obviously, we can eliminate V using the first equation giving \(V=G+\lambda U\).

Next, we solve (3.7b) for v. Together with the Dirichlet boundary condition \(v=V\) at \(z=0\), we obtain a unique solution \(v= {\textsf {V}}_\lambda ^\varepsilon (f,V)\). By classical elliptic regularity theory (see [14]), for all \(\lambda >0\) the linear operator \({\textsf {V}}_\lambda ^\varepsilon \) maps \({\mathrm H}^s(\Omega )\times {\mathrm H}^{s+3/2}(\Sigma )\) boundedly into \({\mathrm H}^{s+2}(\Omega )\). Because we have \(V \in {\mathrm H}^1(\Sigma )\) and \(f \in {\mathrm L}^2(\Omega )\), we treat the two inhomogeneities separately, namely \({\textsf {V}}_\lambda ^\varepsilon (\cdot , 0): {\mathrm L}^2(\Omega )\rightarrow Y^\varepsilon (\Omega ) \subset {\mathrm H}^2(\Omega )\) and \({\textsf {V}}_\lambda ^\varepsilon (0,\cdot ): {\mathrm H}^1(\Sigma ) \rightarrow X_\lambda ^\varepsilon ( \Omega ) \subset {\mathrm H}^{3/2}(\Omega )\). As the equation for U uses \({\partial _z v|_{z=0}}\), we define the bounded, linear operators

$$\begin{aligned} M_\lambda ^\varepsilon :\left\{ \begin{array}{ccc} {\mathrm L}^2(\Omega ) &{}\rightarrow &{} H^{1/2}(\Sigma ),\\ f&{} \mapsto &{} \partial _z {\textsf {V}}_\lambda ^\varepsilon (f,0)|_{z=0}, \end{array}\right. \quad \ \ N_\lambda ^\varepsilon :\left\{ \begin{array}{ccc} {\mathrm H}^1(\Sigma ) &{}\rightarrow &{} {\mathrm L}^2(\Sigma ),\\ V &{} \mapsto &{} \partial _z {\textsf {V}}_\lambda ^\varepsilon (0,V)|_{z=0}, \end{array}\right. \end{aligned}$$
(3.8)

where \(N^\varepsilon _\lambda \) is a (scaled version of the) Dirichlet-to-Neumann operator that maps \({\mathrm H}^{s+1}(\Sigma ) \) continuously into \({\mathrm H}^s(\Sigma )\) for all \(s \in \mathbb {R}\). This can be seen either by the abstract result in [4, Thm. 1.5(c)] or, in our special case with \(\Omega =\Sigma \times {]{-}\infty ,0[}\), by observing the explicit representation \(N_\lambda ^\varepsilon = (\lambda {-} \varepsilon ^2\Delta _x )^{1/2}\) with \(\lambda ,\varepsilon >0\). Choosing \(s=0\), we obtain the case stated in (3.8). Hence, it remains to solve an equation for U alone, namely

(3.9)

Step 3. \( (A_\varepsilon {-}\lambda I)^{-1} F \in D(A_\varepsilon )\). Using \(F=(G,H,f)\in {\mathbf{H}}\) and the mapping properties of \(M_\lambda ^\varepsilon \) and \(N_\lambda ^\varepsilon \), we see that the right-hand side in (3.9) lies in \({\mathrm L}^2(\Sigma )\). Moreover, for \(\lambda >0 \), the operator on the left-hand side generates a bounded and coercive bilinear form on \({\mathrm H}^1(\Sigma )\), because \({\int _\Sigma U N_\lambda ^\varepsilon U \;\!\mathrm {d}x = \int _\Omega \big ( \varepsilon ^2|\nabla _x v|^2 + (\partial _z v)^2 + \lambda v^2\big ) \;\!\mathrm {d}z \;\!\mathrm {d}x \ge 0}\), where \(v= {\textsf {V}}^\varepsilon _\lambda (0,U)\). This is of course the same calculation as in Step 1. Thus, the Lax–Milgram theorem provides a unique solution \(U \in {\mathrm H}^1(\Sigma )\), which by classical linear regularity lies even in \({\mathrm H}^2(\Sigma )\). From \(V=G+\lambda U\), we obtain \(V \in {\mathrm H}^1(\Sigma )\). Finally, we obtain \(v= {\textsf {V}}_\lambda ^\varepsilon (f,V) \in X^\varepsilon _\lambda (\Omega ) \cup Y^\varepsilon (\Omega )\).

Thus, we are done, if the identity (3.2) is established. For this, we take any \(W\in {\mathrm H}^1(\Sigma )\) and consider \(v_\lambda :={\textsf {V}}_\lambda ^\varepsilon (0,W) \in X_\lambda ^\varepsilon (\Omega )\) and \(v_1:= {\textsf {V}}^\varepsilon _1 (0,W) \in X^\varepsilon _1(\Omega )\). By the definition of \({\textsf {V}}_\lambda ^\varepsilon (0,\cdot )\), we see that the difference \(w:=v_\lambda -v_1\) satisfies the linear PDE

$$\begin{aligned} \varepsilon ^2 \Delta _x w+ \partial _z^2 w - w = (1{-}\lambda ) v_\lambda \in {\mathrm H}^{3/2}(\Omega ), \quad w|_{z=0} = 0. \end{aligned}$$

Hence, we conclude \(w \in Y^\varepsilon (\Omega )\), which implies \(v_\lambda = v_1 + w \in X^\varepsilon _1(\Omega ) \cup Y^\varepsilon (\Omega )\) as desired.

Step 4. The case \(\varepsilon =0\). The a priori estimate in Step 1 works for this case, too. The elimination of V and v works similarly, but now with the simplification that \(N_\lambda \) is explicitly given, namely \(N^0_\lambda = \sqrt{\lambda }I\). Together with \(X_\lambda ^0(\Omega )\cup Y^0(\Omega ) = X_1^0(\Omega )\cup Y^0(\Omega ) \) (see above) we see that all results in Steps 1 to 3 also hold for \(\varepsilon =0\).

Step 5. Generation of a semigroup. Steps 1 to 4 show that for all \(\varepsilon \ge 0\) the shifted operator \(A_\varepsilon {-}\frac{1}{2} I\) generates a contraction semigroup on \({\mathbf{H}}\) (cf. [19, Thm. 3.1]) with \(\Vert \mathrm {e}^{t(A_\varepsilon -\frac{1}{2}I)}\Vert _{{\mathbf{H}}\rightarrow {\mathbf{H}}} \le 1\). Thus, \(A_\varepsilon \) generates a semigroup satisfying \(\Vert \mathrm {e}^{t A_\varepsilon }\Vert _{{\mathbf{H}}\rightarrow {\mathbf{H}}} \le \mathrm {e}^{t/2}\).

Step 6. Growth rates for the semigroup. From (3.6), we know that the semigroups \(\mathrm {e}^{t A_\varepsilon }\) satisfy the growth estimate \( \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}\mathrm {e}^{t A_\varepsilon } w \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_\alpha \le \mathrm {e}^{t \alpha /2} \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}w \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_\alpha \) for all \(\alpha \ge 0\). Setting \({\underline{\alpha }}=\min \{1,\alpha \}\) and \(\overline{\alpha }= \max \{1,\alpha \}\) and using the equivalence between \(|\cdot |_\alpha \) and \(\Vert \cdot \Vert _{\mathrm{H^1}} = |\cdot |_1\), we obtain

$$\begin{aligned} \Vert \mathrm {e}^{t A_\varepsilon } w\Vert _{\mathbf{H}}&\le \frac{1}{{\underline{\alpha }}} \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}\mathrm {e}^{t A_\varepsilon } w \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_\alpha \le \mathrm {e}^{t\alpha /2}\,\frac{1}{{\underline{\alpha }}}\, \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}w \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_\alpha \\&\le \mathrm {e}^{t\alpha /2}\,\frac{\overline{\alpha }}{{\underline{\alpha }}} \, \Vert w\Vert _{{\mathbf{H}}} = \max \{ \alpha , 1/\alpha \} \mathrm {e}^{t\alpha /2} \Vert w\Vert _{{\mathbf{H}}}. \end{aligned}$$

Optimizing with respect to \(\alpha >0\) yields the bound \(\mathrm {e}^{t/2}\) for \(t\in [0,2]\) and \(t\,\mathrm {e}/2\) for \(t\ge 2\), which implies the final result \( \Vert \mathrm {e}^{t A_\varepsilon } \Vert _{{\mathbf{H}}\rightarrow {\mathbf{H}}} \le (1+t\,\mathrm {e}/2)\).

The final statement concerning \({\mathbb {E}}_0\) follows by setting \(\alpha =0\), observing \({\mathbb {E}}_0(w)=\frac{1}{2}\varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}w\varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_0^2\), and the contraction property \( \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}\mathrm {e}^{t A_\varepsilon } w \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_0 \le \mathrm {e}^{0\cdot t} \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}w \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_0 = \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}w\varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_0\). Hence, Theorem 3.1 is established. \(\square \)

3.2 Convergence of semigroups

The next result proves the convergence of the resolvents \((A_\varepsilon {-}\lambda I)^{-1}\) as operators from \({\mathbf{H}}\) into itself in the strong operator topology. The critical point is to understand the convergence of the Dirichlet-to-Neumann operators \(N_\lambda ^\varepsilon \) to the limiting operator \(N_\lambda ^0\), see (3.8).

Proposition 3.2

(Strong convergence of resolvents). For all \(\lambda >0\) and all \(F \in {\mathbf{H}}\), we have the strong convergence \((A_\varepsilon {-}\lambda I)^{-1} F \rightarrow (A_0{-}\lambda I)^{-1} F\).

Proof

Throughout the proof, \(\lambda >0\) is fixed.

Step 1. Reduction to F in a dense subset \({\mathbf{Z}}\) of \({\mathbf{H}}\). Let \({\mathbf{Z}}\subset {\mathbf{H}}\) be given such that \({\mathbf{Z}}\) is dense in \({\mathbf{H}}\) and that for all \(F\in {\mathbf{Z}}\) we have \((A_\varepsilon {-} \lambda I)^{-1}F \rightarrow (A_0 {-} \lambda I)^{-1}F\) as \(\varepsilon \rightarrow 0^+\).

For an arbitrary \(F\in {\mathbf{H}}\), we consider \(F_n \in {\mathbf{Z}}\) with \(F_n \rightarrow F\) in \({\mathbf{H}}\) as \( n \rightarrow \infty \). By Step 1 in the proof of Theorem 3.1, we know that the resolvents \((A_\varepsilon {-} \lambda I)^{-1}\) are uniformly bounded by \(C_\lambda \) with respect to \(\varepsilon >0\). Hence, we have

$$\begin{aligned}&\big \Vert (A_\varepsilon {-} \lambda I)^{-1}F - (A_0 {-} \lambda I)^{-1}F \big \Vert _{\mathbf{H}}\\&\le \big \Vert (A_\varepsilon {-} \lambda I)^{-1}( F {-}F_n) \big \Vert _{\mathbf{H}}+\big \Vert (A_\varepsilon {-} \lambda I)^{-1}F_n - (A_0 {-} \lambda I)^{-1}F_n \big \Vert _{\mathbf{H}}\\&\quad +\big \Vert (A_0 {-} \lambda I)^{-1} (F_n {-} F) \big \Vert _{\mathbf{H}}\\&\le C_\lambda \Vert F{-}F_n\Vert _{\mathbf{H}}+ \big \Vert (A_\varepsilon {-} \lambda I)^{-1}F_n - (A_0 {-} \lambda I)^{-1}F_n \big \Vert _{\mathbf{H}}+ C_\lambda \Vert F{-}F_n\Vert _{\mathbf{H}}. \end{aligned}$$

Thus, for a given \(\delta >0\) we can make the difference small by first choosing n so big that \(C_\lambda \Vert F{-}F_n\Vert _{\mathbf{H}}< \delta /3\) and then choosing \(\varepsilon _0 >0\) so small that the middle term is less than \(\delta /3\) for all \(\varepsilon \in {]0,\varepsilon _0[}\) as well. Thus, for all \(F \in {\mathbf{H}}\) we have the convergence \((A_\varepsilon {-} \lambda I)^{-1}F \rightarrow (A_0 {-} \lambda I)^{-1}F\).

Step 2. Higher regularity for smooth right-hand sides F. We use that the system is translation invariant in the domain \(\Sigma \). Thus, if the partial derivatives lie in \({\mathbf{H}}\), then the solutions \(w^\varepsilon = (A_\varepsilon {-}\lambda I)^{-1} F\) have an additional derivative in \(x_j\) direction as well and satisfy the a priori estimate . Thus, for F in the dense subset

$$\begin{aligned} {\mathbf{Z}}= \big \{\, F\in {\mathbf{H}} \, \big | \, \nabla _x F \in {\mathbf{H}} \,\big \} \end{aligned}$$
(3.10)

the improved estimate \(\Vert w^\varepsilon \Vert _{\mathbf{Z}}\le C_\lambda \Vert F\Vert _{\mathbf{Z}}\) holds with \(\Vert F\Vert _{\mathbf{Z}}=\Vert F\Vert _{\mathbf{H}}{+} \Vert \nabla _x F\Vert _{\mathbf{H}}\).

Step 3. Convergence for \(\varepsilon \rightarrow 0^+\). We now assume \(F\in {\mathbf{Z}}\) and compare \(w^\varepsilon =(U^\varepsilon ,V^\varepsilon ,v^\varepsilon )= (A_\varepsilon {-}\lambda I)^{-1} F\) with \(w^0=(U^0,V^0,v^0) = (A_0{-}\lambda I)^{-1}F\). As in Step 1 of the proof of Theorem 3.1, we estimate the difference \(w^\varepsilon - w^0\) as follows (choosing \(\alpha >0\) with \(\alpha <2\lambda \)):

$$\begin{aligned}&\big (\lambda {-}\frac{\alpha }{2}\big ) \big \langle \!\big \langle w^\varepsilon {-}w^0, w^\varepsilon {-}w^0 \big \rangle \!\big \rangle _\alpha \le - \big \langle \!\big \langle (A_0{-}\lambda I) (w^\varepsilon {-}w^0), w^\varepsilon {-}w^0 \big \rangle \!\big \rangle _\alpha \\&\overset{*}{=} \big \langle \!\big \langle (0,0,\varepsilon ^2 \Delta _x v^\varepsilon )^\top , w^\varepsilon {-}w^0 \big \rangle \!\big \rangle _\alpha = - \int _\Omega \varepsilon ^2 \nabla _x v^\varepsilon ( \nabla v^\varepsilon {-}\nabla v^0) \;\!\mathrm {d}z \;\!\mathrm {d}x \\&\le \widehat{C}_\alpha \varepsilon ^2 \Vert w^\varepsilon \Vert _{\mathbf{Z}}(\Vert w^\varepsilon \Vert _{\mathbf{Z}}{+}\Vert w^0\Vert _{\mathbf{Z}}) \le 2\varepsilon ^2 \widehat{C}_\alpha C_\lambda \Vert F\Vert _{\mathbf{Z}}^2 \ \rightarrow 0 \ \text { \ as } \varepsilon \rightarrow 0^+\,. \end{aligned}$$

In the identity \(\overset{*}{=}\), we used the cancellation arising from \((A_0{-}\lambda I)w^0=F\) and

$$\begin{aligned} (A_0{-}\lambda I) w^\varepsilon = (A_\varepsilon {-}\lambda I) w^\varepsilon + (A_0{-}A_\varepsilon ) w^\varepsilon = F -(0,0,\varepsilon ^2 \Delta _x v^\varepsilon )^\top . \end{aligned}$$

By the equivalence of the \({\mathbf{H}}\) norm and the norm induced by \(\big \langle \!\big \langle \cdot , \cdot \big \rangle \!\big \rangle _\alpha \), we conclude \(\Vert w^\varepsilon {-}w^0\Vert _{\mathbf{H}}\) \( \le C \,\varepsilon \), and Proposition 3.2 is proved. \(\square \)

Theorem 3.1 and Proposition 3.2 are the basis for the following result that states that the strongly continuous semigroups \((\mathrm {e}^{t A_\varepsilon })_{t\ge 0}\) on \({\mathbf{H}}\) converge as \(\varepsilon \rightarrow 0^+\) in the strong operator topology. Indeed, the proof of the first part is a direct consequence of the Trotter–Kato theory, see [19, Sec. 3.3], while the second part uses explicit estimates.

Theorem 3.3

(Strong convergence of the solutions). Consider the operators \(A_\varepsilon \) defined in Theorem 3.1 and the induced semigroups \( \mathrm {e}^{t A_\varepsilon }: {\mathbf{H}}\rightarrow {\mathbf{H}}\) for \(t\ge 0\). Then, for all initial conditions \(w_0 \in {\mathbf{H}}\), the solutions \(w^\varepsilon : {[0,\infty [} \rightarrow {\mathbf{H}},\ t\mapsto w^\varepsilon (t)=\mathrm {e}^{t A_\varepsilon } w_0\) satisfy for all \(t\ge 0\) the convergence \(w^\varepsilon (t) \rightarrow w^0(t)\) as \(\varepsilon \rightarrow 0\).

Moreover, for initial conditions with additional derivatives in x-direction, namely \(w_0 \in {\mathbf{Z}}\) (cf. (3.10)) we have the quantitative error estimate

$$\begin{aligned} \Vert w^\varepsilon (t) - w^0(t)\Vert _{\mathbf{H}}\le \varepsilon \,\sqrt{t}\, (2.3{+}t)^2 \,\Vert w_0\Vert _{\mathbf{Z}}\quad \text {for all }t\ge 0. \end{aligned}$$
(3.11)

Proof

It remains to show (3.11). For this, we set \(\delta = w^\varepsilon -w^0\) and perform a simple energy estimate, where we use that \(w^\varepsilon =(U^\varepsilon ,V^\varepsilon ,v^\varepsilon )\) and \(w^0=(U^0,V^0,v^0)\) are sufficiently smooth solutions of (3.1), because we have the extra regularity of \(w_0\in {\mathbf{Z}}\). We employ the norms \(\varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}\cdot \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_\alpha \) from (3.5) and find

where we used \(\Vert \nabla _x w^\varepsilon \Vert _{\mathbf{H}}\le \Vert \mathrm {e}^{t A_\varepsilon }\nabla _x w_0\Vert _{\mathbf{H}}\le \Vert \mathrm {e}^{t A_\varepsilon }\Vert _{{\mathbf{H}}\rightarrow {\mathbf{H}}} \Vert w_0\Vert _{\mathbf{Z}}\) and the growth estimate from Theorem 3.1. With \(\delta (0)=w_0{-}w_0=0\), Gronwall’s lemma yields

$$\begin{aligned} \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}\delta (t) \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_\alpha ^2 \le \varepsilon ^2 \int _0^t 2\, \mathrm {e}^{\alpha (t{-}s)} (1{+}s\,\mathrm {e}/2)^2 \;\!\mathrm {d}s \,\Vert w_0\Vert _{\mathbf{Z}}^2. \end{aligned}$$

For \(t\in [0,2]\), we choose \(\alpha =1\) and obtain

$$\begin{aligned} \Vert \delta (t)\Vert _{\mathbf{H}}^2 = \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}\delta (t) \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_1^2 \le C_* t\,\varepsilon ^2 \Vert w_0\Vert _{\mathbf{Z}}^2 \text { for }t\in [0,2] \end{aligned}$$

where \(C_*=\int _0^2 \mathrm {e}^{2-s}(1{+}s\mathrm {e}/2)^2 \;\!\mathrm {d}s\approx 27.14...\le 2.3^4\). For \(t\ge 2\), let \(\alpha =2/t\le 1\) to obtain

$$\begin{aligned} \Vert \delta (t)\Vert _{\mathbf{H}}^2&\le \frac{1}{\alpha ^2} \varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}\delta (t)\varvec{|}\!\!\,\varvec{|}\!\!\,\varvec{|}_\alpha ^2 \le \frac{t^2}{4} \int _0^t 2\,\mathrm {e}^{2(1-s/t)} (1{+}s\mathrm {e}/2)^2 \;\!\mathrm {d}s \, \varepsilon \Vert w_0\Vert _{\mathbf{Z}}^2\\&=\frac{t^3}{32}\Big ((\mathrm {e}^2{-}5)\mathrm {e}^2t^2+4(\mathrm {e}^2{-}3)\mathrm {e}t+8(\mathrm {e}^2{-}1) \Big )\, \varepsilon ^2 \Vert w_0\Vert _{\mathbf{Z}}^2\\&\le \frac{t}{6}\big ( 1+ t\,\mathrm {e}/2)^4 \, \varepsilon ^2 \Vert w_0\Vert _{\mathbf{Z}}^2. \end{aligned}$$

Combining this with the result for \(t\in [0,2]\) and \((\mathrm {e}/2)^4\le 6\), we find \(\Vert \delta (t)\Vert _{\mathbf{H}}^2 \le t\, (2.3{+}t)^4 \,\varepsilon ^2 \Vert w_0\Vert _{\mathbf{Z}}^2\) for all \(t\ge 0\), which is the desired result (3.11). \(\square \)

4 Energy and dissipation functionals

We now show that the fractionally damped wave equation (2.4) carries a natural energy–dissipation structure. This is done in two different ways. First, we reduce the natural energy–dissipation structure of the limiting system (3.1) with \(\varepsilon =0\) by eliminating the diffusion equation. For this, we first study the one-dimensional diffusion equation on the half line \({]{-}\infty ,0[}\) in detail. Second, we show by a direct calculation that the energy–dissipation structure extends to a more general class of fractionally damped wave equations, where the time derivative of order 3/2 is replaced by order \(1{+}\alpha \) with \(\alpha \in {]0,1[}\).

4.1 Diffusion equation on the half line

We consider the following initial boundary value problem on \(I:={]{-}\infty , 0[}\):

(4.1)

We always assume the compatibility condition \(v_0(0)=\varphi (0)\).

Using the one-dimensional heat kernel \(H(t,y) = (4\pi t)^{-1/2} \exp \!\big ({-}y^2/(4t)\big )\) and the reflection principle, the influence of \(v_0\) is described via \(K_\mathrm {Dir}(t,z,y)= H(t,z{-}y)-H(t,z{+}y)\) giving homogeneous Dirichlet data via \( K_\mathrm {Dir}(t,0,y)=0\):

$$\begin{aligned} v(t,z)= \int _I K_\mathrm {Dir} (t,z,y) v_0(y) \;\!\mathrm {d}y. \end{aligned}$$

To obtain the influence of the inhomogeneous Dirichlet data at \(z=0\), we set \(v_0=0\) and make the ansatz \(v(t,z)= w(t,z)+\varphi (t)\) such that w has to satisfy

$$\begin{aligned} w_t = w_{zz} -\dot{\varphi }(t), \qquad w(t,0)=0, \quad w(0,z)=-\varphi (0)=0. \end{aligned}$$

With Duhamel’s principle (variation-of-constants formula), we obtain

$$\begin{aligned} w(t,z)= -\int _0^t \int _I K_{\mathrm D}(t{-}s,z,y) \dot{\varphi }(s)\;\!\mathrm {d}y \;\!\mathrm {d}s. \end{aligned}$$

Setting \(G(y):=\int _{-\infty }^y H(1,\eta ) \;\!\mathrm {d}\eta \) (such that \(G(-\infty )=0\) and \(G(\infty )=1\)), we obtain

$$\begin{aligned} w(t,z)= \int _0^t \dot{\varphi }(s) \big (G(z/\sqrt{t{-}s}) - G({-}z/\sqrt{t{-} s})\big ) \;\!\mathrm {d}s. \end{aligned}$$

Putting both cases together, the full solution formula for (4.1) reads

$$\begin{aligned} v(t,z) = \int _I K_\mathrm {Dir}(t,z,y) v_0(y) \;\!\mathrm {d}y +\varphi (t) +\int _0^t \dot{\varphi }(s) \Big (G\big (\tfrac{\displaystyle z}{\sqrt{t{-}s}}\big ) - G\big (\tfrac{\displaystyle -z}{\sqrt{t{-}s}}\big )\Big ) \;\!\mathrm {d}s. \end{aligned}$$

For the analysis related to the fractionally damped wave equation, we consider only the case \(v_0 \equiv 0\), which implies \(\varphi (0)=0\) as well by continuity of the boundary-initial data. Doing integration by parts for the time integral and using

we arrive, for the case \(v_0\equiv 0\), at the relation

(4.2)

Using , doing another integration by parts, and using \(\varphi (0)=0\) again, we find

(4.3)

In particular, evaluating at \(z=0\), where \(H(t,0)=1/\sqrt{4 \pi t}\), we find

(4.4)

According to the definition (1.6), the boundary derivative is the fractional Caputo derivative of order 1/2 of \(\varphi \), i.e., .

We now derive an energy–dissipation balance for the diffusion equation by rewriting the natural \({\mathrm L}^2\) integrals in terms of the boundary value \(\varphi \). The starting point is the classical relation

(4.5)

For solutions v of (4.1) with \(v_0\equiv 0\) and \(\varphi (0)=0\), we can rewrite this energy–dissipation balance totally in terms of \(\varphi \) by using the following result.

Proposition 4.1

Assume that v is given via (4.2) and by (4.3), then we can express twice the energy \(\int _I v^2 \;\!\mathrm {d}z\) and the dissipation via

$$\begin{aligned} \int _I v^2 \;\!\mathrm {d}z \quad&= \int _0^t \int _0^t M_0(t{-}r,t{-}s)\,\varphi (r)\varphi (s) \;\!\mathrm {d}r \;\!\mathrm {d}s \quad \text {and} \end{aligned}$$
(4.6a)
(4.6b)

where

$$\begin{aligned} M_j(r,s)&= \int _I K_j(r,z)K_j(s,z) \;\!\mathrm {d}z = \frac{2^j}{\sqrt{4\pi }\,(r{+}s)^{3/2-j}}. \end{aligned}$$
(4.6c)

In particular, \(M_j(r,s)=\widetilde{M}_j(r{+}s)\) and .

Proof

The relations (4.6a) and (4.6b) with \(M_j(r,s)= \int _I K_j(r,z)K_j(s,z) \;\!\mathrm {d}z \) follow simply by the definitions. Thus, it remains to establish the explicit formulas for \(M_0\) and \(M_1\) by exploiting the structure of \(K_0= 2 \partial _z H = -\frac{z}{2t}\, 2H \) and \(K_1 = 2H\). We obtain

$$\begin{aligned} K_j(r,z)K_j(s,z)&= 4 \Big ( \frac{z^2}{4rs}\Big )^{1-j}\frac{1}{4\pi \sqrt{rs}} \exp \Big ( {-}\frac{z^2}{4r} - \frac{z^2}{4s}\Big )\\&= \Big ( \frac{z^2}{4rs}\Big )^{1-j} \frac{1}{\pi \sqrt{rs}} \exp \Big ( {-}\frac{r{+}s}{4rs}\,z^2\Big ). \end{aligned}$$

An explicit integration with \(\int _0^\infty \mathrm {e}^{-a^2 x^2} \;\!\mathrm {d}x = \sqrt{\pi }/(2a)\) and \(\int _0^\infty x^2\mathrm {e}^{-a^2 x^2} \;\!\mathrm {d}x = \sqrt{\pi }/(4a^3)\) yields the stated formulas for \(M_0\) and \(M_1\). \(\square \)

It is a surprising fact that \(M_0\) and \(M_1\) depend only on the sum \(r{+}s\) rather than on r and s individually. The relations in (4.6) allow us to rewrite the energy–dissipation balance (4.5) in terms of \(\varphi \) alone. We obtain the identity

$$\begin{aligned}&\frac{{\mathrm d}}{{\mathrm d}t} \int _0^t\int _0^t \frac{1}{2} M_0(t{-}r,t{-}s)\, \varphi (s) \varphi (r) \;\!\mathrm {d}r \;\!\mathrm {d}s\nonumber \\&= \varphi (t) \int _0^t \frac{\dot{\varphi }(\tau )}{\sqrt{\pi (t{-}\tau )}} \;\!\mathrm {d}\tau - \int _0^t\int _0^t M_1(t{-}r,t{-}s)\, \dot{\varphi }(s) \dot{\varphi }(r) \;\!\mathrm {d}r \;\!\mathrm {d}s\,. \end{aligned}$$
(4.7)

4.2 An energy–dissipation relation for fractional derivatives

Here we show that the identity (4.7) can be derived in an independent way, not using the Dirichlet-to-Neumann map for the one-dimensional diffusion equation. We even generalize the result to the case of general fractional derivatives \({}^{\mathrm C}{\mathrm D}^\alpha \varphi \), where (4.7) is the special case \(\alpha =1/2\). For general \(\alpha \in {]0,1[}\), we set

$$\begin{aligned} N_0^\alpha (r,s)= \frac{\alpha }{\Gamma (1{-}\alpha ) \,(r{+}s)^{1+\alpha }} \quad {\text {and}} \quad N_1^\alpha (r,s)= \frac{1}{\Gamma (1{-}\alpha )\,(r{+}s)^{\alpha }}. \end{aligned}$$
(4.8)

With this, we obtain the following result.

Proposition 4.2

For all \(\varphi \in {\mathrm C}^1({[0,\infty [})\) with \(\varphi (0)=0\), we have the identity

$$\begin{aligned}&\frac{{\mathrm d}}{{\mathrm d}t} \int _0^t\!\!\int _0^t \frac{1}{2} N^\alpha _0(t{-}r,t{-}s)\, \varphi (s) \varphi (r) \;\!\mathrm {d}r \;\!\mathrm {d}s\\&= \varphi (t)\, {}^{\mathrm C}{\mathrm D}^\alpha \varphi (t) - \int _0^t\!\!\int _0^t N^\alpha _1(t{-}r,t{-}s)\, \dot{\varphi }(s) \dot{\varphi }(r) \;\!\mathrm {d}r \;\!\mathrm {d}s. \end{aligned}$$

Proof

We set \(r=t{-}\rho \) and \(s=t{-}\sigma \) and obtain

$$\begin{aligned}&\frac{{\mathrm d}}{{\mathrm d}t} \int _0^t\int _0^t \frac{1}{2} N^\alpha _0(t{-}r,t{-}s)\, \varphi (s) \varphi (r) \;\!\mathrm {d}s \;\!\mathrm {d}r \\&= \frac{{\mathrm d}}{{\mathrm d}t} \int _0^t\int _0^t \frac{1}{2} N^\alpha _0(\rho ,\sigma )\, \varphi (t{-}\sigma ) \varphi (t{-}\rho ) \;\!\mathrm {d}\sigma \;\!\mathrm {d}\rho \\&\overset{1}{=}\int _{\rho =0}^t \int _{\sigma =0}^t \frac{1}{2} N^\alpha _0(\rho ,\sigma ) \big ( \varphi (t{-}\sigma ) \dot{\varphi }(t{-}\rho ) + \big ( \dot{\varphi }(t{-}\sigma ) \varphi (t{-}\rho ) \big ) \;\!\mathrm {d}\sigma \;\!\mathrm {d}\rho .\\&\overset{2}{=} \int _{\rho =0}^t \int _{\sigma =0}^t N^\alpha _0(\rho ,\sigma ) \varphi (t{-}\sigma ) \dot{\varphi }(t{-}\rho ) \;\!\mathrm {d}\sigma \;\!\mathrm {d}\rho . \end{aligned}$$

Here \(\overset{1}{=}\) uses \(\varphi (0)=0\) such that the boundary terms arising from \(\frac{{\mathrm d}}{{\mathrm d}t} \int _0^t g(\tau )\;\!\mathrm {d}\tau = g(t)\) vanish. In \(\overset{2}{=} \), we simply use the symmetry \(N^\alpha _0(r,s)=N^\alpha _0(s,r)\).

In the next step, we perform an integration by parts with respect to \(\sigma \in [0,t]\) and use the fundamental relation \(N^\alpha _0(\rho ,\sigma ) = - \partial _\sigma N^\alpha _1(\rho ,\sigma )\). Hence, we continue

$$\begin{aligned}&= \int _{\rho =0}^t \Big \{ \big [{-}N^\alpha _1(\rho ,\sigma )\varphi (t{-}\sigma )\big ]\Big |_0^t -\int _0^t \big ({-}N^\alpha _1(\rho ,\sigma )\big ) \big ({-}\dot{\varphi }(t{-}\sigma )\big ) \;\!\mathrm {d}\sigma \Big \} \dot{\varphi }(t{-}\rho ) \;\!\mathrm {d}\sigma \;\!\mathrm {d}\rho \\&=\int _{\rho =0}^t N^\alpha _1(\rho ,0)\varphi (t)\dot{\varphi }(t{-}\rho ) \;\!\mathrm {d}\sigma \;\!\mathrm {d}\rho -\int _0^t \int _0^t N^\alpha _1(\rho ,\sigma )\dot{\varphi }(t{-}\sigma ) \dot{\varphi }(t{-}\rho ) \;\!\mathrm {d}\sigma \;\!\mathrm {d}\rho , \end{aligned}$$

where we again have used \(\varphi (0){=}0\). The definition of the function \(N^\alpha _1\) gives \(N^\alpha _1(t{-}\tau ,0)= 1/ \big (\Gamma (1{-}\alpha ) (t{-}\tau )^\alpha \big ) \), such that the first term is indeed equal to \(\varphi \;({}^{\mathrm C}{\mathrm D}^\alpha \varphi )\). With this, the result is established. \(\square \)

We emphasize that the above result does not need the exact form of \(N_0^\alpha \) and \(N_1^\alpha \) as given in (4.8). We only exploited the relations

$$\begin{aligned} N^\alpha _0(r,s)= N^\alpha _0(s,r), \quad N^\alpha _0(r,s)= - \partial _s N^\alpha _1(r,s), \quad N^\alpha _1(r,0)=1/\big ( \Gamma (1{-}\alpha )\,r^\alpha \big ). \end{aligned}$$

Clearly, there are many more functions satisfying these conditions. However, we also want positive semi-definiteness of the kernels \(N^\alpha _j\), i.e.,

$$\begin{aligned} \int _0^t \int _0^t N^\alpha _j(r,s) \psi (s)\psi (r) \;\!\mathrm {d}s \;\!\mathrm {d}r \ge 0 \quad \text {for all }\psi \in {\mathrm C}^0({[0,\infty [}),\ t>0, {\text { and }}j\in \{0,1\}. \end{aligned}$$

For general \(N^\alpha _j\), this positive semi-definiteness is a significant restriction, but for our chosen cases it can be established as follows:

$$\begin{aligned} 0&\le \int _{y=0}^\infty \Big (\int _{r=0}^t \frac{y^{\alpha -1/2}}{r^\alpha } \mathrm {e}^{-y^2/r} \psi (r)\;\!\mathrm {d}r \Big )^2 \;\!\mathrm {d}y\\&= \int _{y=0}^\infty \int _{r=0}^t\int _{s=0}^t \frac{y^{\alpha -1/2}}{r^\alpha } \mathrm {e}^{-y^2/r} \psi (r) \frac{y^{\alpha -1/2}}{s^\alpha } \mathrm {e}^{-y^2/s} \psi (s) \;\!\mathrm {d}s \;\!\mathrm {d}r \;\!\mathrm {d}y\\&= \int _{r=0}^t\int _{s=0}^t \int _{y=0}^\infty \frac{y^{2\alpha -1}}{(rs)^\alpha } \mathrm {e}^{- y^2(r{+}s)/(rs)} \;\!\mathrm {d}y \; \psi (r)\psi (s) \;\!\mathrm {d}s \;\!\mathrm {d}r. \end{aligned}$$

Using \(\int _0^\infty |y|^{2\alpha -1} \mathrm {e}^{-by^2} \;\!\mathrm {d}y = \Gamma (\alpha )/(2 b^\alpha )\), we obtain the desired result

$$\begin{aligned} 0 \le \int _{r=0}^t\int _{s=0}^t \frac{\Gamma (\alpha )}{(r{+}s)^\alpha } \,\psi (r)\psi (s) \;\!\mathrm {d}s \;\!\mathrm {d}r \quad {\text {for all }}\;\psi \in {\mathrm C}([0,t]), \end{aligned}$$

which holds for all \(\alpha \ge 0\).

4.3 Energetics for the fractionally damped wave equation

To derive the physically relevant energy–dissipation balance for the fractionally damped wave equation

$$\begin{aligned} \ddot{U}(t,x) + \int _0^t \frac{\ddot{U}(s,x)}{\sqrt{\pi (t{-}s)}} \;\!\mathrm {d}s \ = \ \Delta U(t,x) \end{aligned}$$
(4.9)

we use the limiting system (2.2). The latter is a classical system of partial differential equations, and it is easy to write down the physically motivated energy functional \({\mathbb {E}}\) and the corresponding dissipation function \({\mathcal D}\).

The total energy \({\mathbb {E}}\) is the sum of the kinetic and potential energy in the membrane \(\Sigma \) plus the kinetic energy in the lower half space \(\Omega = \Sigma \times {]{-}\infty ,0[}\), where we consider v as the horizontal component of a shear flow. This leads to \({\mathbb {E}}_0\) as defined in (3.4), and along the solutions of (2.2) the energy–dissipation balance takes the form

$$\begin{aligned} \frac{{\mathrm d}}{{\mathrm d}t} {\mathbb {E}}_0(U(t),\dot{U}(t), v(t)) = -{\mathcal D}(U,\dot{U},v):= \int _\Omega \frac{1}{2} (\partial _zv)^2 \;\!\mathrm {d}z \;\!\mathrm {d}x. \end{aligned}$$

This shows that the only dissipation occurs by the (shear) viscosity of the fluid in the lower half space \(\Omega \).

As explained in Sects. 2.4 and 4.1, we can eliminate v via (using \(\varphi (t,x)=\dot{U}(t,x,0)\) and assuming \(v(0,x,z)=0\) of all \((x,z)\in \Omega \))

$$\begin{aligned} v(t,x,z)=\int _0^t 2 \partial _z H(t{-}\tau ,z) \dot{U}(\tau ,x) \;\!\mathrm {d}\tau , \quad {\text {where }} H(t,z)=\frac{ \mathrm {e}^{-z^2/(4t)}}{\sqrt{4\pi t}}. \end{aligned}$$

This allows us to eliminate \(\partial _z v(t,x,0)\) via (4.4) and we obtain the fractionally damped wave equation (4.9).

Using the formulas derived in Proposition 4.1, we obtain the reduced energy function \({\mathbf{E}}\) and the reduced dissipation function \({\mathbf{D}}\) in the form

$$\begin{aligned} {\mathbf{E}}(U(t),[\dot{U}]_{[0,t]})&:= \int _\Sigma \Big \{\frac{1}{2}\dot{U}(t)^2 + \frac{1}{2}|\nabla _xU(t)|^2 \end{aligned}$$
(4.10a)
$$\begin{aligned}&\qquad \quad \ + \int _0^t\!\!\int _0^t \frac{1}{4\sqrt{\pi }\,(2t{-}r{-}s)^{3/2}} \, \dot{U}(r,x) \dot{U}(s,x) \;\!\mathrm {d}r\;\!\mathrm {d}s \Big \} \;\!\mathrm {d}x,\nonumber \\ {\mathbf{D}}(U(t),[\dot{U}]_{[0,t]})&:= \int _\Sigma \Big \{ \int _0^t\!\!\int _0^t \frac{1}{\sqrt{\pi }\,(2t{-}r{-}s)^{1/2}}\, \ddot{U}(r,x) \ddot{U}(s,x) \;\!\mathrm {d}r\;\!\mathrm {d}s \Big \}\;\!\mathrm {d}x. \end{aligned}$$
(4.10b)

Clearly, along solutions of the fractionally damped wave equation (4.9) we have the reduced energy–dissipation balance

$$\begin{aligned} \frac{{\mathrm d}}{{\mathrm d}t} {\mathbf{E}}(U(t),[\dot{U}]_{[0,t]}) = - {\mathbf{D}}(U(t),[\dot{U}]_{[0,t]}). \end{aligned}$$
(4.11)

Of course, it is possible to check this identity directly without any reference to the limiting system (2.2) involving the hidden state variable v. For this, we do the standard argument for energy conservation for the wave equation plus the calculation in the proof of Proposition 4.2 for the parts non-local in time.

In the related works [8, 23, 25], other energy functionals were constructed for equations with fractional time derivatives. However, the approach there is quite different and is less inspired by the true energy and dissipation hidden in the eliminated state variable v.

Indeed, we may generalize the energy–dissipation balance (4.11) to the case of fractional damping of order \(1{+}\alpha \in {]1,2[}\). We consider (4.9) as a special case of the equation

$$\begin{aligned} \ddot{U} + {}^{\mathrm C}{\mathrm D}^\alpha \dot{U} = \Delta U \quad {\text { on }} \Sigma . \end{aligned}$$
(4.12)

Taking into account the calculations in Sect. 4.2, we define the energy \({\mathbf{E}}^\alpha \) and the dissipation function \({\mathbf{D}}^\alpha \) via

$$\begin{aligned} {\mathbf{E}}^\alpha (U(t),[\dot{U}]_{[0,t]})&:= \int _\Sigma \Big \{\frac{1}{2}\dot{U}(t)^2 + \frac{1}{2}|\nabla _xU(t)|^2 \nonumber \\&\qquad + \int _0^t\!\!\int _0^t\!\! \frac{\alpha /2}{\Gamma (1{-}\alpha )(2t{-}r{-}s)^{1{+}\alpha }} \, \dot{U}(r,x) \dot{U}(s,x) \;\!\mathrm {d}r\;\!\mathrm {d}s \Big \} \;\!\mathrm {d}x, \end{aligned}$$
(4.13a)
$$\begin{aligned} {\mathbf{D}}^\alpha (U(t),[\dot{U}]_{[0,t]})&:= \int _\Sigma \!\!\, \Big \{ \! \int _0^t\!\!\int _0^t \!\!\frac{1}{\Gamma (1{-}\alpha )(2t{-}r{-}s)^{\alpha }}\, \ddot{U}(r,x) \ddot{U}(s,x) \;\!\mathrm {d}r\;\!\mathrm {d}s \Big \}\;\!\mathrm {d}x. \end{aligned}$$
(4.13b)

Clearly, for sufficiently smooth solutions of the fractionally damped wave equation (4.12) we have the reduced energy–dissipation balance

$$\begin{aligned} \frac{{\mathrm d}}{{\mathrm d}t} {\mathbf{E}}^\alpha (U(t),[\dot{U}]_{[0,t]}) = - {\mathbf{D}}^\alpha (U(t),[\dot{U}]_{[0,t]}). \end{aligned}$$
(4.14)

5 Conclusion and outlook

In this work, we have shown that the fractionally damped wave equation can be obtained as a scaling limit from a bulk-interface coupling between a wave equation for a membrane and a viscous fluid motion in the adjacent half space. The coupling is such that the natural mechanical energies act as a Liapunov function. We have identified the physical scaling parameters like the equivalent membrane thickness \(\ell _\mathrm {thick}=\rho _{\mathrm{memb}}/\rho _{\mathrm{bulk}}\) for the vertical scaling and the effective travel length \(\ell _\mathrm {trav}(t_*) =\rho _{\mathrm{memb}}^{3/2}\kappa ^{1/2}/(\mu \rho _{\mathrm{bulk}})\) for the horizontal scaling. Thus, taking the limit \(\varepsilon \rightarrow 0\) in the critical parameter

$$\begin{aligned} \varepsilon = \frac{\ell _\mathrm {thick}}{\ell _\mathrm {trav}(t_*)} = \frac{\mu }{\sqrt{\rho _{\mathrm{memb}}\, \kappa \,}} \end{aligned}$$

leads to the appearance of the fractionally damped wave equation.

The first main outcome of the mathematical analysis is that the system is stable uniformly with respect to \(\varepsilon \) and that it converges strongly in the natural energy space \({\mathbf{H}}\) in the sense of linear semigroup theory. For initial data with higher horizontal regularity, a convergence rate could be derived. Thus, the fractional time derivative of order 3/2 appears naturally as a consequence of the Dirichlet-to-Neumann map of a one-dimensional parabolic equation on the half line.

The second outcome of our approach is the energy–dissipation structure for the fractionally damped wave equation which is derived by integrating out the “hidden states” v in the fluid layer in the full mechanical energy–dissipation structure of the coupled system of partial differential equations. As expected, we obtain quadratic functionals for the reduced energy and the reduced dissipation function that are non-local in time, thus keeping track of information stored in the hidden state variable v. It is surprising that both quadratic functionals obtained have memory kernels \(M_j(t{-}r,t{-}s)\) that depend only on the sum \((t{-}r)+(t{-}s)\). It is certainly important to understand where this special structure comes from and how it relates to more general energy–dissipation structures as introduced in [8, 23, 25].

A major restriction occurs through our assumption \(v(0,x,z)=0\) for a.a. \((x,z) \in \Omega \), which implies \(\dot{U}(0,x)=0\). We expect that this assumption can be avoided by suitably generalizing the Caputo derivative and by extending the memory kernel to negative time, thus allowing for some pre-initial conditions. This will be the content of further research.

This work is understood as a first step to understand the principles behind damping based on fractional time derivatives. In subsequent works, we plan to extend the analysis to a more physical model, namely that of a true membrane over a viscous incompressible fluid governed by the Navier–Stokes equations. The approach based on partial differential equations developed here, will then allow us to study the full vector-valued case \(v(t,x,z)\in \mathbb {R}^d\) including the associated nonlinearities. It will be interesting to see under what conditions the relevant scalings in the nonlinear setting will be the same as in the linear theory in [12, 13]. Moreover, it will be critical to see the occurrence of fractional damping, which relies on the linearity of the Dirichlet-to-Neumann map of the parabolic equation on the vertical half line.