In this paper, we develop a Taylor expansion (TE) based fast multipole method (FMM)for low frequency 3D Helmholtz Green’s function in layered media. Two forms ofTaylor expansions, with either non-symmetric or symmetric derivatives of layered mediaGreen’s functions, are used for the implementations of the proposed TE-FMM. In theimplementation with non-symmetric derivatives, an algorithm based on discrete compleximage approximations and recurrence formulas is shown to be very efficient and accuratein computing the high order derivatives. Meanwhile, the implementation based onsymmetric derivatives is more robust and pre-computed tables for the high orderderivatives in translation operators are used. Numerical tests in layered media havevalidated the accuracy and O(N) complexity of the proposed algorithms.