# An energy approach to space-time Galerkin BEM for wave propagation problems

A. Aimi, M. Diligenti, C. Guardasoni, I. Mazzieri, S. Panizzi

► To cite this version:

A. Aimi, M. Diligenti, C. Guardasoni, I. Mazzieri, S. Panizzi. An energy approach to space-time Galerkin BEM for wave propagation problems. International Journal for Numerical Methods in Engineering, 2009, 80 (9), 10.1002/nme.2660. hal-01352570

HAL Id: hal-01352570

https://hal.science/hal-01352570

Submitted on 8 Aug 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An energy approach to space-time Galerkin BEM for wave propagation problems

A. Aimi¹, M. Diligenti¹, C. Guardasoni², I. Mazzieri³ and S. Panizzi¹

¹Department of Mathematics, Università di Parma, V.le G.P. Usberti 53/A, Parma, Italy

²Department of Mathematics, Università degli Studi di Milano, V. Saldini 50, Milan, Italy

³Department of Mathematics, Politecnico di Milano, P. za Leonardo da Vinci 32, Milan, Italy

In this paper we consider Dirichlet or Neumann wave propagation problems reformulated in terms of boundary integral equations with retarded potential. Starting from a natural energy identity, a space-time weak formulation for 1D integral problems is briefly introduced, and continuity and coerciveness properties of the related bilinear form are proved. Then, a theoretical analysis of an extension of the introduced formulation for 2D problems is proposed, pointing out the novelty with respect to existing literature results. At last, various numerical simulations will be presented and discussed, showing unconditional stability of the space-time Galerkin boundary element method applied to the energetic weak problem.

KEY WORDS: wave propagation; energy identity; boundary integral equation; weak formulation; Galerkin boundary element method

## 1. INTRODUCTION

Time-dependent problems that are frequently modeled by hyperbolic partial differential equations can be dealt with the boundary integral equations (BIEs) method. The ideal situation is when we have a homogeneous partial differential equation with constant coefficients, the initial conditions vanish and the data are given only on the boundary of a domain independent of time. In this situation the transformation of the problem to a BIE follows the same well-known method for elliptic boundary value problems. Causality condition and time invariance imply that the integral equations are of Volterra type in time variable and of convolution type in time, respectively. Boundary element methods (BEMs) have been successfully applied to many such problems from
fields as electromagnetic wave propagation, computation of transient acoustic wave, linear elasto-dynamics, fluid dynamics, etc. [1–7]. Frequently claimed advantages over domain approaches are the dimensionality reduction, the easy implicit enforcement for radiation conditions at infinity, the reduction of an unbounded exterior domain to a bounded boundary, the high accuracy achievable and simple pre- and post-processing for input and output data.

In principle, both the frequency-domain and time-domain BEM can be used for hyperbolic boundary value problems. Most earlier contributions concerned direct formulations of BEM in the frequency domain, often using the Laplace or Fourier transforms and addressing wave propagation problems. After this transformation a standard boundary integral method for an elliptic (Helmoltz) problem is applied and then the transformation back to time domain employs special methods for the inversion of Laplace or Fourier transforms.

On the other side, the consideration of the time-domain (transient) problem yields directly the unknown time-dependent quantities. In this case, the representation formula in terms of single layer and double layer potentials uses the fundamental solution of the hyperbolic partial differential equation and jump relations, giving rise to retarded BIEs. Usual numerical discretization procedures include collocation techniques and Laplace–Fourier methods coupled with Galerkin boundary elements in space. The convolution quadrature method for the time discretization has been developed in [8–10]. It provides a straightforward way to obtain a stable time stepping scheme using the Laplace transform of the kernel function.

The application of Galerkin boundary elements in both space and time has been implemented by several authors, but in this direction only the weak formulation due to Bamberger and Ha Duong [11–14] furnishes genuine convergence results. During the last 20 years, the Bamberger and Ha Duong method has been successfully applied to many problems in transient wave propagations (see e.g. [15–17]). The technique they use to find the weak formulation and to prove stability results may be summarized in the following steps: Fourier–Laplace transform in time variable; uniform estimates with respect to complex frequencies of the corresponding Helmholtz problem; application of the Paley–Wiener theorem and Parseval identity. In particular, this final step provides a space–time weak formulation closely related to the energy functional of the wave equation, whose associated quadratic form turns out to be coercive with respect to a suitable weighted Sobolev norm. The only drawback of the method is that passage to complex frequencies leads to stability constants that grow exponentially in time, as stated in [18]. We refer to the surveys [18, 19] for a more complete bibliography on the subject.

In this paper, we consider a Dirichlet or Neumann problem in $\mathbb{R}^n$, $n=1,2$, for a temporally homogeneous (normalized) scalar wave equation outside the obstacle $\Gamma$ in the time interval $[0, T]$, reformulated as a BIE with retarded potential. We avoid the passage to complex frequencies by simply exploiting the well-known energy-flux relation:

$$ \frac{\partial}{\partial t} \left( \frac{1}{2} u_t^2 + \frac{1}{2} |\nabla u|^2 \right) - \nabla_x \cdot (u_t \nabla u) = 0 \quad (1) $$

satisfied by any solution $u$ of the (real valued) wave equation. Under, for instance, Dirichlet boundary conditions, (1) yields the identity

$$ \mathcal{E}(T, u) := \frac{1}{2} \int_{\mathbb{R}^n \setminus \Gamma} \left( \left( \frac{\partial u(\mathbf{x}, T)}{\partial t} \right)^2 + |\nabla u(\mathbf{x}, T)|^2 \right) d\mathbf{x} = \int_{\Gamma \times [0, T]} \frac{\partial u}{\partial t}(\mathbf{x}, t) \left[ \frac{\partial u}{\partial \mathbf{n}} \right](\mathbf{x}, t) dt d\gamma_{\mathbf{x}} \quad (2) $$
where $[\partial u/\partial n]$ represents both the jump of the normal derivative of $u$ along $\Gamma$ and the unknown density function $\varphi$ of the single layer BIE associated with the differential problem. From (2), we obtain a natural quadratic form in the unknown density and we derive a suitable space-time weak formulation of the integral problem.

We remark that the idea of introducing a space-time weak formulation for the transient wave problem based on the energy identity (2) is not new. In fact, in [17] it was already exploited to get a satisfactory stability result for the acoustic wave equation with the aid of an absorbing boundary condition. Unfortunately, the case of the Dirichlet or Neumann problem, considered here, does not lead to any natural coerciveness property, as we shall see in Section 3.

Starting from the consideration that for suitable geometries of $\Gamma$, and thanks to the finite speed of propagation property, the square root of $\mathcal{E}(T, u)$ defines a norm of the density $\varphi$, special attention is devoted to the coerciveness properties of the energy functional: in Section 2, devoted to the simple case of 1D problems, we have a precise estimate in $L^2(\Sigma_T)$ and consequently, unconditionally stable schemes with well-behaved stability constants even for large times; in Section 3 we analyze, via Fourier transform, the 2D extension in the case of a flat obstacle. We prove that the energetic quadratic form cannot be coercive with respect to any Sobolev norm (Theorem 3.1) and we outline two possible strategies to overcome the stability problem: (1) with the aid of a constraint on the oscillations in the space variable of the unknown density, we provide (Theorem 3.2) a coerciveness estimate in $L^2(\Sigma_T)$; (2) we show how the introduction of the Hilbert transform with respect to time would provide, at least theoretically, an unconditionally stable weak formulation (Theorem 3.3). Then in Section 4 Galerkin BEM discretization is introduced and in Section 5 some interesting numerical integration issues are addressed.

At last numerical simulations for 1D and 2D wave propagation problems will be presented and discussed. We will compare, when possible, results obtained with the energetic Galerkin BEM with analogous literature results. Instabilities phenomena, which arise starting from classical $L^2$ weak formulation of the BIEs both in the case of 1D and 2D problems, as shown, respectively, in [15, 20], are never present in the energetic procedure, which appears to be unconditionally stable.

We finally note that the present analysis has its natural counterpart in the context of 3D problems: theoretical results of Section 3 can be extended without modifications to any spatial dimension $n \ge 3$; 3D simulations, focused on applied seismology, are scheduled by the authors as next work.

## 2. 1D WAVE EQUATION

### 2.1. Dirichlet problem and energetic weak formulation of the related BIE

We consider the retarded potential representation of solutions of the Dirichlet problem for a 1D wave equation, from which we will deduce a suitable boundary integral reformulation of the differential problem.

Let $\Omega = (0, L) \subset \mathbb{R}$ and let $u(x, t)$ be the solution to the wave problem

$$u_{tt} - u_{xx} = 0, \quad x \in \mathbb{R} \setminus \{0, L\}, \quad t \in (0, T) \qquad (3)$$

$$u(x, 0) = u_t(x, 0) = 0, \quad x \in \mathbb{R} \setminus \{0, L\} \qquad (4)$$

$$u(x, t) = g(x, t), \quad (x, t) \in \Sigma_T := \{0, L\} \times [0, T] \qquad (5)$$
where $g(x, t)$ is a given function. Note that $u$ is considered as the solution on the whole $\mathbb{R}$, not only in $\Omega$. Whenever necessary we shall distinguish the internal solution $u^-$, i.e. for $x \in \Omega$, from the external one $u^+$, i.e. for $x \in \mathbb{R} \setminus \Omega$.

In order to rewrite problem (3)–(5) as a BIE, we need to recall the expression of the forward fundamental solution $G(x, t)$ of the wave operator:

$$G(x, t) = \frac{1}{2} H[t - |x|] = \frac{1}{2} H[t](H[x + t] - H[x - t]) \quad (6)$$

where $H[t]$ is the Heaviside step function. Using the fundamental solution (6) we obtain the single layer representation formula for the solutions of the Equation (3), in terms of the unknown density vector-valued function $\varphi(t) = (\varphi(0, t), \varphi(L, t))^{\top} =: (\varphi_0(t), \varphi_L(t))^{\top}$. In fact, for $t \in \mathbb{R}$, $x \in \mathbb{R} \setminus \{0, L\}$, it holds:

$$u(x, t) = (V\varphi)(x, t) = \sum_{\xi=0, L} \int_{-\infty}^{+\infty} G(x - \xi, t - \tau) \varphi_{\xi}(\tau) \, d\tau = \frac{1}{2} \sum_{\xi=0, L} \int_{-\infty}^{t - |x - \xi|} \varphi_{\xi}(\tau) \, d\tau \quad (7)$$

It is well known that

$$\varphi = \left[ \frac{\partial u}{\partial v} \right] := \frac{\partial u^-}{\partial v} - \frac{\partial u^+}{\partial v} \quad \text{where} \quad \frac{\partial}{\partial v} = \begin{cases} -\frac{\partial}{\partial x}, & x=0 \\ \frac{\partial}{\partial x}, & x=L \end{cases} \quad (8)$$

Since problem (3)–(5) is formulated on the time interval $[0, T]$, in order to keep our notations as simple as possible, hereafter we shall consider functions $\varphi$ defined on the whole real line but having support only in the fixed time interval $[0, T]$. From formula (7), taking the limits as $x \to 0$ and $x \to L$, we obtain the following system of BIEs at the endpoints of the interval $(0, L)$:

$$(V\varphi)(0, t) = \frac{1}{2} \left[ \int_{0}^{t} \varphi_{0}(\tau) \, d\tau + \int_{0}^{t-L} \varphi_{L}(\tau) \, d\tau \right] = g_{0}(t), \quad t \in [0, T] \quad (9)$$

$$(V\varphi)(L, t) = \frac{1}{2} \left[ \int_{0}^{t-L} \varphi_{0}(\tau) \, d\tau + \int_{0}^{t} \varphi_{L}(\tau) \, d\tau \right] = g_{L}(t),$$

The system (9) can be written with the compact notation

$$V\varphi = g \quad (10)$$

where $g(t) = (g_0(t), g_L(t))^{\top}$ represents the given boundary function in $x=0$ and $x=L$.

In order to formulate the operator equation (10) in a suitable functional framework, we assume $V$ as defined in $L^2(\Sigma_T) := L^2(0, T) \times L^2(0, T)$. With this choice, from (9) it is easy to verify that the range of $V$ lies in $H^1_{\{0\}}(\Sigma_T)$, the space of $H^1(0, T) \times H^1(0, T)$ functions vanishing for $t=0$. Hence hereafter, for 1D problems, we shall consider

$$V: L^2(\Sigma_T) \to H^1_{\{0\}}(\Sigma_T)$$

A classical way to introduce a weak formulation for (10) is to project the BIE using $L^2(\Sigma_T)$ scalar product. Unfortunately, the resulting bilinear form is not coercive; further, in this way, it is implicit that Equation (10) must be understood considering the compact operator $V: L^2(\Sigma_T) \to L^2(\Sigma_T)$, which obviously cannot have continuous inverse. As a consequence, it is not surprising that classical
$L^2$ weak problem gives rise to instability phenomena in the discretization phase, as it has been shown in [20]. An alternative approach is suggested by the well-known conservation law satisfied by the (real-valued) solutions to the d'Alembert equation:

$$0=u_t(u_{tt}-u_{xx})=\frac{\partial}{\partial t}(\frac{1}{2}u_t^2+\frac{1}{2}u_x^2)-\frac{\partial}{\partial x}(u_t u_x)$$

Integrating with respect to space-time in $\mathbb{R} \times (0, T)$ and taking into account that $u$ and $u_t$ vanish for $t=0$, we get the energy identity

$$\frac{1}{2} \int_{-\infty}^{+\infty} (u_t^2 + u_x^2) \, dx \Big|_{t=T} = \int_{0}^{T} u_t \cdot \left[ \frac{\partial u}{\partial v} \right] \, dt = \int_{0}^{T} (V\varphi)_t \cdot \varphi \, dt =: \mathcal{E}(T, \varphi) \quad (11)$$

where the dot denotes the scalar product in $\mathbb{R}^2$.

The quadratic form appearing in the last term of (11) leads to a natural space-time weak formulation of the corresponding BIE (10) with robust theoretical properties. In fact, the main advantage of this approach is that the quadratic form given by the energy, i.e.

$$\mathcal{E}(T, \varphi) = \langle (V\varphi)_t, \varphi \rangle_{L^2(\Sigma_T)}$$

is, in the 1D case, both continuous and coercive in the appropriate spaces, i.e. exactly the functional spaces where the Dirichlet problem is well-posed [21].

In order to derive continuity and coerciveness properties of the total energy $\mathcal{E}(T, \varphi)$, we concentrate our analysis on the operator $A: L^2(\Sigma_T) \to L^2(\Sigma_T)$, defined as

$$A\varphi(t) := (V\varphi)_t(t) = \begin{bmatrix} (V\varphi)_t(0,t) \\ (V\varphi)_t(L,t) \end{bmatrix} = \frac{1}{2} \begin{bmatrix} \varphi_0(t) + H[t-L]\varphi_L(t-L) \\ H[t-L]\varphi_0(t-L) + \varphi_L(t) \end{bmatrix}, \quad t \in [0,T] \quad (12)$$

By an application of the Cauchy-Schwarz inequality, we have immediately that $A$ is a continuous operator. We state this property in the following.

### Proposition 2.1

For any given time $T$, the operator $A: L^2(\Sigma_T) \to L^2(\Sigma_T)$, defined in (12) is bounded, with norm $\|A\| \le 1$.

More interesting are the positivity properties of the quadratic form associated with the operator $A$. Instead of the non-symmetric operator $A$, we shall consider its symmetric part

$$A_s = \frac{A + A^*}{2}$$

where the asterisk denotes the adjoint of an operator. In fact, having introduced the bilinear form $a_\epsilon(\varphi, \psi) : L^2(\Sigma_T) \times L^2(\Sigma_T) \to \mathbb{R}$ defined by

$$a_{\epsilon}(\varphi, \psi) := \langle A\varphi, \psi \rangle_{L^2(\Sigma_T)} = \int_{0}^{T} (V\varphi)_{t}(0, t)\psi_{0}(t) \, dt + \int_{0}^{T} (V\varphi)_{t}(L, t)\psi_{L}(t) \, dt \quad (13)$$

for every $\varphi \in L^2(\Sigma_T)$, it holds

$$a_{\epsilon}(\varphi, \varphi) = \langle A_s \varphi, \varphi \rangle_{L^2(\Sigma_T)} \quad (14)$$
The main result of this section is stated in the following

### Theorem 2.1

For every $T>0$, there exists a positive constant $c(T)$ such that

$$a_{\varepsilon}(\varphi, \varphi) \geq c(T) \|\varphi\|_{L^2(\Sigma_T)}^2, \quad \varphi \in L^2(\Sigma_T) \qquad (15)$$

Moreover, let $N$ be the least positive integer such that $T \le NL$; then one has the explicit bound:

$$c(T) \ge \sin^2 \left( \frac{\pi}{2(N+1)} \right) \qquad (16)$$

The proof of this result is given in Appendix.

At this point we can write down the energetic weak problem related to the BIE (10), which admits a unique, stable solution:

given $g \in H_{\{0\}}^1(\Sigma_T)$, find $\varphi \in L^2(\Sigma_T)$ such that

$$a_{\varepsilon}(\varphi, \psi) = \langle g_t, \psi \rangle_{L^2(\Sigma_T)} \quad \forall \psi \in L^2(\Sigma_T) \qquad (17)$$

## 2.2. Neumann problem

Similar considerations as those developed in the previous subsection can be done for the wave problem with Neumann boundary condition:

$$u_{tt} - u_{xx} = 0, \quad x \in \mathbb{R} \setminus \{0, L\}, \quad t \in (0, T) \qquad (18)$$

$$u(x, 0) = u_t(x, 0) = 0, \quad x \in \mathbb{R} \setminus \{0, L\} \qquad (19)$$

$$\frac{\partial u}{\partial v}(x, t) = g(x, t), \quad (x, t) \in \Sigma_T := \{0, L\} \times [0, T] \qquad (20)$$

where $g(x, t)$ is a given function. In this case we have the double layer representation of the solution $u(x, t)$ through the unknown retarded potential $\phi$:

$$u(x, t) = (K \phi)(x, t) = \sum_{\xi=0, L} \int_{-\infty}^{+\infty} \frac{\partial}{\partial v} G(x - \xi, t - \tau) \phi(\xi, \tau) \, d\tau \qquad (21)$$

where $G(\cdot, \cdot)$ is given in (6). After a straightforward calculation we obtain the more explicit formula

$$u(x, t) = \frac{1}{2} \frac{x}{|x|} \phi_0(t - |x|) - \frac{1}{2} \frac{x - L}{|x - L|} \phi_L(t - |x|)$$

from which

$$[u](0, t) := u^{-}(0, t) - u^{+}(0, t) = \phi_{0}(t), \quad [u](L, t) := u^{-}(L, t) - u^{+}(L, t) = \phi_{L}(t)$$

By deriving formula (21), taking the limits as $x \to 0$ and $x \to L$ and using Neumann datum $g(x, t)$, we obtain the boundary equations for the unknown potential $\phi$ at the endpoints of the interval $(0, L)$:

$$D\phi := \frac{\partial}{\partial v} K\phi = g \qquad (22)$$
that is

$$ (D\phi)(0,t) = \frac{1}{2}[\phi_{0,t}(t) - \phi_{L,t}(t-L)] = g_0(t) \\ (D\phi)(L,t) = -\frac{1}{2}[\phi_{0,t}(t-L) - \phi_{L,t}(t)] = g_L(t) $$

In order to derive a weak formulation of Equation (22), our starting point is again the energy identity, which for the Neumann boundary condition assumes the form

$$ \int_{0}^{T} \frac{\partial u}{\partial v} \cdot [u]_{t} \, dt = \int_{0}^{T} D\phi \cdot \phi_{t} \, dt =: \mathcal{E}(T, \phi) $$

Thus, having defined the bilinear form

$$ \tilde{a}_{\varepsilon}: H_{\{0\}}^{1}(\Sigma_{T}) \times H_{\{0\}}^{1}(\Sigma_{T}) \to \mathbb{R}, \quad \tilde{a}_{\varepsilon}(\phi, \eta) := \langle D\phi, \eta_{t} \rangle_{L^{2}(\Sigma_{T})} $$

the coerciveness of $\tilde{a}_{\varepsilon}(\cdot, \cdot)$ follows at once from the observation that

$$ \tilde{a}_{\varepsilon}(\phi, \phi) = \int_{0}^{T} A \tilde{\phi}_{t} \cdot \tilde{\phi}_{t} \, dt $$

where

$$ \tilde{\phi}(t) = (\phi_0(t), -\phi_L(t))^\top $$

and $A$ is the operator defined in (12). Then from Theorem 2.1 we have the following

**Theorem 2.2**

For every $T>0$, there exists a positive constant $c(T)$ such that

$$ \tilde{a}_{\varepsilon}(\phi, \phi) \geq c(T) \|\phi_t\|_{L^2(\Sigma_T)}^2, \quad \phi \in H_{\{0\}}^{1}(\Sigma_T) $$

The constant $c(T)$ is bounded from below as in (16).

We can conclude that for 1D problems, Dirichlet or Neumann boundary conditions are similar for what concerns the analysis of the energetic weak formulations related to the corresponding BIEs.

## 3. 2D WAVE EQUATION

### 3.1. Dirichlet problem and its energetic weak formulation

Here, we will consider a Dirichlet problem for the wave equation exterior to an open arc $\Gamma \subset \mathbb{R}^2$:

$$ u_{tt} - \Delta u = 0, \quad \mathbf{x} \in \mathbb{R}^2 \setminus \Gamma, \quad t \in (0, T) \qquad (23) $$

$$ u(\mathbf{x}, 0) = u_t(\mathbf{x}, 0) = 0, \quad \mathbf{x} \in \mathbb{R}^2 \setminus \Gamma \qquad (24) $$

$$ u(\mathbf{x}, t) = g(\mathbf{x}, t), \quad (\mathbf{x}, t) \in \Sigma_T := \Gamma \times [0, T] \qquad (25) $$
In this case the boundary datum $g(\mathbf{x}, t)$ represents the value of the excitation field over $\Gamma$. As we have seen for 1D problems, we can consider the single-layer representation of the solution of (23)–(25):

$$u(\mathbf{x}, t) = \int_{\Gamma} \int_{0}^{t} G(r, t-\tau) \varphi(\xi, \tau) \, d\tau \, d\gamma_{\xi}, \quad \mathbf{x} \in \mathbb{R}^{2} \setminus \Gamma, \ t \in (0, T) \qquad (26)$$

where $r = \|\mathbf{r}\|_2 = \|\mathbf{x} - \xi\|_2$, $\varphi = [\partial u / \partial \mathbf{n}]$ is the jump of the normal derivative of $u$ along $\Gamma$ and

$$G(r, t-\tau) = \frac{1}{2\pi} \frac{H[t-\tau-r]}{[(t-\tau)^2 - r^2]^{1/2}} \qquad (27)$$

is the forward fundamental solution of the 2D wave operator. With a limiting process for $\mathbf{x}$ tending to $\Gamma$ and using the assigned Dirichlet boundary condition we obtain the space-time BIE

$$\int_{\Gamma} \int_{0}^{t} G(r, t-\tau) \varphi(\xi, \tau) \, d\tau \, d\gamma_{\xi} = g(\mathbf{x}, t), \quad \mathbf{x} \in \Gamma, \ t \in (0, T) \qquad (28)$$

which can be written with the compact notation

$$V\varphi = g \qquad (29)$$

The energetic weak formulation of problem (29) is defined as in (17), with

$$a_{\varepsilon}(\varphi, \psi) = \langle (V\varphi)_t, \psi \rangle_{L^2(\Sigma_T)} = \int_{\Gamma} \int_{0}^{T} (V\varphi)_t(\mathbf{x}, t)\psi(\mathbf{x}, t) \, dt \, d\gamma_{\mathbf{x}} \qquad (30)$$

where $\psi$ is a suitable test function, belonging to the same functional space of $\varphi$. With an integration by parts and (24), we can also write

$$a_{\varepsilon}(\varphi, \psi) = \int_{\Gamma} (V\varphi)(\mathbf{x}, T)\psi(\mathbf{x}, T) \, d\gamma_{\mathbf{x}} - \int_{\Gamma} \int_{0}^{T} (V\varphi)(\mathbf{x}, t)\psi_{t}(\mathbf{x}, t) \, dt \, d\gamma_{\mathbf{x}} \qquad (31)$$

In order to apply the Fourier transform to the analysis of the bilinear form (30), we assume that the obstacle $\Gamma$ is flat, that is $\Gamma = \{(x, 0) : x \in [0, L]\} \subset \mathbb{R}$. Then, for a given regular function $\varphi$ having support in $\Sigma_T$, we reinterpret the function $u$ in formula (26), as the solution of the jump problem:

$$u_{tt} - u_{xx} - u_{yy} = 0, \quad (x, y) \in \mathbb{R}^2 \setminus \Gamma, \ t > 0 \qquad (32)$$

$$\left[ \frac{\partial u}{\partial y} \right] (x, 0, t) = \varphi(x, t), \quad 0 \le x \le L, \ t > 0 \qquad (33)$$

$$[u](x, 0, t) = 0, \quad 0 \le x \le L, \ t > 0 \qquad (34)$$

$$u(x, y, 0) = u_t(x, y, 0) = 0, \quad (x, y) \in \mathbb{R}^2 \qquad (35)$$

For any $f \in L^2(\Sigma_T)$, we define the Fourier transform with respect to a single variable, either $x$ or $t$, as follows:

$$\hat{f}(\xi, t) := \mathcal{F}_x(f)(\xi, t) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{+\infty} e^{-ix\xi} f(x, t) \, dx$$

$$\hat{f}(x, \omega) := \mathcal{F}_t(f)(x, \omega) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{+\infty} e^{-it\omega} f(x, t) \, dt$$
while

$$ \tilde{f}(\xi, \omega) := \mathcal{F}(f)(\xi, \omega) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} e^{-i(t\omega + \xi x)} f(x, t) \, dt \, dx \quad (36) $$

will denote the Fourier transform with respect to both variables. As it is well-known by the Paley–Wiener theorem, $\tilde{f}(\omega, \xi)$ is an entire analytical function of exponential type.

We assume for the moment that $\varphi \in C_0^\infty(\Sigma_T)$ and we perform the transformation of the operator $\varphi \mapsto (V\varphi)_t$ into a multiplication operator in the frequency variables. Having set

$$ \hat{u}(\xi, y, t) := \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{+\infty} e^{-ix\xi} u(x, y, t) \, dx, \quad y \neq 0, \ t > 0 $$

we obtain the following homogeneous Cauchy problem with a jump condition for $\hat{u}_y$ in $y=0$, for the Klein–Gordon equation with mass $\xi^2$:

$$ \hat{u}_{tt} - \hat{u}_{yy} - \xi^2 \hat{u} = 0, \quad y \neq 0, \ t > 0 \quad (37) $$

$$ \left[ \frac{\partial \hat{u}}{\partial y} \right](\xi, 0, t) = \hat{\phi}(\xi, t), \quad t > 0 \quad (38) $$

$$ [\hat{u}](\xi, 0, t) = 0, \quad t > 0 \quad (39) $$

$$ \hat{u}(\xi, y, 0) = \hat{u}_t(\xi, y, 0) = 0, \quad y \neq 0 \quad (40) $$

In order to represent $\hat{u}$ as a single layer integral, we need the expression of the fundamental solution of the Klein–Gordon equation, which is given by (see [22])

$$ K(\xi, y, t) = \frac{1}{2} H[t - |y|] J_0(|\xi|\sqrt{t^2 - y^2}) \quad (41) $$

where $J_0$ represents the Bessel function of order 0. Then we have

$$ \hat{u}(\xi, y, t) = \frac{1}{2} \int_{0}^{t-|y|} J_0(|\xi|\sqrt{(t-\tau)^2 - y^2}) \hat{\phi}(\xi, \tau) \, d\tau, \quad \xi \in \mathbb{R}, \ y \neq 0, \ t > 0 \quad (42) $$

and by taking the limit $y \to 0$,

$$ \hat{u}(\xi, 0, t) = \frac{1}{2} \int_{0}^{t} J_{0}(|\xi|(t-\tau)) \hat{\phi}(\xi, \tau) \, d\tau, \quad \xi \in \mathbb{R}, \ t > 0 \quad (43) $$

Since

$$ \frac{\partial \hat{u}}{\partial t} = \frac{\partial}{\partial t} \mathcal{F}_x[u] = \mathcal{F}_x \left[ \frac{\partial u}{\partial t} \right] = \mathcal{F}_x[(V\varphi)_t] $$

$J_0(0) = 1$ and $J'_0(t) = -J_1(t)$, deriving in (43) with respect to $t$, we get

$$ \mathcal{F}_x[(V\varphi)_t] = \frac{\partial \hat{u}}{\partial t}(\xi, 0, t) = \frac{1}{2}\hat{\phi}(\xi, t) - \frac{|\xi|}{2} \int_0^t J_1(|\xi|(t-\tau))\hat{\phi}(\xi, \tau) \, d\tau \quad (44) $$
We need the following result (see [23]) on the Fourier transform of $H[t]J_1(t)$:

$$G(\omega) := \sqrt{2\pi} \mathcal{F}[H[t]J_1(t)](\omega) = 1 + \frac{i\omega\chi_{\{|\omega|<1\}}}{\sqrt{1-\omega^2}} - \frac{|\omega|\chi_{\{|\omega|>1\}}}{\sqrt{\omega^2-1}}$$

where $\chi_\mathcal{L}$ is the characteristic function of the set $\mathcal{L}$. Therefore

$$\sqrt{2\pi} \mathcal{F}_t[H[t]J_1(|\xi|t)](\omega) = \frac{1}{|\xi|} G\left(\frac{\omega}{|\xi|}\right) = \frac{1}{|\xi|} \left[ 1 + \frac{i\omega\chi_{\{|\omega|<|\xi|\}}}{\sqrt{|\xi|^2 - \omega^2}} - \frac{|\omega|\chi_{\{|\omega|>|\xi|\}}}{\sqrt{\omega^2 - |\xi|^2}} \right]$$

Finally, from (44), we obtain the desired expression of the Fourier transform of $(V\varphi)_t$, when $\varphi \in C_0^\infty(\Sigma_T)$:

$$\begin{aligned} \mathcal{F}[(V\varphi)_t](\xi, \omega) &= \frac{1}{2}\tilde{\varphi}(\xi, \omega) - \frac{1}{2}G\left(\frac{\omega}{|\xi|}\right)\tilde{\varphi}(\xi, \omega) \\ &= \frac{1}{2}\frac{|\omega|\chi_{\{|\omega|>|\xi|\}}}{\sqrt{\omega^2 - |\xi|^2}}\tilde{\varphi}(\xi, \omega) - \frac{1}{2}\frac{i\omega\chi_{\{|\omega|<|\xi|\}}}{\sqrt{|\xi|^2 - \omega^2}}\tilde{\varphi}(\xi, \omega) \end{aligned} \quad (45)$$

Owing to the Parseval identity, for $\psi \in \mathcal{S}(\mathbb{R}^2)$ and $\varphi \in C_0^\infty(\Sigma_T)$, we get the following representation of the bilinear form

$$\begin{aligned} a_{\varepsilon}(\varphi, \psi) &= \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} (V\varphi)_t(x, t) \psi(x, t) \, dt \, dx \\ &= \frac{1}{2} \int_{-\infty}^{+\infty} \mathrm{d}\xi \int_{|\omega|>|\xi|} \frac{|\omega|}{\sqrt{\omega^2 - \xi^2}} \tilde{\varphi}(\xi, \omega) \tilde{\psi}(\xi, \omega) \, \mathrm{d}\omega \\ &\quad - \frac{i}{2} \int_{-\infty}^{+\infty} \mathrm{d}\xi \int_{|\omega|<|\xi|} \frac{\omega}{\sqrt{\xi^2 - \omega^2}} \tilde{\varphi}(\xi, \omega) \tilde{\psi}(\xi, \omega) \, \mathrm{d}\omega \end{aligned} \quad (46)$$

Of course, when the support of the function $\psi$ lies in $\Sigma_T$, the double integral on the left-hand side of (46) is restricted to $\Sigma_T$. In particular, for $\psi = \varphi$, from the property $\mathcal{E}(T, \varphi) = \Re(a_\varepsilon(\varphi, \varphi))$, or remembering that for real-valued functions $\tilde{\varphi}(\xi, -\omega) = \tilde{\varphi}(\xi, \omega)$, we have

### Proposition 3.1

For any $\varphi \in C_0^\infty(\Sigma_T)$, the energy at time $T$ of the solution $u$ to the problem (32)–(35), is given by the following formula

$$\mathcal{E}(T, \varphi) = a_{\varepsilon}(\varphi, \varphi) = \int_{\Sigma_T} (V\varphi)_t(x, t) \varphi(x, t) \, dt \, dx = \frac{1}{2} \int_{-\infty}^{+\infty} \mathrm{d}\xi \int_{|\omega|>|\xi|} \frac{|\omega||\tilde{\varphi}(\xi, \omega)|^2}{\sqrt{\omega^2 - \xi^2}} \mathrm{d}\omega \quad (47)$$

We define $\mathcal{V}(a_\varepsilon) \subset L^2(\Sigma_T)$ the set of functions $\varphi, \psi$ such that $|a_\varepsilon(\varphi, \psi)|<+\infty$. In order to characterize $\mathcal{V}(a_\varepsilon)$, we assume for the moment that $\varphi, \psi$ are complex-valued functions. In this case, also the imaginary part of $a_\varepsilon(\varphi, \varphi)$ does not vanish, thus

$$a_{\varepsilon}(\varphi, \varphi) = \frac{1}{2} \int_{-\infty}^{+\infty} \mathrm{d}\xi \int_{|\omega|>|\xi|} \frac{|\omega||\tilde{\varphi}(\xi, \omega)|^2}{\sqrt{\xi^2 - \omega^2}} \mathrm{d}\omega - \frac{i}{2} \int_{-\infty}^{+\infty} \mathrm{d}\xi \int_{|\omega|<|\xi|} \frac{\omega|\tilde{\varphi}(\xi, \omega)|^2}{\sqrt{\omega^2 - \xi^2}} \mathrm{d}\omega \quad (48)$$
It follows that

$$ \mathcal{V}(a_{\varepsilon}) = \left\{ \varphi \in L^2(\Sigma_T) : \int_{\mathbb{R}^2} \frac{|\omega| |\tilde{\phi}(\xi, \omega)|^2}{\sqrt{|\xi^2 - \omega^2|}} d\omega d\xi < +\infty \right\} $$

Now we can state the following

**Proposition 3.2**

The space $L^2((0, T); H^{1/4}(\Gamma))$, that is the closure of $C_0^\infty(\Sigma_T)$ with respect to the norm

$$ \|\varphi\|_{L^2((0,T); H^{1/4}(\Gamma))}^2 = \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} (1 + |\xi|)^{1/2} |\tilde{\phi}(\xi, \omega)|^2 d\xi d\omega $$

is contained in $\mathcal{V}(a_\varepsilon)$. Moreover, there exists a positive constant $C$ such that, for every $\varphi$, $\psi \in L^2((0, T); H^{1/4}(\Gamma))$,

$$ a_{\varepsilon}(\varphi, \psi) \le C(1+T) \|\varphi\|_{L^2((0,T); H^{1/4}(\Gamma))} \|\psi\|_{L^2((0,T); H^{1/4}(\Gamma))} $$

*Proof*

Obviously we have

$$ \int_{\mathbb{R}^2} \frac{|\omega| |\tilde{\phi}(\xi, \omega)|^2}{\sqrt{|\xi^2 - \omega^2|}} d\omega d\xi = \int_{\{|\omega| > |\xi|\}} \frac{|\omega| |\tilde{\phi}(\xi, \omega)|^2}{\sqrt{\omega^2 - \xi^2}} d\omega d\xi + \int_{\{|\omega| < |\xi|\}} \frac{|\omega| |\tilde{\phi}(\xi, \omega)|^2}{\sqrt{\xi^2 - \omega^2}} d\omega d\xi \quad (49) $$

Let us estimate the first addendum. We have

$$ \int_{-\infty}^{+\infty} d\xi \int_{|\omega| > |\xi|} \frac{|\omega| |\tilde{\phi}(\xi, \omega)|^2}{\sqrt{\omega^2 - \xi^2}} d\omega = \int_{-\infty}^{+\infty} d\xi \int_{|\xi| < |\omega| < |\xi| + 1} \frac{|\omega| |\tilde{\phi}(\xi, \omega)|^2}{\sqrt{\omega^2 - \xi^2}} d\omega \\ + \int_{-\infty}^{+\infty} d\xi \int_{|\omega| > |\xi| + 1} \frac{|\omega| |\tilde{\phi}(\xi, \omega)|^2}{\sqrt{\omega^2 - \xi^2}} d\omega =: I_1^a + I_2^a $$

From the definition of $\tilde{\phi}$ and thanks to the Cauchy–Schwarz inequality, we get

$$ |\tilde{\phi}(\xi, \omega)|^2 \le \frac{T}{2\pi} \int_0^T |\hat{\phi}(\xi, t)|^2 dt $$

Therefore, by a simple integration and thanks to the Parseval identity, we have

$$ \begin{aligned} I_1^a &\le \frac{T}{2\pi} \int_{-\infty}^{+\infty} \left( \int_0^T |\hat{\phi}(\xi, t)|^2 dt \right) d\xi \int_{|\xi| < |\omega| < |\xi| + 1} \frac{|\omega|}{\sqrt{\omega^2 - \xi^2}} d\omega \\ &= \frac{T}{\pi} \int_{-\infty}^{+\infty} \left( \int_0^T |\hat{\phi}(\xi, t)|^2 dt \right) \sqrt{2|\xi| + 1} d\xi \le \frac{\sqrt{2T}}{\pi} \|\varphi\|_{L^2((0,T); H^{1/4}(\Gamma))}^2 \end{aligned} $$

On the other hand, from the inequality

$$ \frac{|\omega|}{\sqrt{\omega^2 - \xi^2}} \le \sqrt{|\xi| + 1} \quad (|\omega| > |\xi| + 1) $$
we have the estimate

$$I_2^a \le \int_{-\infty}^{+\infty} \int_{|\omega| > |\xi| + 1} \sqrt{|\xi| + 1} |\tilde{\phi}(\xi, \omega)|^2 d\omega d\xi \le \|\varphi\|_{L^2((0,T); H^{1/4}(\Gamma))}^2$$

The estimate for the second term in (49) is quite similar. Also in this case we may split the integral as follows:

$$\int_{\{|\xi| > |\omega|\}} \frac{|\omega| |\tilde{\phi}(\xi, \omega)|^2}{\sqrt{\xi^2 - \omega^2}} d\xi d\omega = \int_{|\xi| > 1} d\xi \int_{|\omega| < |\xi| - 1} \frac{|\omega| |\tilde{\phi}(\xi, \omega)|^2}{\sqrt{\xi^2 - \omega^2}} d\xi d\omega \\ + \int_{-\infty}^{+\infty} \int_{|\xi| - 1 < |\omega| < |\xi|} \frac{|\omega| |\tilde{\phi}(\xi, \omega)|^2}{\sqrt{\xi^2 - \omega^2}} d\xi d\omega =: I_1^b + I_2^b$$

By an argument similar to the previous one, we get

$$I_1^b \le \int_{|\xi| > 1} d\xi \int_{|\omega| < |\xi| - 1} \frac{(|\xi| - 1) |\tilde{\phi}(\xi, \omega)|^2}{\sqrt{2|\xi| - 1}} d\omega \le \|\varphi\|_{L^2((0,T); H^{1/4}(\Gamma))}^2$$

On the other hand

$$I_2^b \le \frac{T}{2\pi} \int_{-\infty}^{+\infty} \left( \int_0^T |\hat{\phi}(\xi, t)|^2 dt \right) d\xi \int_{|\xi| - 1 < |\omega| < |\xi|} \frac{|\omega|}{\sqrt{\xi^2 - \omega^2}} d\omega \\ = \frac{T}{\pi} \int_{-\infty}^{+\infty} \mu(|\xi|) \left( \int_0^T |\hat{\phi}(\xi, t)|^2 dt \right) d\xi$$

where $\mu(|\xi|) = |\xi|$ for $|\xi| < 1$ and $\mu(|\xi|) = \sqrt{2|\xi| + 1}$ for $|\xi| > 1$. Thus

$$I_1^b \le \frac{\sqrt{2T}}{\pi} \|\varphi\|_{L^2((0,T); H^{1/4}(\Gamma))}^2 \quad \square$$

<figure><img src="image_2.png" /></figure>

Now we come to the question of the coerciveness of the energy functional. Two simple considerations can be drawn by the fact that the domain of integration in formula (47) is the cone $C := \{(\xi, \omega) : |\omega| > |\xi|\}$. The first one is that the energy is a strictly positive functional. This follows immediately from the fact that $\tilde{\phi}$ is an entire analytic function; therefore, it cannot vanish in $C$ unless it vanishes on the whole $\mathbb{R}^2$. The second and more relevant consideration is that most of the information regarding oscillations in the space variable is not taken into account by the energy. This crucial remark is made more precise by the following theorem, where we show that the quadratic form $a_\epsilon(\varphi, \varphi)$ cannot be coercive with respect to any Sobolev norm.

### Theorem 3.1

There exists a sequence $(\varphi_n)_{n \in \mathbb{N}}$ of nonvanishing functions in $C_0^\infty(\Sigma_T)$, such that, for any $s \in \mathbb{R}$ we have

$$\lim_{n \to \infty} \frac{a_\epsilon(\varphi_n, \varphi_n)}{\|\varphi_n\|_s^2} = 0 \qquad (50)$$

where $\|\cdot\|_s^2$ stands for the norm in $H^s(\mathbb{R}^2)$.
Proof

It suffices to show (50) for negative Sobolev exponents. The definition of the sequence $(\varphi_n)_{n \in \mathbb{N}}$ is very simple: let $\varphi \in C_0^\infty(\Sigma_T)$ be a nonvanishing function. We set

$$ \varphi_n(x, t) = \cos(nx)\varphi(x, t) \quad (51) $$

We shall show that

(i) $\forall k \in \mathbb{N}^+$ $a_\varphi(\varphi_n, \varphi_n) = O(1/n^{2k})$ for $n \to +\infty$,

(ii) $\forall s \ge 0 \exists c > 0$ such that $\| \varphi_n \|_{-s}^2 \ge c/n^{2s}$.

It is clear that (i) and (ii) yield (50). From the definition of $\varphi_n$, we have

$$ \tilde{\varphi}_n(\xi, \omega) = \frac{1}{2} \tilde{\varphi}(\xi - n, \omega) + \frac{1}{2} \tilde{\varphi}(\xi + n, \omega) $$

then

$$ |\tilde{\varphi}_n(\xi, \omega)|^2 \le \frac{1}{2} |\tilde{\varphi}(\xi - n, \omega)|^2 + \frac{1}{2} |\tilde{\varphi}(\xi + n, \omega)|^2 $$

From (47), it follows that

$$ \begin{aligned} a_\varphi(\varphi_n, \varphi_n) ={}& \frac{1}{2} \int_{-\infty}^{+\infty} d\xi \int_{|\omega|>|\xi|} \frac{|\omega||\tilde{\varphi}_n(\xi, \omega)|^2}{\sqrt{\omega^2 - \xi^2}} d\omega \le \frac{1}{4} \int_{-\infty}^{+\infty} d\xi \int_{|\omega|>|\xi|} \frac{|\omega||\tilde{\varphi}(\xi - n, \omega)|^2}{\sqrt{\omega^2 - \xi^2}} d\omega \\ & + \frac{1}{4} \int_{-\infty}^{+\infty} d\xi \int_{|\omega|>|\xi|} \frac{|\omega||\tilde{\varphi}(\xi + n, \omega)|^2}{\sqrt{\omega^2 - \xi^2}} d\omega \end{aligned} \quad (52) $$

Let us estimate the second term in (52), the estimate of the first one being very similar. By a simple change of variable, we get

$$ \int_{-\infty}^{+\infty} d\xi \int_{|\omega|>|\xi|} \frac{|\omega||\tilde{\varphi}(\xi+n, \omega)|^2}{\sqrt{\omega^2 - \xi^2}} d\omega = \int_{-\infty}^{+\infty} d\xi \int_{|\omega|>|\xi-n|} \frac{|\omega||\tilde{\varphi}(\xi, \omega)|^2}{\sqrt{\omega^2 - (\xi-n)^2}} d\omega \quad (53) $$

Since $\varphi \in C_0^\infty(\mathbb{R}^2)$, its Fourier transform $\tilde{\varphi}$ is rapidly decreasing. Therefore, for every $k \in \mathbb{N}^+$ there exists $C_k > 0$ such that

$$ |\tilde{\varphi}(\xi, \omega)|^2 \le \frac{C_k}{(1 + \omega^2 + \xi^2)^{k+2}} $$

Note that if $(\xi, \omega) \in \mathscr{C}_n := \{(\xi, \omega) : |\omega| > |\xi - n|\}$ then $\omega^2 + \xi^2 \ge n^2/2$; it follows that for every $k \in \mathbb{N}^+$ there exists $C_k > 0$ such that

$$ |\tilde{\varphi}(\xi, \omega)|^2 \le \frac{C_k}{n^{2k} [1 + \omega^2 + \xi^2]^2} \quad \forall (\xi, \omega) \in \mathscr{C}_n $$

From this last inequality and owing to (53) we have

$$ \int_{\mathscr{C}_n} \frac{|\omega||\tilde{\varphi}(\xi, \omega)|^2}{\sqrt{\omega^2 - (\xi - n)^2}} d\omega d\xi \le \frac{C_k}{n^{2k}} \int_{\mathscr{C}_n} \frac{|\omega|}{\sqrt{\omega^2 - (\xi - n)^2}} \frac{1}{[1 + \omega^2 + \xi^2]^2} d\omega d\xi \quad (54) $$

This concludes the proof of (i) since, as it is not difficult to show, the last integral in (54) is finite.
We now turn to the property (ii). Having set $B = \{\omega^2 + \xi^2 \le 1\}$, by a simple change of variable, we obtain

$$
\begin{align*}
\| \varphi_n \|_{-s}^2 &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac{|\tilde{\varphi}_n(\xi, \omega)|^2}{[1+\omega^2+\xi^2]^s} \, d\omega \, d\xi \\
&= \frac{1}{4} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac{|\tilde{\varphi}(\xi-n, \omega) + \tilde{\varphi}(\xi+n, \omega)|^2}{[1+\omega^2+\xi^2]^s} \, d\omega \, d\xi \\
&\geq \frac{1}{4} \int_{\{\omega^2+(\xi-n)^2 \leq 1\}} \frac{|\tilde{\varphi}(\xi-n, \omega) + \tilde{\varphi}(\xi+n, \omega)|^2}{[1+\omega^2+\xi^2]^s} \, d\omega \, d\xi \\
&= \frac{1}{4} \int_B \frac{|\tilde{\varphi}(\xi, \omega) + \tilde{\varphi}(\xi+2n, \omega)|^2}{[1+\omega^2+(\xi+n)^2]^s} \, d\omega \, d\xi \\
&\geq \frac{1}{4(n-1)^{2s}} \int_B |\tilde{\varphi}(\xi, \omega) + \tilde{\varphi}(\xi+2n, \omega)|^2 \, d\omega \, d\xi \tag{55}
\end{align*}
$$

We remark that

$$
\begin{align*}
|\tilde{\varphi}(\xi, \omega) + \tilde{\varphi}(\xi + 2n, \omega)|^2 &= |\tilde{\varphi}(\xi, \omega)|^2 + |\tilde{\varphi}(\xi + 2n, \omega)|^2 + 2\Re[\tilde{\varphi}(\xi, \omega)\tilde{\varphi}(\xi + 2n, \omega)] \\
&\geq |\tilde{\varphi}(\xi, \omega)|^2 - 2|\tilde{\varphi}(\xi, \omega)||\tilde{\varphi}(\xi + 2n, \omega)| \\
&\geq \frac{1}{2}|\tilde{\varphi}(\xi, \omega)|^2 - 2|\tilde{\varphi}(\xi + 2n, \omega)|
\end{align*}
$$

Since for every $h>0$ there exists $C_h>0$ such that

$$
|\tilde{\varphi}(\xi + 2n, \omega)| \le \frac{C_h}{[1 + (\xi + 2n)^2]^h} \le \frac{C_h}{n^{2h}} \quad (|\xi| \le 1)
$$

it follows that

$$
\int_B |\tilde{\varphi}_n(\xi + 2n, \omega)|^2 d\omega d\xi \to 0 \quad \text{for } n \to +\infty
$$

Thus, for $n$ sufficiently large,

$$
\int_B |\tilde{\varphi}(\xi, \omega) + \tilde{\varphi}(\xi + 2n, \omega)|^2 d\omega d\xi \ge \frac{1}{4} \int_B |\tilde{\varphi}(\xi, \omega)|^2 d\omega d\xi \quad (56)
$$

Owing to (55) and (56), we conclude that there exists $C>0$ such that

$$
\| \varphi_n \|_{-s}^2 \ge \frac{C}{n^{2s}}
$$

which proves (ii).

<figure><img src="image_3.png" /></figure>

The main consequence of Theorem 3.1 is that in order to obtain satisfactory *a priori* bounds on the energetic Galerkin approximated solutions of problem (29), we need to complement the information coming from the quadratic form $a_g(\varphi, \varphi)$ with some other arising from alternative arguments.

In this respect, we present two possible strategies:

(1) to obtain a constraint on the oscillations in the space variable of the approximating solutions. In fact, as we shall see below, under a suitable constraint, the energetic bilinear form is coercive with respect to the $L^2(\Sigma_T)$ norm.
(2) to fill the gap in the frequency space that we have in formula (47), considering the energetic bilinear form modified by an additional term which takes into account also the skew-symmetric part of the operator $\varphi \mapsto (V\varphi)_t$. We shall briefly discuss this point at the end of this subsection.

In order to investigate the first point, let us fix $P>0$ and consider the following maximum problem:

$$ \lambda_0 = \max \left\{ \int_{-P}^{P} |\hat{f}(\omega)|^2 d\omega : f \in L^2(0, T), \int_{0}^{T} |f(t)|^2 dt = 1 \right\} $$

It is not difficult to see that $\lambda_0 = \lambda_0(PT) < 1$, $\lambda_0(\rho)$ is an increasing function and that $\lim_{\rho \to 0^+} \lambda_0(\rho) = 0$, $\lim_{\rho \to +\infty} \lambda_0(\rho) = 1$. The number $\lambda_0(PT)$ is the first eigenvalue of the so called time-band limited operator which has been studied indeph in the sixties by Slepian, Landau and Pollack (see e.g. [24]). Let us fix $\xi$ and set $P = |\xi|$, $f = \hat{\phi}(\xi, \cdot)$. We get

$$ \begin{aligned} \int_{|\omega| > |\xi|} |\tilde{\phi}(\xi, \omega)|^2 d\omega &= \int_{-\infty}^{\infty} |\tilde{\phi}(\xi, \omega)|^2 d\omega - \int_{-|\xi|}^{|\xi|} |\tilde{\phi}(\xi, \omega)|^2 d\omega \\ &\geq (1 - \lambda_0(|\xi|T)) \int_{-\infty}^{\infty} |\tilde{\phi}(\xi, \omega)|^2 d\omega \end{aligned} \quad (57) $$

Since $|\omega| > \sqrt{\omega^2 - \xi^2}$ for $|\omega| > |\xi|$, we obviously have from (47),

$$ a_{\varepsilon}(\varphi, \varphi) \geq \frac{1}{2} \int_{-\infty}^{+\infty} d\xi \int_{|\omega| > |\xi|} |\tilde{\phi}(\xi, \omega)|^2 d\omega \quad (58) $$

thus, integrating with respect to $\xi$ in (58) and owing to (57),

$$ a_{\varepsilon}(\varphi, \varphi) \geq \frac{1}{2} \int_{-\infty}^{+\infty} d\xi \int_{-\infty}^{+\infty} (1 - \lambda_0(|\xi|T)) |\tilde{\phi}(\xi, \omega)|^2 d\omega $$

Unfortunately, when $\rho \to +\infty$, the function $1 - \lambda_0(\rho)$ tends to zero exponentially. More precisely, Fuchs in [25] has shown that

$$ 1 - \lambda_0(\rho) \sim 4 \sqrt{\frac{\pi}{2}} \sqrt{\rho} e^{-\rho}, \quad \rho \to +\infty $$

Therefore, by the properties of $\lambda_0(\rho)$, we can only conclude that there exists a constant $C>0$ such that

$$ 1 - \lambda_0(\rho) \geq C \sqrt{1 + \rho} e^{-\rho} \quad (\rho > 0) \quad (59) $$

As a consequence, we obtain the following very poor coerciveness estimate

$$ a_{\varepsilon}(\varphi, \varphi) \geq C \int_{-\infty}^{+\infty} d\xi \int_{-\infty}^{+\infty} \sqrt{1 + |\xi|T} e^{-|\xi|T} |\tilde{\phi}(\xi, \omega)|^2 d\omega $$

Now, what can we say if the frequencies of the function $\varphi$ are not spread over the whole $\mathbb{R}^2$ and, for instance, $\tilde{\phi}$ has 'mass' concentrated on the strip $|\xi| \le R$, for some $R>0$?
We try to give a quantitative answer to this question in the following way. Let us fix $R>0$ and consider the following closed set in $L^2(\Sigma_T)$:

$$K_R := \left\{ \tilde{\varphi} \in L^2(\Sigma_T) : \int_{-\infty}^{+\infty} d\xi \int_{|\xi| \le R} |\tilde{\varphi}(\xi, \omega)|^2 d\omega \ge \frac{1}{2} \int_{-\infty}^{+\infty} d\xi \int_{-\infty}^{+\infty} |\tilde{\varphi}(\xi, \omega)|^2 d\omega \right\}$$

From (57), (58), and the above-mentioned properties of $\lambda_0$, we get

$$\begin{aligned} a_{\varepsilon}(\varphi, \varphi) &\ge \int_{-R}^{R} d\xi \int_{|\omega|>|\xi|} |\tilde{\varphi}(\xi, \omega)|^2 d\omega \ge \int_{-R}^{R} (1-\lambda_0(|\xi|T)) d\xi \int_{-\infty}^{+\infty} |\tilde{\varphi}(\xi, \omega)|^2 d\omega \\ &\ge (1-\lambda_0(RT)) \int_{-R}^{R} d\xi \int_{-\infty}^{+\infty} |\tilde{\varphi}(\xi, \omega)|^2 d\omega \end{aligned} \quad (60)$$

Thus, owing to (59) and the Parseval identity, we have just proved the following:

**Theorem 3.2**

For every $\varphi \in K_R \cap \mathcal{V}(a_{\varepsilon})$, we have

$$a_{\varepsilon}(\varphi, \varphi) \ge \frac{C}{2} \sqrt{1+RT} e^{-RT} \|\varphi\|_{L^2(\Sigma_T)}^2 \quad (61)$$

We can summarize as follows: if the approximating solutions $\varphi_h$ satisfy the constraint $\varphi_h \in K_R$ for some $R>0$, from the inequality (61) one gets a stability estimate in $L^2(\Sigma_T)$.

As an application of this argument, let us fix an integer $n \ge 1$, set $\Delta x = L/n$ and consider the functions

$$\varphi(x, t) = \sum_{k=0}^{n-1} f_k(t) \psi_k(x)$$

where $\psi_k(x) = H[(k+1)\Delta x - x] - H[k\Delta x - x]$ and $f_k(t) \in L^2(0, T)$, $k=0, \dots, n-1$. We denote by $\mathcal{V}_n$ the (infinite-dimensional) space of such functions and remark that $\mathcal{V}_n \subset L^2((0, T), H^{1/4}(\Gamma)) \subset \mathcal{V}(a_{\varepsilon})$.

**Lemma 3.1**

We have the inclusion $\mathcal{V}_n \subset K_R$, provided $R \ge 8L/(\pi \Delta x^2)$.

**Proof**

By computing the Fourier transform of $\varphi \in \mathcal{V}_n$, we get

$$\tilde{\varphi}(\xi, \omega) = \sqrt{\frac{2}{\pi}} e^{-i(\Delta x \xi / 2)} \frac{\sin\left(\frac{\Delta x \xi}{2}\right)}{\xi} \sum_{k=0}^{n-1} e^{-ik\Delta x \xi} \hat{f}_k(\omega)$$

thus

$$|\tilde{\varphi}(\xi, \omega)|^2 = \frac{2 \sin^2\left(\frac{\Delta x \xi}{2}\right)}{\xi^2} \left| \sum_{k=0}^{n-1} e^{-ik\Delta x \xi} \hat{f}_k(\omega) \right|^2 \le \frac{2n}{\pi} \frac{\sin^2\left(\frac{\Delta x \xi}{2}\right)}{\xi^2} \sum_{k=0}^{n-1} |\hat{f}_k(\omega)|^2$$
Since, for $\varphi \in \mathcal{V}_n$, $\|\varphi\|_{L^2(\Sigma_T)}^2 = \Delta x \sum_{k=0}^{n-1} \int_0^T |f_k(t)|^2 dt = \Delta x \sum_{k=0}^{n-1} \int_{-\infty}^{+\infty} |\hat{f}_k(\omega)|^2 d\omega$, we have

$$ \int_{-\infty}^{+\infty} d\xi \int_{|\xi|>R} |\tilde{\phi}(\xi, \omega)|^2 d\omega \le \frac{2n}{\pi \Delta x} \left( \int_{|\xi|>R} \frac{\sin^2\left(\frac{\Delta x \xi}{2}\right)}{\xi^2} d\xi \right) \|\varphi\|_{L^2(\Sigma_T)}^2 $$

The conclusion follows by the elementary inequality:

$$ \int_{|\xi|>R} \frac{\sin^2\left(\frac{\Delta x \xi}{2}\right)}{\xi^2} d\xi \le \frac{2}{R} $$

<figure><img src="image_4.png" /></figure>

Putting together the estimate (61) and Lemma 3.1 we obtain the following:

**Proposition 3.3**

Every $\varphi \in \mathcal{V}_n$ satisfies inequality (61) with $R=8L/(\pi \Delta x^2)$.

The theoretical value of the coerciveness constant in (61) is very small but we believe that it can be improved by a subtler analysis.

Finally, we sketch a few considerations on the second strategy, point (2) above. In general terms, we have the following situation: let A be a linear operator (bounded or not) on a Hilbert space; let us consider the equation

$$ A\varphi = A_s\varphi + A_{ss}\varphi = g \qquad (62) $$

where $A_s$ and $A_{ss}$ stand for the symmetric and the skew-symmetric part of A. Multiplying by $\varphi$ in (62), we obviously get

$$ a(\varphi, \varphi) = \langle A\varphi, \varphi \rangle = \langle A_s\varphi, \varphi \rangle = \langle g, \varphi \rangle $$

Now, if the quadratic form $a(\varphi, \varphi)$ enjoys some reasonable coerciveness property, one applies the usual arguments to prove the convergence of the Galerkin approximants, and eventually is allowed to neglect the skew-symmetric component of the operator. Unfortunately, this is not the case for the BIE treated in this paper. Thus, we pose the question: how can we resume information from the skew-symmetric part of the operator?

A possible answer could be that of choosing a suitable unitary operator $J$ and considering an auxiliary bilinear form $\bar{a}(\varphi, \psi) = \langle A\varphi, J\psi \rangle$, in order to modify the weak formulation of the problem introducing $b(\varphi, \psi) = a(\varphi, \psi) + \bar{a}(\varphi, \psi)$. In our case, a natural choice is provided by the Hilbert transform $\mathcal{H}$ with respect to the time variable. For our purposes, it is convenient to define the Hilbert transform as a Fourier multiplier:

$$ \mathcal{H}: L^2(\mathbb{R}^2) \to L^2(\mathbb{R}^2), \quad \mathcal{F}[\mathcal{H}](\xi, \omega) := \operatorname{sign}(\omega) \tilde{\phi}(\xi, \omega) $$

As it is well-known, one has

$$ \mathcal{H}^{-1} = \mathcal{H}^* = -\mathcal{H}, \quad \|\mathcal{H}\varphi\|_{H^s(\mathbb{R}^2)} = \|\varphi\|_{H^s(\mathbb{R}^2)} \quad \text{with } s \in \mathbb{R} \qquad (63) $$
Moreover, as a consequence of formula (46) we have for $\mathcal{H}\varphi \in \mathcal{S}(\mathbb{R}^2)$,

$$ \bar{a}(\varphi, \varphi) = \int_{0}^{+\infty} \int_{0}^{L} (V\varphi)_{t}(x, t) (\mathcal{H}\varphi)(x, t) \, dx \, dt = -\frac{1}{2} \int_{-\infty}^{+\infty} d\xi \int_{|\omega|<|\xi|} \frac{|\omega| |\tilde{\varphi}(\xi, \omega)|^2}{\sqrt{\xi^2 - \omega^2}} \, d\omega \quad (64) $$

From (46) and (48) we conclude that in our case the modified quadratic form would read

$$ \begin{aligned} b(\varphi, \varphi) &:= \int_{0}^{+\infty} \int_{0}^{L} (V\varphi)_{t}(x, t)(\varphi(x, t) - (\mathcal{H}\varphi)(x, t)) \, dx \, dt \\ &= \frac{1}{2} \int_{-\infty}^{+\infty} d\xi \int_{-\infty}^{+\infty} \frac{|\omega| |\tilde{\varphi}(\xi, \omega)|^2}{\sqrt{|\xi^2 - \omega^2|}} \, d\omega \end{aligned} \quad (65) $$

Regarding the coerciveness of the quadratic form $b(\varphi, \varphi)$, it holds, for instance

$$ 2b(\varphi, \varphi) \ge \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} v(|\omega|, |\xi|) |\tilde{\varphi}(\xi, \omega)|^2 d\xi d\omega, \quad v(|\omega|, |\xi|) := \frac{|\omega| \wedge |\xi|}{|\xi|} $$

Starting from here, we can obtain the coerciveness with respect to $L^2((0, T); H^{-1/2}(\Gamma))$ norm. In fact, for every $0<r\le 1$, evidently one has $v(|\omega|, |\xi|)>r/(|\xi|+1)$ for $|\omega|>r$; hence

$$ 2b(\varphi, \varphi) \ge r \int_{-\infty}^{+\infty} \int_{|\omega|>r} \frac{|\tilde{\varphi}(\xi, \omega)|^2}{|\xi|+1} d\xi d\omega \quad (66) $$

On the other side, remembering inequalities on band limited functions, one has, for every $\xi \in \mathbb{R}$,

$$ \int_{|\omega|>r} |\tilde{\varphi}(\xi, \omega)|^2 d\omega \ge (1-\lambda_0(rT)) \int_{-\infty}^{+\infty} |\tilde{\varphi}(\xi, \omega)|^2 d\omega $$

in particular choosing $r=1 \wedge (1/T)$, and applying (66), we have

$$ 2b(\varphi, \varphi) \ge [1 \wedge (1/T)](1 - \lambda_0(1)) \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} \frac{|\tilde{\varphi}(\xi, \omega)|^2}{|\xi|+1} d\xi d\omega $$

Hence, we can conclude with the following

### Theorem 3.3

There exists a positive constant $C$ such that, for every $\varphi \in C_0^\infty(\Sigma_T)$

$$ b(\varphi, \varphi) \ge \frac{C}{1+T} \|\varphi\|_{L^2((0,T); H^{-1/2}(\Gamma))}^2 \quad (67) $$

We note that, unless $\varphi=0$, it is impossible to have both $\varphi$ and $\mathcal{H}\varphi$ supported on the finite time interval $[0, T]$. However, thanks to (63), one can define, denoting by $\mathcal{I}$ the identity operator, the new bilinear form:

$$ b(\varphi, \psi) = \langle (\mathcal{I} + \mathcal{H})(V\varphi)_t, \psi \rangle_{L^2(\mathbb{R}^2)} $$

from which, one can derive, for test functions supported in $\Sigma_T$, the modified energetic weak formulation of the problem (29)

$$ \langle (\mathcal{I} + \mathcal{H})(V\varphi)_t, \psi \rangle_{L^2(\Sigma_T)} = \langle (\mathcal{I} + \mathcal{H})g_t, \psi \rangle_{L^2(\Sigma_T)} $$
We hope that starting from these suggestions, one can achieve a relevant progress toward a complete theory for the convergence of the energetic Galerkin BEM applied to the BIE (29).

### Remark

In the past years and in different contexts, several authors [13, 26, 27] have dealt with the properties of the Neumann pseudo-differential operator

$$ V: \left[ \frac{\partial u}{\partial n} \right] \longmapsto u, \quad (x, t) \in \Gamma \times (0, +\infty) \qquad (68) $$

related to the wave equation in the case of a flat boundary $\Gamma$. Owing to an obvious division problem in the symbol of the operator, any analysis has been done by using a Fourier-Laplace transform with nonvanishing imaginary part $\sigma$ in the phase variable $\tau+i\sigma$. In particular in Reference [13], devoted to the transient BIE for the acoustic equation, Ha Duong performs a detailed analysis of the symbol and obtains, under the restriction $\sigma \ge \sigma_0 > 0$, optimal results in terms of regularity and stability of the associated bilinear form. As shown in the same paper, this turns out to be equivalent to a coerciveness property of the functional

$$ \int_{0}^{T} \mathcal{E}(s, u) \, ds $$

On the other side, it should be noted that, as far as we know, in any numerical implementation of energy-related Galerkin methods for the transient wave equation, the parameter $\sigma$ has always been set equal to 0.

In the present paper we have chosen not to address the analysis of the Neumann operator (68), for which we refer to [13, 14], but to restrict our attention to the properties of the associated quadratic form $\mathcal{E}(T, u)$ (the energy at a fixed time $T$). Thus, the main difference with the cited papers is that we have assumed $\sigma=0$ showing that even in this case it is possible to obtain some stability results, as stated in Theorems 3.2 and 3.3.

## 3.2. Neumann problem

For completeness we introduce also a 2D Neumann model problem, for which we will give numerical results in the last section. Let us consider the scattering problem by a crack in an unbounded elastic isotropic medium. Let $\Omega=\mathbb{R}^2\setminus\Gamma$ be this medium and $\Gamma$ the crack represented by an open arc. Let $\Gamma^-$ and $\Gamma^+$ denote the lower and upper faces of the crack and $n$ the normal unit vector to $\Gamma$, oriented from $\Gamma^-$ to $\Gamma^+$. As usual, the total displacement field can be represented as the sum of the incident field (the wave propagating without the crack) and the scattered field. In a 3D elastic isotropic medium, there are three plane waves propagating in a fixed direction: the P wave, the SH wave and the SV wave. The 2D antiplane problem corresponds to an incident SH wave, when all quantities are independent of the third component $z$ (in particular, the crack has to be invariant with respect to $z$).

The scattered wave satisfies the following Neumann problem for the wave operator (without loss of generality we will consider a dimensionless problem which can be obtained after an appropriate
scaling of the units):

$$u_{tt} - \Delta u = 0, \quad \mathbf{x} \in \mathbb{R}^2 \setminus \Gamma, \quad t \in (0, T) \qquad (69)$$

$$u(\mathbf{x}, 0) = u_t(\mathbf{x}, 0) = 0, \quad \mathbf{x} \in \mathbb{R}^2 \setminus \Gamma \qquad (70)$$

$$\frac{\partial u}{\partial \mathbf{n}}(\mathbf{x}, t) = g(\mathbf{x}, t), \quad (\mathbf{x}, t) \in \Sigma_T := \Gamma \times [0, T] \qquad (71)$$

In (69)–(71) the unknown function $u$ stands for the third component of the displacement field and $g$ is the datum, which is the opposite of the normal derivative of the incident wave along $\Gamma$, i.e. $g = -\partial u^I / \partial \mathbf{n}$.

Let us consider the double-layer representation of the solution of (69)–(71):

$$u(\mathbf{x}, t) = \int_{\Gamma} \int_{0}^{t} \frac{\partial}{\partial \mathbf{n}_{\xi}} G(r, t-\tau) \phi(\xi, \tau) \, d\tau \, d\gamma_{\xi}, \quad \mathbf{x} \in \mathbb{R}^2 \setminus \Gamma, \quad t \in (0, T) \qquad (72)$$

where $G$ is given in (27) and $\phi = [u]$ is the jump of $u$ along $\Gamma$. Taking the normal derivative with respect to $\mathbf{x}$ of (72), with a limiting process for $\mathbf{x}$ tending to $\Gamma$ and using the assigned Neumann boundary condition, we obtained the space-time hypersingular BIE

$$\int_{\Gamma} \int_{0}^{t} \frac{\partial^2}{\partial \mathbf{n}_{\mathbf{x}} \partial \mathbf{n}_{\xi}} G(r, t-\tau) \phi(\xi, \tau) \, d\tau \, d\gamma_{\xi} = g(\mathbf{x}, t), \quad \mathbf{x} \in \Gamma, \quad t \in (0, T) \qquad (73)$$

which can be written with the compact notation

$$D\phi = g \qquad (74)$$

The energetic weak problem related to (74) will be of the form:

$$\tilde{a}_{\varepsilon}(\phi, \eta) = \langle g, \eta_t \rangle_{L^2(\Sigma_T)}$$

where

$$\tilde{a}_{\varepsilon}(\phi, \eta) := \langle D\phi, \eta_t \rangle_{L^2(\Sigma_T)} = \int_{\Gamma} \int_{0}^{T} (D\phi)(\mathbf{x}, t) \eta_t(\mathbf{x}, t) \, dt \, d\gamma_{\mathbf{x}}$$

and $\eta$ is a suitable test function belonging to the same functional space of $\phi$. The hypersingular integral operator $D$ can be equivalently expressed in the following way:

$$D\phi(\mathbf{x}, t) = \int_{\Gamma} \frac{\partial^2 r}{\partial \mathbf{n}_{\mathbf{x}} \partial \mathbf{n}_{\xi}} \int_{0}^{t} G(r, t-\tau) \left[ \phi_{t}(\xi, \tau) + \frac{\phi(\xi, \tau)}{t-\tau+r} \right] d\tau d\gamma_{\xi} \\ + \int_{\Gamma} \frac{\partial r}{\partial \mathbf{n}_{\mathbf{x}}} \frac{\partial r}{\partial \mathbf{n}_{\xi}} \int_{0}^{t} G(r, t-\tau) \left[ \phi_{tt}(\xi, \tau) + 2 \frac{\phi_{t}(\xi, \tau)}{t-\tau+r} + 3 \frac{\phi(\xi, \tau)}{(t-\tau+r)^2} \right] d\tau d\gamma_{\xi} \qquad (75)$$
4. GALERKIN BEM DISCRETIZATION

In the case of 1D problems, for instance with Dirichlet boundary conditions, we introduce, for the discretization phase, a uniform decomposition of the time interval $[0, T]$ with time step $\Delta t = T/N_{\Delta t}$, $N_{\Delta t} \in \mathbb{N}^{+}$, generated by the $N_{\Delta t}+1$ instants

$$t_k = k\Delta t, \quad k=0, \dots, N_{\Delta t}$$

Denoting by $\mathbb{P}_d, d \ge 0$, the space of polynomials of degree less than or equal $d$, we consider the standard finite element space

$$X_{\Delta t}^d = \{v_{\Delta t} \in L^2(0, T) : v_{\Delta t|_{[t_k, t_{k+1}]} \in \mathbb{P}_d, d \ge 0, k = 0, \dots, N_{\Delta t} - 1\}$$

Then, considering the finite dimensional space $W_{\Delta t}^d = X_{\Delta t}^d \times X_{\Delta t}^d \subset L^2(\Sigma_T)$, we can write down the discrete form of the energetic weak problem. Referring to (17) we have:

given $g \in H_{\{0\}}^1(\Sigma_T)$, find $\varphi_{\Delta t} \in W_{\Delta t}^d$ such that

$$a_{\varepsilon}(\varphi_{\Delta t}, \psi_{\Delta t}) = \langle g_t, \psi_{\Delta t} \rangle_{L^2(\Sigma_T)} \quad \forall \psi_{\Delta t} \in W_{\Delta t}^d \quad (76)$$

Denoting by $\{v_h\}$ a basis for $X_{\Delta t}^d$, the unknown function $\varphi_{\Delta t}$ can be expressed as:

$$\varphi_{0,\Delta t}(t) = \sum_{h} \varphi_{0,h} v_h(t), \quad \varphi_{L,\Delta t}(t) = \sum_{h} \varphi_{L,h} v_h(t)$$

and the discrete problem can be equivalently written as the linear system

$$A_{\varepsilon}x_{\varepsilon} = b_{\varepsilon} \quad (77)$$

in the unknowns the coefficients $\varphi_{0,h}$ and $\varphi_{L,h}$.

Remark

The energetic weak formulation for 1D problems has been introduced as a one-shot analysis till final time $T$. In fact, the linear system $A_{\varepsilon}x_{\varepsilon} = b_{\varepsilon}$ has a vector solution which allows to recover the whole time histories $\varphi_{\Delta t}(0, t)$ and $\varphi_{\Delta t}(L, t)$ till the final instant of the observation time interval. If instead we search the approximate solution step by step, for instance using constant or linear basis functions, i.e. we analyze the problem on each time interval of the type $(0, h\Delta t)$, for $h = 1, \dots, N_{\Delta t}$, we will have to solve $N_{\Delta t}$ linear systems of order two, with the same coefficient matrix, in the unknown vector $(\varphi_{0,h}\varphi_{L,h})^{\top}$, using the approximate solutions obtained in the previous steps to update the right-hand sides. This time marching procedure is formally equivalent to solve a global linear system having a lower block triangular matrix, with coincident diagonal blocks, using a block forward substitution algorithm.

For exterior 2D problems, we consider for a crack $\Gamma$ of length $L$, a boundary mesh constituted by $M$ straight elements $\{e_1, \dots, e_M\}$, with length$(e_i) \le \Delta x$, $e_i \cap e_j = \emptyset$ and such that $\bigcup_{i=1}^M \bar{e}_i$ coincides with $\bar{\Gamma}$ if the crack is (piece-wise) linear, or is a suitable approximation of $\bar{\Gamma}$, otherwise. The functional background compels one to choose spatially shape functions belonging to $L^2(\Gamma)$ for Dirichlet problems and to $H_0^1(\Gamma)$ for Neumann problems. Hence we use standard piecewise polynomial boundary element functions $w_j(\mathbf{x})$, $j = 1, \dots, M_{\Delta x}$, suitably defined in relation to the introduced mesh over $\Gamma$.
For time discretization we consider a uniform decomposition of the time interval $[0, T]$ with time step $\Delta t = T/N_{\Delta t}$, $N_{\Delta t} \in \mathbb{N}^{+}$, generated by the $N_{\Delta t} + 1$ instants

$$t_k = k\Delta t, \quad k = 0, \dots, N_{\Delta t}$$

as done for the 1D case, and we choose temporally piecewise constant shape functions for Dirichlet problems and piecewise linear shape functions for Neumann problems, although higher degree shape functions can be used. Note that, for this particular choice, our shape functions, denoted with $v_k(t)$, $k=0, \dots, N_{\Delta t}-1$, will be defined as

$$v_k(t) = H[t - t_k] - H[t - t_{k+1}]$$

for Dirichlet problems, or as

$$v_k(t) = R(t - t_k) - 2R(t - t_{k+1}) + R(t - t_{k+2})$$

where $R(t - t_k) = ((t - t_k)/\Delta t)H[t - t_k]$ is the ramp function, for Neumann problems. Hence, the unknown approximate solution of the problem at hand will be expressed as

$$\sum_{k=0}^{N_{\Delta t}-1} \sum_{j=1}^{M_{\Delta x}} \alpha_j^{(k)} w_j(\mathbf{x}) v_k(t) \quad (78)$$

The Galerkin BEM discretization coming from energetic weak formulation produces the linear system

$$A_{\mathcal{E}}x_{\mathcal{E}} = b_{\mathcal{E}} \quad (79)$$

In the case of a Dirichlet problem, the matrix elements, after a double analytic integration in the time variables, are of the form

$$\sum_{\alpha, \beta=0}^{1} (-1)^{\alpha+\beta} \int_{\Gamma} w_{i}(\mathbf{x}) \int_{\Gamma} \mathcal{B}(r, t_{h+\alpha}, t_{k+\beta}) w_{j}(\xi) \, d\gamma_{\xi} \, d\gamma_{\mathbf{x}} \quad (80)$$

where, having set $\Delta_{hk} = t_h - t_k$

$$\mathcal{B}(r, t_h, t_k) = -\frac{1}{2\pi} H[\Delta_{hk} - r] [\log(\Delta_{hk} + \sqrt{\Delta_{hk}^2 - r^2}) - \log r] \quad (81)$$

In the case of a Neumann problem, the matrix elements, after a double analytic integration in the time variables, are of the form

$$\sum_{\alpha, \beta, \delta=0}^{1} (-1)^{\alpha+\beta+\delta} \int_{\Gamma} w_{i}(\mathbf{x}) \int_{\Gamma} \mathcal{C}(r, t_{h+\alpha}, t_{k+\beta+\delta}) w_{j}(\xi) \, d\gamma_{\xi} \, d\gamma_{\mathbf{x}} \quad (82)$$

where

$$\mathcal{C}(r, t_h, t_k) = \frac{H[\Delta_{hk} - r]}{2\pi\Delta t} \left\{ \frac{\mathbf{r} \cdot \mathbf{n}_x \mathbf{r} \cdot \mathbf{n}_\xi}{r^2} \frac{\Delta_{hk} \sqrt{\Delta_{hk}^2 - r^2}}{r^2} + \frac{(\mathbf{n}_x \cdot \mathbf{n}_\xi)}{2} \left[ \log(\Delta_{hk} + \sqrt{\Delta_{hk}^2 - r^2}) - \log r - \frac{\Delta_{hk} \sqrt{\Delta_{hk}^2 - r^2}}{r^2} \right] \right\} \quad (83)$$
Anyway, the above elements depend on the difference $t_h - t_k$ and in particular they vanish if $t_h \le t_k$. Hence, matrix $A_\epsilon$ has a block lower triangular Toeplitz structure. Each block has dimension $M_{\Delta x}$. If we indicate with $A_\epsilon^{(\ell)}$ the block obtained when $t_h - t_k = (\ell+1)\Delta t$, $\ell=0, \dots, N_{\Delta t}-1$, the linear system can be written as

$$ \begin{pmatrix} A_{\epsilon}^{(0)} & 0 & 0 & \dots & 0 \\ A_{\epsilon}^{(1)} & A_{\epsilon}^{(0)} & 0 & \dots & 0 \\ A_{\epsilon}^{(2)} & A_{\epsilon}^{(1)} & A_{\epsilon}^{(0)} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ A_{\epsilon}^{(N_{\Delta t}-1)} & A_{\epsilon}^{(N_{\Delta t}-2)} & \dots & A_{\epsilon}^{(1)} & A_{\epsilon}^{(0)} \end{pmatrix} \begin{pmatrix} \alpha^{(0)} \\ \alpha^{(1)} \\ \alpha^{(2)} \\ \vdots \\ \alpha^{(N_{\Delta t}-1)} \end{pmatrix} = \begin{pmatrix} b_{\epsilon}^{(0)} \\ b_{\epsilon}^{(1)} \\ b_{\epsilon}^{(2)} \\ \vdots \\ b_{\epsilon}^{(N_{\Delta t}-1)} \end{pmatrix} \quad (84) $$

where

$$ \alpha^{(\ell)} = (\alpha_j^{(\ell)}) \quad \text{and} \quad b_{\epsilon}^{(\ell)} = (b_{\epsilon,j}^{(\ell)}) \quad \text{with } \ell=0, \dots, N_{\Delta t}-1, \quad j=1, \dots, M_{\Delta x} \quad (85) $$

The solution of (84) is obtained with a block forward substitution, i.e. at every time instant $t_\ell = (\ell+1)\Delta t$, $\ell=0, \dots, N_{\Delta t}-1$, we solve a reduced linear system of the type:

$$ A_{\epsilon}^{(0)} \alpha^{(\ell)} = b_{\epsilon}^{(\ell)} - (A_{\epsilon}^{(1)} \alpha^{(\ell-1)} + \dots + A_{\epsilon}^{(\ell)} \alpha^{(0)}) \quad (86) $$

Procedure (86) is a time-marching technique, where the only matrix to be inverted is the positive-definite $A_\epsilon^{(0)}$ diagonal block, while all the other blocks are used to update at every time step the right-hand side. Owing to this procedure we can construct and store only the blocks $A_\epsilon^{(0)}, \dots, A_\epsilon^{(N_{\Delta t}-1)}$ with a considerable reduction in computational cost and memory requirement.

Further, it can be proved that, for a fixed space discretization and for vanishing $\Delta t$, matrices $(\Delta t)^{-1} A_\epsilon^{(0)}$ for Dirichlet problems and $\Delta t A_\epsilon^{(0)}$ for Neumann problems, tend to the mass matrix $M^{(0)}$ of order $M_{\Delta x}$, with elements

$$ M_{ij}^{(0)} = \frac{1}{2} \int_{\Gamma} w_i(\mathbf{x}) w_j(\mathbf{x}) \, d\gamma_{\mathbf{x}} $$

## 5. NUMERICAL INTEGRATION ISSUES

Weakly singular and hypersingular double integrals in space variables (see (80), (82)) are efficiently calculated with numerical quadrature schemes widely used for BIEs related to elliptic problems [28], coupled with a suitable regularization technique, after a careful subdivision of the integration domain due to the presence of the Heaviside function.

Let us go into some details. In the sequel, we will refer to one of the double integrals in (82), i.e.

$$ \int_{\Gamma} w_i(\mathbf{x}) \int_{\Gamma} \mathcal{C}(r, t_h, t_k) w_j(\xi) \, d\gamma_{\xi} \, d\gamma_{\mathbf{x}} \quad (87) $$
<figure><img src="image_5.png" /><figcaption>Figure 1. Double integration domain (coincident elements) for different values of $\Delta_{hk}$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 2. Behavior of outer integrand function derivative for different values of $\Delta_{hk}$.</figcaption></figure>

Further, using the standard element by element technique, the evaluation of every double integral of the form (87) is reduced to the assembling of local contributions of the type

$$ \int_{e_i} \tilde{w}_i^{(d_i)}(\mathbf{x}) \int_{e_j} \mathcal{C}(r, t_h, t_k) \tilde{w}_j^{(d_j)}(\xi) \, d\gamma_{\xi} \, d\gamma_{\mathbf{x}} \quad (88) $$

where $\tilde{w}_i^{(d_i)}(\mathbf{x})$ defines one of the local lagrangian basis function in the space variable of degree $d_i$ defined over the element $e_i$ of the boundary mesh.

Looking at the kernel $\mathcal{C}(r, t_h, t_k)$ in (83), we observe space singularities of type $\log r$ or $O(r^{-2})$ as $r \to 0$, which are typical of 2D static weakly singular and hypersingular kernels. We remark that when the kernel is hypersingular and $e_i \equiv e_j$ we have to define both the inner and the outer integrals as Hadamard finite parts, while if $e_i$ and $e_j$ are consecutive, only the outer integral must be understood in the finite part sense. The correct interpretation of double integrals is the key point for any efficient numerical approach based on element by element technique (see [28]).

Further, in (83) we find the Heaviside function $H[\Delta_{hk}-r]$ and the square root $\sqrt{\Delta_{hk}^2-r^2}$, which are responsible for different types of troubles, that we will illustrate for the case of coincident elements $e_i \equiv e_j$ of length $2l_i$, where $r = |s-z|$ in the local variables of integration.
<figure><img src="image_7.png" /><figcaption>Figure 3. Behavior of $\partial/\partial z\sqrt{\Delta_{hk}^2 - |s-z|^2}$ for $\Delta_{hk}=0.15, 0.05, 0.025$.</figcaption></figure>

In Figure 1 we show the double integration domain represented by the intersection between the square $[0, 2l_i] \times [0, 2l_i]$ and the strip $|s-z|<\Delta_{hk}$ where the Heaviside function is not trivial, for different values of $\Delta_{hk}$.

The numerical quadrature in the outer variable of integration $s$ has been optimally performed subdividing, when necessary, the outer interval of integration, as shown in the same figure. In fact, the derivative with respect to $s$ of the outer integrand function, after the inner integration, presents jump in correspondence to the subdivision points, as depicted in Figure 2 for the function

$$ \frac{d}{ds} \int_{0}^{2l_i} H[\Delta_{hk} - |s-z|] \log(\Delta_{hk} + \sqrt{\Delta_{hk}^2 - |s-z|^2}) \tilde{w}_j^{(0)}(z) dz $$

Without this subdivision, one should use a lot of quadrature nodes to achieve the single precision accuracy.

Then, after subdividing, when necessary, the outer interval of integration in suitable subdomains, the inner numerical integration is still difficult, due to the presence of the square root function in the kernel. In fact, the argument of $\sqrt{\Delta_{hk}^2 - |s-z|^2}$ is always positive but it can assume very small values and in the limit for the argument tending to zero the derivative of the square root with respect to the inner variable of integration $z$ becomes unbounded. This behavior happens along the oblique boundary of the double integration domain, as shown in Figure 3 for different values of $\Delta_{hk}$, and produces a bad performance, for instance in the evaluation of the integral

$$ \int_{0}^{2l_i} \tilde{w}_{i}^{(0)}(s) \int_{0}^{2l_i} H[\Delta_{hk} - |s-z|] \log(\Delta_{hk} + \sqrt{\Delta_{hk}^2 - |s-z|^2}) \tilde{w}_{j}^{(0)}(z) dz ds \quad (89) $$

even of the classical Gauss-Legendre quadrature formula, in the sense that one should use a lot of quadrature nodes to achieve the single precision accuracy. To overcome this difficulty, we have considered the regularization procedure introduced in [29], which suitably pushes the Gaussian nodes towards the end-points of the inner interval of integration and modify the Gaussian weights in order to regularize integrand functions with mild boundary 'singularities'.

In Figure 4, we show the computational costs of Gaussian quadrature formula and regularization procedure in relation to the achievement of the single precision accuracy in the evaluation of the double integral (89) for $\Delta_{hk}=0.15, 0.05, 0.025$.

The numerical treatment of weak singularities and hypersingularities, respectively of type $\log r$ or $O(r^{-2})$ as $r \to 0$, has been operated through the quadrature schemes proposed in [28] and widely used in the context of Galerkin BEM related to BIEs coming from elliptic problems. In Figure 5
<figure><img src="image_8.png" /><figcaption>Figure 4. Computational cost of Gaussian quadrature and regularization procedure, for $\Delta_{hk}=0.15, 0.05, 0.025$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 5. Computational cost of HFP quadrature coupled with regularization procedure, for $\Delta_{hk}=0.15, 0.05, 0.025$.</figcaption></figure>

we show the performance of the Hadamard Finite Part quadrature formula coupled, if necessary, with the regularization technique cited above, with respect to the single precision accuracy in the numerical evaluation of the hypersingular integral

$$ \int_{0}^{2l_i} \tilde{w}_{i}^{(1)}(s) \int_{0}^{2l_i} H[\Delta_{hk} - |s-z|] \frac{\Delta_{hk} \sqrt{\Delta_{hk}^2 - |s-z|^2}}{|s-z|^2} \tilde{w}_{j}^{(1)}(z) dz ds $$

for $\Delta_{hk}=0.15, 0.05, 0.025$.

The above considerations have to be done, in the element-by-element technique, also for the other geometrical disposition of the double integration mesh elements $e_i, e_j$. In particular, if we consider contiguous elements $e_i, e_{i+1}$, with length respectively $2l_i, 2l_{i+1}$, forming an angle $\theta$, with $0<\theta<\pi$, the double integration domain is constituted by the intersection between the rectangle $[0, 2l_i] \times [0, 2l_{i+1}]$ and the 2D domain $r^2 - \Delta_{hk}^2 < 0$, that for not aligned elements is an ellipsis centered in the unique singularity point $(2l_i, 0)$. The directions of the two axes of the ellipsis are $(1, 1)$ and $(-1, 1)$ with semi-length, respectively $(1+\cos(\theta))^{-1/2}$ and $(1-\cos(\theta))^{-1/2}$. In Figure 6, we show various types of intersections, i.e. double integration domains, for different values of $\Delta_{hk}$ and different angles $\theta$ between contiguous elements.
<figure><img src="image_10.png" /><figcaption>Figure 6. Double integration domain (contiguous elements) for different values of $\Delta_{hk}$ and different angles $\theta$.</figcaption></figure>

A specific illustration of the efficient numerical integration schemes we have used for the discretization of weakly singular and hypersingular BIEs related to wave propagation problems will be the subject of a forthcoming paper.

## 6. NUMERICAL RESULTS

As first 1D test problem, let us consider a 1D domain of length $L = \pi/2$, subject to the Dirichlet boundary conditions:

$$u(0, t) = 0, \quad u(L, t) = \sin^2(t) H[t] H\left[\frac{\pi}{2} - t\right] + H\left[t - \frac{\pi}{2}\right]$$

The observation time interval is $[0, 6L]$. For the energetic Galerkin BEM discretization, we have considered time steps of the type $\Delta t = \pi/(2N_{\Delta t})$, $N_{\Delta t} \in \mathbb{N}^+$, such that the time $\pi/2$ required by the wave, travelling with unitary speed, to cover the distance between the two end-points of the domain is a multiple of $\Delta t$. Traction in $x=0$ and $x=L$ has been approximated by linear shape functions. The numerical solution obtained solving (77), with $\Delta t = \pi/32$, is reported in Figure 7. On the other hand, if the observation time interval is $[0, 10]$ it is natural to choose the time step $\Delta t = 0.1$, even if the time $\pi/2$ required by the wave to cover the distance between the two end-points of the domain is not a multiple of it. With this choice, classical Galerkin BEM coming from $L^2$ weak formulation would produce huge instability phenomena [20]. This does not happen with the energetic Galerkin BEM: in Figure 8 we show the numerical solution obtained solving (77), using $\Delta t = 0.1$ and linear shape functions.

The energetic procedure, applied to various 1D wave propagation problems [20], appears to be unconditionally stable. In fact, even though the stability constant $1/c(T)$, with $c(T)$ given as in Theorem 2.1, has an asymptotic behavior of the type $O(T^2)$ when $T \to \infty$, the approximate
<figure><img src="image_11.png" /><figcaption>Figure 7. Piecewise linear approximations $\varphi_{\pi/32}(0, t)$, $\varphi_{\pi/32}(L, t)$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 8. Piecewise linear approximations $\varphi_{0.1}(0, t)$, $\varphi_{0.1}(L, t)$.</figcaption></figure>

solution remains stable even for large times. In order to show this particular feature, we consider a second Dirichlet test problem on the time interval $[0, 100]$, related to a 1D domain of length $L=1$, subject to the boundary conditions:

$$u(0, t) = 0, \quad u(L, t) = H[t]$$
<figure><img src="image_13.png" /><figcaption>Figure 9. Numerical stability of the energetic approach for large times for the second 1D test problem.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 10. Numerical stability of the energetic approach for large times for the third 1D test problem.</figcaption></figure>

In Figure 9 we show the numerical solution of the test problem, obtained starting from energetic weak formulation using $\Delta t=0.1$ and constant shape functions, till the final time $T=100$. It tries to reproduce the Dirac functions appearing at $x=0$ and $x=1$ with unitary time period.

Further, Figure 10 presents the numerical solution referred to another Dirichlet test problem on the time interval $[0, 100]$, related to the previous 1D domain subject to the boundary conditions:

$$u(0,t)=0, \quad u(L,t)=H[t]\sin(t)$$
The new challenging features of the datum $u(L, t)$ with respect to the previous one are: the fact that it never collapses to a constant for $t>\bar{t}>0$ and it presents a jump in the first derivative for $t=0$, which is reproduced in the solution after reflection at time intervals of unitary length. The approximation $\varphi_{\Delta t}$ is carried out starting from the energetic weak formulation, using $\Delta t=0.1$ and constant shape functions in the discretization phase. In this case, the analytical solution is

$$ \varphi_0(t) = -2 \sum_{k=0}^{49} H[t-2k-1] \cos(t-2k-1), \quad \varphi_L(t) = 2 \sum_{k=0}^{50} H[t-2k] \cos(t-2k) $$

and it was possible to evaluate the approximation error, which turned out to be

$$ \max_{t \in (0, 100)} |\varphi(t) - \varphi_{\Delta t}(t)| \approx 8.98 \times 10^{-4} $$

The interested reader is referred to [30] for many other numerical results related to 1D wave propagation analysis in layered media always obtained with the energetic Galerkin BEM approach.

In the following, we will present numerical results obtained for 2D both Dirichlet and Neumann problems, starting from the proposed energetic weak formulation.

As first 2D test problem, we consider (23)–(25) with $\Gamma = \{(x, 0), x \in [0, 1]\}$ and Dirichlet boundary datum

$$ g(x, t) = H[t-kx]f(t-kx) \quad \text{where} \quad f(z) = \begin{cases} \sin^2\left(\frac{\omega z}{2}\right) & \text{if } 0 \le z \le \frac{\pi}{\omega} \\ 1 & \text{if } z \ge \frac{\pi}{\omega} \end{cases} \quad (90) $$

with $\omega=8\pi$, $k=\cos\theta$ and $\theta \in (0, \pi)$. For this problem we have chosen a uniform decomposition of $\Gamma$ constituted by 40 nods ($\Delta x=0.025$) and a subdivision of the time interval $[0, 2]$ operated by 160 instants ($\Delta t=0.0125$). In Figure 11 we present the numerical solution, that in the sequel will be indicated with $\varphi(x, t)$, obtained in some points of $\Gamma$ for different values of the angle $\theta$ and considering spatial constant shape and test functions. When $\theta=\pi/2$, the excitation field is uniform on $\Gamma$ at initial time, hence symmetric points of $\Gamma$ behave in the same way; further, note that points $x=\frac{1}{4}$, $x=\frac{3}{4}$ and $x=\frac{1}{2}$ are excited till the time instant $\pi/\omega=\frac{1}{8}$, but while the first two are influenced at $t=\frac{1}{4}$ by the wave coming from the endpoints of $\Gamma$ and travelling with unitary velocity, the midpoint of $\Gamma$ is influenced at $t=\frac{1}{2}$. When $\theta=\pi/4$, the excitation field is not uniform on $\Gamma$ at initial time; hence, differently from the previous case, symmetric points of $\Gamma$ do not behave in the same way at the beginning of the simulation.

Now we consider, on the same domain of the previous test problem, the Dirichlet datum shown in Figure 12, and we fix the observation time interval $[0, 10]$. As uniform temporal and spatial discretization steps we use $\Delta t=0.1$ and $\Delta x=0.0125$, respectively, and we adopt spatial constant shape and test functions. In Figure 13 we show the time history of the density $\varphi$ obtained in the point $x=\frac{1}{2}$. As one can note, it has the same form of the derivative of the boundary datum and hence it vanishes for long times. In Figure 14 we present a section of the solution $u(x, y, t)$ of the problem (23)–(25), in the same point $x=\frac{1}{2}$, for $t \in [0, 10]$: as one can observe, the wave travelling away from the boundary $\Gamma$ assumes the same structure of the Dirichlet boundary datum but with diminishing intensity. Figure 15 shows the time-dependent behavior of the solution $u(x, t)$ outside $\Gamma$.
<figure><img src="image_15.png" /><figcaption>Figure 11. Density $\varphi(x, t)$ in some points of $\Gamma$ for different values of $\theta$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 12. Dirichlet boundary condition.</figcaption></figure>

Further, let us consider problem (23)–(25) assigned on the semi-circular arc

$$ \Gamma = \{ \mathbf{x} \in \mathbb{R}^2 : \mathbf{x} = (\cos \alpha, \sin \alpha), \alpha \in [0, \pi] \} $$

depending on the clockwise angle $\alpha$, and with Dirichlet boundary datum

$$ g(\alpha, t) = H[t] f(t) \cos \alpha $$

where $f$ has been already given in (90) with $\omega=8\pi$. The observation time interval is $[0, 10]$. As uniform temporal discretization step we use $\Delta t=0.1$ and $\Gamma$ is uniformly approximated by 40
<figure><img src="image_17.png" /><figcaption>Figure 13. Density $\varphi(\frac{1}{2}, t)$ obtained for $\Delta t = 0.1$ and $\Delta x = 0.0125$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 14. Section of the solution $u(x, y, t)$ in $x = \frac{1}{2}$ in some time instants.</figcaption></figure>

straight boundary elements where we adopt constant shape and test functions. In Figure 16, we show the approximate solution $\varphi(\alpha, 10)$ at the final instant of analysis. Since the Dirichlet datum becomes independent of time, the transient solution on $\Gamma$ tends to the stationary one $\varphi_{\infty}(\alpha)$, which
<figure><img src="image_19.png" /><figcaption>Figure 15. Solution $u(x, t)$ outside $\Gamma$ in some time instants.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 16. Densities $\varphi(\alpha, 10)$, obtained on $\Gamma$ at final time $T=10$, and $\varphi_\infty(\alpha)$.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 17. Graph of the time function $||\varphi(\cdot, t) - \varphi_{\infty}(\cdot)||_{L^1(\Gamma)}$.</figcaption></figure>

is the solution of the BIE related to the following problem:

$$ 
\begin{array}{ll}
-\Delta u_{\infty} = 0 \text{ in } \mathbb{R}^2 \setminus \Gamma, & u(\mathbf{x}) = O(1) \text{ for } \|\mathbf{x}\|_2 \to \infty \\
u_{\infty} = \cos \alpha & \text{on } \Gamma
\end{array}
\qquad (91) $$

Also $\varphi_{\infty}(\alpha)$ has been reported in Figure 16 and, as one can observe, the two curves overlap. In Figure 17 we have reported the graph of $||\varphi(\cdot, t) - \varphi_{\infty}(\cdot)||_{L^1(\Gamma)}$, which is a function of time.

### Remark

In Figure 18, we show the typical behavior of the minimum eigenvalue of the symmetric part of matrices $A_{\varepsilon}$, which is related to the coercivity constant of the energetic bilinear form, for fixed $\Delta x$ and vanishing $\Delta t$; the graphs are presented for Dirichlet problems outside rectilinear obstacles of length $L=0.2$, $L=0.4$, respectively, having fixed $\Delta x=0.1$ for the spatial discretization and space-time piecewise constant shape and test functions; in these cases, the final observation time is $T=1$. The obtained curves are in agreement with theoretical result stated in Proposition 3.3.

Now, we present numerical results related to 2D Neumann exterior problem (69)–(71), for $\Gamma = \{(x, 0) : x \in [0, 1]\}$. The incident wave $u^I(\mathbf{x}, t)$ is a plane wave propagating in direction $\mathbf{k}$ with unitary amplitude:

$$ u^I(\mathbf{x}, t) = f(t - \mathbf{k} \cdot \mathbf{x}) \quad \text{with} \quad \mathbf{k} = (\cos \theta, \sin \theta) \qquad (92) $$

Hence, the Neumann datum on $\Gamma$ in (71) will be:

$$ g(x, t) = - \left. \frac{\partial}{\partial \mathbf{n}_{\mathbf{x}}} f(t - \mathbf{k} \cdot \mathbf{x}) \right|_{\Gamma} \qquad (93) $$

We show the results obtained for two different functions, that have been chosen for the known asymptotic behavior of the solution, which allows us to validate the approximate solution, in the sequel indicated with $\phi(x, t)$.
<figure><img src="image_22.png" /><figcaption>Figure 18. Typical behavior of the minimum eigenvalue of $A_ε$ symmetric part.</figcaption></figure>

Plane harmonic wave: we consider a wave which becomes harmonic after a fixed time (see [15]):

$$f(t) = \begin{cases} 0 & \text{if } t < 0 \\ \frac{1}{2}(1 - \cos \omega t) & \text{if } 0 \le t \le \frac{\pi}{\omega} \\ \sin \left(\frac{\omega t}{2}\right) & \text{if } t \ge \frac{\pi}{\omega} \end{cases} \qquad (94)$$

where $\omega$ represents the frequency. In this case the solution has to also become harmonic, with the same period as the incident wave, i.e. $P=2\pi/\tilde{\omega}$, where $\tilde{\omega}=\omega/2$. The fixed circular frequency $\omega=8\pi$ is such that the wavelength $\lambda=2\pi/\omega$ is equal to a quarter the crack length.

We choose a uniform decomposition of the crack $\Gamma$ in 20 subintervals ($\Delta x=0.05$) and we decompose the observation time interval [0, 10] in 400 equal parts ($\Delta t=0.025$). For this numerical simulation, we choose spatial linear shape and test functions.

In Figure 19 we show the time harmonic behavior for $\theta=\pi/2$ of the crack opening displacement (COD) $\phi$ at $x=0.4$, obtained starting from the energetic weak formulation. Note that the solution becomes immediately not trivial since the incident wave strikes the whole crack simultaneously. In Figure 20 we present the approximated COD at instants 2, 4, 5, 10; the four curves overlap each other since the period is $P=0.5$.

In Figure 21 we show the time harmonic behavior for $\theta=\pi/3$ of the crack opening displacement $\phi$ at $x=0.4$. Note that the COD is zero till the time instant $t^*=\cos(\theta)x=x/2$, since the incident wave, differently from the previous case, does not invest the whole crack simultaneously.

In order to verify that the period of $\phi$ coincides with the period of the incident wave, we show in Figure 22 the approximate solution $\phi$ on $\Gamma$ in time instants separated by multiples of the time period.

Plane linear wave: let us consider another example, always taken from [15], where the Neumann boundary conditions comes from this choice of $f$:

$$f(t) = tH[t]$$
<figure><img src="image_23.png" /><figcaption>Figure 19. Density $\phi(0.4, t)$ for $\theta = \pi/2$.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 20. COD at $t=2, 4, 5, 10$ for $\theta = \pi/2$.</figcaption></figure>

In this case, the Neumann datum (93) tends to the constant value $g_{\theta} = \sin\theta$ when $t$ tends to infinity.
The solution $u$ tends to the solution of the static problem

$$ 
\begin{align} 
-\Delta u_{\infty} &= 0 \quad \text{in } \mathbb{R}^2 \setminus \Gamma, \quad u(\mathbf{x}) = O(\|\mathbf{x}\|_2^{-1}) \quad \text{for } \|\mathbf{x}\|_2 \to \infty \tag{95} \\ 
\frac{\partial u_{\infty}}{\partial \mathbf{n}} &= g_{\theta} \quad \text{on } \Gamma \notag 
\end{align} 
$$
<figure><img src="image_25.png" /><figcaption>Figure 21. Density $\phi(x, t)$ calculated in $x=0.4$ for $\theta=\pi/3$, with a zoom.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 22. Solution $\phi(x, t)$ in different time instants separated by $5P$.</figcaption></figure>

and the associated jump $\phi_{\theta}^{\infty}(x)=[u_{\infty}]$ across $\Gamma$ can be computed explicitly:

$$ \phi_{\theta}^{\infty}(x) = \sin\theta\sqrt{x(1-x)} $$

Hence, we can compare the solution $\phi(x, t)$ with the static solution $\phi_{\theta}^{\infty}(x)$.
<figure><img src="image_27.png" /><figcaption>Figure 23. Density $\phi(x, t)$ evaluated in some points of $\Gamma$ for $\theta = \pi/3$.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 24. Solution $\phi(x, t)$ in different time instants.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>Figure 25. Density $\phi(x, t)$ evaluated in some points of $\Gamma$ for $\theta = \pi/4$.</figcaption></figure>

Let us consider the crack $\Gamma$ decomposed in 20 subintervals ($\Delta x = 0.05$), the final observation time $T = 10$, the time step $\Delta t = 0.025$ and spatial linear shape and test functions. Also in this case, the numerical solution obtained for large times and reported in Figure 23 for $\theta = \pi/3$ is in perfect agreement with the corresponding one reported in [15]. Note that points of $\Gamma$, symmetric with respect to $x = 0.5$, behaves in different ways at the beginning of the simulation, since the incident wave does not strike simultaneously the crack; then they assume a symmetric behavior for sufficiently large times. Finally, to verify that the whole numerical solution $\phi(x, t)$ stabilizes itself, we show in Figure 24 the COD on $\Gamma$ in different time instants $t \ge 4$.

In Figure 25, we show analogous numerical results obtained for $\theta = \pi/4$, and in Figure 26 we present the approximate COD for $T = 4, 5, 10$ together with the analytical solution of the corresponding static problem: the four curves overlap each other. Numerical results are in agreement with those reported in [15], but differently from what is said in that paper, they seem to be independent of the ratio $\Delta t/\Delta x$.

At last, to complete numerical simulations, we present some results involving the total displacement field $u(x, t)$ obtained by the superposition of the incident wave $u^I(x, t)$ and the reflected and diffracted waves caused by the presence of a crack $\Gamma \subset \mathbb{R}^2$. The temporal profile of the incident wave, that strikes the crack and from which we have deduced the Neumann datum on $\Gamma$, is shown in Figure 27 and it is similar to that one considered in [31, 32].

In the first simulation we deal with the crack $\Gamma = \{(x, 0), x \in [-1, 1]\}$ struck perpendicularly by the incident wave. The observation time interval is $[0, 4]$. For the discretization, we fix a uniform subdivision of $\Gamma$ in 40 elements ($\Delta x = 0.05$), the time step $\Delta t = 0.1$ and spatial linear shape and test functions.

In Figure 28 we show the total recovered displacement in a square around the crack for different time instants. These results show how the plane wave reaches the crack and how the diffraction caused at the edges of the crack degenerates the wavefront. The effect of the diffraction on the upper half of the square creates a shadow. On the other hand, diffraction can be observed on the lower half of the square, too. Note that at the beginning of the simulation, the reflected wave on the upper half of the square cancel out with the incident wave. As time increases, the wavefront recovers and the scattering effect caused by the crack on the plane wave diminishes.
<figure><img src="image_30.png" /><figcaption>Figure 26. COD at $t=4, 5, 10$ compared with static solution for $\theta=\pi/4$.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Figure 27. Temporal profile of the incident wave in the last two examples.</figcaption></figure>

In the second simulation we consider the semi-circular arc

$$ \Gamma = \{ \mathbf{x} \in \mathbb{R}^2 : \mathbf{x} = (\cos \alpha, \sin \alpha), \alpha \in [0, \pi] \} $$

depending on the clockwise angle $\alpha$, struck by the plane wave with an incident angle of amplitude $\pi/3$. The observation time interval is $[0, 4]$. As uniform temporal discretization step we use $\Delta t = 0.1$ and $\Gamma$ is uniformly approximated by 40 straight boundary elements where we adopt spatial linear
<figure><img src="image_32.png" /><figcaption>Figure 28. Total recovered displacement $u(x, t)$ around the straight crack.</figcaption></figure>

shape and test functions. In Figure 29 we show the total recovered displacement in a square around the crack for different time instants. We are able to identify at $t=0.1$ the incident wave, at $t=1.0$ the reflected wave, at $t=1.6$ the first circular diffracted wave generated at the left edge of the crack, at $t=2.5$ the first diffracted wave generated at the right edge of the crack.
<figure><img src="image_33.png" /><figcaption>Figure 29. Total recovered displacement $u(x, t)$ around the semicircular crack.</figcaption></figure>

# APPENDIX A

## Proof of Theorem 2.1

Let us introduce the modified anticipated and retarded shift operators:

$$ (S_{+}f)(t) := H[T-L-t]f(t+L), \quad (S_{-}f)(t) := H[t-L]f(t-L) $$
Denoting by $\mathcal{R}$ the reflection matrix $[0, 1; 1, 0]$, from a straightforward calculation we obtain

$$A^* \varphi(t) = \frac{1}{2} [\varphi(t) + S_+ \mathcal{R} \varphi(t)]$$

hence the symmetric part of $A$ has the following expression:

$$2A_s \varphi(t) = \varphi(t) + \frac{1}{2}(S_- + S_+) \mathcal{R} \varphi(t), \quad t \in [0, T] \qquad (A1)$$

There is a clear analogy between the formulation (A1) of the operator $A_s$ and the usual second-order finite differences scheme. In what follows we proceed to make clear this analogy and eventually we shall reduce the action of $A_s$ to that of a matrix operator closely connected to finite differences. First of all we may assume $T=NL$ for some positive integer $N$, otherwise we would trivially extend $\varphi$ up to $NL$ where $N = \lceil T/L \rceil$. Then we reformulate (A1) in the following equivalent way:

$$2(S_+^k A_s \varphi)(t) = S_+^k \varphi(t) + \frac{1}{2}(S_- + S_+) \mathcal{R} S_+^k \varphi(t), \quad t \in [0, L], \quad k=0, 1, \dots, N-1 \qquad (A2)$$

Let us define, for $k=0, 1, \dots, N-1$,

$$\Phi_k(t) := (\Phi_{0,k}(t), \Phi_{L,k}(t))^{\top} := \mathcal{R}^k S_+^k \varphi(t), \quad t \in [0, L]$$

and

$$\tilde{\Phi}(t) := (\tilde{\Phi}_0(t), \tilde{\Phi}_L(t)), \quad t \in [0, L]$$

where

$$\tilde{\Phi}_0(t) := (\Phi_{0,k}(t))_{k=0,1,\dots,N-1} = (\varphi_0(t), S_+\varphi_L(t), S_+^2\varphi_0(t), \dots, S_+^{N-1}\varphi_{m_1L}(t)), \quad t \in [0, L]$$

$$\tilde{\Phi}_L(t) := (\Phi_{L,k}(t))_{k=0,1,\dots,N-1} = (\varphi_L(t), S_+\varphi_0(t), S_+^2\varphi_L(t), \dots, S_+^{N-1}\varphi_{m_2L}(t)), \quad t \in [0, L]$$

with $m_1 = [(N-1) - 2\lfloor(N-1)/2\rfloor]$ and $m_2 = m_1 + (-1)^{N-1}$. Similarly, for the left-hand side in (A2), we set

$$\Psi_k(t) := (\Psi_{0,k}(t), \Psi_{L,k}(t))^{\top} := 2(\mathcal{R}^k S_+^k A_s \varphi)(t) \quad \text{and} \quad \tilde{\Psi}(t) := (\tilde{\Psi}_0(t), \tilde{\Psi}_L(t)), \quad t \in [0, L]$$

Finally, we observe that multiplying both sides of formula (A2) by $\mathcal{R}^k$, we obtain $(\Phi_{-1}(t) \equiv \Phi_N(t) \equiv 0)$:

$$\Psi_k(t) = \frac{1}{2}\Phi_{k-1}(t) + \Phi_k(t) + \frac{1}{2}\Phi_{k+1}(t), \quad t \in [0, L], \quad k=0, 1, \dots, N-1$$

Thus, each component of $\Psi(t)$ can be expressed as

$$\tilde{\Psi}_0(t) = \mathcal{M}\tilde{\Phi}_0(t), \quad \tilde{\Psi}_L(t) = \mathcal{M}\tilde{\Phi}_L(t)$$

where $\mathcal{M} = \text{tridiag}[\frac{1}{2}, 1, \frac{1}{2}]$ is a tridiagonal matrix of order $N$. Note that the only difference between the matrix $\mathcal{M}$ and the usual finite difference matrix is the sign of the elements equal to $\frac{1}{2}$. In fact, the two matrices are similar through the diagonal matrix which alternates 1 and $-1$ along the principal diagonal. It is well known that the $N \times N$ finite difference matrix is positive definite and its spectrum is given by the $N$ eigenvalues [33]

$$\omega_k^2 = 2 \sin^2 \left( \frac{k\pi}{2(N+1)} \right), \quad k=1, \dots, N$$
Remembering (14), the conclusion follows from the following identity:

$$ \int_{0}^{T} |\varphi(t)|^2 dt = \sum_{k=0}^{N-1} \int_{0}^{L} |\mathcal{R}^k S_{+}^k \varphi(t)|^2 dt = \sum_{k=0}^{N-1} \int_{0}^{L} |\Phi_k(t)|^2 dt = \int_{0}^{L} |\tilde{\Phi}(t)|^2 dt $$

and from the inequality

$$ \begin{aligned} 2\langle A_s \varphi, \varphi \rangle_{L^2(\Sigma_T)} &= \sum_{k=0}^{N-1} \int_0^L (\mathcal{R}^k S_+^k A_s \varphi)(t) \cdot (\mathcal{R}^k S_+^k \varphi)(t) dt = \sum_{k=0}^{N-1} \int_0^L \Psi_k(t) \cdot \Phi_k(t) dt \\ &= \int_0^L \mathcal{M} \tilde{\Phi}_0(t) \cdot \tilde{\Phi}_0(t) dt + \int_0^L \mathcal{M} \tilde{\Phi}_L(t) \cdot \tilde{\Phi}_L(t) dt \\ &\geq \omega_1^2 \int_0^L |\tilde{\Phi}(t)|^2 dt = \omega_1^2 \|\varphi\|_{L^2(\Sigma_T)}^2 \end{aligned} \quad \square $$

REFERENCES

1. Antes H. A boundary element procedure for transient wave propagations in two-dimensional isotropic media. *Finite Elements in Analysis and Design* 1985; **1**:313–322.

2. Bachelot A, Bounhoure L, Pujols A. Time dependent integral method for Maxwell's system. *Mathematical and Numerical Aspects of Wave Propagation*. SIAM: Philadelphia, 1995; 151–159.

3. Bachelot A, Bounhoure L, Pujols A. Couplage éléments finis-potentiels retardés pour la diffraction électromagnétique par un obstacle hétérogène. *Numerische Mathematik* 2001; **89**:257–306.

4. Chudinovich I. Boundary equations in dynamic problems of the theory of elasticity. *Acta Applicandae Mathematicae* 2001; **65**:169–183.

5. Frangi A. Elastodynamics by BEM: a new direct formulation. *International Journal for Numerical Methods in Engineering* 1999; **45**:721–740.

6. Gaul L, Schanz M. A comparative study of three boundary element approaches to calculate the transient response of viscoelastic solids with unbounded domains. *Computer Methods in Applied Mechanics and Engineering* 1999; **179**(1–2):111–123.

7. Maier G, Diligenti M, Carini A. A variational approach to boundary element elastodynamic analysis and extension to multidomain problem. *Computer Methods in Applied Mechanics and Engineering* 1991; **92**:193–213.

8. Lubich C. Convolution quadrature and discretized operational calculus. I. *Numerische Mathematik* 1988; **52**: 129–145.

9. Lubich C. Convolution quadrature and discretized operational calculus. II. *Numerische Mathematik* 1988; **52**: 413–425.

10. Lubich C. On the multistep time discretization of linear initial-boundary value problems and their boundary integral equations. *Numerische Mathematik* 1994; **67**:365–389.

11. Bamberger A, Ha Duong T. Formulation variationnelle espace-temps pour le calcul par potentiel retardé de la diffraction d'une onde acoustique. I. *Mathematical Methods in the Applied Sciences* 1986; **8**:405–435.

12. Bamberger A, Ha Duong T. Formulation variationnelle pour le calcul de la diffraction d'une onde acoustique par une surface rigide. *Mathematical Methods in the Applied Sciences* 1986; **8**:598–608.

13. Ha Duong T. On the transient acoustic scattering by a flat object. *Japan Journal of Applied Mathematics* 1990; **7**:489–513.

14. Ha Duong T. On retarded potential boundary integral equations and their discretization. In *Topics in Computational Wave Propagation. Direct and Inverse Problems*, Davies P, Duncan D, Martin P, Rynne B (eds). Springer: Berlin, 2003; 301–336.

15. Becache E. A variational boundary integral equation method for an elastodynamic antiplane crack. *International Journal for Numerical Methods in Engineering* 1993; **36**:969–984.

16. Becache E, Ha Duong T. A space-time variational formulation for the boundary integral equation in a 2D elastic crack problem. *Mathematical Modelling and Numerical Analysis* 1994; **28**:141–176.
17. Ha Duong T, Ludwig B, Terrasse I. A Galerkin BEM for transient acoustic scattering by an absorbing obstacle. *International Journal for Numerical Methods in Engineering* 2003; **57**:1845–1882.

18. Costabel M. Time-dependent problems with the boundary integral method. In *Encyclopedia of Computational Mechanics*, Chapter 25, Stein E, de Borst R, Hughes TJR (eds). Wiley: New York, 2004.

19. Costabel M. Developments in boundary element methods for time-dependent problems. In *Problems and Methods in Mathematical Physics*, Jentsch L, Tröltzsch F (eds), vol. 134. B. G. Tuebner: Leipzig, 1994; 17–32.

20. Aimi A, Diligenti M, Guardasoni C, Panizzi S. A space–time energetic formulation for wave propagation analysis by BEMs. *Rivista di Matematica della Università di Parma* 2008; **8**(7):171–207.

21. Lasieka I, Lions JL, Triggiani R. Non homogeneous boundary value problems for second order hyperbolic operators. *Journal de Mathématiques Pures et Appliquées* 1986; **65**:149–192.

22. Duffy DG. *Green's Functions with Applications*. Chapman & Hall, CRC: London, Boca Raton, 2001.

23. Davies B. *Integral Transforms and their Applications*. Springer: New York, 2002.

24. Slepian D. Some comments on Fourier analysis. *Uncertainty and Modelling SIAM Review* 1983; **25**(3):379–393.

25. Fuchs WHJ. On the eigenvalues of an integral equation arising in the theory of band-limited signals. *Journal of Mathematical Analysis and Applications* 1964; **9**:317–330.

26. Lebeau G, Schatzman M. A wave problem in a half-space with a unilateral constraint at the boundary. *Journal of Differential Equations* 1984; **53**:309–361.

27. Miyatake S. Neumann operator for wave equation in a half space and microlocal orders of singularities along the boundary. *Seminaire sur les Equations aux Derivees Partielles*, vol. XVI, 8. Ecole Polytech: Palaiseau, 1993; 1992–1993.

28. Aimi A, Diligenti M, Monegato G. New numerical integration schemes for applications of Galerkin BEM to 2D problems. *International Journal for Numerical Methods in Engineering* 1997; **40**:1977–1999.

29. Monegato G, Scuderi L. Numerical integration of functions with boundary singularities. *Journal of Computational and Applied Mathematics* 1999; **112**:201–214.

30. Aimi A, Diligenti M. A new space–time energetic formulation for wave propagation analysis in layered media by BEMs. *International Journal for Numerical Methods in Engineering* 2008; **75**:1102–1132.

31. Sanchez-Sesma FJ, Iturraran-Viveros U. Scattering and diffraction of SH waves by a finite crack: an analytical solution. *Geophysical Journal International* 2001; **145**:749–758.

32. Iturraran-Viveros U, Vai R, Sanchez-Sesma FJ. Scattering of elastic waves by a 2-D crack using the indirect boundary element method (IBEM). *Geophysical Journal International* 2005; **162**:927–934.

33. Ciarlet PG. *Introduction à l'analyse Numérique Matricielle et à l'optimisation*. Masson: Paris, 1985.