Overcoming Computational Bottlenecks in Quantum Hydrodynamics: A Volume-Based Integral Formalism

TL;DR

The paper addresses computational bottlenecks in simulating quantum corrections to plasmonic responses by introducing a volume-integral-equation (VIE) approach for the self-consistent hydrodynamic Drude model (SC-HDM). By exploiting spherical symmetry and symmetry-inspired basis functions, the method yields an effectively 1D radial discretization, enabling efficient treatment of nonlocality and electron spill-out while remaining modular and extendable to more sophisticated models. Key contributions include a solver that recovers classical local and nonlocal results, accurate reproduction of SC-HDM-specific features such as the Bennett resonance, and a practical route to extract Feibelman parameters ($d_ot$) for integration with surface-response models. The framework promises a scalable, benchmark-enabled pathway for quantum hydrodynamic modeling of nanoparticles and complex geometries, with potential to feed semiclassical models like SRM and reduce reliance on expensive TD-DFT calculations.

Abstract

Mesoscopic models of the optical response of metals have emerged as fundamental building blocks in quantum plasmonics, in principle overcoming the computational bottlenecks of ab initio techniques by implementing aspects of the atomistic description of the metal in otherwise classical calculations. Nonetheless, even these approaches are eventually hindered by demanding computations due to sophisticated material response. Here, this issue is addressed for the advanced Self-Consistent Hydrodynamic Drude Model (SC-HDM), which captures both nonlocal electron dynamics and electron spill-out, through a Volume Integral Equation (VIE) method. Adopting an IE-based method shifts perspective from the commonly employed Differential Equation (DE)-based ones, demonstrating significant computational efficiency. The VIE approach is a valuable methodological scaffold: It addresses SC-HDM and simpler models, but can also be adapted to more advanced ones. For spherical nanoparticles (NPs), using the inherent symmetries, similar performance for three increasingly complicated material models is achieved, breaking the taboo that increased sophistication in material response requires taxing simulations. Mesoscopic material-response functions can be readily extracted from the VIE implementation, thus circumventing the need for lengthy microscopic calculations. This method opens a new way of modeling quantum hydrodynamic NPs and will serve as essential benchmarking tool for recipes addressing more complicated geometries.

Figures (4)

• Figure 1: (a) A local spherical NP of radius $R$ is illuminated by a plane wave described by electric field $\mathbf{E}$, while the direction of propagation is denoted by the wave vector $\mathbf{k}$. Notice the well-defined boundary, due to the hard-wall assumption. (b) A spherical NP described by SC-HDM. (c) Zoom-in to depict better the spill-out. A soft shell wraps around the NP's hard ionic core, which is described by the jellium model and terminates at $r=R$. The shell, representing the decaying electron tail beyond the geometric boundary, formally extends to infinity, but practically is assigned a thickness $s$, the spill-out allowance. (d) Evolution of the treatment of the SC-HDM NP, from a DE-based perspective to an IE-based one and finally one exploiting the symmetries of the spherical geometry. In the first, careful meshes must be designed. The NP terminates at a distance $r=R+s$ from the origin. This nanosphere must be very finely discretized, in order to capture the wavelength of the matter waves of electrons (in the order of a few Angström Mortensen21). Beyond the nanosphere extends a dielectric medium of zero electron ground-state and induced electron density. This medium is physically infinite; in DE-based methods a finite box is required. An adaptive mesh, dense in the near-field region (in order to capture its fast variations) and sparser far from it is typically used. To minimize the error of this abrupt truncation, Perfectly Matched Layers (PMLs) are employed, which must be also properly meshed. Our VIE does not enter into such discussions, discretizing only the nanosphere and encapsulating the infinite background in the Green's dyad. However, even this meshing redundant in light of the symmetries of the problem: The use of basis and testing functions such as the ones of Equations (\ref{['eq:theory_basis']}) and (\ref{['eq:theory_testing']}) enables significant analytical progress and leads to the possibility of 1D discretization, as suggested in the bottom schematic (the two concentric circles represent the jellium edge and the truncation of the NP, from in to out.
• Figure 2: Validation of our computational scheme through OpenSANS and quasistatic theory. (a) Extinction spectra of Na spherical NPs, normalized to the geometrical cross-section $\pi R^2$, within the LRA description. The dashed line indicates the LRA quasistatic result $\omega_p/\sqrt{3}$. (b) Relative error between the extinction spectra of the VIE approach and OpenSANS regarding LRA. (c) Extinction spectra of Na spherical NPs, normalized to the geometrical cross-section $\pi R^2$, using the HDM description. The dashed line indicates the LRA quasistatic result $\omega_p/\sqrt{3}$. (d) Relative error between the extinction spectra of the VIE approach and OpenSANS regarding HDM. (e) Comparison of the dipolar LSP as a function of the NP radius between OpenSANS (crosses), the VIE approach (circles) for the LRA and the HDM with the relevant quasistatic results (LRA is represented by the red line and HDM by the blue one).
• Figure 3: Comparison of results between LRA, HDM, and SC-HDM. (a) Normalized extinction cross sections within LRA (blue line), HDM (orange line), and SC-HDM (black line) Dipolar LSPs are marked with a D, while the Bennett resonance is marked with a B. (b) The Feibelman parameter $d_\perp$ retrieved from the SC-HDM calculation. (c) Induced charge density for the dipolar LSP resonance of HDM (orange D in panel (a)). (d) Induced charge density for the dipolar LSP resonance of SC-HDM (black D in panel (a)). (e) Induced charge density for the Bennett resonance of SC-HDM (black B in panel (a)). The blue circle is a guide to the eye which corresponds to the end of the nanosphere for HDM and the jellium edge for SC-HDM. (f) Simulation time for each model (averaged over $20$ runs).
• Figure : Numerical modeling of mesoscopic material response models that capture the quantum dynamics of electrons does not have to come with discouraging computational bottlenecks. The main message is that through a shift in perspective in modeling towards Integral Equation methods and exploiting symmetry-based arguments, it is possible to capture complicated phenomena such as nonlocality and electron spill-out using a personal laptop.