  In this paper, a generalized finite element method (GFEM) with optimal local
approximation spaces for solving high-frequency heterogeneous Helmholtz
problems is systematically studied. The local spaces are built from selected
eigenvectors of carefully designed local eigenvalue problems defined on
generalized harmonic spaces. At both continuous and discrete levels, $(i)$
wavenumber explicit and nearly exponential decay rates for local and global
approximation errors are obtained without any assumption on the size of
subdomains; $(ii)$ a quasi-optimal convergence of the method is established by
assuming that the size of subdomains is $O(1/k)$ ($k$ is the wavenumber). A
novel resonance effect between the wavenumber and the dimension of local spaces
on the decay of error with respect to the oversampling size is implied by the
analysis. Furthermore, for fixed dimensions of local spaces, the discrete local
errors are proved to converge as $h\rightarrow 0$ ($h$ denoting the mesh size)
towards the continuous local errors. The method at the continuous level extends
the plane wave partition of unity method [I. Babuska and J. M. Melenk,
Int.\;J.\;Numer.\;Methods Eng., 40 (1997), pp.~727--758] to the
heterogeneous-coefficients case, and at the discrete level, it delivers an
efficient non-iterative domain decomposition method for solving discrete
Helmholtz problems resulting from standard FE discretizations. Numerical
results are provided to confirm the theoretical analysis and to validate the
proposed method.
