An accurate modeling of the optical interactions in metallic nanostructures with subnanometer features requires an accurate description of quantum effects at the scale of billions of atoms. At such scale, first-principle methods are not computationally viable. Quantum hydrodynamic theory (QHT) has emerged as a powerful method that includes nonlocal contributions of the kinetic energy and the spatial dependence of the electron density, and it can predict both plasmon energy and spill-out effects in large metal nanoparticles. In this paper, we introduce a hybridizable discontinuous Galerkin method for solving Maxwell's equations coupled with a QHT model in order to account for quantum effects in three-dimensional metallic nanostructures. The coupled system of Maxwell's equations and QHT model is not only nonlinear but also multi-scale due to the interaction between the micrometer electromagnetic waves and the nanometer cavities of metallic nanostructures. We present extensive numerical experiments to validate the QHT model and demonstrate the capability of the HDG method to provide accurate solutions in the presence of strong nonlinearities and multiple length scales. These results offer a possibility to enhance nonlinear optical effects or to harness quantum mechanical electron tunneling by engineering metallic nanostructures at the quantum level. •An HDG method to simulate the equilibrium equation of quantum hydrodynamic theory.•An HDG method to simulate the first-order linear system of quantum hydrodynamic theory coupled with Maxwell's equations.•2D simulation of metal nanowire and full 3D simulation of a periodic nanocoax array accounting for quantum effects.•Comparison of the optical response using quantum effects to that of local and nonlocal electron models for both applications.