Space-time stochastic Galerkin boundary elements for acoustic scattering problems

Acoustic emission or scattering problems naturally involve uncertainties about the sound sources or boundary conditions. This article initiates the study of time domain boundary elements for such stochastic boundary problems for the acoustic wave equation. We present a space-time stochastic Galerkin boundary element method which is applied to sound-hard, sound-soft and absorbing scatterers. Uncertainties in both the sources and the boundary conditions are considered using a polynomial chaos expansion. The numerical experiments illustrate the performance and convergence of the proposed method in model problems and present an application to a problem from traffic noise.


Introduction
This article introduces a boundary element method for the quantification of uncertainties in timedependent sound emission and scattering problems.In these problems the nature of the sound source and the reflection or absorbing properties of the scatterer is often difficult to measure, for example, in traffic noise to assess the sound emission of vehicles [9].
The quantification of uncertainties is by now a key consideration in the assessment of numerical models in engineering [12,21,43,44].We refer to [7,8,10,11,15,20] for relevant works on partial differential equations involving random variables, in time-dependent or time-independent situations.The survey [31] discusses modern stochastic collocation and Galerkin approaches, especially for timeindependent problems.Specifically for the wave equation in a bounded domain, reference [40] considers a stochastic collocation method for problems with an uncertain wave velocity, but deterministic initial and boundary conditions.Extensions to elastodynamics are considered in [41].Also geometric uncertainties of the domain have recently been considered in the frequency domain [34,35].
In this article we propose and analyze a stochastic space-time Galerkin boundary element method for the acoustic wave equation in unbounded, exterior domains with uncertain boundary conditions, based on a polynomial chaos expansion [49,50,39] of the random variables.Our approach is based on recent advances for time domain boundary elements, which reduce the exterior wave equation to an integral equation on the boundary.Such methods are particularly relevant for emission or scattering problems, which involve wave propagation over large distances.Both space-time Galerkin and convolution quadrature methods have been developed, see [13,32,45] for an overview.Recent developments in the directions of the current work include stable formulations [6,38,47], efficient discretizations [4,16,23,26,27,42], compression of the dense matrices [1,46] as well as complex coupled and interface problems [2,3,5,17,24,29,30,37].
To be specific, the method presented here reduces stochastic boundary value problems for the exterior acoustic wave equation to integral equations on the boundary of the computational domain.Using a polynomial chaos expansion of the random variables, a high-dimensional integral equation is obtained which is then discretized by a Galerkin method in space and time.The resulting highdimensional, but sparse, space-time system is solved in a marching-on-in-time time-stepping scheme [48].The approach is presented for stochastic Dirichlet, Neumann and acoustic boundary conditions, with uncertain sources and acoustic parameters.We develop the theoretical background to establish stability and convergence of the proposed discretization of the stochastic boundary integral equations.The numerical experiments in three space dimensions illustrate the performance and accuracy of this method for increasing stochastic, respectively space-time degrees of freedom.An application to a problem from traffic noise is presented.
Note that the size of the linear algebraic system for stochastic Galerkin finite elements in three space dimensions remains challenging, in spite of extensive research (see p. 566 in [31] for a review of relevant works).Our use of boundary element methods for the wave equation crucially leads to a sparse system of equations in two space dimensions.This article is organized as follows: In Chapter 2, we introduce stochastic Dirichlet, Neumann and acoustic boundary value problems for the wave equation.We introduce the relevant boundary integral operators for the wave equation and formulate the boundary value problems as integral equations on the boundary.Furthermore we discuss abstract stochastic space-time Galerkin methods for the stochastic Dirichlet, stochastic Neumann and stochastic acoustic problems.Their stability follows from the weak coercivity of the formulation.In Chapter 3, we use a polynomial chaos expansion and discretize the Galerkin formulation using orthonormal basis functions in the stochastic variables.The discretization in space and time uses a tensor-product ansatz, leading to a marching-on-in time scheme.Chapter 4 presents the numerical experiments for uncertain sources and acoustic parameters, both in model problems and for a standard problem in traffic noise.Details of the discretization and the theoretical framework of Sobolev spaces are presented in appendices.
Notation: In this article we denote the time derivative of a function f by ∂ t f or also by ḟ .We further write f ≲ g provided there exists a constant C such that f ≤ Cg.If the constant C is allowed to depend on a parameter σ, we write f ≲ σ g.

Stochastic boundary problems and integral formulations
We consider the exterior wave equation where ν is the outer unit normal vector to Γ, x ∈ Γ.The variable ξ ranges over the stochastic space Ξ, with a probability measure π(ξ).
We recast the boundary value problems as time dependent boundary integral equations in a suitable scale of function spaces.We first consider the wave equation ( 1) with stochastic Dirichlet (SD), respectively Neumann conditions (SN), before acoustic boundary conditions (SA) are considered.

Dirichlet and Neumann boundary conditions
For the Dirichlet problem (SD), we make an ansatz for the solution using the single layer potential in time domain, almost everywhere for x ∈ Γ and π-a.s. for ξ ∈ Ξ.Here, G is a fundamental solution of the wave equation (1), As u| Γ = Vϕ, we obtain an equivalent formulation of the wave equation ( 1) with stochastic Dirichlet boundary conditions (SD), u| Γ = f , as a time dependent integral equation Given a solution for the density ϕ in (5), the solution to the wave equation is obtained using the single layer ansatz (2).
An analogous formulation for the Neumann problem uses a double layer potential ansatz for the solution u: where ψ again is a causal function, ψ(τ, y, ξ) = 0 for τ < 0. The resulting integral formulation for the wave equation ( 1) with stochastic Neumann boundary conditions (SN), −∂ ν u = f , is the hypersingular equation Here, Given a solution for the density ψ in (7), the solution to the wave equation is obtained using the double layer ansatz (6).
The numerical experiments will quantify errors in space, time and stochastic variables using natural energy and L 2 norms.More generally, space-time anisotropic Sobolev spaces provide a convenient mathematical setting to study time dependent boundary integral equations like ( 5) and ( 7), which we review in Appendix A. The Sobolev space L 2 w (Ξ; is used to measure the error of the numerical solution in the stochastic variables in an L 2 -norm, the spatial L 2 -error of the s-th derivative, and the temporal L 2 -error of the r + s-th derivative.
Using these spaces, we define the bilinear form bilinear form associated to the Dirichlet, respectively Neumann problem by Here d σ t = e −2σt dt with σ > 0.
The weak formulation of the single-layer equation ( 5) then reads, given data Similarly, the weak formulation of the hypersingular equation ( 7) reads as follows, for data Below we will consider Galerkin discretizations of ( 11) and ( 12) using piecewise polynomial, conforming approximations in the subspaces The discretization is discussed in detail in Section 3. The resulting Galerkin boundary element method of the single-layer equation ( 11) is given by: Find ϕ h,∆t,ξ ∈ V such that for all φ h,∆t,ξ ∈ V B s D (ϕ h,∆t,ξ , φ h,∆t,ξ ) = ( ḟ , φ h,∆t,ξ ) .
The following results detail the theoretical numerical analysis of the Galerkin approximations (13), resp.(14), including the stability and convergence of the numerical methods.Readers interested in the implementation and performance may refer to Sections 3 and 4.
For the analysis, we first note a stochastic extension of the well-known mapping properties known for the layer potentials between Sobolev spaces related to the energy, see e.g.[32].A detailed exposition of the mathematical background of time domain integral equations and their discretizations is available in the monograph by Sayas [45], including methods based on convolution quadrature.See [13,32] for more concise introductions.
Theorem 2.1.The single layer and hypersingular operators are continuous for r ∈ R: Proof.Because the integral operators V and W do not depend on the stochastic variable, the proof follows from the mapping properties without stochastic variables, see [32,25].
and ∥ϕ∥ 2 and ∥ψ∥ 2 Proof.Again, the proof of Theorem 2.2 follows immediately from the results without stochastic variables, see e.g.[32] or in the specific notation used here, Proposition 3.1 of [27] for a) and Theorem 2.1 b) of [28] for b).
As a first consequence, we obtain the stability of the proposed methods: Then there exists a unique solution ϕ ∈ L 2 w (Ξ, H r σ (R + , H − 1 2 (Γ))) of (11) and a unique solution ϕ h,∆t,ξ ∈ V of (13), and Then there exists a unique solution ψ ∈ L 2 w (Ξ, H r+1 σ (R + , H 1 2 (Γ))) of (12) and a unique solution ψ h,∆t,ξ ∈ W of (14), and Proof.The proof of this theorem follows immediately from Theorem 6 of [23].In particular, measurability follows from the continuous dependence of the solution operator.
From the continuity and (weak) coercivity of the bilinear forms, it is now standard [32,25,28] to conclude the convergence of the Galerkin solutions in (13), (14).Details are omitted here.

Acoustic boundary conditions
For the acoustic boundary problem (SA), we assume that the coefficient α(t, x, ξ) in the acoustic boundary condition (SA) is bounded from above and below, almost everywhere for t ≥ 0, x ∈ Γ, π-a.s.ξ ∈ Ξ.
In addition to the single layer operator V and the hypersingular operator W, for the boundary integral formulation of this problem we require the double layer operator K and the adjoint double layer operator K ′ for x ∈ Γ, t > 0: For the derivation and the analysis of the relevant boundary integral formulation, following [32,33] we introduce the auxiliary interior problem with the boundary condition to the solution u i of the auxiliary problem in X , Then the following representation formula holds where φ = u i − u and p = ∂u i ∂ν − ∂u ∂ν denote the jumps of the Dirichlet, resp.Neumann traces on Γ. Letting x ∈ R n \ Γ tend to Γ, one obtains the following classical relations on Γ: We may substitute these into the acoustic boundary conditions (SA), (SA i ) to obtain a system of boundary integral equations on Γ: Specifically, for scattering problems with an incoming wave u inc one uses For the time interval [0, T ], where T is finite or T = ∞, we introduce the relevant bilinear form , given by The bilinear form now allows us to reduce the wave equation ( 1) with acoustic boundary conditions (SA) to a boundary integral equation on Γ.The weak form is given by Find for all (ψ, q) ∈ L 2 w (Ξ; Here, where In the deterministic case, this formulation has been considered in [32,33,27]. The Galerkin discretization of ( 22) using piecewise polynomial, conforming approximations in the subspace The discretization is discussed in detail in Section 3.
The following results detail the theoretical numerical analysis of the Galerkin approximations (24), including the stability and convergence of the numerical methods.Readers interested in the implementation and performance may again refer to Sections 3 and 4.
We first note the following mapping properties of K and K ′ : Theorem 2.4.The following operators are continuous for r ∈ R: Proof.Because the integral operators V and W do not depend on the stochastic variable, the proof follows from the mapping properties without stochastic variables, see [32,25].
The following theorem summarizes the continuity and coercivity properties of a s T .As above, the continuity follows from the mapping properties in Theorem 2.1 and Theorem 2.4. .
Proof.The upper estimate follows from the mapping properties of the integral operators in Theorem 2.1 and Theorem 2.4.The coercivity follows from Equation (64) of [32], 2 dx is the total energy at time T .While only timeindependent α were considered in [32], the proof of these formulas holds verbatim also when α depends on time.
We obtain the stability of the proposed method: Then the weak form (22) of the acoustic problem and its discretization (24) )), resp.(φ h,∆t,ξ , p h,∆t,ξ ) ∈ Z, which depend continuously on the boundary conditions.
As in Subsection 2.1, the continuity and (weak) coercivity of the bilinear form allow to conclude the convergence of the Galerkin solutions in (24).Details are omitted here.

Stochastic space time boundary element method
This section establishes the space-time and stochastic discretization of ( 11), (12) and (22).We first consider the stochastic discretization which uses the Stochastic Galerkin (SG) method.The resulting SG system is the solved in space-time by a Petrov-Galerkin method.

Stochastic Galerkin discretization
In order to derive the SG discretization we expand the solution of ( 11), (12), or (22) into a generalized Fourier series using a suitable orthonormal basis.To do so, we assume that the stochastic space be an orthonormal basis w.r.t. to the probability measure π m , i.e. for all i, j ∈ N 0 we have Furthermore, we define the corresponding polynomial approximation space for a fixed polynomial degree of J ∈ N 0 as follows, P J (Ξ) := {p : Ξ → R | p is a polynomial of degree J}.We let N N 0 := {κ = (κ m ) m∈N , κ m ∈ N 0 , for all m ∈ N}, and for J ∈ N 0 we define the truncated multi-index set By we denote the corresponding multivariate orthonormal polynomials and we define the tensor-polynomial approximation space of polynomials of degree J ∈ N 0 , Following [18,50], the solution ϕ of ( 11) and ( 12) can be written as Moreover, we truncate the infinite series in (28) at polynomial degree J ∈ N 0 and write, using the index set (26), The deterministic Fourier modes ϕ κ = ϕ κ (t, x) in ( 28) are defined by From the Fourier modes (30) we immediately extract expectation and variance of ϕ via Using the ansatz (29) and testing against ψ(t, x, ξ) κ := ψ(t, x)Ψ κ (ξ) yields, in combination with the the orthonormality relation (25), the following SG system of (11) for all κ ∈ J .Using the same ansatz and test functions in stochastics yield the following SG System for (22): For our computations in Section 4, we choose Ξ = [−1, 1] n+1 with Legendre polynomials of degree J as basis functions.

Space-time discretization of the Stochastic Galerkin system of the stochastic Dirichlet problem
For the discretization in a space-time setting, we use a tensor product ansatz, dividing the space from the time components.In space, we mesh the hypersurface Γ into triangular faces Γ i with diameter h i and h = max i h i .Here, we introduce the space of piecewise polynomial functions V p h with degree p ≥ 0. Further the corresponding basis {φ p j } is continuous fpr p ≥ 1.In time, we decompose R + into subinervals . ., and (∆t) = sup n (∆t) n .Here, we denote with {β q j } a corresponding basis of the space V q (∆t) of piecewise polynomial functions of degree q ≥ 0. The basis functions are continuous and vanish at t = 0 for q ≥ 1.Therefore, we divide the space-time cylinder with the basis functions {β q j φ p k }.Altogether the κ−th (space-time dependent) stochastic mode in ( 31) is given by ϕ κ h,∆t ∈ V p,q h,∆t for all κ ∈ J. Finally the fully discrete numerical approximaion of (11) becomes This leads to the following discrete variational formulation of ( 13): Find ϕ J h,∆t ∈ V such that for all l ∈ J and ψ l h,∆t ∈ V p,q h,∆t For our computations in Section 4, we set a uniform space mesh with diameter h and an equidistant timemesh with length ∆t.Furthermore we choose piecewise constant ansatz and test functions in space and time, i.e. ϕ κ h,∆t , ψ l h,∆t ∈ V 0,0 h,∆t .This leads to a lower triangular space-time Toeplitz matrix system.As usual, we solve this space-time system by backsubstitution, leading to a marching-onin-time (MOT) time stepping scheme [48].However, we have to note the stochastic contribution.Benefitting from the orthonormal polynomial basis (see (25)), we get a block diagonal matrix system in the stochastic space time framework (see (39)), where we apply the MOT scheme (J + 1) n times (see (40)).For more details about the computation, see the Appendix B.1.

Space-time discretization of the Stochastic Galerkin system of the stochastic acoustic problem
The stochastic space-time discretization spaces remain the same as in the subsections before.Further we compuate the stochastic acoustic problem on a uniform mesh with equidistant (∆t) and Ξ = [−1, 1] n+1 with Legendre polynomials of degree J.For simplicity, we detail the implementation for a time-independent absorption coefficient α, which we take to be of the form α(x, ξ) = α 1 (x)α 2 (ξ).In space-time, we choose the following ansatz and test functions: h,∆t piecewise linear in space and time • Ansatz function p κ h,∆t ∈ V 0,1 h,∆t piecewise constant in space and piecewise linear in time h,∆t piecewise linear in space and piecewise constant in time • Test function q l h,∆t ∈ V 0,0 h,∆t piecewise constant in space and time The Galerkin discretization (24) contains the retarded potential operators, namely the single layer potential V, the double layer potential K, the adjoint double layer potential K ′ and the hypersingular integral operator W. For the precise discretization of these retarded operators with this choice of ansatz and test functions in space and time, we refer to the Appendix in [30].This leads again to a space-time system for which we apply the MOT scheme.Unfortunately, since (24) depends on α as well as 1  α , we lose the stochastic block space-time diagonal system.At least we are able to separate the stochastic space from the time, leading to a lower triangle block stochastic space Toeplitz matrix, where we solve this system via an MOT scheme with stochastic space matrices.
Precisely, for every time step we solve: with ⊗ the Kronecker-Product and p o resp.φ o contains stochastic entries where each stochastic entry contains entries in space.Further and I const resp.I lin the identity with constant ansatz and test functions in space resp.linear ansatz and test functions in space.V k , K k , K ′ k , W k at timestep k + 1, k ∈ N 0 are space matrices for the corresponding single layer, double layer, adjoint double layer and hypersingular integral operators.For more detals on these retarded potential matrices, see [30,22] and for more details on how to construct (35) see the Appendix B.2.
4 Numerical experiments 4.1 Numerical example for the single layer potential with stochastic Dirichlet data In this numerical example we consider the stochastic Dirichlet Problem (11) on the unit sphere Γ = S 2 (see Figure 1 for a 2nd and 4th refinement level) with time [0, 2] and Ξ = [−1, 1] 3 .The stochastic right hand side with k = (0.2, 0.2, 0.2) is given via The spatial and temporal, uniform refinements are prescribed in Table   1 To quantify the error in the numerical approximation we consider the energy norm given by Upon using the orthonormality relation ( 25), the expression in ( 36) can be simplified to To obtain a corresponding benchmark energy we either consider a sufficient fine stochastic resolution or perform an extrapolation of the energy with respect to the space-time degrees of freedom (DOF).We are first interested in the convergence of the numerical approximation with respect to the number of orthonormal polynomials, i.e. an increasing SG polynomial degree, for a fixed space-time mesh.In order to compute a benchmark energy we compute the numerical approximation for the corresponding refinement level with a SG polynomial degree of 12.The mean and variance of the benchmark solution is depicted in Figure 2 for T = 2.
In Figure 3a) we plot the relative error between the approximated energy and the benchmark energy versus the SG polynomial degree in one stochastic dimension, where we fix the space-time mesh to be either the fourth or fifth refinement level, cf.Table 1.We can see clearly that the numerical error decays spectrally, which we expect from approximating a smooth solution with orthonormal polynomials.In Figure 3b) we also plot the convergence of the approximated energy with respect to increasing space-time DOFs (cf.Table 1).In order to compute the benchmark energy, we fix the SG polynomial degree to fourth and five and perform an extrapolation of the energy with regard to the space-time DOF.To estimate the rate of convergence we use a linear least squares fit using all data points.We compute a slope of roughly one-half, which corresponds to a rate of one half in terms of the space-time DOFs and in terms of the mesh width h to a rate of 3/2.
Here ∥x∥ is the euclidean norm of x, ⃗ n x the exterior normal vector of the face Γ i , where x lies and H the Heaviside function.Since we are using different ansatz and test functions we are not able to compute the energy of the solution, but rather quantify the error in the L 2 ((0, T ) × Γ, L 2 (Ξ))-norm, i.e. we consider the difference Similarly to the stochastic Dirichlet problem we compute the benchmark norm either by considering a sufficiently fine stochastic resolution or we perform an extrapolation of the norm with respect to the space-time DOFs.The corresponding space-time refinement levels are detailed in Table 2.
Similarly to the previous experiment for the stochastic Dirichlet problem we consider the convergence of the relative error with respect to the number of orthonormal polynomials, i.e the SG polynomial degree, for a fixed space-time mesh or with increasing refinement levels (see Table 2) for a fixed SG polynomial degree.In Figure 5 we present the mean and variance of p at time T = 1 for the highest refinement level of the cube, c.f. Table 2 for SG polynomial degree 4.
At the edges of the cube, where the mean has the lowest value, the variance seems to have the highest value.This indicates, that the stochastic coefficients influences the strength of the edge singularities.
Further Figure 6a) shows the convergence of (p, φ) each for the relative L 2 ((0, T ) × Γ; L 2 (Ξ))-error with increasing SG polynomial degree for a fixed space-time mesh.
With the increase of each SG polynomial degree, we almost see an error reduction of 10 −2 for φ until SG polynomial degree 4. Afterwards a clear reduction ins't visible anymore, whereas for p the convergence rate reduces with the increase of each SG polynomial degree until 6.Since ansatz and test functions aren't conform, this leads to a more unknown behaviour.
Figure 6b) depicts the convergence when increasing the space-time refinement level.Here we consider three different fixed SG polynomial degrees.Computing the rate of convergence in terms of DOFs reveals a rate of 2/3 for both components, the solution on Γ, φ and it's normal derivative on Γ, p.This corresponds to a rate of 2 in terms of the mesh width h.In our next numerical example, we let α ∼ U(1.5, 2.5) on the top of the cube Γ 1 := [−1, 1] 2 × {1} and and α ∼ U(3, 4) on Γ 2 := Γ \ Γ 1 .In Figure 7 the mean and variance of p are presented at time T = 1 for the highest refinement level and SG polynomial degree 4. The plot clearly differs from Figure 5, where the variance only has entries on the edges.Similar to Figure 5, we observe the highest values on the edges for the mean as well as for the variance.
In Figure 8, we plot an analoguous result to the Figure 8.The first plot shows a reduction of the error of 10 −3 until SG polynomial degree 2 for φ and 10 −2 until SG polynomial degree 3 for p.Then there is a very slow convergence until degree 6.For the second plot, we again observe a convergence rate of −2/3 in a fixed space-time DOF for φ.For p we can't tell a clear rate, but it tends to be −2/3, too, except for the 3rd to the 4th refinement, like in the corresponding Figure before.
In our final example we consider the stochastic coefficient of the impedance function α ∼ U(0.1, 0.9)   In Figure 9, we plot again the mean and the variance of p with the same settings as in the corresponding plots before.We observe that indeed edge singularitutes of the mean is independent of the variance.As in Figure 7 the interior of the boundary is zero on the whole boundary except the top and the edges.For the mean the change of α at the top causes that a high value switches into another edge.
In Figure 10a) correspondingly to the other Figures 7 and 5, we see an error reduction of 10 −2 until degree 2 for p and φ.Afterwards they decay spectrally.For a fixed SG polynomial degree in Figure 10b), we see for φ a rate of 1/4 and for p a rate of 1.
Altogether, we conclude that not only the convergence rates, also the edge and corner singularities depend highly on the impedance α despite choosing nonconforming ansatz and test functions.

Application to traffic noise
We finally indicate the feasibility and relevance of the developed techniques for problems in traffic noise.We consider the sound impact of a grown slick 205/55R16 passenger car tyre, of diameter around 0.6m at 2 bar pressure and subject to 3415N axle load at 50km/h on a flat street with an ISO 10844 surface.The domain is here given by the complement of the tyre in the upper half space, with the origin in the centre of the tyre, and its boundary Γ is shown in Figure 11.The sound pressure evolves according to a wave equation with wave speed c = 343m/s.Following [9], the sound emission corresponds to Neumann boundary condition (SN) with right hand side f = −ρ ∂ 2 un ∂t 2 , which has been obtained from simulations of the tyre displacement u n normal to the surface.For simplicity, we here consider f to be deterministic.On the street surface acoustic boundary conditions (SA) with right hand side f = 0 are imposed.The coefficient α depends strongly on the characteristics of the street surface [14], and it is therefore of interest to study the sound emission for a realistic range α ∼ U(1, 10) of its values.We extend α by 0 to the surface of the tyre.
As in [9,19], we use a single layer ansatz for the sound pressure, leading to the second-kind boundary integral equation (− The sound pressure is recovered from ϕ in a postprocessing step.Note that this formulation (38) is commonly used by engineers because, unlike (22), it only involves the two operators V and K ′ .However, the mathematical analysis even of deterministic second-kind boundary integral equations remains widely open for time-dependent problems.
The truncated boundary Γ is triangulated by the mesh in Figure 11 with 14092 elements.We set ∆t = 0.01/343s ≃ 2.915 × 10 −5 s and T = 1/343s ≃ 2.915 × 10 −3 s.We discretize the weak formulation of (38) on this mesh by piecewise constant ansatz and test functions in space and time, and stochastic polynomial degree 1.Details are provided in Appendix B. 3.
In Figure 12, we plot the mean and the variance of the density ϕ.As the source term f is nonzero only on the tyre, the density vanishes on the street outside a neighborhood of size cT of the tyre.After postprocessing the density ϕ using the single layer potential, Figure 13 depicts the sound pressure in a point (0.125m, 0m, −0.31m) near the street surface close to the tyre.
The results in the time domain also provide an efficient way to obtain information for a wide band of frequencies.The A-weighted sound pressure level gives an approximation to the human perception of noise [36].Figure 14 shows the A-weighted mean sound pressure level in dB of the acoustic wave emitted from the tyre, averaged over 17 points at a distance of 0.4m from the center, as well as one standard deviation.It is obtained by a discrete fast Fourier transform in time from the sound pressure.
Simulations as above of the sound impact of a passenger car tyre provide the computional basis for investigations in traffic noise [9].For realistic problems, both sources and the properties of the street surface are uncertain, and the proposed methods allow to quantify the mean sound pressure level and its variance as key ingredients for further analyses.

Conclusions
In this article we have introduced and studied a stochastic space-time Galerkin boundary element method for the acoustic wave equation.Stochastic Dirichlet, Neumann and acoustic boundary conditions are considered, with uncertain sources and acoustic parameters.The approach builds on recent advances for time domain boundary integral equations.
After a polynomial chaos expansion of the random variables, we obtained a high-dimensional spacetime boundary integral equation.This integral equation is discretized using a Galerkin method based on tensor products of piecewise polynomial basis functions in space, respectively in time, leading to a high-dimensional marching-on-in-time time-stepping scheme [48].After presenting the formulation and its theoretical background, we discussed the detailed discretization in the stochastic and spacetime variables.The numerical experiments for the wave equation in 3d empirically confirmed the performance and the convergence rates for increasing stochastic, respectively space-time degrees of freedom.They also confirmed the applicability of the proposed methods for applications to traffic noise.
They are Hilbert spaces, and we note that the basic case r = s = 0 is the weighted L 2 -space with scalar product ∞ 0 e −2σt Γ uvds x dt.Because Γ is Lipschitz, like in the case of standard Sobolev spaces these spaces are independent of the choice of α i and φ i when |s| ≤ 1.
Further we set σ = 0, then we get the system: Using the linearity of V, we can separate the integral into the stochastic and space-time part.Furthermore, since else for k = 1, . . ., k, considering the stochastic integrals we only need to look at the case i = j, where we nummerate i = i Now let us consider the remaining space time integral.
where χ E l−m is the indicator function with the lightcones we get the space time system Altogether we get the stochastic space time system: Now we are able to solve this system by applying forward substitution (marching on in time scheme) (J + 1) n times.For every l and i = 0, . . ., (J + 1) n − 1 solve

B.2 Discretization of the stochastic acoustic problem
In this subsection, we write more details about the computation of the stochastic acoustic problem, in particular (21).The coefficients α, α −1 depend on random variables.Therefore the acoustic boundary conditions the system of equations leads no longer into a block diagonal system.However, the implementation of the stochastic part is independent of the implementation of the integral operators V, K, K ′ , W and may be done separately.
For the retarded potential operators, we get similar results as before.However we have to consider the identity part with care.
We separate (41) by exploiting the linearity of the integral operators.The stochastic terms of the third and fourth terms can be handled as in the case for V the section before.Since the first term and the second term depend on α, which may depend on the stochastic and the space, we have to be careful.Consider the first and second term.Using the ansatz and test functions, we get for the time integral:   For the computation of the boundary integral operators we refer to the Appendix of [29] and [22].Altogether, we get a marching on in time scheme for the stochastic space system.Therefore for every time step l, solve the stochastic space system: , where V l−1 , K l−1 , K ′ l−1 , W l−1 are the corresponding matrices for the l-th time step.
Here again we use Legendre polynomials and piecewise constant ansatz functions in time resp. in space.
For the test function, we again use qh,∆t = γ m ∆t (t)η r h (x)Ψ j (ξ), where γ m ∆t (t) and η r h (x) are piecewise constant functions in time, resp. in space, with m = 1, . . ., N o , r = 1, . . ., N s and j = 0, . . ., (J +1) n −1.Therefore, we get:  For the computation of the single layer operator, we refer to the computation above in B.1 and for the adjoint boundary layer operator, we refer to the Appendix of [29] and to [22].As in the main text of the article, we obtain a marching-on-in-time scheme to solve the following system for the solution in each time step.If α only depends on the stochastic variable, the individual operators factorize into tensor products a polyhedral domain or screen, with X ⊂ R d , where we focus on the cases d = 2 and d = 3.On the boundary Γ = ∂X we either prescribe stochastic Dirichlet (SD), stochastic Neumann (SN) or stochastic acoustic boundary conditions (SA), involving a random right hand side f involving the Heaviside function H for d = 2 and the Dirac distribution δ for d = 3.The density ϕ is assumed to be causal: ϕ(τ, y, ξ) = 0 for τ < 0.Taking the trace of the ansatz (2) on Γ, we obtain the single layer operator V Vϕ(t, x, ξ) = R + Γ G(t − τ, x, y) ϕ(τ, y, ξ) dτ ds y .

Figure 4 :
Figure 4: 2nd and 4th spatial refinement level of the unit cube Section 4.2.

Figure 5 :
Figure 5: Plot of a) Mean and b) variance of the first component of the solution at time T = 1, refinement level 5 and SG polynomial degree 4 and α ∼ U(0.1, 1.9) on Γ.

Figure 7 : 3 Figure 8 :
Figure 7: Plot of a) Mean and b) variance of the first component of the solution at time T = 1 for refinement level 5 and SG polynomial degree 4 and α ∼ U(1.5, 2.5) on Γ 1 and α ∼ U(3, 4) on Γ 2 .

Figure 9 :
Figure 9: Plot of a) Mean and b) variance of the first component of the solution at time T = 1, refinement level 5 and SG polynomial degree 4 and α ∼ U(0.1, 0.9) on Γ 1 and α ∼ U(3, 4) on Γ 2 .

Figure 12 :
Figure 12: Plot of a) the mean and b) the variance of the density at time T = 2.915 × 10 −3 s.

Figure 14 :
Figure 14: A-weighted mean sound pressure and confidence intervall as function of frequency, averaged over 17 points.

B. 1
Discretization of the Dirichlet-Problem Let ξ = (ξ 1 , . . ., ξ n ), x = (x 1 , x 2 , x 3 ) with N o the amount of timesteps, N s the amount of elements and J the polynomial degree of Legendre polynomials, which are uniformly distributed at [−1, 1].Then we set the ansatz functions φ

Table 1 :
Space-time mesh hierarchy for the numerical experiment from Section 4.1.2nd and 4th refinement level depicted in Figure . The CFL (Courant-Friedrichs-Levy) rate is fixed at 0.605.