This article aims at a proper definition and resolution of the parabolic Anderson model on Heisenberg groups $\mathbf{H}_{n}$. This stochastic PDE is understood in a pathwise (Stratonovich) sense. We consider a noise which is smoother than white noise in time, with a spatial covariance function generated by negative powers $(-Δ)^{-α}$ of the sub-Laplacian on $\mathbf{H}_{n}$. We give optimal conditions on the covariance function so that the stochastic PDE is solvable. A large portion of the article is dedicated to a detailed definition of weighted Besov spaces on $\mathbf{H}_{n}$. This definition, related paraproducts and heat flow smoothing properties, forms a necessary step in the resolution of our main equation. It also appears to be new and of independent interest. It relies on a recent approach, called projective, to Fourier transforms on $\mathbf{H}_{n}$.