  In this paper, we will introduce a new heterogeneous fast multipole method
(H-FMM) for 2-D Helmholtz equation in layered media. To illustrate the main
algorithm ideas, we focus on the case of two and three layers in this work. The
key compression step in the H-FMM is based on a fact that the multipole
expansion for the sources of the free-space Green's function can be used also
to compress the far field of the sources of the layered-media or domain Green's
function, and a similar result exists for the translation operators for the
multipole and local expansions. The mathematical error analysis is shown
rigorously by an image representation of the Sommerfeld spectral form of the
domain Green's function. As a result, in the H-FMM algorithm, both the
"multipole-to-multipole" and "local-to-local" translation operators are the
same as those in the free-space case, allowing easy adaptation of existing
free-space FMM. All the spatially variant information of the domain Green's
function are collected into the "multipole-to-local" translations and therefore
the FMM becomes "heterogeneous". The compressed representation further reduces
the cost of evaluating the domain Green's function when computing the local
direct interactions. Preliminary numerical experiments are presented to
demonstrate the efficiency and accuracy of the algorithm with much improved
performance over some existing methods for inhomogeneous media. Furthermore, we
also show that, due to the equivalence between the complex line image
representation and Sommerfeld integral representation of layered media Green's
function, the new algorithm can be generalized to multi-layered media with
minor modification where details for compression formulas, translation
operators, and bookkeeping strategies will be addressed in a subsequent paper.
