  An accurate description of the optical response of subwavelength metallic
particles and nanogap structures is a key problem of plasmonics. Quantum
hydrodynamic theory (QHT) has emerged as a powerful method to calculate the
optical response of metallic nanoparticles (NPs) since it takes into account
nonlocality and spill-out effects. Nevertheless, the absorption spectra of
metallic NPs obtained with conventional QHT, i.e., incorporating Thomas-Fermi
(TF) and von Weizs\"acker (vW) kinetic energy (KE) contributions, can be
affected by several spurious resonances at energies higher than the main
localized surface plasmon (LSP). These peaks are not present in reference
time-dependent density-functional-theory (TD-DFT) spectra, where, instead, only
a broad shoulder exists. Moreover, we show here that these peaks incorrectly
reduce the LSP peak intensity and have a strong dependence on the simulation
domain size so that a proper calculation of QHT absorption spectra can be
problematic. In this article, we introduce a more general QHT method accounting
for KE contributions depending on the Laplacian of the electronic density
($q$), thus, beyond the gradient-only dependence of the TFvW functional. We
show that employing a KE functional with a term proportional to $q^2$ results
in an absorption spectrum free of spurious peaks, with LSP resonance of correct
intensity and numerically stable Bennett state. Finally, we present a novel
Laplacian-level KE energy functional that is very accurate for the description
of the optical properties of NPs with different sizes as well as for dimers.
