AbstractWe propose four approaches to solve time‐harmonic Maxwell's equations in ℝ3 through the finite element method (FEM) in a bounded region encompassing parameter inhomogeneities, coupled with the multiple multipole program (MMP) in the unbounded complement. MMP belongs to the class of methods of auxiliary sources and of Trefftz methods, as it employs point sources that spawn exact solutions of the homogeneous equations. Each of these sources is anchored at a point that, if singular, is placed outside the respective domain of approximation. In the MMP domain, we assume that material parameters are piecewise constant, which induces a partition into one unbounded subdomain and other bounded, but possibly very large, subdomains, each requiring its own MMP trial space. Hence, in addition to the transmission conditions between the FEM and MMP domains, one also has to impose transmission conditions connecting the MMP subdomains. Coupling approaches arise from seeking stationary points of Lagrangian functionals that both enforce the variational form of the equations in the FEM domain and match the different trial functions across subdomain interfaces. We discuss the following approaches: Least‐squares‐based coupling using techniques from PDE‐constrained optimization. Discontinuous Galerkin coupling between the meshed FEM domain and the single‐entity MMP subdomains. Multi‐field variational formulation in the spirit of mortar FEMs. Coupling through the Dirichlet‐to‐Neumann operator. We compare these approaches in a series of numerical experiments with different geometries and material parameters, including examples that exhibit triple‐point singularities and infinite layered media.