with k > 2 and m x,y,z ∈ N 0 arise in various contexts of quantum chemistry. For instance, the operator representing the polarization energy, Ĥpol = -1 2 f A f , which depends quadratically on the electric field operator f and some effective polarizability tensor A, is of this type. For a single atom having a dipole polarizability α interacting with a single electron, the polarization potential takes the form V pol (r) = -1 2 αr -4 . Integrals of such polarization operators are found for example in core polarization potentials (CPP) that are used to model core-valence electron correlation in atoms and molecules and the direct reaction field method (DRF) which is a polarizable embedding scheme for quantum mechanics/molecular mechanics (QM/MM) simulations. DRF requires matrix elements of polarization operators between basis functions on different atomic centers. Historically these integrals were approximated by a Taylor expansion of the electric fields around an arbitrary reference position to avoid the exact evaluation. However, since the choice of the expansion center is not a differentiable operation, analytic gradients of the polarization energy with respect to nuclear coordinates were not easily accessible. Quadrupole and higher-order polarizability corrections also require matrix element for a range of different values of m x,y,z and k. In general, terms of the form in Eqn. appear in multipole expansions of the electrostatic potential.