Quantum adiabatic optimization has long been expected to outperform classical methods in solving NP-type problems. While this has been proven in certain experiments, its main applications still reside in academic problems where the size of the system to be solved would not represent an obstacle to any modern desktop computer. Here we develop a systematic procedure to find the global optima of the modularity function to discover community structure in complex networks solely relying on pure annealers rather than hybrid solutions. We bypass the one-hot encoding constraints by hierarchically and recursively encoding binary instances of the problem that can be solved without the need to guess the exact penalties for the Lagrange multipliers. We study the variability, and robustness of the annealing process as a function of network size, directness of connections, topology, and the resolution of the communities. We show how our approach produces meaningful and at least equally optimal solutions to state-of-the-art community detection algorithms while maintaining tractable computing times. Lastly, due to its recursive nature, the annealing process returns intermediate subdivisions thus offering interpretable rather than black-box solutions. These \textit{dendrograms} can be used to unveil normal and pathological hidden hierarchies in brain networks hence opening the door to clinical workflows. Overall, this represents a first step towards an applicable practice-oriented usage of pure quantum annealing potentially bridging two segregated communities in modern science and engineering; that of network science and quantum computing.