  The quantum mechanical position operators, and their products, are not
well-defined in systems obeying periodic boundary conditions. Here we extend
the work of Resta who developed a formalism to calculate the electronic
polarization as an expectation value of a many-body operator, to include higher
multipole moments, e.g., quadrupole and octupole. We define $n$-th order
multipole operators whose expectation values can be used to calculate the
$n$-th multipole moment when all of the lower moments are vanishing (modulo a
quantum). We show that changes in our operators are tied to flows of $n-1$-st
multipole currents, and encode the adiabatic evolution of the system in the
presence of an $n-1$-st gradient of the electric field. Finally, we test our
operators on a set of tightbinding models to show that they correctly determine
the phase diagrams of topological quadrupole and octupole models, capture an
adiabatic quadrupole pump, and distinguish a bulk quadrupole moment from other
mechanisms that generate corner charges.
