1 IntroductionWave propagation simulations, governed by the Helmholtz equation, in bounded heterogeneous and unbounded homogenous media are fundamental for numerous applications [13, 33, 39].Finite element methods (FEM) are efficient for simulating the Helmholtz equation in a bounded heterogeneous medium, say, Ω0 Rm (m = 2, 3). The standard (non-coercive) variational formulation of the variable coefficient Helmholtz equation in H1(Ω0) [33] has been widely used for developing and analyzing the sign-indefinite FEM, see for example [3,7,12,27,29, 40]. The open problem of developing a coercive variational formulation for the heterogeneous Helmholtz model was solved recently in [28], and an associated preconditioned sign-definite high-order FEM was also established using direct and domain decomposition methods in [28]. For a large class of applications the wave propagation occurs in the bounded heterogeneous medium and also in its complement, Rm Ω0, the exterior unbounded homogeneous medium. Using the fundamental solution, the constant coefficient Helmholtz equation exterior to Ω0 can be reformulated as an integral equation (IE) on the boundary of Ω0. Algorithms for simulating the boundary IE (BIE) are known as boundary element methods (BEM). Several coercive and non-coercive BIE reformulations [13,39] of the exterior Helmholtz model have been used to develop algorithms for the exterior homogeneous Helmholtz models, see for example the acoustic BEM survey articles [11, 34], respectively, by mathematical and engineering researchers, eachwith over 400 references.The exterior wave propagation BEM models lead to dense complex algebraic systems, and the standard variational formulation based interior wave FEM models lead to sparse complex systems with their eigenvalues in the left half of the complex plane [26, 38]. Developing efficient preconditioned iterative solvers for such systems has also dominated research activities over the last two decades [19], in conjunction with efficient implementations using multigrid and domain decomposition techniques, see [25, 27] and references therein.For applications that require solving both the interior heterogenous and exterior homoge- neous problems, various couplings of the FEM and BEM algorithms with appropriate con- ditions on polygonal interfaces have also been investigated in the literature [5, 6, 32]. The review article [43] describes some theoretical validations of the coupling approaches considered in the earlier literature and delicate choices of the coupling interface. The coupling methods in [5,6,31,32,43] lead to very large algebraic systems with both dense and sparse structures. For wave propagation models, given the complexity involved in even separately solving the FEM and BEM algebraic systems, it is efficient to avoid large combined dense and sparse structured systems arising from the coupling methods in [5, 6, 31, 32, 43].Such complicated-structured coupled large-scale systems can be avoided, for the Helmholtz PDE interior and exterior problems, using the approach proposed in [35] and recently further explored in [24] using high-order elements for a class of applications with complex hetero- geneous structures. The FEM-BEM algorithms in [24, 35] are based on the idea of using a non-overlapping smooth interface to couple the interior and exterior solutions. As described in [24, Section 6], there are several open mathematical analysis problems remain to be solved in the coupling and FEM-BEM framework of [24, 35].The choice of smooth interface in the FEM-BEM algorithms of [24, 35] is crucial because the methods require solving several interior and exterior wave problems to setup the interface