We introduce a new method for analysing the Bose-Hubbard model for an array of bosons with nearest neighbor interactions. It is based on a number-theoretic implementation of the creation and annihilation operators that constitute the model. One of the advantages of this approach is that it facilitates computation with arbitrary accuracy, enabling nearly perfect numerical experimentation. In particular, we provide a rigorous computer assisted proof of quantum phase transitions in finite systems of this type. Furthermore, we investigate properties of the infinite array via harmonic analysis on the multiplicative group of positive rationals. This furnishes an isomorphism that recasts the underlying Fock space as an infinite tensor product of Hecke spaces, i.e., spaces of square-integrable periodic functions that are a superposition of non-negative frequency harmonics. Under this isomorphism, the number-theoretic creation and annihilation operators are mapped into the Kastrup model of the harmonic oscillator on the circle. It also enables us to highlight a kinship of the model at hand with an array of spin moments with a local anisotropy field. This identifies an interesting physical system that can be mapped into the model at hand.