On the space-time discretization of variational retarded potential boundary integral equations

This paper discusses the practical development of space-time boundary element methods for the wave equation in three spatial dimensions. The employed trial spaces stem from simplex meshes of the lateral boundary of the space-time cylinder. This approach conforms genuinely to the distinguished structure of the solution operators of the wave equation, so-called retarded potentials. Since the numerical evaluation of the arising integrals is intricate, the bulk of this work is constituted by ideas about quadrature techniques for retarded layer potentials and associated energetic bilinear forms. Finally, we glimpse at algorithmic aspects regarding the efficient implementation of retarded potentials in the space-time setting. The proposed methods are verified by means of numerical experiments, which illustrate their capacity.


Introduction
The philosophy of space-time methods is to consider space and time as components of space-time rather than disconnected entities. Space-time finite elements are based on meshes of the d + 1-dimensional space-time domain, where d ∈ N denotes the number of spatial dimensions. Especially over the course of the last decade, space-time finite element methods have achieved remarkable progress [1,2,3,4]. Advantages of this methodology are the natural treatment of non-stationary domains [5,6], adaptivity [7,8] and efficient parallelization techniques [9,10]. In the context of hyperbolic problems, space-time approaches facilitate locally explicit solution strategies exploiting causality and finite speed propagation [11,12,13].
While boundary integral equations (BIEs) have proven to be a compelling device for exterior scattering problems or transparent boundary conditions [14], the development of genuine space-time boundary element methods (BEMs) is in its infancy. Typical time domain BEMs are based on semi-discretization. In particular, they employ trial functions which are the product of separate functions in space and time [15,16,17,18,19]. An earlier attempt at relinquishing this product structure is due to Frangi [20], who exploits "causal" shape functions to discretize BIEs of the wave equation for d = 2. These functions can be interpreted as a predecessor to trial functions defined on unstructured space-time meshes. By giving up the usual product structure, however, one is confronted with more complicated integrals. The evaluation of these integrals is a major obstacle, stalling the practical development of space-time BEMs.
In the context of BIEs of parabolic problems, Tausch and collaborators [21,22] are actively developing quadrature techniques for these integrals. In [23], we proposed a tentative space-time BEM for the wave equation for d = 3. Integral formulations of the wave equation, especially for odd d ≥ 3, are of extraordinary structure, reverberating through their name retarded potential boundary integral equations (RPBIEs). The space-time methodology is particularly apt for treating the distinguished nature of RPBIEs. Therefore, this paper is intended to advance our earlier work.
The novelty of this paper lies in the utilization of space-time boundary elements to discretize variational formulations of RPBIEs. Although the mathematical analysis of Galerkin methods for RPBIEs is yet incomplete [24], they have already been applied successfully [25,26]. The integral operators acting on the surface density w : Σ → R (Σ is the space-time boundary) are of the form k(x, y)w(y)dS(y).
Here, x, y ∈ R 4 are points in space-time with spatial components x, y ∈ R 3 and k : R 3 × R 3 → R is the integral kernel. The set Q(x) is the intersection of Σ and a quadratic hypersurface, namely the backward light cone, which depends on x. Energetic bilinear forms with trial and test functions w, v read [27,28] k(x, y)w(y)∂ t v(y)dS(y)dS(x).
From here on, we refer to these integrals as inner (integral operator T k ) and outer (Galerkin testing). The perhaps most successful quadrature techniques for BIEs of elliptic problems treat both integrals together as one high-dimensional integral [29,30]. While this approach has compelling advantages, the design of such high-dimensional quadrature methods for hyperbolic problems is complicated due to the nonlinear behavior of x → Q(x). This is the reason why typical quadrature schemes employed in classical semidiscretizations of RPBIEs treat these integrals separately [31,32]. The present paper stays in line with these approaches in the sense that the inner and outer integral are treated individually. On the one hand, an alternative to the quadrature scheme for the inner integral we developed in [23] is proposed. On the other hand, a suitable formula for the outer integral and a tentative numerical integration method are discussed.
The paper is organized as follows. In Section 2, we exhibit the model initial-boundary value problem, two related RPBIEs, and their variational formulations. Section 3 discusses space-time boundary elements and quadrature techniques for RPBIEs. An algorithm which aims at the efficient implementation of retarded potentials is presented in Section 3.3. The purpose of Section 4 is to verify the proposed schemes via numerical experiments. Section 5 provides a brief conclusion of this work.

Retarded potential boundary integral equations
Let R n , n ∈ N be equipped with the usual Euclidean inner product ·, · and induced norm · . The unit sphere is denoted S n−1 := {x ∈ R n : x = 1} and we abbreviate S := S 2 . Consider a bounded open domain Ω − ⊂ R 3 whose exterior is denoted Ω + := R 3 \ Ω − . The Lipschitz boundary Γ := ∂ Ω − is equipped with the unit outward normal vector field ν Γ : Γ → S. Let d ∈ {+, −} and dist Γ : R 3 → R be the signed distance function of Γ defined by dist Γ : x → d inf y∈Γ x − y for x ∈ Ω d . Throughout this work, time coordinates are defined as geometrized time, i.e, the product of ordinary time and wave velocity, see [23,Section 2]. Let T > 0 be the simulation end time and Q d := (0, T ) × Ω d be the space-time cylinder with lateral boundary Σ := (0, T ) × Γ. To simplify notation, we introduce the fixed decomposition of points in space-time R 4 x := (t, x), R 4 y := (τ, y), with times t, τ ∈ R and spatial components x, y ∈ R 3 . Since Ω d is stationary, dist Σ : [0, T ] × R 3 → R is the time-invariant signed distance function of Σ given by dist Σ : x → dist Γ (x). Moreover, the space-time normal vector field ν Σ : Σ → S 3 has vanishing time component .
Let the (Lipschitz continuous) function φ Ξ : R 4 → R be defined by where the gradient ∇ is split into the time derivative ∂ t and the spatial gradient ∇ x . The three-dimensional hypersurface Ξ(x) := {y ∈ R 4 : φ Ξ (x − y) = 0} is the backward light cone with apex at x, see [23, Fig.  1].

Integral form of the wave equation
Let := ∂ 2 t − ∆ x be the d'Alembertian and u : Q d → R be subject to the homogeneous wave equation As a model problem consider Dirichlet boundary conditions with given datum g : where γ d 0 denotes the trace operator, see [33]. The normal derivative of u is denoted by γ d 1 u and it holds In this paper, we focus on boundary integral representations of the solution of (2)-(4). The involved integral operators employ the forward fundamental solution G of the d'Alembertian in three spatial dimensions [34,Operator 51] where δ 0 denotes the Dirac delta function. We tacitly exploit δ 0 (t) = δ 0 (−t) for any t ∈ R.
Definition 1. Let x ∈ Q − ∪ Q + . Define for sufficiently smooth w : Σ → R the retarded single layer potential and for sufficiently smooth v : Σ → R the retarded double layer potential with the kernel functions k i : In Section 2.2, the formulas of Definition 1 are revisited, abolishing the Dirac delta functions. It holds S w = 0 and D v = 0 in Q − ∪ Q + for any (admissible) w and v, respectively. Furthermore, Kirchhoff's formula u = d D γ d 0 u − S γ d 1 u represents solutions of (2) and (3) uniquely by their Cauchy data (γ d 0 u, γ d 1 u), see [35,Section 3.5]. Application of the trace induces the retarded single layer and double layer boundary integral operators where the latter formula with the factor 1/2 holds almost everywhere on Σ. Note that the integral representations of V and K are given by (5) and (6), respectively (for x ∈ Σ). In this work, we examine method unknown surface density BIE to be solved solution of (2)-(4) Table 1: Two approaches to solve (2)-(4) via BIEs. The indirect method is based on the ansatz u = S w, whose trace yields the BIE γ d 0 S w = g. In contrast, the direct approach employs Kirchhoff's formula u = d D g − S γ d 1 u , describing the wave field in terms of its Cauchy data. The yet unknown Neumann trace γ d 1 u is determined via the trace of Kirchhoff's formula.
two different approaches to solve (2)-(4) by means of BIEs. They are explained in Table 1 to provide a concise overview.
In order to enable space-time Galerkin discretizations, suitable variational formulations of the RPBIEs in Table 1 are required. Finding compelling space-time bilinear forms of RPBIEs has been the goal of multiple research efforts [36,27,24]. However, computable bilinear forms proven to be coercive in the same (Sobolev space) norm in which they are bounded are elusive, see [27, Theorems 3.1 and 3.3], [24,Corollary 4.6], and [15,Theorem 3]. Therefore, we resort to well-established bilinear forms which are of manageable complexity (nevertheless nontrivial) and supported by experience. The chosen setting is as in [15,Theorem 3] with a weight of ω I = 0, which corresponds to [27,Equations (29)- (31)]. Let H 0 , H 1 be appropriate Sobolev spaces and b V : The sole purpose of the abstract spaces H 0 and H 1 is to distinguish qualitative properties of the input densities of b V and b K in the discretization process of Section 3. For given Dirichlet datum g the functional b I (g, ·) : H 0 → R is defined by w → Σ g∂ t w dS, where the subscript I denotes the identity map. A variational formulation of the indirect approach in Table 1 is: For the direct method in Table 1 we use the formulation: As shown in [24, Propositions 3.4 and 3.7], a bilinear form similar to b V is positive definite iff T is sufficiently small, however, its induced norm is not equivalent to the H s (Σ)-norm for any s ∈ R [27, Theorem 3.1]. Still, promising numerical evidence is reported in [16,37].

Retarded layer potential integrals from the light cone's perspective
In this segment, we recast the integral operators of Definition 1 to a natural representation in the spacetime context. In [23], we employ local parametrizations of the space-time boundary Σ to derive a suitable formula for retarded potentials. In this work, we seek a representation in terms of the light cone Ξ instead.
Theorem 2 (Coarea formula). Let m, n ∈ N with m > n, f : R m → R n be Lipschitz continuous, J f (x) be its n-dimensional Jacobian at x ∈ R m , and g : R m → R be integrable. It holds A proof of the coarea formula can be found in [38,Theorem 3.2.12]. Assuming that f and g are as in Theorem 2 with the addition that the function R m → R, This, in combination with the sifting property of δ 0 , leads to if the integral on the right hand side exists. A formula similar to (10) for n = 1 can be found in [39,Theorem 6.1.5]. Let f := dist Σ and g : R 4 → R be such that γ 0 g := γ − 0 g = γ + 0 g holds. In this case, (10) yields where we used γ 0 ∇ dist Σ = ν Σ . We turn our attention to the operators in Definition 1 and introduce an operator that unifies the integral formulas of S, D, V, and K. Let x ∈ R 4 and k : R 3 × R 3 → R be a kernel function as in (7). For sufficiently smooth f : R 4 → R with bounded support we define the retarded Newtonian potential N k by Applying (10) and (1) to (12) yields k(x, y) f (y)dS(y).
For sufficiently smooth w : Σ → R we define analogously the retarded layer potential T k by which models the operators in Definition 1 via S = T k 1 and D = T k 3 + T k 2 ∂ t . For given w : Σ → R consider an extension w : R 4 → R such that w = γ − 0 w = γ + 0 w holds. Insertion of f := wδ 0 • dist Σ in (12) in conjunction with (11) yields the identity N k ( wδ 0 • dist Σ ) = T k w, where T k w is as in (14). Application of (13) leads to the desired representation By recasting the potentials of Definition 1 to the form (15) we have yet traded the Dirac delta on Ξ(x) for a Dirac delta on Σ. We incorporate a parametrization of Ξ(x) to obtain a computationally sensible formula. In the following, SO(3) denotes the special orthogonal matrix group in three spatial dimensions. For given x ∈ R 4 , r 0 > 0, and R ∈ SO(3) define ψ x : P → R 4 by where e S : [0, 2π) × [0, π] → S is defined by e S : (ϕ, θ ) → cos ϕ sin θ sin ϕ sin θ cos θ .

Space-time discretization and numerical evaluation of retarded layer potentials
As already indicated, the novelty of the proposed method lies in the utilization of space-time boundary element spaces as in [23]. The space-time boundary Σ is represented by , a mesh composed of N ∈ N open nonoverlapping tetrahedrons σ . We refer to the subsets σ ⊂ Σ as panels (not elements), see [30, Section 1.2] and [41, Section 2.3]. The mesh size is denoted h := max σ ∈Σ N diam σ . Simplex space-time meshes are constructed via the algorithm outlined in [42] and we resort to lowest order trial spaces.
Definition 4. Let P n (σ ) be the space of polynomials of order up to n ∈ N 0 := N ∪ {0} in the tetrahedron σ . Define the (discontinuous) space of indicator functions and the space of (continuous) hat functions by is used to approximate functions in H 1 . Consequently, the discretized version of (8) reads: In (9) the integral operator K acts on the given Dirichlet data g. In such cases it is common practice in BEMs for elliptic problems to approximate the data, see, e.g., [43,Chapter 12]. To this end, we employ the L 2 (Σ)-orthogonal projection Q 1 h : Assuming g ∈ L 2 (Σ) holds, the discretization of (9) with w h ≈ γ d 1 u reads: The following sections are concerned with the numerical evaluation of the involved operators.

A quadrature method for the "inner integral"
In this section, we devise a numerical integration scheme for (17) tailored to tetrahedral panels.
Definition 5 (Inner integral (17)). Let x ∈ R 4 be arbitrary but fixed and σ be a tetrahedron embedded in R 4 with normal vector ν ∈ S 3 . The unit outward conormal vectors of the four triangular faces of σ are denoted ν i ∈ S 3 , i = 1, . . . , 4 and satisfy ν, ν i = 0. Let k : R 3 × R 3 → R be as in (7), w : σ → R be analytic, and ψ x : P → Ξ(x) be as in Definition 3. Define the integral kernel k ψ : P → R and the integral I by The notation introduced in Definition 5 is employed throughout the remainder of this section. We denote the tangent hyperplane of the panel The condition ψ x (ζ ) ∈ T σ with ψ x as in Definition 3 is equivalent to where the spatial and time components of ν are denoted by ν x ∈ R 3 and ν t = 0, respectively.
Definition 6. For given normal vector ν x ∈ S let R ∈ SO(3) be such that R ν x = 0 0 1 holds. Furthermore, for given apex x ∈ R 4 and panel σ ∈ Σ N let r 0 := sup y∈σ (t −τ) be the largest time separation between x and σ .
Note that sup y∈σ (t −τ) ≤ 0 implies Ξ(x)∩σ ∈ {∅, {x}} and, therefore, I = 0. As a consequence, the potential conflict between the assumption r 0 > 0 in Definition 3 and r 0 as in Definition 6 is of no practical significance. An explicit formula for R in Definition 6 is provided in [40,Remark 3.12]. Inserting e S from Definition 3 and R from Definition 6 into (21) leads to with ρ 0 := x − x σ , ν /r 0 . The partial derivatives of the level set function in (22) We restrict these derivatives to the solution of (22), namely ρ 0 = ρ cos θ , and obtain the maps ρ → −ρ 0 /ρ as well as ρ → ρ 2 − ρ 2 0 . While the magnitude of the first derivative is monotonically decreasing, the latter is increasing. This shows that for sufficiently small ρ the ρ-direction is suitable for parametrizing the solution of (22), while for large ρ the θ -direction becomes the better choice, see Fig. 1a. The point where the partial derivatives are of equal magnitude is given by Define the two domains and the parametrizations i : D i → P, i = 1, 2 by 1 : This facilitates a parametrization of the integral in Definition 5 via 1 and 2 which involves integrals in the implicitly defined subsets Definition 6 enables the use of a finite patch D 2 in (24): the following choice is sufficient to capture the entire panel σ  Fig. 1a curves beneath the line θ = π/2 are solutions of (22) for ρ 0 > 0, while those above correspond to ρ 0 < 0. The black dots indicate the points (ρ eq , θ eq ), where the tangent to the curve has unit slope. Figure 1b illustrates the piecewise parametrization of ψ −1 x (Ξ(x) ∩ T σ ) in terms of 1 and 2 for ρ 0 = 1/5. Depiction of the ϕ-component is omitted because the curve is merely extruded into the third direction.
Lemma 7. Let k ψ be as in Definition 5 with kernel function k i , i = 1, 2, 3 as in (7) and R be as in Definition 6. The integral kernel k ψ • j , j = 1, 2 is smooth in D j .
There exist several procedures for evaluating integrals like (24) accurately, see, e.g., [44,45,46]. The algorithm employed in this paper is a combination of quadtree subdivision and exact parametrizations of the zero level set. In a nutshell, it attempts to identify the shape of the subset of B i that lies in a quadtree cell among a few predefined scenarios. For these admissible cases, exact parametrizations of the relevant subset are constructed and the transformed integrals are approximated accurately by standard tensor-Gauss quadrature rules. If this case identification fails, the algorithm resorts to subdivision. Our approach is based on the method proposed in [47], however, in contrast to the cited source, the zero level set is parametrized by means of the ideal transformation discussed in [48]. A similar approach is elaborated in [49]. We denote the depth of the quadtree by r max ∈ N 0 and n G ∈ N is the number of Gaussian quadrature points per direction. For each admissible quadtree cell at most 2n 2 G quadrature points are employed. The reader is referred to [40, Section 3.8.1] for details regarding the implementation.

A quadrature method for the "outer integral"
In order to evaluate the bilinear forms in (18) and (19), integrals of the form Σ w∂ t v h dS have to be computed, where v h ∈ S 0 h (Σ N ) and w is either in S 1 h (Σ N ) or defined through the action of T k . To this end, we consider a fixed panel σ ⊂ Σ with Lipschitz boundary ∂ σ . The unit outward conormal vector field ν ∂ σ : ∂ σ → S 3 satisfies ν ∂ σ (x), ν Σ (x) = 0 for any x ∈ ∂ σ for which it exists. Let v ∈ C 1 (Σ) and define which jumps only across ∂ σ . Note that v σ represents a basis function of S 0 where ν ∂ σ ,t is the time component of the unit outward conormal vector. The application of the cited theorem is justified because time is a tangential coordinate on Σ. For w ∈ C(Σ) and where v| ∂ σ denotes the trace of the restriction v| σ to ∂ σ . Let the bilinear form if A w is continuous across ∂ σ . Note that the use of the lowest order test space causes integrals on σ to vanish in (28). Since ∂ σ is composed of four 2-simplices, we only require a quadrature technique for triangles. In classical time domain discretization schemes, certain singularities of functions induced by retarded layer potentials have been studied, leading to carefully developed quadrature schemes [50,51]. Such an analysis in the space-time context could unveil the regularity of the function x → T k w(x) for different kernels k and regularity classes of w. These investigations might drive the design of tailored quadrature schemes for (28). While such comprehensive surveys lie beyond the scope of this work, the occurrence of singularities in the function x → T k w(x) is hinted in A by virtue of an example. Due to the lack of smoothness, we suggest to apply composite midpoint rules in order to evaluate (28) for A = T k , see Fig. 2. Figure 2: Illustration of quadrature rules employed for (28); for m Q ∈ N the reference triangle is subdivided into m 2 Q congruent triangles. Due to the lack of smoothness of the integrand, the midpoint rule is applied in each triangle. The triangles and midpoints, indicated by black dots, are displayed for m Q = 1, 2, 3, 5 (from left to right).

An algorithm for computing the set of lit panels efficiently
As discussed in Section 2.2, retarded layer potentials evaluated at x ∈ R 4 integrate along Ξ(x) ∩ Σ. Given a mesh Σ N , the set of panels lit by the backward light cone is denoted Both Σ and Ξ(x) are three-dimensional hypersurfaces, implying that Ξ(x) ∩ Σ is two-dimensional, unless it degenerates. The dimensions of these sets suggest that although |Σ N | = N holds, we expect |Σ Ξ N (x)| = O(N 2/3 ) as N → ∞ for a sequence of quasiuniform meshes. Given a density function w : Σ → R, our goal is to implement the evaluation of the linear retarded layer potential of Section 2.2 where w σ is defined as in (26). In an approach we are inclined to label "naive", Σ Ξ N (x) is constructed by considering each panel σ ∈ Σ N individually and verifying if Ξ(x)∩σ is nonempty. Clearly, this procedure involves O(N) operations (the maximum amount of operations necessary to verify Ξ(x) ∩ σ = ∅ is independent of N), spoiling the O(N 2/3 ) behavior dictated by the cardinality of Σ Ξ N (x). In the subsequent paragraphs, we exhibit a straightforward algorithm for constructing Σ Ξ N (x) more efficiently. Assume we were given a set X N (x) ⊂ Σ N such that Σ Ξ N (x) ⊂ X N (x) and |X N (x)| ≤ C|Σ Ξ N (x)| held for some C ≥ 1 independent of N. We proceed naively on X N (x) by checking every σ ∈ X N (x) if Ξ(x) ∩ σ is nonempty and if so, it is a member of Σ Ξ N (x). By assumption |X N (x)| < C|Σ Ξ N (x)| = O(N 2/3 ) holds, hence the computational cost of this approach is dictated by the amount of operations necessary to set up X N (x). The algorithm for assembling X N (x) is based on a hierarchical organization of the panels.
any v ∈ leaves(T) has two disjoint successors whose union is v.
The vertices are called clusters and we identify V with T, i.e., we write v ∈ T instead of v ∈ V .
The construction of the cluster tree T is performed as discussed in [52, Section 1.4.1.1, Equation (1.21)], which involves O(N log(N)) operations for quasiuniform meshes [52,Theorem 1.27]. In essence, T depends on Σ N and n min only, therefore, it is set up once and used for every evaluation point. For each v ∈ T we find a bounding ball . Let x ∈ R 4 be given and B r (y) be the closed ball of radius r > 0 around y ∈ R 4 . It holds min z∈B r (y) where d r : [0, ∞) → r, r √ 2 is defined by Proof. We employ the decomposition z := (s, z) with s ∈ R and z ∈ R 3 . It holds z ∈ B r (y) iff z − y ≤ r 2 − (s − τ) 2 holds. We expand (29) and the triangle inequality yields The maximum is attained at for any z ∈ B r (y). The bound is sharp, because it holds Considering the lower bound, we apply the reverse triangle inequality to (29) where we used s − τ ≥ −|s − τ| ≥ − r 2 − z − y 2 . We abbreviate β := z − y and declare f r : The sharpness of this bound is confirmed by We turn our attention to the case x − y < r/ √ 2. Since f r (β ) < 0 holds iff β < x − y holds, the minimum of f r is located at β = x − y . Insertion of β = z − y = x − y in (30) yields confirms the sharpness of the stated bound for x − y < r/ √ 2 and the proof is complete.
Theorem 9 is applied in line 2 of Algorithm 1. There is no root of z → φ Ξ (x − z) for z ∈ B r (y), i.e., the bounding ball of v ∈ T is not lit by Ξ(x), iff either its minimum value is positive or its maximum value is negative. As an initialization, set L(x) := ∅ and call ApproximateLitLeaves(L(x), root T). Once the algorithm concludes, set X N (x) := v∈L(x) i∈v σ i . In line 1 of Algorithm 1, the routine GetBoundingSphere(v) returns a precomputed bounding sphere that encloses all σ i , i ∈ v. In our implementation, we use a slightly modified version of the algorithm laid out in [53], which computes a nonminimal bounding sphere. Although the cited source exhibits the algorithm explicitly for R 3 , its extension to R n , n ∈ N is obvious. Remark. The proposed algorithm is based on concepts typically encountered in fast BEMs. Nevertheless, this approach does not constitute a traditional "fast method"; it implements evaluation procedures of exact (apart from quadrature) retarded potential integral operators efficiently. The necessity for such implementational tricks, even outside the realm of fast methods, arises because retarded potentials are not classically global operators, but their integrals are supported on (subsets of) the hypersurface Ξ(x). (L(x), v). Given the current iterate of L(x) ⊂ leaves T and v ∈ T.

Experiment 1: computation of lit panels
The first experiment investigates the performance of the method discussed in Section 3.3, which computes the set of lit panels Σ Ξ N (x) efficiently. Two computational domains are considered, namely the unit cube Ω − = − As a preliminary test, the cardinalities of the set of lit panels |Σ Ξ N (x)| and the proxy set |X N (x)| are investigated. The evaluation point x C is chosen and three upper bounds for the size of leaf-level clusters n min ∈ {1, 5, 50} are employed. Results of this study are displayed in Fig. 3. On the one hand, the conjectured O(N 2/3 ) behavior of |Σ Ξ N (x C )| can be observed. On the other hand, the results suggest the existence of a constant C > 1 such that |X N (x C )| < C|Σ Ξ N (x C )| holds, which is the key assumption in Section 3.3. Furthermore, |Σ Ξ N (x C )| < |X N (x C )| holds in all considered cases, even for n min = 1. A second example is considered, which aims at demonstrating the increase in performance achieved by the proposed technique. The naive approach (check every σ ∈ Σ N if σ ∈ Σ Ξ N (x) holds) is compared to the procedure outlined in Section 3.3: Both approaches construct the same set Σ Ξ N (x), however, the elapsed times differ. Let t N be the execution time of the naive approach and t n be the time required to perform above list of three steps. Again, the subscript n := n min ∈ {1, 50} is the maximum size of leaf-level clusters. All execution times (provided in ordinary time, seconds) reported in Fig. 4 are minimum values of five consecutive runs of the stated procedures. The displayed results suggest that t N behaves like O(N), while t n features an O(N 2/3 ) behavior. The three solid lines in Figs. 4a and 4b overlap, indicating that t N depends little on the actual position of x. However, t n depends heavily on x, at least in Fig. 4a. This shows that the proposed  algorithm can yield particularly large reductions of the execution time for points x such that |Σ Ξ N (x)| is small. The difference between t 50 and t 1 is noteworthy, suggesting that the extreme choice n min = 1 is advantageous if Σ Ξ N (x) is computed for sufficiently many evaluation points x. Finally, note that all exhibited elapsed times are obtained on an Intel ® Core™ i7-8700 desktop machine with a clock speed of 3.2 GHz. The absolute values of the execution times in Fig. 4 are of little significance because the implementation is both single-threaded and immature. Nevertheless, the presented results reveal the improvement of the asymptotic behavior due to the method proposed in Section 3.3.

Experiment 2: quadrature method
The goal of this section is to verify the quadrature scheme of Section 3.1. To this end, we revisit the experiments carried out in [23,Section 4.1] and compare the results. The experimental setup is recapped for the sake of completeness. As a computational domain we consider the unit cube Ω − = − 1 2 , 1 2 3 with T = 5. Given a function u : Q + → R subject to (2) and (3), we evaluate the function where J : Γ → [0, 1] is the solid angle, see [43,Equation (6.11)]. All integral operators in (31) are approximated by the quadrature method introduced in Section 3.1. This is the only relevant source of the error u − u because u = u would hold if all integral operators were evaluated exactly (note that the exact Cauchy data are used in (31)). The chosen solution u of (2) and (3) is a spherical wave function where y S ∈ Ω − is set to y S := −0.1 −0.2 −0.3 and µ : R → R is given by This choice of µ is smooth µ ∈ C ∞ (R) and causal, i.e., µ(t) = 0 holds for all t ≤ 0. In the first example, a mesh Σ N of N = 180 panels is considered. The evaluation point is set to The relative error measure is evaluated for d ∈ {0, 0.1, 1, 3}. The quadrature scheme discussed at the end of Section 3.1 has two main input parameters, namely the number of quadrature points per direction n G ∈ N and the depth of the quadtree r max ∈ N 0 . We consider r max ∈ {10, 20} and study the convergence with respect to n G . Results of this experiment are exhibited in Fig. 5a for r max = 10 and in Fig. 5b for r max = 20. Clearly, e d decays rapidly as n G is increased. However, convergence ceases once the error falls below a certain threshold, which depends on r max . The existence of such a threshold suggests that certain quadtree cells fit no admissible scenario even after r max steps of subdivision. These cells are treated by low-order approximations and, therefore, convergence with respect to n G is capped. Nevertheless, for r max = 20 the achievable error e d ≈ 1 × 10 −12 is already rather close to machine epsilon. It is noteworthy that the case d = 0, which involves weakly singular kernel functions, is handled just as well as the cases with d > 0. This behavior is due to the employed transformations, which regularize the integrand, see Lemma 7. Finally, it is emphasized that the results for r max = 20 are quite comparable to the data provided in our earlier work [23, Figure 3(a)].
We consider a further test in order to support the capacity of the quadrature scheme for weakly singular integral kernels. As in [23, Section 4.1] a different mesh of the computational domain Σ N , consisting of N = 288 panels, is employed. The examined relative error measure is given by is the set of centroids of the panels in Σ N . Results of this convergence study are displayed in Fig. 5c for r max ∈ {7, 14}. Again, we observe that e Σ decays swiftly as n G is increased, until it falls below a certain magnitude, which depends on r max . For r max = 7 convergence ceases at e Σ ≈ 1 × 10 −5 , while r max = 14 leads to an error threshold more than three orders of magnitude smaller. Overall, the results depicted in Fig. 5c are quite similar to the ones displayed in [23, Figure 3(b)]. We conclude that the quadrature approach laid out in Section 3.1 is, albeit immature, indeed capable of computing highly accurate pointwise evaluations of retarded layer potentials in the space-time setting. c) error e Σ averaged over 288 points Figure 5: Convergence study of the quadrature scheme discussed in Section 3.1; r max and n G denote the quadtree depth and the number of quadrature points per direction, respectively. The total number of quadrature points behaves like O(n 2 G ).

Experiment 3: space-time BEMs
The final experiment is intended to verify the space-time BEMs discussed in Section 3 and to illustrate their capacity. In all following tests, the parameters for the inner quadrature are set to (r max , n G ) := (7, 8), while the outer quadrature employs m Q := 3 points per direction, see Fig. 2.
The first test investigates the indirect BEM (18) and we employ the exact solutions of V w = g derived in [54] for spherical scatterers Γ = S. Denote by Y m n : S → C the spherical harmonic function of degree n ∈ N 0 and order m ∈ Z such that −n ≤ m ≤ n holds. Let g : Σ → R be defined by g : , where Re denotes the real part and g 0 : R → R reads In this case, the solution w : Σ → R of V w = g is given by w : Note that the minimum in e opt is attained by the L 2 (Σ)-orthogonal projection of w onto S 0 h (Σ N ). A convergence study is displayed in Fig. 6. Both e opt and e BEM exhibit first-order convergence with respect to h. Therefore, the BEM approximation seems to satisfy a quasioptimality principle in L 2 (Σ), or in other words, there seems to exist a T -dependent constant C(T ) > 1 such that e BEM ≤ C(T )e opt holds.
The final example investigates the performance of the direct BEM (19). The computational domain is set to Ω − = − 1 2 , 1 2 3 and the employed reference solution is given by (32), where µ ∈ C 2 (R) reads The discretized RPBIE (19) is solved for w h ≈ γ + 1 u and the relative error measures of (33) are computed. Additionally, the error in the wave field u − u h is studied, where u h := D Q 1 h g − S w h is given by the discretized Kirchhoff's formula. We consider 26 evaluation points x i := (T, x i ) for i = 1, . . . , 26, where each x i ∈ R 3 lies on the boundary of the cube − 3 5 , 3 5 3 . The following relative error measure is reported   This provides further evidence that the BEM solution seems to satisfy a quasioptimality principle. Figure 7b indicates that the examined mesh sizes still lie in the preasymptotic regime of e Q . On average, we observe that the pointwise error in the wave field converges quadratically with respect to h. If the theory of BIEs for elliptic problems indeed carried over to hyperbolic RPBIEs, we could indeed expect second-order convergence, see [43,Equation (12.21)].
In both tests the proposed space-time Galerkin BEMs for RPBIEs yield optimal convergence rates in the L 2 (Σ)-norm. The results are tremendously better than the ones obtained by the space-time collocation BEM we developed in [23]. Furthermore, the provided evidence confirms that the quadrature scheme for space-time bilinear forms discussed in Section 3.2 yields sufficiently accurate matrix entries such that the overall convergence of the BEM solution is not spoiled (even for few quadrature points m Q = 3).

Conclusion
This paper presents a discretization scheme for variational integral equations of the wave equation based on space-time boundary elements. We derive an integral formula for retarded layer potentials which fits the space-time setting exceptionally well. A carefully constructed parametrization of the light cone simplifies these integrals greatly for piecewise flat boundary meshes. This enables the application of existing quadrature schemes developed for implicitly defined domains. Since retarded layer potentials induce non-smooth functions, the evaluation of related space-time bilinear forms is accomplished via a tentative low-order approach. Nevertheless, numerical evidence suggests that the proposed methods provide sufficiently accurate evaluations of retarded layer potentials and associated Galerkin matrix entries. In all examined tests, the error of the Galerkin approximation converges quasioptimally in the L 2 (Σ)norm. Furthermore, the efficient computation of the set of panels lit by the light cone is addressed. The proposed algorithm is based on a hierarchical structure of the mesh and facilitates computations of the set of lit panels in (nigh) optimal complexity. Although several ideas and techniques described in this work are still in an early stage of development, numerical experiments indicate their potential.