{'fraction_non_alphanumeric': 0.04535806055189364, 'fraction_numerical': 0.02182996313261088, 'mean_word_length': 3.683549602344077, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian identification of sound sources with the Helmholtz equation

In this work we discuss the problem of identifying sound sources from pressure measurements with a Bayesian approach. The acoustics are modelled by the Helmholtz equation and the goal is to get information about the number, strength and position of the sound sources, under the assumption that measurements of the acoustic pressure are noisy. We propose a problem specific prior distribution of the number, the amplitudes and positions of the sound sources and algorithms to compute an approximation of the associated posterior. We also discuss a finite element discretization of the Helmholtz equation for the practical computation and prove convergence rates of the resulting discretized posterior to the true posterior. The theoretical results are illustrated by numerical experiments, which indicate that the proven rates are sharp.


Introduction
Throughout the paper, unless noted otherwise, all functions and measures are complex valued. Let D ⊆ R d , d = 2, 3, be an bounded convex polygonal/polyhedral domain. Furthermore let the boundary be decomposed in the form ∂D =: Γ = Γ N ∪ Γ Z , with H d−1 (Γ Z ) > 0, where Γ N := J j=1 Γ j can be written as the union of some faces of Γ and Γ Z := Γ \Γ N . H d−1 denotes here the d − 1-dimensional Hausdorff measure. We model the acoustic pressure y u with angular frequency ζ and speed of sound c as the solution of The Helmholtz equation describes a stationary wave. It appears in the modelling of acoustics or electromagnetism (see for example (Dautray and Lions, 2012, Chapter 2, §8, Section 7) or (Kosljakov et al., 1964, Chapters VI, XXV) for a more physics based approach). To be more specific, g ∈ H 1 2 (Γ N ) is a complex function which models the amplitude of a primary sound source on a part of the boundary and τ (u) is a complex linear combination, possibly even a series, of Dirac measures with support in D. The Dirac measures model sound sources, which we want to identify. ρ describes the density of the fluid and γ ζ ∈ C is the wall impedance given by The frequency-dependent material constants α ζ > 0 and β ζ > 0 are related to the viscous and elastic response of the isolating material. However, in what follows, we will assume only that β ζ = 0. The boundary condition on Γ Z allows the modelling of an absorbing viscoelastic material covering the boundary walls. The boundary conditions on Γ N models external influence on the acoustic pressure.
Our aim is to deduce the distribution of the number, locations and amplitudes of the sound sources τ (u) from pressure measurements y := (y u (z j )) m j=1 + η ∈ C m at finitely many locations in the domain. The pressure measurements are taken at distinct points (z j ) m j=1 in D. We define the observation operator ( 3) where 1 κ is a suitable subset chosen from a sequence space of amplitudes α and positions x. Later in Section 3.2, we analyse under which assumptions and restrictions G is actually well-defined.
The measurement noise is denoted by η ∈ C m . We analyse the following inverse problem with a Bayesian approach, see for example (Stuart, 2010): In particular, we give precise meaning to η, y and u as random variables. We study (4) with a Bayesian approach which means that we are going to propose a problem specific prior, i.e. a distribution, which models prior knowledge of the distribution of the number, amplitudes and positions of the sound sources. We are then going to deduce the posterior distribution of those quantities, which incorporates our prior knowledge and observations from measurements. For the general principle we refer to (Stuart, 2010) or (Sullivan, 2015) for example. This work is structured as follows. In Section 2 we give a precise meaning to (1) and prove existence and regularity of solutions to (1). We proceed to approximate (1) by finite elements and prove a priori estimates for the error between the exact acoustic pressure and its finite element approximation. Moreover, we show that the corresponding discretization of observation operator G h is well-defined and meaningful.
In Section 3, we specify the prior and posterior distribution of the number, amplitudes and positions of the unknown sound sources.
Based on the theoretical stochastic framework we propose a sequential Monte Carlo method in Section 4 to solve (4). In particular, we propose an algorithm to approximate the posterior distribution of the number, the amplitudes and the positions of the sources.
Section 5 is reserved for numerical results which illustrate the theoretical framework.
Before we continue we would like to put this work into perspective. In (Bermúdez et al., 2004;Pieper et al., 2018), deterministic models are considered to recover number, amplitude and positions of the sound sources from measurements of the acoustic pressure. Both papers rely on techniques of optimal control theory. While the cost functional in (Bermúdez et al., 2004) is smooth, a non-smooth and sparsity promoting cost functional is considered in (Pieper et al., 2018). Morover, (Bermúdez et al., 2004) includes a discussion of finite element approximations.
Concerning Bayesian inverse problem as (4), prior distributions with similar structures as in this work can be found in computational statistics: The focus in (Green, 1995) is to simulate distributions with state spaces of differing dimensions. In particular, (Green, 1995) considers a measure π on measurable space (Q, K) where Q = k∈K {k} × Q k , Q k ⊆ R k , K ⊆ N. Examples where such distributions can be observed are Bayesian mixture modelling (Richardson and Green, 1997) or non-linear classification and regression (Denison et al., 2002). In our Bayesian inverse problem (4) we will work with a prior distribution which has its support in 1 κ . Elements in 1 κ with finitely many amplitudes α different from zero can be identified with elements of Q. However, a model with Q as in (Green, 1995) in our setting would only allow for a finite number of amplitudes α and positions x. This is not the case for 1 κ . Thus we are in a more general setting.
More application-oriented research on the location of sound sources with a Bayesian approach can be found in (Asano et al., 2013). Here, the Gaussian prior and the finite dimensional observation operator are based on a frequency decomposition of the signals. In particular, the authors use a Markov chain Monte Carlo method and compare the simulated results to real experiments. Due to the different structure of the observation operator in (Asano et al., 2013), the model does not consider the Helmholtz equation or its discretization. Further practical results are the subject of (Nakamura et al., 2009), where a similar problem with a dynamically moving robot is studied.

Preliminaries
For convenience of the reader we recap some basic results on Sobolev spaces and some auxiliary results. We also introduce the notation for the rest of the work. Most definitions are not given in full generality, but can be easily generalized. Moreover, we define several solution concepts to the Helmholtz equation in this section.
The absolute value of a complex number is denoted by | · |. We write · for the Euclidean norm in R d and denote by dist(x, y) the Euclidean distance of x and y. We also define the distances dist(x, y).
for all complex valued functions v ∈ C ∞ (D) with compact support in D. We write ∂ i f := h. Higher derivatives are defined analogously. For a multi-index α ∈ N d 0 we define We also write ∇f T = (∂ x1 f, . . . , ∂ xn f ), provided all weak derivatives of first order exist. The notion of a weak derivative and the following definitions can obviously be adapted to subdomains of D.
Definition 2 Let p ∈ [1, ∞] and k ∈ N 0 be given. We define the k-th Sobolev space as a generalization of C k (D) for k ∈ N: For k = 0 we write W k,p (D) := L p (D). Another special case is p = 2, then we write H k (D) := W k,2 (D).
We define the natural norm on W k,p (D) by For the next theorem we need the notion of a Lipschitz domain. A Lipschitz domain is a domain D in R d such that locally the boundary of D can be written as the graph of a Lipschitz map. A precise definition is given in (Wloka, 1987, Definition 2.4). There is said to have the N 0,1 property. We emphasize that any domain appearing in this paper, for example D or D \ B κ (x) for x ∈ D, κ ≥ 0, are Lipschitz domains.
Theorem 1 (Sobolev Embedding Theorem) For each Lipschitz domain E ⊂ R d we have the continuous embedding More precisely, there exists a C > 0 such that Proof The result can be found in many works including Sobolev spaces. See for example (Wloka, 1987, Theorem 6.2) and note that d ∈ {2, 3}.
In general, functions in H k (D) are not well-defined on the boundary, but their boundary values can still be defined in the sense of traces.
Theorem 2 (Trace Theorem) There exists a linear and continuous map γ tr : Proof Note that D is a Lipschitz domain and thus (Wloka, 1987, Theorem 8.7) applies.

Definition 3 We define
In the following, we introduce different solution concepts for the Helmholtz equation.
Definition 4 The function y ∈ C 2 (D)∩C(D) is called classical solution to the Helmholtz equation if it satisfies for all v ∈ H 1 (D).
Definition 6 (Very weak solution) The function y ∈ L 2 (D) is called very weak solution to the Helmholtz Both notions of weak and very weak solutions are derived from partial integration of smooth functions. In case of a very weak solution of the Helmholtz equation with u = ∞ =1 α δ x the first integral on the right has to be understood as The very weak solution to the Helmholtz equation has singularities at sources x ∈ D if u = ∞ =1 α δ x and α = 0. This makes it complicated to estimate the norm of the very weak solution on the entire domain D. Therefore, we restrict the possible source locations x to a source domain and the given measurement points z ∈ D to a measurement domain and demand a positive distance of the domains to each other.

Definition 7 (Source and measurement domain)
Let D ⊆ D be given. For κ > 0 the source domain D κ is defined as The measurement domain is defined as With restriction to these domains we can derive estimates of the norms · H 2 (Mκ) , · W 2,∞ (Mκ) or point evaluations of the form y u (z) for z ∈ M κ with sources x ∈ D κ . We further allow for an infinite amount of sources.
We denote this Banach space by 1 . We define the set 1 κ as the restriction of 1 to the source domain Moreover, we introduce the mapping τ which maps intensities and sources to a series of dirac measures The general assumption is now that D κ and 1 κ are nonempty. The definition of 1 κ requires that 0 ∈ D κ , otherwise 1 κ is empty. This assumption can be bypassed by an appropriate shift of the domain. Thus we assume w.l.o.g. that 0 ∈ D κ . We will show error estimates in the · 1 .
In this section we study the very weak solution y u of the Helmholtz equation. We proceed as in (Bermúdez et al., 2004, Section 3.2) to prove existence, uniqueness and norm estimates with respect to the data. Furthermore, we are interested in sensitivity properties of y u in u. Later, we want u to follow some prior probability distribution. This requires that all appearing constants in this section are independent of u. After the results for the infinite dimensional problem have been established, we continue to prove similar results for solutions y h u of a discretized problem and show convergence rates.
In the following, we will repeat the main steps from (Bermúdez et al., 2004, section 3.2) to show existence and regularity of the very weak solution. We prove both at the same time using a decomposition argument.

Proposition 1 The Helmholtz equation
has a unique very weak solution G x for any x ∈ D κ . It satisfies The constant C κ > 0 depends on κ but not on x.
Proof We introduce the fundamental solution of the Helmholtz equation in the whole space, i.e.
This equation has to be understood in the sense of (tem- We introduce the space S (R d ) as in (Zeidler, 1990(Zeidler, , pp. 1058(Zeidler, -1061. For functions v : R d → C we consider the seminorms With this definition the solution of (6) is given by Here, Y 0 denotes the zero-order second-kind Bessel function (Dautray and Lions, 2012, Chapter. 2, §8.6 Proposition 27). We also refer to (Bowman, 1958, Chapter 1) for some background on Bessel functions. In particular, there holds Φ x ∈ C ∞ (R \ {x}) and Φ x | D ∈ L 2 (D).
We introduce p x as the weak solution of According to (Bermúdez et al., 2004, Theorem 3.3) and the remarks before and after (Bermúdez et al., 2004, (3.10)) the solution p x exists, is unique and p x ∈ H 2 (D).

Moreover, it satisfies
This implies that G x = p x + Φ x | D is the very weak solution to (5). In particular for any κ > 0 we have Note that Bowman, 1958, Chapter 1). Moreover, for a constant C > 0 which is independent of κ and x. Because dist(x, Γ ) > κ, the right hand side in (8) can be bounded similarly using the same arguments. The proof for d = 2 is analogous. The estimate in the W 2,∞ (M κ )-norm follows from (Bermúdez et al., 2004, Lemma 3.5) with D 0 := M κ .

Remark 1
We have to 'cut out' discs around x ∈ D in Proposition 1 so that the expressions in (9) remain finite. In particular, the fundamental solution Φ x in R d in the proof of Proposition 1 has a singularity at x and Φ x ∈ H 2 (D).
Proposition 2 Let x ∈ D κ be given. Then for any y ∈ D κ with x − y < 1 2 κ there holds with a constant C κ depending on κ but not on x and y.
Proof Observe that with p x , p y , Φ x , Φ y from the proof of Proposition 1. Let z ∈ M κ . Note that for each ξ in the line [x, y] we have: By the proof of the previous we have for 1 ≤ k ∈ N and ξ ∈ [x, y] For for any 1 ≤ k ∈ N we have by the mean value theorem We conclude Using the definition of p x and equation (8) the expression p x − p y H 2 (D\Bκ(x)) can be bounded by .
Once again using the mean value theorem and bounding derivatives we get the bound which proves the claim.

Inhomogeneous Neumann Data and Multiple Sound Sources
In this section we present an existence and uniqueness result for the very weak solution of the Helmholtz equation (1) for g = 0. For a given sequence u = (α , x ) ∞ =1 ∈ 1 κ we allow the right hand side to be a series of Diracs given by τ (u) = ∞ =1 α δ x . Theorem 3 Let u ∈ 1 κ and g ∈ H 1 2 (Γ N ) be given. Then the problem is well-defined in the very weak sense. It has a unique very weak solution y u ∈ H 2 (M κ )∩L 2 (D) which satisfies with a constant C κ > 0 depending on κ but not on u.
Proof Let u = (α , x ) ∞ =1 ∈ 1 κ and g ∈ H 1 2 (Γ N ) be arbitrary. First we show well-definedness of the very weak formulation of (11). For a suitable testfunction v ∈ H 2 (D) this formulation is given by The only critical term is the one containing the series of Diracs. Since v ∈ H 2 (D) ⊂ C(D) by Theorem 1, this term can be bounded by This shows well-definedness of the very weak formulation. The next step is to show existence, uniqueness and the norm bound of the very weak solution. We use linearity of equation (11) to split up the solution in the form y u = y u,0 + y 0,g where y u,0 and y 0,g shall satisfy The unique very weak solution of the first problem is given by y u,0 = ∞ =1 α G x according to Proposition 1. We apply the bound in the same proposition and obtain a bound of the very weak solution of the form where the constant C κ > 0 does not depend on u.
We apply (Bermúdez et al., 2004, Theorem 3.3) to obtain existence and uniqueness of a weak solution y 0,g for the second equation of (12), which satisfies We combine the results for y u,0 and y 0,g to conclude the theorem.
Theorem 4 Let g ∈ H 1 2 (Γ N ) and u, v ∈ 1 κ be given such that u−v 1 < 1 2 κ. Let y u and y v be the very weak solutions of the Helmholtz equation (1) with the same Neumann data g but with different right hand sides τ (u) respectively τ (v). Then the solutions satisfy be given with u − v 1 < 1 2 κ. We apply decomposition (12) to show that the difference of the solutions satisfies This allows us to assume g = 0.
A straight forward computation shows We use the assumption u − v 1 < 1 2 κ which implies that for any ∈ N Applying these inequalities and Proposition 1 together with Proposition 2 for every ∈ N we further estimate which proves the claim.

Finite Element Spaces
We follow (Bermúdez et al., 2004) and discretize the Helmholtz equation (1) by piecewise linear finite elements. To this we end consider a family of triangula- The corners of the triangles are called nodes of D. We assume that no hanging nodes or hanging edges exist. More precisely, we suppose that for all h > 0 and any T, T ∈ T h the intersection T ∩ T is either empty, a single point, a common edge of T and T or a common facet of T and T . The last case is only relevant for d = 3. The domain D is polygonal/polyhedral which implies the existence of such a triangulation. For each triangulation T h we write h T for the diameter of a triangle T ∈ T h . The mesh size of T h is given by h := max T ∈T h h T . We denote by ρ T the diameter of the largest ball contained in T ∈ T h . We make the following assumption for the remaining part of the paper: This or similar conditions are called shape regularity and quasi-uniformity, as this for example prevents our triangles from becoming too acute or too flat. We note that these conditions imply that there exists a constant C > 0 such that We associate with each triangulation T h the finite element space V h which consists of real valued functions globally continuous in D and linear on each element T ∈ T h :

Galerkin Approximations
Definition 9 (Discrete solution) The function y h ∈ V h is called discrete solution, or Galerkin Approximaton, to the Helmholtz equation for all v h ∈ V h . We prove existence and uniquness of solutions to (14) in Theorem 5.
The following proposition allows us to work with L 2 (D) functions instead of diracs in the discrete setting. For such a construction also see (Scott, 1973/74, Theorem 1 and after Lemma 3).
Proof Let x ∈ D and define the functional j x as follows Because V h is finite dimensional it is a Banach space together with the norm · L 2 (D) . Thus the Riesz representation theorem implies the existence of a δ By construction and the equivalence of norms on V h we have Note that this C depends on h, but not on x.
For the Lipschitz continuity let x 1 , x 2 ∈ D. We have that δ x1,h , δ x2,h ∈ V h and both are Lipschitz continuous with constants ∇δ x1,h L ∞ (D) and ∇δ x2,h L ∞ (D) , see for example (Ambrosio et al., 2000, Proposition 2.13). This yields Again using the equivalency of norms we can use the bound from before: ∇δ x1,h L ∞ (D) ≤ C δ x1,h ≤ C, where C does not depend on x, but on h. The same estimate holds for δ x2,h concluding the proof.
with C independent of f and g.
Proof Let f ∈ L 2 (D), g ∈ H 1 2 (Γ N ) be given. Then has a unique weak solution y f by (Bermúdez et al., 2004, Theorem 3.3). It satisfies with a constant C independent of f or g. We introduce the sesquilinear form a satisfies Gårding's inequality by (Bermúdez et al., 2004, Lemma 3.2), which means that there exist C 1 , C 2 > 0 such that Now one can simply use the same ideas used in the real case in (Schatz, 1974).
A direct consequence of Proposition 4 is the existence and uniqueness of a discrete Green's function G Corollary 1 There exists a h 0 > 0 such that for all h ∈ (0, h 0 ] and all x ∈ D κ the discrete formulation Remark 2 The existence of G x h strongly depends on the mesh size h. The upper bound h 0 for the mesh size depends on the coefficients of Gårding's inequality used in the existence proof (above equation (7) in (Schatz, 1974)). In (Bermúdez et al., 2004, Lemma 3.2) those coefficients are stated explicitly. The bound h 0 is proportional to 1 + ζ 2 /c 2 −1 which guarantees existence of G x h . Hence, for a large ratios ζ c one has to work with fine meshes to compute reasonable solutions.
Furthermore, the solution is given by where y g,h is the discrete solution of the Helmholtz equation with zero forcing term and Neumann data g.
Proof By Lemma 1 and Proposition 4 it is easy to see that exists and is well defined as the series converges in Linearity of the integral now shows that is indeed a solution to (18).

Pointwise error estimates
The proofs in this section follow the proofs in (Bermúdez et al., 2004, Section 4). The authors explicitly state though that they to not analyse the dependence of the appearing constants with respect to x − z , where x ∈ D κ denotes the position of a Dirac and z ∈ M κ an arbitrary measurement point. First we restate a slight variant of (Schatz and Wahlbin, 1977, Corollary 5.1).
The constants C, C do not depend on h, f, f h , D 1 and D 2 .
where h 0 and C κ do not depend on x or z.
Proof Let z ∈ M κ and let x ∈ D κ . Lemma 1 shows that there exists a h 0 independent of x such that G x h exists. We want to apply Proposition 5 to the real and imaginary part of G Here C does not depend on the used balls or on the particular z, but only on κ. Therefore we choose D 1 = B 1 8 κ (z) and D 2 = B 1 4 κ (z) and apply Proposition 1 to obtain . Linearity of real part and imaginary part in the weak formulation yields, after short computation, that for This proves (20) according to Proposition 5. Again by Proposition 1 we obtain with C κ depending on κ but not on x or z. This implies is estimated as in the proof of (Bermúdez et al., 2004, Theorem 4.4), which is based on the proof of (Schatz and Wahlbin, 1977, Theorem 6.1) and (Bermúdez et al., 2004, Lemma 3.5). Tracking the constants in both proofs shows the independence of x and z.
Theorem 6 Let z ∈ M κ , u ∈ 1 κ and g ∈ H 1 2 (Γ N ). Then there exists a h 0 > 0 such that for all h ∈ (0, h 0 ] the discrete solution y u,h ∈ V h satisfies with a constant C κ > 0 independent of u, h and z. Proof We are allowed to evaluate y u and y u,h at z because y u ∈ H 2 (M κ ) ⊂ C(M κ ) and y u,h ∈ V h . For the error estimate observe We apply Lemma 6 to bound the term under the supremum uniformly in x and obtain the desired rate of | ln h|h 2 for the first term. (Bermúdez et al., 2004, Theorem 4.2) gives the estimate for the second term. It is stated in such a way that the minimal mesh size may depend on g, but checking its proof and the reference to (Schatz and Wahlbin, 1977, Theorem 5.1) in the proof of (Bermúdez et al., 2004, Theorem 4.2) shows that it is independent of g.
the discrete solution of the Helmholtz equation with right hand side f = δ x1,h − δ x2,h and Neumann data g = 0. We use Proposition 4 to get with C independent of x 1 and x 2 . Proposition 3 now shows the result.
Corollary 2 Let z ∈ M κ . Then there exists a h 0 > 0 such that for h ∈ (0, h 0 ] and g ∈ H 1 2 (Γ N ) the solution mapping y (·),h (z) : 1 κ → C is continuous. More precisely there holds Proof For x 1 , x 2 ∈ D κ we use the equivalence of norms on V h and Proposition 7 to show Lemma 1 shows that |G x1 h (z)| ≤ C with C independent of x 1 . Hence a similar computation as in the proof of Corollary 4 replacing the Greens function with its discrete counterpart and replacing · H 2 (Mκ) with | · (z)| shows the result.

The Inverse Problem
Up to this point all definitions and all the analysis was carried out in a deterministic setting. Now let us recall problem (4). Our aim is to analyse the following inverse problem with a Bayesian approach: The observed data y ⊆ C m is fixed. We also recall the definition of the observation operator, namely G is well defined by Corollary 3 and the embedding H 2 (M κ ) ⊂ C(M κ ). In the Bayesian setting we assume that u follows some prior distribution µ 0 , which we will introduce in Section 3.1. We then want to incorporate the knowledge from the measured data y ∈ C m to obtain the posterior distribution µ y . Let us assume that the measurement noise η is multivariate complex normal distributed such that η has a probability density proportional to where · Σ is a norm on C induced by a scalar product.
Here, Σ ∈ C 2m×2m is a Hermitian, positive definite complex matrix of the form Σ = Γ C CΓ . In particular, η is multivariate complex normal distributed with expectation value zero, positive definite variance matrix Γ ∈ C m×m and Hermitian relation matrix C ∈ C m×m . Moreover, z 2 Σ is defined as Accordingly, the random vector Re(η) Im(η) is multivariate normal distributed in R 2m with mean zero and covariance matrix .
We will shown in Section 3.3 that our particular noise assumption implies that the posterior distribution is given by for F ∈ F, where F is a Sigma algebra on 1 κ . The constant Λ is a normalization constant given by The term exp(− 1 2 y − G(u) 2 Σ ) can be viewed as penalization for deviating too far away from the observed data. Thus the posterior combines both prior knowledge of the solution and the measurements favouring those u ∈ 1 κ which closely predict the observation. Before we continue we recall some basic definitions and tools.
Definition 10 (Radon-Nikodym Deriviative) Let (X, F, ν) be a measure space. Let f : X → R ≥0 be a ν-measurable function. Define the measure µ on X by Then f is called Radon-Nikodym derivative of µ with respect to ν and is denoted by If two measures µ 1 and µ 2 have a Radon-Nikodym derivative with respect to a third measure ν, then we are able to define a metric between them, see (Stuart, 2010, Section 6.7).
Definition 11 (Hellinger Distance) Let µ 1 , µ 2 , ν be measures on X such that µ 1 , µ 2 have a Radon-Nikodym Derivative with respect to ν. Then the Hellinger distance between µ 1 and µ 2 is defined as The Hellinger distance tells us how well two measures agree. Using Gaussian noise we obtain exponential terms in the Radon-Nikodym derivatives. In order to estimate those, we require an simple lemma, which immediately follows from the mean value theorem.

The Prior
Recall that the Bayesian approach needs prior knowledge in form of a probability distribution which fits to problem (4). As we want to recover the number of sound sources, their amplitudes, and their positions in D κ we require a prior distribution which is suits these requirements. Let ( 1 κ , F) be the measurable set where F is the Borel σ-algebra associated with the open sets generated by the norm · 1 on the open set 1 κ . For n ∈ N 0 we define the sets 1 n,κ := {(α , x ) ∞ =1 ∈ 1 κ | α k = 0, x k = 0, ∀k > n} ⊆ 1 .
Using this notation we are able to construct a specific probability measure on ( 1 κ , F).
This implies that samples from µ 0 are µ 0 -almost surely in 1 n,κ for some n.
Proof The proof is straight forward.
The motivation to choose such a measure follows from property (25). This ensures that the forcing term for (1) consists µ 0 -almost surely of a finite number of Diracs, that is because u ∈ 1 n,κ for some n ∈ N 0 . A practical approach to construct such a measure is now to express u in terms of random variables n, α, x such that the Diracs on the right hand side of the Helmholtz equation are given by The function q is the probability mass function of the random variable n, and given a particular n we are able to define the probability measure µ 0 n for the intensities (α n 1 , ..., α n n ) ∈ C n and positions (x n 1 , ..., x n n ) ∈ D n κ . We now state a simple to prove proposition from which we can deduce whether u ∼ µ 0 has moments.
Proposition 8 Assume that µ 0 is constructed as in Theorem 7 and that for any n ∈ N the measure µ 0 n have p-th moment. Let id denote the identity on 1 κ . Then there holds and if the last series is finite, then µ 0 has p-th moment.

Properties of the observation operator and potential
Now we investigate in the observation operator G. It is important to see that the Helmholtz equation only influences the posterior through the observation operator G. Therefore, we derive suitable properties of G from the underlying PDE model.

Proposition 9
The observation operator G is µ 0 -measurable and satisfies .
Proof Because µ 0 is a measure on ( 1 κ , F) and F is the Borel F-Algebra we can apply Corollary 4 for every point z ∈ M κ to show that G : 1 κ → C m is continuous and in particular µ 0 -measurable. Again Corollary 3 applied to every observation point z ∈ M κ yields , which shows the bound and G ∈ L p ( 1 κ ; µ 0 ) if µ 0 has p-th moment.
Definition 12 Define the potential Ψ as follows: We assume that · Σ is derived from an inner product (·, ·) Σ on C m .
Because the norm · Σ is defined on a finite dimensional space, where all norms are equivalent, we exclusively use the · Σ -norm w.l.o.g.
Given the regularity of the observation operator G we are able to deduce several properties of the potential Ψ .
Proof Property (i) follows directly from the definition of the potential Ψ . To prove property (ii) apply Markov's inequality to the nonnegative random variable G 2 Σ ∈ L 1 ( 1 κ ; µ 0 ) (see Proposition 9) to get Thus we are able to define X as the complementary set of events and conclude µ 0 (X) ≥ 1 2 . Using the triangle and Young's inequality we obtain for u ∈ X, y ∈ C m which proves the claim using Proposition 9 to estimate G(u) Σ .

Posterior
In the following, we define the posterior by means of a Radon-Nikodym derivative: Definition 13 Let the observation y ∈ C m be given. The posterior density with respect to the prior measure is defined by where Λ(y) is a normalization constant given by Λ(y) = 1 κ exp(−Ψ (u, y))dµ 0 (u).
Next we prove that the posterior measure µ y is welldefined.
Theorem 8 Let y ∈ C m . Assume that µ 0 has second moment. Then the posterior defined by (27) is a welldefined probability measure on 1 κ . Moreover there exist constants C 1 , C 2 > 0 independent of y such that the normalization constant satisfies Proof We follow (Stuart, 2010, Theorem 4.1). First notice that the potential Ψ (·, y) is µ 0 -measurable. To prove the lower bound of Λ(y) we apply Lemma 2 (ii) to get an appropriate set X ⊆ 1 κ with µ 0 (X) > 0 such that This implies Λ(y) > 0. For the upper bound we use Lemma 2 (i) and that µ 0 is a probability measure to immediately get This proves well-definedness of the posterior.
In the next theorem we show that the posterior measure is stable with respect to small variations in y in the Hellinger distance d Hell .
Proof First let us define the function This function satisfies Here we used Lemma 1 for the first inequality, Lemma 2 (v) for the second inequality and that µ 0 has second moment for the last inequality. The normalization constants as functions of y are Lipschitz continuous. Indeed, by Lemma 2 (i), (v) and again Lemma 1 we get The function x → x − 1 2 is locally Lipschitz for x > 0 and applying Theorem 8 so that With this result in mind take a look at two times the square of the Hellinger distance: Here, the last inequality follows from the previous estimates together with f ≤ 1 and Λ(y) ≥ C.
Theorem 8 has shown that the posterior is in fact welldefined. Often it is desirable to show existence of moments of it.
Theorem 10 Let y ∈ C m . Assume that µ 0 has p-th moment with p ≥ 2 and let (X, · X ) be a Banach space. Furthermore, let f : 1 κ → X be a µ 0 -measurable function such that for µ 0 -almost every u ∈ 1 κ there holds f (u) X ≤ C u p 1 . Then f ∈ L 1 ( 1 κ ; µ y ). Moreover, µ y has p-th moment.
Proof There holds where we have used that Ψ ≥ 0. This shows f ∈ L 1 ( 1 κ ; µ y ). Observe that the posterior µ y has p-th moment since we can choose f = id, where id is the identity on 1 κ .
Theorem 11 Assume that µ 0 has second moment and let (X, · X ) be a Banach space. Let f : 1 κ → X satisfy f ∈ L 2 ( 1 κ ; µ y1 ) and f ∈ L 2 ( 1 κ ; µ y2 ). Let r > 0 and y 1 , y 2 ∈ B r (0) ⊂ C m . Then there is a C > 0 such that In particular, computing moments is stable under small perturbations in the measurements y.

Sampling
In this section, we derive a method to sample from the posterior. First, we introduce a discretized observation operator to obtain a computable discretized posterior measure and show several properties of it. We proceed to apply a Sequential Monte Carlo method to generate samples from the posterior. Finally, we derive a method to decrease the complexity of evaluating the observation operator multiple times.

Discrete Approximation of the Posterior Measure
In practice we are not able to compute the solution of the Helmholtz equation exactly. As a consequence we have to replace the observation operator G by an approximation G h , where h is a discretization parameter. In our setting h is the mesh diameter. The discrete observation operator is then given by The next lemma states several properties of the discrete observation operator and its relation to G.
Lemma 3 There exists a h 0 > 0 such that for every h ∈ (0, h 0 ] and every u ∈ 1 κ the discrete observation operator satisfies with constants C independent of u and h. Furthermore, G h is µ 0 -measurable and if µ 0 has p-th moment then G h ∈ L p ( 1 κ ; µ 0 ).
Proof From Theorem 6 we have with a constant C > 0 independent of u and h. This estimate applied to the measurement points (z j ) m j=1 ∈ M m κ shows (29). By the triangle inequality we have which implies (28) using the a-priori estimate, h ≤ h 0 and Proposition 9. G h is µ 0 -measurable because it is continuous by Corollary 2. In particular, An application of (28) under the assumption that µ 0 has p-th moment shows G h ∈ L p ( 1 κ ; µ 0 ).
We further have to work with a discretized potential Ψ h . Let us fix some arbitrary measurement observation y ∈ C m in (4) and define the discrete potential The next lemma states some essential properties of Ψ h .
Lemma 4 There exist C, h 0 > 0 such that for every h ∈ (0, h 0 ] the discrete potential satisfies for all u ∈ 1 κ , y ∈ C m . .
Furthermore, Lemma 2 is valid for Ψ h instead of Ψ h with constants and sets independent of h.
Proof We compute Now the a-priori error estimate in Lemma 3 together with G(u) ≤ C u 1 and Young's inequality imply .
That Lemma 2 holds with Ψ replaced by Ψ h with constants and sets independent of h is shown by a straight forward computation following the proof of Lemma 2. We omit most of it for brevity. The only slight difficulty is to define the set X ⊆ 1 κ in (ii) independently of h. Therefore, let and conclude µ 0 (X h ) ≥ 1 2 for any h ∈ (0, h 0 ] as shown in the proof of (ii) of Lemma 2. Using Young's inequality we estimate sup h∈(0,h0] where we used h ≤ h 0 , that µ 0 has second moment and G ∈ L 2 ( 1 κ ; µ 0 ). Hence, setting implies X h ⊆ X and therefore µ 0 (X) ≥ 1 2 . Now we apply Young's inequality to obtain for all u ∈ X, y ∈ C m . This proves the statement.
We apply Lemma 4 and Lemma 3 to show that Lemma 2 holds if we replace G and Ψ by their discrete counterparts G h and Ψ h . We emphasize that all the estimates in Lemma 2 are valid uniformly in h as long as h ∈ (0, h 0 ] for a suitably small h 0 . Then applying Theorem 8 ensures that the discrete posterior µ y h is well-defined, i.e. µ y h is defined by the Radon-Nikodym derivative with respect to the prior Here, Λ h (y) is a normalization constant defined by We summarize this result in the next theorem.
Theorem 12 Let y ∈ C m . Assume that µ 0 has a finite second moment. Then there exists h 0 > 0 such that for any h ∈ (0, h 0 ] the discrete posterior measure µ y h is welldefined. There exist constants C 1 , C 2 > 0 independent of h such that the normalization constants satisfy In particular (30) yields the well-definedness of µ y h .
Proof The proof is analogous to the proof of Theorem 8 using the fact that Lemma 2 holds for the discrete potential Ψ h with constants and sets independent of h.
The next theorem states that under certain conditions the rate of convergence in h of the observation operator carries over to the rate of convergence of µ y h to µ y in the Hellinger distance.
Theorem 13 Let y ∈ C m . Assume that µ 0 has a finite fourth moment. Then there exists a h 0 > 0 such that for any h ∈ (0, h 0 ] we have d Hell (µ y , µ y h ) ≤ C| ln h|h 2 .
C does not depend on h, but on y.
Proof Proving this theorem is similar to proving Theorem 9. We define the functions Using the properties of the exponential function from Lemma 1 together with Lemma 4 yields where we have used hat µ 0 has fourth moment in the last inequality. For the normalization constants observe that We are able to bound Λ(y) and Λ h (y) away from zero independently of h using Theorem 8 and Theorem 12 which shows that because x → x − 1 2 is locally Lipschitz for x > 0. A straight forward computation yields the result: Here we have used the previous estimates, f h ≤ 1 for any h and Λ(y) > 0.
Theorem 14 Let y ∈ C m . Assume that µ 0 has a finite fourth moment. Let (X, · X ) be a Banach space and let f : 1 κ → X have second moments with respect to both µ y and µ y h . Then there holds Proof As in the proof of Theorem 11 this follows from (Stuart, 2010, Lemma 6.37) using Theorem 13.

Sequential Monte Carlo
Throughout Section 4.2 let y ∈ C m be a given observation. In the following section we will use the Sequential Monte Carlo Method (SMC) from (Dashti and Stuart, 2015) to draw samples from the posterior measure. We also derive an error estimate.
Let J ∈ N and for j ∈ {0, ..., J} define a sequence of measures µ j << µ 0 by Note that µ 0 is equal to the prior µ 0 and µ J equal to the posterior measure µ y . Our goal is to approximate µ J sequentially using information of each µ j to construct the next approximation µ j+1 . One idea behind SMC is the approximation of each measure µ j by a weighted sum of Dirac measures with u (n) j ∈ 1 κ and weights w (n) j ≥ 0 that sum up to 1. We define the operator Here P( 1 κ , F) denotes the space of probability measures on ( 1 κ , F). We also define We remark that the operator L satisfies µ j+1 = Lµ j , j = 1, ..., J.
We further define P j : 1 κ ×F → [0, 1] as µ j -invariant Markov kernel. That means The idea of the kernel is to redraw samples in each iteration of the algorithm to better approximate µ j+1 . (See 4. in Algorithm 1.) This allows us to define the discrete measures according to (Dashti and Stuart, 2015, Algorithm 5.12.) as follows Algorithm 1: SMC-Algorithm N , n = 1, · · · , N and define µ N j by (32). 4. Resampling: drawû In Algorithm 1 we see the SMC-Algorithm described in (Dashti and Stuart, 2015, Section 5.3.).
Let us generalize (Dashti and Stuart, 2015, Theorem 5.13) in such a way that it applies to our setting.
Theorem 15 For every measurable and bounded function f the measure µ N J satisfies where E SMC is the expectation with respect to the randomness in the SMC algorithm.
Proof We first prove a variant of (Dashti and Stuart, 2015, Lemma 5.17) using similar techniques. Define g := exp(−J −1 Ψ ) and let f with f ∞ ≤ 1 be given. Then from the proof of (Dashti and Stuart, 2015, Lemma 5.17) and defining η j = S N P j µ N j we conclude A quick calculation shows Hence we obtain We define the distance d op for probability measures ν, η on ( 1 κ , F) as follows Now we show that (Dashti and Stuart, 2015, Theorem 5.13) holds with weaker assumptions. First we use (35) and the triangle inequality to conclude A straight forward continuation similar to the proof of (Dashti and Stuart, 2015, Theorem 5.13) leads to the statement Here we remark that both (Dashti and Stuart, 2015, Lemma 5.15) and (Dashti and Stuart, 2015, Lemma 5.16) needed for the proof of (Dashti and Stuart, 2015, Theorem 5.13) still hold. From (36) we conclude that (33) holds for f ∞ ≤ 1 and the statement for a measurable bounded function f follows by a scaling argument.
The prior source positions are uniformly distributed . We in particular assume that the random variables α k j and x k j are independent. For the observational noise we consider η = (η l ) 3 l=1 with η 1 , η 2 , η 3 all uniformly and identically distributed with mean The domain D is uniformly triangulated with triangles of diameter h. Throughout our experiments we choose the Markov kernel which is clearly invariant with respect to µ j (or any other measure).

The SMC Method in Action
We now present results on the source positions of samples from the posterior. We consider the SMC as introduced in Section 4 with J = 1 and N = 448000. Since we then only have one tempering step for j = J = 1 we drop the j from the index. We choose a uniform mesh of size h = √ 2 · 2 −5 for all experiments in this section. Algorithm 1 delivers the following data k (n) ∈ N 0 , κ for n = 1, . . . , N and weights w (1) , . . . , w (N ) . We define the function P emp (x) as an approximation to the probability that a source is located at x ∈ D in the following sense ≈ P µ y (u has a source in B ε (x)). Here ε = 0.04 and K ε is a smooth cutoff function approximating the indicator function 1 Bε(0) . In Figure 1 (a) we can see the function P emp decently recovers the true source positions x exact,1 and x exact,l .
We further want to analyze positions of pairs. In particular we would like to find out if it happens that for a sample two positions are close to the same x exact,l (l = 1 or l = 2). That would mean we do not identify two distinct sound sources. Given that a sample has two sources, i.e. k = 2, and one is in the set Q we look for the probability that the second source is close to some x ∈ D but not in Q. Formally for the index set I(Q) = {n ∈ {1, ..., N }|k (n) = 2, x K ε x − x (n) l ≈ P µ y (u has a source in B ε (x) \ Q | u has a source in Q, k = 2).
The results for this function are given in Figure 1 (b) and (c). We conclude that if one of the sources is located near x exact,1 then the other near x exact,2 and vice versa. We also would like to address the notion of a maximum a posteriori estimator (MAP). We define the empirical MAP-index and the empirical MAP-index given k sources as follows In our experiments both quantities are well-defined and unique. The empirical MAP-estimator is the n MAP -th sample of the last (and only) iteration j = J = 1 of the SMC. Analogously we define the empirical MAPestimator given k sources. Further results of this experiment are given in Table 1 and visualized in Figure  2. Note that Table 1 contains only results for k = 2, 3 number of sound sources. The occurrence probability for k = 2, 3 sound sources is below 1% and therefore not mentioned.

Convergence in Mean Square Error for functions
Theorem 15 shows linear convergence in the Mean Square Error (MSE) for bounded measurable functions f integrated over the SMC posterior µ N J . We are interested testing the convergence of different functions f 1 , ..., f 5 listed in table 2. These functions allow us to extract certain information from the posterior measure. We chose z prediction := (1/2, 1/4) ∈ D to predict the pressure at a point distinct from the measurement points z 1 , ..., z 3 . For f 5 we chose time t = 1. We remark that Theorem 15 shows convergence results only for f 2 . Nevertheless we observe convergence for the other functions in figure 3. For this experiment we fix the mesh size to h ref = √ 2 · 2 −5 and only vary the number of samples. We chose N ref = 448.000 and use µ N ref h ref as reference measure. We used j=J=1 as tempering step. For the expectation integral "E SMC " we used a standard Monte Carlo method with 50 runs.

Convergence in the Hellinger Distance
We also verify Theorem 13 numerically. It states convergence of the discretized to the true posterior with

Convergence in h for functions
We want to verify Theorem 14 which states the convergence of the error for functions f with second moments with respect to both µ y and µ y h . We use the functions f 1 , ..., f 5 defined in Table 2 and a reference measure µ y ref : For all experiments we choose N ref = 20.000 samples from the SMC sampler with temperature J = 1. For this experiment we do not use the resampling step of the SMC sampler. For the other measures we vary the mesh parameter h = √ 2 · 2 −k for k = 2, ..., 6 and average the approximation error e h in (38) over 20 runs. The results are seen in figure 5 and verify the stated convergence for the functions f 1 , ..., f 5 .

Conclusion
In this paper we studied source identification for the Helmholtz equation using the Bayesian approach. We derived suitable properties and regularity for the solution of the Helmholtz equation with a series of Diracs as forcing term. We defined a prior probability distribution on the number of sources, their position and intensity and proved well-posedness of the Bayesian inverse problem. We have shown that the discretized so-lution operator for the Helmholtz equation leads to a discretized posterior which converges to the true posterior in the Hellinger distance. We used the Sequential Monte Carlo method to sample from the posterior. In the last section, we applied these techniques to verify our results numerically.