The formulation of density-functional expansion methods is extended to treat the second and higher-order terms involving the response density and spin densities with an arbitrary single-center auxiliary basis. The two-center atomic orbital products are represented by the auxiliary functions centered about those two atoms, and the mapping coefficients are determined from a local constrained variational procedure. This two-center variational procedure allows the mapping coefficients to be pretabulated and splined as a function of internuclear separation for efficient look up. The splines of mapping coefficients have a range no longer than that of the overlap integrals, and the auxiliary density appears as a single point-multipole expansion to all nonoverlapping atoms, thus allowing for the trivial implementation of a linear-scaling algorithm. The method is tested using Gaussian multipole expansions, and the effect of angular and radial completeness is explored. Several auxiliary basis sets are parametrized and compared to an auxiliary basis analogous to that used in the self-consistent-charge density-functional tight-binding model, and the method is demonstrated to greatly improve the representation of the density response with respect to a reference expansion model that does not use an auxiliary basis.
