In this paper, a novel physics-aware harmonic state estimation (HSE) method is developed for power distribution systems with low-observability. The only available measurements are the harmonic synchrophasors from a small number of Harmonic Phasor Measurement Units (H-PMUs) on a power distribution feeder. The proposed HSE method provides a novel and practical application of H-PMUs, which are an emerging class of smart grid sensors, in order to address an important and challenging problem in power distribution system monitoring. More importantly, the proposed HSE method can address the highly challenging scenario, in which, not only the network has a low-observability condition, but also the number and the location(s) of the harmonic sources are unknown. The proposed HSE method captures physics-based sparsity pattern in the analysis of harmonic phasors, and then it incorporates such sparsity pattern in the formulation of a novel mixed-integer quadratic programming (MIQP) optimization which can be solved by standard solvers. The method that is proposed in this paper to address low-observability is innovative and it is specific to HSE. It does not have any similar counterpart in the literature, whether for harmonic state estimation for in ordinary (i.e., not harmonic) distribution system state estimation. The effectiveness of the proposed method is tested and confirmed in multiple case studies and compared with the existing methods.