and implementation details. In Section 4 we demonstrate the efficiency of the FEM-BEM algorithm for simulating wave propagation in two distinct classes of (smooth and non-smooth) heterogenous media.2 Decomposition framework and well-posedness analysisLet Ω0 Rm, m = 2, 3, be a bounded domain. The ratio of the speed of wave propagation inside the heterogeneous (and not necessarily connected) region Ω0 and on its free-space exterior Ωc := Rm Ω0 is described through a refractive index function n that we assume in this articleto be piecewise smooth with 1 − n having compact support in Ω0 (i.e, n|Ωc ≡ 1).The main focus of this article is to study the wave propagation in Rm, induced by theimpinging of an incident wave uinc, say, a plane wave with wavenumber k > 0. More precisely, the continuous wave propagation model is to find the total field u(:= us + uinc) ∈ H1 (Rm) that satisfies the Helmholtz equation and the Sommerfeld RC:∆u + k2n2 u = 0, in Rm, (2.1) . ∂ us − ikus = o(|r| m+1 ), as |r| → ∞.It is well known that (2.1) is uniquely solvable [36]. (Later in this section, we introduce the classical Sobolev spaces Hs, for s ≥ 0, with appropriate norms.) A decomposition frameworkThe heterogeneous-homogeneous model problem (2.1) is decomposed by introducing two artifi- cial curves/surfaces Γ and Σ with interior Ω1 and Ω2 respectively satisfying Ω0 Ω1 Ω1 Ω2. We assume from now on that Γ is smooth and Σ is a polygonal/polyhedral boundary. A sketch of the different domains is displayed in Figure 1. Henceforth, Ωc := Rm Ωi, i = 0, 1, 2.We introduce the following decomposed heterogeneous and homogeneous media auxiliarymodels:For a given function f inp H1/2(Σ), we seek a propagating wave field w so that w and its trace γΣw on the boundary Σ satisfy ∆w + k2n2 w = 0, in Ω2, (2.2) Throughout the article, we assume that this interior problem is uniquely solvable. We introduce the following operator notation for the heterogeneous auxiliary model: For any Lipschitz m- or (m 1)-dimensional (domain or manifold) D Ω2, we define the solution operator KDΣ associated with the auxiliary model (2.2) asKDΣf inp := w|D. (2.3)Two cases will be of particular interest for us: KΩ Σf inp, which is nothing but w satisfy-2 Σing (2.2), and KΓΣf inp = γΓw, the trace of the solution w of (2.2) on Γ⊂ Ω2.