A finite field method for calculating spherical tensor molecular polarizability tensors αlm;l′m′ = ∂Δlm/∂ϕl′m′* by numerical derivatives of induced molecular multipole Δlm with respect to gradients of electrostatic potential ϕl′m′* is described for arbitrary multipole ranks l and l′. Inter-conversion formulae for transforming multipole moments and polarizability tensors between spherical and traceless Cartesian tensor conventions are derived. As an example, molecular polarizability tensors up to the hexadecapole-hexadecapole level are calculated for water at the HF, B3LYP, MP2, and CCSD levels. In addition, inter-molecular electrostatic and polarization energies calculated by molecular multipoles and polarizability tensors are compared to ab initio reference values calculated by the Reduced Variation Space (RVS) method for several randomly oriented small molecule dimers separated by a large distance. It is discussed how higher order molecular polarizability tensors can be used as a tool for testing and developing new polarization models for future force fields.