Laplacian-level quantum hydrodynamic theory for plasmonics

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 since it takes into account nonlocality and spill-out effects. Nevertheless, the absorption spectra of metallic particles from QHT is affected, at energies higher than the main plasmon peak, by several additional peaks, which are instead largely damped in reference time-dependent density-functional theory calculations. Moreover, we show here that these peaks have a strong dependence on the simulation domain-size so that the numerical convergence of QHT calculations is problematic. In this article, we introduce a QHT method accounting for kinetic energy contributions depending on the Laplacian of the electronic density, thus beyond the gradient-only dependence of conventional QHT. In this way, only the main plasmon peak is obtained, with a numerically stable absorption spectrum and more accurate intensities of the plasmon peak. Thus, the Laplacian-level QHT represents a novel, efficient and accurate platform to study plasmonic systems.

Another approach would be to treat the electron system semi-classically: a fluid, characterized by the macroscopic local quantities such as the electron density n (r, t) and the electron velocity field v (r, t) [25][26][27][28] but at the same time considering quantum effects through energy functionals of the electron density fluctuations. This approach is known as hydrodynamic theory (HT). The HT is part of a larger class of methods based on the orbitalfree (OF) [29][30][31] description of quantum electronic systems dating back to the works of Thomas [32] and Fermi [33]. Although the interest in OF-DFT methods has gradually decreased in favor of Kohn-Sham (KS) orbitalbased methods, the last decades have witnessed a reinvigorated interest due to the ideal scaling of computational resources with respect to the size of the electronic system offered by the OF-DFT approach [34]. Most of the research efforts in this field, however, have been devoted * cristian.ciraci@iit.it to static properties [35][36][37][38] and, more recently, also to response properties with the time-dependent OF-DFT [39][40][41][42]. In both cases, the central quantity that controls the accuracy of these methods is the noninteracting kinetic energy (KE) functional.
The most simple KE functional is the Thomas-Fermi (TF) functional, which accounts for the Pauli exclusion principle for a homogeneous system of noninteracting electrons [29] and it yields the electron quantum pressure p (r, t) ∝ n (r, t) 5/3 [43] that accounts for the nonlocal electron response. It has been demonstrated that TF-HT is able to provide surprisingly accurate predictions that match well experiments with noble metal NPs, such as Au [44] and Ag [45], both qualitatively and quantitatively. Nevertheless, for alkali metals or aluminum, the TF-HT predicts a blueshift of the localized surface plasmon resonance with respect to the classical Mie resonance [46], in contradiction with the redshift from the experiments [47] and TD-DFT calculations [48]. The origin of this difference lies in neglecting the spill-out of the plasmon-induced charges at the NP surface [46]. In fact, the TF-HT approach employs (with some recent exception [49]) a spatially uniform electronic density inside the NP and zero outside (i.e., hard-wall boundary) [15].
To properly address spill-out effects, the spatial dependence of electron density as well as a correction to the KE functional, in order to describe the density variation effects, must be introduced. The simplest functional that depends on the gradient of the density is the von Weizsäcker (vW) functional [29,50]. The TF-HT with a fraction (λ, with 0 < λ ≤ 1) of the vW correction (i.e., the TFλvW KE functional) is usually referred to in the literature as the QHT since the vW functional does not have a classical counterpart. The QHT has been largely used in plasma physics [51][52][53][54][55][56], and more recently, for plasmonic response properties of metal NPs of different geometries [57][58][59][60][61][62][63], as well as for surfaces [64][65][66] and strongly coupled plasmonic structures [62,67,68]. It has been shown that the QHT can predict plasmon resonance, spill-out, and retardation effects in noble and simple metal NPs, matching very well with TD-DFT calculations [59,62]. There are also other works on the development of the QHT that consider the viscous contribution of electron fluid [67,69] and formulation of HT for nonlinear phenomena [70][71][72].
However, it is important to highlight that the QHT results depend on the approximation made for the KE functional (e.g., λ parameter) as well as on the electronic density, which is an input quantity. The input electronic density can be obtained from a preceding OF-DFT calculation using the same KE functional used for the response, i.e., the self-consistent QHT approach [58]. Other approaches can use the exact KS density [59] or, more efficiently, from a model density [57,59].
Although QHT can describe different quantum effects relevant in plasmonics, the QHT is not unaffected by drawbacks: i) Various QHT works [58][59][60][61]66], show the presence of additional resonances above the main plasmon peak and below the plasma frequency (ω p ): these resonances originate from the spatial variation of the electronic density as first pointed out by Bennet [73]. Using the hardwall boundary, no Bennet states are present for NPs in QHT [74][75][76] nor in TF-HT [76][77][78][79]: in these approaches, several peaks (volume plasmons) occur due to nonlocality, but only at frequencies larger than ω p . Such Bennet states are thus peculiar to hydrodynamic models with nonuniform density. Instead, in TD-DFT calculations of large jellium spheres, only a shoulder above the main plasmon peak (around 4 eV) is present [16] and it has been attributed to the interaction between single-particle transitions and surface modes [80][81][82][83]. To better illustrate this point (often overlooked in the literature), we report in Fig. 1 (a) a direct comparison of the absorption spectra of a Na jellium nanosphere between reference TD-DFT and QHT theory. TD-DFT can be considered as a reference for QHT because the latter can be directly derived from TD-DFT equations (and for the two-electron case the methods coincide) [56,67,84]. Although the energy position of the main peak, the localized surface plasmon (LSP), is very well reproduced, additional peaks are not present in the reference TD-DFT spectrum. Results in Fig. 1 (a) represent the current state-of-the-art of QHT calculations: clearly, the presence of the other peaks strongly limit the QHT accuracy and applicability.
ii) For energies above a critical frequency ω c [59,65], the induced density (i.e., the first-order change of the electronic density due to the excitation) have both an oscillating and exponentially decaying behavior, as shown in Fig. 1 (b), which is problematic to treat numerically. We will show in this paper, that all the energy position of all the peaks above the critical frequency strongly changes with the computational domain size so that a for a Na jellium sphere with radius R = 2.167, /nm as obtained from TD-DFT and QHT (using KS density and λ = 1); (b) Induced charge density n1 in atomic units (a.u.) at different energies. The inset in (a) schematically demonstrates the nanosphere interaction with the incident plane wave. See the Appendix A for definitions and details on the absorption spectra calculation.
numerically converging QHT spectrum is challenging to obtain.
iii) The TFλvW functional is known to be quite a rough approximation to the exact KE, and different limitations of this functional have been shown in different contexts, e.g., lack of dynamical corrections [40,55,65,85] and incorrect response for homogeneous electron gas [30,86,87]. Thus, the great accuracy of QHT calculations with the TFλvW functional obtained in some cases should be related to some error cancellation and, therefore, cannot have general validity.
In order to overcome these limitations, in this article, we complement the QHT approach with the recently developed Laplacian-level KE functionals [87], that performed well for semiconductors and metals in the framework of OF-DFT. In particular, the Pauli-Gaussian second-order and Laplacian (PGSL) functional has an improved Lindhard response [87], which is a fundamen-tal property for the description of metallic systems. We derive the Laplacian-level QHT linear-response equations in the frequency-domain: this is a completely novel implementation for the QHT, so far only limited to the TFλvW KE functionals. Note that Laplacianlevel KE functionals are much simpler than fully nonlocal functional based on the Lindhard response in the reciprocal space [30,85,86] and can be easily applied to finite systems [87]. We demonstrate that in the QHT-PGSL approach, only the main plasmon peak appears in the absorption spectrum, which is stable to the changes of computational domain size. In fact, in QHT-PGSL, the induced density decays in the same way for all frequencies, and no critical frequencies exist anymore. Moreover, as we will show later, in QHT-PGSL, the intensity of the main plasmon peak is more accurate than the conventional QHT approach.
The article is organized as follows: in Section II we introduce the equations governing the Laplacian-level QHT theory, which contains also the conventional QHT approach as a special case; in Section III we discuss theoretically the properties of the induced density in the tail region in spherical systems, showing that the QHT-PGSL has an unexpected and completely different behavior with respect to the conventional QHT; in Section IV, we provide numerical details of our implementation which can efficiently describe system with spherical and cylindrical symmetry; in Section V, we compare in details the absorption spectra of Na jellium nanosphere from TD-DFT, QHT and QHT-PGSL, showing their different dependence on the computational domain size as well as their oscillator strength; in Section VI, we describe the numerical results of the induced density decay for Na jellium nanospheres which confirms the theoretical prediction of Section III; in Section VII, we study the energy and oscillator strength as a function of particle size; in Section VIII, we present results for spherical dimers. Finally, conclusion and future perspectives are drawn in Section IX.

II. PAULI-GAUSSIAN AND LAPLACIAN-LEVEL FUNCTIONALS IN QUANTUM HYDRODYNAMICS
The linearized QHT response [28,88] is governed by the following equations [58,59] for the electric field E and polarization vector P: where c is the speed of light, ε 0 and µ 0 are the vacuum permittivity and permeability, m e and e are the electron mass and charge (in absolute value), γ is the phenomenological damping rate and ω p (r) = e 2 n 0 (r) / (m e 0 ) is the plasma frequency with n 0 (r) being the ground-state (equilibrium) electron density. δG [n] δn 1 is the first-order term for the potential associated to the energy functional G [n] given by: where E LDA XC [n] is the exchange-correlation (XC) energy functional in local density approximation (LDA), while T s is the noninteracting KE functional.
• α = 0, β = 0. It corresponds to the case where the QHT is improved with the addition of the PGα functional. We will refer to this case as QHT-PGα.
• α = 0, β = 0. This is the more complex case in which the Laplacian-level correction Lβ is included in the energy functional. This case will be referred to as QHT-PGαLβ.
In order to calculate the potential, we take the functional derivative of T s [n] [89] and obtain: where the subscripts i = n, w, q denote the corresponding partial derivatives. The detailed derivation of Eq. (7) is given in Sec. I of supplementary materials (SM) [90]. The first-order term of the potential, δTs δn 1 , is obtained using a perturbation approach where the perturbed density is taken as n = n 0 +n 1 , with n 1 = 1 e ∇·P being the electron density perturbation. After some tedious algebra and neglecting higher-order terms, we obtain the following expression for the linear potential (we refer to SM Sec. I and II for the full derivation [90]): where: δT TFvW δT PGα where the superscript (0) indicates that the function is evaluated at n = n 0 . First-order term are the contribution associated with the PG and L corrections respectively.
The parameters α and β can be determined in a nonempirical way by imposing exact asymptotic solutions. In particular, we set α = 40/27 in order to satisfy second-order gradient expansion [87,92] and use PGS for PG40/27. Moreover, we follow the results of Ref. 87 and fix β = 0.25 such that the overall correction functional PGSL0.25 accurately reproduces the linear response function of a noninteracting homogeneous electron gas at both small and large wavevectors [87]. For brevity, we will use the acronym PGSL for PGSL0.25.

III. ASYMPTOTIC ANALYSIS
We take the divergence of Eq. (1b), and we use the quasistatic approximation (so that ε 0 ∇·E = ∇·P = en 1 ), obtaining: To obtain the asymptotic form of Eq. (10) we assume that [59]: where κ > 0 is the decay constant of the ground-state density, and νκ is the decay constant of the (dipole excited) induced density. The right-end side (RHS) of Eq. (10) is asymptotically vanishing and it decays as where d 1 is the dipole moment of n 1 (see Ref. 59). For the left-end side (LHS) of Eq. (10) we obtain (after some algerba, see SM, Sec. III [90]): while the terms with E LDA XC and T PGα s decay exponentially faster than n 1 . We underline that Eqs. (14) and (15) represent only the leading terms in the asymptotic region.
The case with β = 0 has been described in Ref. [59]: when ω is higher than a critical excitation energy the asymptotic decay is complex-valued and oscillating. Otherwise, the asymptotic decay is exponential, and ν depends on ω.
When β > 0, we found, interestingly, that the Lβ term gives an exponentially increasing contribution, due to the division by n 2/3 0 , which dominates over the term in Eq. (14) as well as the term at the RHS. Thus, the asymptotic solution does not depend on ω, as in the conventional QHT approach with the TFvW functional, but it is related to the solution of the sixth-degree polynomial in ν in Eq. (15), which are: −1.320, +0.543, +2/3, 1.123, +5/3, +2.987.
Only for those values of ν, the LHS term vanishes asymptotically, as it does the RHS. Some of these solutions are not possible or unstable, i.e., those with ν ≤ 2/3, as the term n 1 /n 2/3 0 will not decay asymptotically. The other three values of ν give the right asymptotic solution, but a high order analytical analysis or a full numerical solution is required to select the actual value of ν. Interestingly, all these solutions have ν > 1, which is another difference with respect to the QHT approach with the TFvW functional [59], where ν < 1.

IV. NUMERICAL IMPLEMENTATION
The system of Eqs. (1) with Eq. (2) and expressions (9) is solved for a plane wave excitation using a commercial implementation of the finite-element method (FEM) [93]. In order to easily compute absorption spectra for spheres and sphere dimers, we have implemented our equations using the 2.5D technique, which significantly reduces the computational time for axisymmetric structures [59,94,95]. A detailed explanation of the FEM implementation can be found in the Appendix B. A completely independent implementation has also been carried out using a finite-difference method for spherical systems in the quasi-static approximation [96]: the results obtained with the two methods are numerically the same. In order to solve the system of Eqs. (1) an ex-pression for ground-state density function n 0 (r) is required. Throughout the article we consider the following two ground-state density functions: i) the exact KS density n KS 0 (r) calculated using a DFT in-house code [59], and ii) a model density defined as [57,59]: normalized with a condition n Mod 0 dV = N e , where N e is the number of electrons. For κ Mod coefficient κ Mod = 1.05/a 0 value is fixed fitted with asymptotic decay of the KS electronic density decay [59]. Fig. 2 shows n KS 0 and n Mod 0 densities for a Na (Wigner-Seitz radius r s = 4 a.u.) jellium nanosphere with N e = 1074 electrons (nanosphere radius R = 2.167 nm). Note that n Mod 0 does not display Friedel oscillations inside the nanosphere volume (surface marked with vertical line), which are instead present in n KS 0 . The inset shows that the asymptotic decay is the same for both cases.

V. ABSORPTION SPECTRA
In Fig. 3 we report the comparison of the normalized absorption cross-section for a Na jellium nanosphere with N e = 1074 electrons as obtained using QHT, QHT-PGS, QHT-PGSL as well as the TD-DFT approaches (see Appendix A for definitions and details).   3 shows that the energy of the localized surface plasmon (LSP) resonance (first main peak) for QHT and QHT-PGS are in good agreement (within 10 meV) with TD-DFT (≈ 3.22 eV), which is broader due to quantumsize effects, while KS/QHT-PGSL and Mod/QHT-PGSL give the LSP peaks at ≈ 3.374 eV and ≈ 3.313 eV, respectively, which are blue-shifted with respect to TD-DFT results (for further analysis of the position of LSP peak, see Sec. VII). As discussed in the Introduction, QHT gives accurate energy of the LSP and predicts additional peaks at higher energies, which are not present in the TD-DFT. Almost the same situation is obtained from QHT-PGS, meaning that even the more general gradient approximation in Eq. (6b) does not solve the problem of additional peaks. On the other hand, the QHT-PGSL absorption spectrum is quite different. The main difference between QHT and QHT-PGSL is not the energy shift of the LSP, but the absence of additional resonances in the latter.
Although the QHT predicts very well the LSP resonance when compared to more sophisticated TD-DFT approaches, the presence of additional peaks is a major shortcoming. These peaks have, in fact, energy higher than ω c ≈ 3.55 eV (see Ref. [59]) and thus can hardly be treated in an efficient numerical scheme. In fact, this is shown in Fig. 4, where QHT normalized absorption cross-section (σ/σ 0 ) for the same jellium nanosphere are calculated for increasing size of the simulation domain. These calculations have been done with a novel finite-difference code for spherical systems [96], which reproduces exactly the FEM results reported in this work but is more accurate in the asymptotic region. Results were obtained with KS (upper panel) and model (lower panel) ground-state densities. Clearly, as the domain size increases, more and more modes appear (and with reduced intensities) in the spectrum without any limit. Thus the absorption spectrum is very sensitive to the domain-size. We note that no previous report in literature has considered the numerical convergence of Bennet states in QHT calculations. With an infinite computational domain size, there should be an infinite number of states with infinitely small peak intensity, i.e., no peaks can be distinguished anymore, and only an unstructured shoulder could be present. This is indeed shown in Fig. 5 where we report the QHT and TD-DFT results for two different computational domain size. We have used a larger broadening for QHT (namely γ=0.2 eV) so that it will give the same intensity at the LSP peak as TD-DFT. While the TD-DFT results are converged with standard domain-size, convergence seems to appear for QHT only with a domain-size of 200 a.u., where no more Bennet peaks can be distinguished and only a shoulder is present. However, this shoulder which starts at ω c is significantly higher (about a factor of 2.5) than the TD-DFT one, which starts later, at about 3.7 eV. Clearly a domain-size of 200 a.u. to obtain a converged absorption spectrum, is not reasonable for any applications in plasmonics, and it has been obtained only with a specialized code for reference calculations [96].
In Fig. 5, we also report the QHT/PGSL* results, where the * indicates that the spectra have been redshifted by 0.15 eV in order to have the same LSP energy position of QHT; the broadening is γ = 0.24 eV so that also the peak intensity is the same. The plot shows that the QHT-PGSL* does not change at all with the computational domain-size, and overall it is much closer to TD-DFT than QHT.
A more quantitative comparison of methods can be done by considering the integrated absorption crosssection: which will converge to (πe 2 )/(2 0 m e c)N e for ω → ∞, where N e is the number of electrons [81,97,98]. The integrated absorption is plotted in Fig. 6, and it shows that I(ω) for QHT and QHT-PGSL* converge to the same value for high energies. However, while the in- tegrated absorption curve for TD-DFT and QHT-PGSL* are very close to each other, the growth in QHT is much slower, meaning that the oscillator strength (i.e., the energy-integrated intensity) in QHT is split in several Bennet modes, whereas the single peak in QHT-PGSL* contains it all. In fact, the integrated absorption for QHT at ω = ω c is about 15% smaller than QHT-PGSL* and TD-DFT. In Sec. VII, a more detailed analysis of oscillator strength and absorption cross-section for different number of electrons is presented. Here, we remark that these features are not limited to spherical NPs, but could happen in other geometries or materials: in fact, for ω c the identical expression was obtained for a jellium sphere [59] and slab [65]. Thus, in general, one could have for LSP ω lsp ω c or even ω lsp > ω c . In such cases, the QHT will not accurately predict the LSP energy.

VI. INDUCED CHARGE DENSITY
As discussed in the Sec. III, the decay of conventional QHT induced densities is frequency-dependent, and solutions are pure exponentially decaying at the metal surface only if the incident plane wave energy is lower than ω c , whereas using the PGSL functional a fixed exponential decay should be obtained.
This can be verified numerically by plotting the computed induced charge density n 1 (associated with the absorption). In Fig. 7 we plot |n 1 | (in the logarithmic scale) as obtained from the KS/QHT and KS/QHT-PGSL for a Na jellium nanosphere with N e = 1074. To have clear comparison of decay rates, the curves for |n 1 | are shifted to have the maximum at z = R and normalized to |n 1 (R)|, while n 0 density is only normalized to n 0 (R).
For the KS/QHT induced density the decay slope shows a clear dependence on the incident energy ω, becoming oscillatory for ω > ω c = 3.55 eV (such This is indeed shown in Fig. 5 where we report the QHT and TD-DFT results for two different computational domain size. We have used a larger broadening for QHT (namely γ=0.2 eV) so that it will give the same intensity at the LSP peak as TD-DFT. While the TD-DFT results are converged with standard domain-size, convergence seems to appear for QHT only with a domain-size of 200 a.u., where no more Bennet peaks can be distinguished and only a shoulder is present. However, this shoulder which starts at ω c is significantly higher (about a factor of 2.5) than the TD-DFT one, which starts later, at about 3.7 eV. Clearly a domain-size of 200 a.u. to obtain a converged absorption spectrum, is not reasonable for any applications in plasmonics, and it has been obtained only with a specialized code for reference calculations [96].
In Fig. 5, we also report the QHT/PGSL* results, where the * indicates that the spectra have been redshifted by 0.15 eV in order to have the same LSP energy position of QHT; the broadening is γ = 0.24 eV so that also the peak intensity is the same. The plot shows that the QHT-PGSL* does not change at all with the computational domain-size, and overall it is much closer to TD-DFT than QHT.
A more quantitative comparison of methods can be done by considering the integrated absorption crosssection: which will converge to (πe 2 )/(2 0 m e c)N e for ω → ∞, where N e is the number of electrons [81,97,98]. The integrated absorption is plotted in Fig. 6, and it shows that I(ω) for QHT and QHT-PGSL* converge to the same value for high energies. However, while the in- This is indeed shown in Fig. 5 where we report the QHT and TD-DFT results for two different computational domain size. We have used a larger broadening for QHT (namely γ=0.2 eV) so that it will give the same intensity at the LSP peak as TD-DFT. While the TD-DFT results are converged with standard domain-size, convergence seems to appear for QHT only with a domain-size of 200 a.u., where no more Bennet peaks can be distinguished and only a shoulder is present. However, this shoulder which starts at ω c is significantly higher (about a factor of 2.5) than the TD-DFT one, which starts later, at about 3.7 eV. Clearly a domain-size of 200 a.u. to obtain a converged absorption spectrum, is not reasonable for any applications in plasmonics, and it has been obtained only with a specialized code for reference calculations [96].
In Fig. 5, we also report the QHT/PGSL* results, where the * indicates that the spectra have been redshifted by 0.15 eV in order to have the same LSP energy position of QHT; the broadening is γ = 0.24 eV so that also the peak intensity is the same. The plot shows that the QHT-PGSL* does not change at all with the computational domain-size, and overall it is much closer to TD-DFT than QHT.
A more quantitative comparison of methods can be done by considering the integrated absorption crosssection: which will converge to (πe 2 )/(2 0 m e c)N e for ω → ∞, where N e is the number of electrons [81,97,98]. The integrated absorption is plotted in Fig. 6, and it shows that I(ω) for QHT and QHT-PGSL* converge to the same value for high energies. However, while the in- tegrated absorption curve for TD-DFT and QHT-PGSL* are very close to each other, the growth in QHT is much slower, meaning that the oscillator strength (i.e., the energy-integrated intensity) in QHT is split in several Bennet modes, whereas the single peak in QHT-PGSL* contains it all. In fact, the integrated absorption for QHT at ω = ω c is about 15% smaller than QHT-PGSL* and TD-DFT. In Sec. VII, a more detailed analysis of oscillator strength and absorption cross-section for different number of electrons is presented. Here, we remark that these features are not limited to spherical NPs, but could happen in other geometries or materials: in fact, for ω c the identical expression was obtained for a jellium sphere [59] and slab [65]. Thus, in general, one could have for LSP ω lsp ω c or even ω lsp > ω c . In such cases, the QHT will not accurately predict the LSP energy.

VI. INDUCED CHARGE DENSITY
As discussed in the Sec. III, the decay of conventional QHT induced densities is frequency-dependent, and solutions are pure exponentially decaying at the metal surface only if the incident plane wave energy is lower than ω c , whereas using the PGSL functional a fixed exponential decay should be obtained.
This can be verified numerically by plotting the computed induced charge density n 1 (associated with the absorption). In Fig. 7 we plot |n 1 | (in the logarithmic scale) as obtained from the KS/QHT and KS/QHT-PGSL for a Na jellium nanosphere with N e = 1074. To have clear comparison of decay rates, the curves for |n 1 | are shifted to have the maximum at z = R and normalized to |n 1 (R)|, while n 0 density is only normalized to n 0 (R).
For the KS/QHT induced density the decay slope shows a clear dependence on the incident energy ω, becoming oscillatory for ω > ω c = 3.55 eV (such tegrated absorption curve for TD-DFT and QHT-PGSL* are very close to each other, the growth in QHT is much slower, meaning that the oscillator strength (i.e., the energy-integrated intensity) in QHT is split in several Bennet modes, whereas the single peak in QHT-PGSL* contains it all. In fact, the integrated absorption for QHT at ω = ω c is about 15% smaller than QHT-PGSL* and TD-DFT.
In Sec. VII, a more detailed analysis of oscillator strength and absorption cross-section for different number of electrons is presented. Here, we remark that these features are not limited to spherical NPs, but could happen in other geometries or materials: in fact, for ω c the identical expression was obtained for a jellium sphere [59] and slab [65]. Thus, in general, one could have for LSP ω lsp ω c or even ω lsp > ω c . In such cases, the QHT will not accurately predict the LSP energy.

VI. INDUCED CHARGE DENSITY
As discussed in the Sec. III, the decay of conventional QHT induced densities is frequency-dependent, and solutions are pure exponentially decaying at the metal surface only if the incident plane wave energy is lower than ω c , whereas using the PGSL functional a fixed exponential decay should be obtained.
This can be verified numerically by plotting the computed induced charge density n 1 (associated with the absorption). In Fig. 7 we plot |n 1 | (in the logarithmic scale) as obtained from the KS/QHT and KS/QHT-PGSL for a Na jellium nanosphere with N e = 1074. To have clear comparison of decay rates, the curves for |n 1 | are shifted to have the maximum at z = R and normalized to |n 1 (R)|, while n 0 density is only normalized to n 0 (R).
For the KS/QHT induced density the decay slope shows a clear dependence on the incident energy ω, becoming oscillatory for ω > ω c = 3.55 eV (such KS/QHT induced densities are not converged with respect to the computational domain size, as discussed in Sec. V).
On the other hand, the KS/QHT-PGSL calculations yield the same slope for all excitation energies, as analytically demonstrated in Sec. III. A numerical fit of the decay gives a value of ν close to +1.12, i.e., the slowest from asymptotically decaying solutions (with ν > 2/3) (17).
It is important to note that the TD-DFT calculations ( Fig. 7 (c)) give qualitatively similar results as the QHT-PGSL. In fact, for TD-DFT we get the same decay slope for the induced density (at least for ω < ω p ). However, as discussed in Sec. III in QHT-PGSL we have ν > 1 while ν < 1 in the conventional QHT, meaning that spill-out effects are somehow smaller in QHT-PSGL. Nonetheless, we need to point out that this feature is peculiar to PGSL, which is one of the few Laplacian-level KE functionals, and PGSL has not been developed for QHT calculations. Thus, other Laplacian-level KE functional can be developed with better features. Another important aspect is the numerical stability of the QHT-PGSL approach: not only the absorption spectra does not depend on the domain-size, but the fact that the decay constant is fixed and independent from the frequency allows to use mixed boundary condition for an exponential decay (i.e.,r · ∇n 1 + νκn 1 = 0), allowing to have converged results even with a very small domainsize.
The impact of the Laplacian-level corrections studied in this article is not limited to the density tail decay but can be easily appreciated in the linear scale as well. In Fig. 8, we report the imaginary part of the induced charge density n 1 , Figs. 8 (a) -(d), and norm of the total (incident + scattered) field for the KS ground-state density (for the model density see Fig. S1 of the SM [90]). The cross-section plane for the map plots in Figs. 8 (b) -(g) is perpendicular to the incident plane wave propagation. We clearly observe that the spill-out is less pronounced for QHT-PGSL, as the imaginary part of n 1 is more squeezed into the nanosphere volume. Consequently, for QHT-PGSL the peak value of Im (n 1 ) and field enhancement are the biggest. This behavior is directly related to the blue shift of LSP resonance obtained using QHT-PGSL with respect to QHT and QHT-PGS shown in Fig. 3.

VII. ELECTRON NUMBER VARIATION EFFECTS
An important aspect in nanoplasmonic systems is the LSP resonance dependence on the NP size [47,48,99]. As already mentioned in the introduction, LSP for Na jellium nanosphere must undergo red shift with respect to Mie theory results. In Fig. 9 (a) (horizontal axis is in the logarithmic scale), we show the LSP resonance energy of various Na jellium nanospheres with the number of electrons N e varying from 254 to 6174 (R = r s N  has some better properties than the TFvW functional [87], yield similar but somehow smaller accuracy. Finally, PGSL overestimates the LSP peak, but an overestimation of 85 meV is quite small, considering that that choice of the exchange-correlation functional can shift the results even more [100][101][102][103]. To describe the accuracy of a given theoretical method for the calculation of the absorption spectra, not only the energy of the LSP peak have to be considered, but also the oscillator strength, f osc , associated to it. The oscillator-strength is readily available in an eigenvalue formulation of QHT [57,104]. Our QHT implementation is frequency-dependent, and, therefore, the oscillator strength is not directly computed, but it can be extracted from the absorption spectra using the fit- ting procedure described in Sec. IV of the SM [90]. The oscillator strength of the LSP peak can also be extracted from the TD-DFT spectra, if the onset of the plasmon shoulder is considered (Sec. IV of the SM [90]). Previous attempts to compute the f osc of the LSP peak are based on the sum-rule approaches [83]. In Fig. 9 (b), we report f osc of the LSP peak, as obtained from TD-DFT, Mod/QHT, Mod/QHT-PGS, and Mod/QHT-PGSL. Fig. 9 (b) shows that for all methods LSP converges to the Mie results for large N e . However, f osc for Mod/QHT and Mod/QHT-PGS is largely underestimated, as the main-plasmon peak is subdivided in different peaks, as previously discussed. On the other hand, the main-peak of Mod/QHT-PGSL contains almost all the oscillator strength. The mean absolute relative error for these systems are 6.7%, 6.4%, and 2.9% for Mod/QHT, Mod/PGS, and Mod/PGSL. Thus PGSL overestimates the energies with respect to QHT but gives a better description of the oscillator strength.
Results using the KS-density (not reported) are similar, but in this case the KS/QHT-PGSL overestimates the TD-DFT even more, as also shown in Fig. 3. This can be traced back to the higher oscillating behavior of the KS density inside the NP, as shown in Fig. 1. Such quantum oscillations induce higher values of the Laplacian and thus higher contributions to the energy. On the other hand, with the model density, both the gradient and the Laplacian are very small inside the NP. An 'exact' KE functional should be able to describe both situations, but this is not the case of the PGSL functional, which has not been optimized for jellium nanosphere nor for the QHT approach.
In any case, for applications involving large systems, only the model density can be used since the calculation of n KS 0 via KS-DFT or OF-DFT would require additional computational effort.

VIII. APPLICATION TO SPHERICAL DIMER
Our FEM implementation allows to calculate absorption spectra for axisymmetric structures. An important example of such a system is a nanosphere dimer. The NP dimer has been widely studied in literature since it supports gap plasmons that can squeeze light down to sub-nanometer volumes, making it an ideal system for exploring quantum and nonlocal phenomena [21,67,[105][106][107]. Here we consider a dimer of Na jellium spheres with 1074 electrons each. In Fig. 10 (a) we present the comparison of the absorption cross-section as calculated from the Mod/QHT, Mod/QHT-PGS and Mod/QHT-PGSL (the cross-section is normalized to the 2σ 0 = πR 2 with R being the radius of a single sphere). The plane wave that excites the structure is polarized along the z-axis, and the input ground-state density is a sum of model densities (18) of two spheres. As we can see, the Mod/QHT and Mod/QHT-PGS give oscillations in the spectrum, which are absent in Mod/QHT-PGSL approach. Our convergence analysis showed that these oscillations, as in the case of the sphere (see Fig. 3), persist with the change of the computational domain. These oscillations should not be confused with the small undulation next to the main plasmon peak that is more clearly visible in gap = 1 nm case ( Fig. 10 (a)). This undulation  comes from the charge-transfer plasmon and gets higher for smaller sizes of the gap [21]. For all considered cases of gap size, Mod/QHT-PGSL gives blueshifted plasmon resonance energy with respect to other methods. The respective values of the plasmon resonance are shown in the map plots of the total field enhancement in Fig. 10.
There we also see that the field gets more enhanced for Mod/QHT-PGSL as for other approaches. Clearly, as Table I shows, more field is concentrated in the gap for Mod/QHT-PGSL. In turn, this means that more energy is moved to the main plasmon peak, hence, it is blue shifted for Mod/QHT-PGSL with respect to Mod/QHT and Mod/PGSL as shown in absorption curves.

IX. CONCLUSIONS AND FUTURE PERSPECTIVES
We have extended the quantum hydrodynamic theory to Laplacian-level kinetic energy functionals. In particular, we studied the PGSL functional, which is accurate for OF-DFT calculations and well reproduced the linear response of the homogeneous electron gas [87]. We investigated in detail Na jellium nanospheres and results are compared to reference TD-DFT calculations. The obtained results are focused around two main findings: 1. QHT and QHT-PGS that are gradient-level functionals of electron density, together with an LSP resonance give additional resonances in the absorption spectrum of Na jellium nanospheres. Well defined additional resonances are not present in TD-DFT nor in QHT with an infinite computational domain-size. In both cases, only a shoulder is present at the high energy side to the plasmon peak, with the TD-DFT result being much smaller and at higher energy than in QHT. On the other hand, complementing QHT with a Laplacian-level functional (QHT-PGSL) of electron density yields only the LSP peak in the absorption spectrum, yielding an overall spectrum and oscillator strength closer to TD-DFT.
2. The theoretical and numerical asymptotic analysis of the induced charge density as obtained from QHT, shows that the decay slope is changing at different energies of incident radiation. Contrarily, QHT-PGSL shows the same decay slope for all energies. Despite the spill-out effects are somehow underestimated in QHT-PGSL, this strongly simplifies the boundary condition, so that converged calculation can be obtained with a very small domainsize.
Our results thus demonstrate that the convergence of the QHT absorption spectra is problematic, and most of the QHT results reported so far are thus not accurate enough for energies above LSP resonance. The Laplacian-level QHT, on the other hand, does not suffer from these problems and can be successfully applied to investigate quantum and spill-out effects in different sys- tems, as shown in Sec. IX, where we report the results for NP dimers.
From a computational point-of-view, QHT-PSGL is close to the widely used TF-HT with hard-wall: in fact, the size of the computational domain can be set very close to the size of the nanoparticles, as discussed in Sec. VI. However, QHT-PGSL solves the two main drawbacks of the TF-HF approach: the neglect of spill-out effects and the blueshift of LSP energies.
QHT-PGSL thus allows an efficient and numerically converged computation of collective excitations in quantum systems: the main drawbacks of QHT-PGSL are related to the too fast decay of the induced density, which then causes reduced spill-out effects and increased plasmon energy. However, we point out that the QHT-PGSL functional, used in this work, has not been developed for QHT, but for bulk properties of bulk metal and semiconductors. We believe the Laplacian-level QHT results can be improved considered more sophisticated functional that could depend, for example, on a different function of q and/or on a product of w and q.
The Laplacian-level QHT is thus a new platform, very promising for the future, as the Laplacian ingredient includes much more degree of freedom in developing accurate KE functionals, than more conventional functionalbased on density-gradient. So far, however, the development of semilocal KE functional focused only on groundstate properties, considering only the total KE energy and the KE potential (i.e., the first functional derivatives). Instead, for the QHT response properties, the KE kernel (i.e., the second functional derivative) is required, but, so far, it has not been considered at all in the semilocal KE functional development [87,[108][109][110][111][112][113]. In addition, it is crucial to understand the role of static and dynamic corrections to the energy functional. Although, here, we have considered only static corrections at the second-order gradient and Laplacian-level, the analysis of dynamic correction represents another important route to explore. Overall, we believe that our current results will help to better understand the role of functional dependence on electron density in plasmonic systems.

Appendix A: Absorption spectrum
In QHT, the absorption cross-section is calculated as with I 0 being the intensity for the incident plane wave with frequency ω. The electric field E and the polarization vector P are obtained solving Eqs. (1a, 1b). Considering the very small size of the of the investigated nanoparticles, only dipole modes are excited. An important parameter for the shape of the absorption spectra is the damping parameter (γ), see Eq. 1b. If not stated differently, in all QHT calculations, we use γ = 66 meV. The normalized absorption cross-section (absorption efficiency) is then obtained by normalizing σ to the geometric cross-section of a nanosphere σ 0 = πR 2 with R being the radius of the nanosphere.
The TD-DFT absorption spectra have been computed with a finite-difference in-house code (with spherical symmetry) introduced in Ref. 59; a radial uniform grid is used to represent KS orbitals and densities. In TD-DFT, no retardation effects are included, and only longitudinal electric field are considered [19]. The absorption crosssection [114][115][116][117], is calculated as where the polarizability is given by with χ (r, r , ω) = δn (r) /δ (eV ext (r )) being the inter-acting density response function [19], which is obtained solving the Dyson equation In Eq. (A4) v coul is the Coulomb interaction, f LDA XC is the adiabatic LDA XC kernel, and χ 0 is the noninteracting density response function, which is computed using Green's function [115] using occupied KS orbitals from the ground-state calculation (again using LDA). The broadening parameter for the Green's function calculations is, if not stated differently, Γ 0 = 33 meV.

Appendix B: FEM implementation
In order to lower the order of derivatives we multiply Eq. (1b) by test functionP and integrate by parts, which give us: where we assumed that the integral on the boundary goes to zero. Even after integration by parts δG[n] δn 1 potential contains derivatives up to the fourth-order of n 1 (see the Exps. (9)), so auxiliary variables should be added to lower the order of differentiation. By introducing two variables F = ∇n 1 and O = ∇ ∇ 2 n 1 = ∇ · F we have only first-order derivatives. Considering axisymmetry of considered structures, we adopt 2.5D technique [59,94,95], and the dependence of E, P, F, and O on az-imuthal coordinate is taken in e −imφ form with m ∈ Z. The dependence on m for test functionsẼ,P,F andÕ is of e imφ form. Thus, instead of a three-dimensional problem, we can have 2m max + 1 problems (with m max being the maximum value for m). Moreover, for the dimensions considered in the current work m max = 0 is enough for the convergence of results. Finally, only one two-dimensional problem needs to be solved. Hence, we come to the following system of equations: where (0) superscript denotes the zero-order coefficients of the v (ρ, φ, z) = m∈Z v (m) (ρ, z) e −imφ vector field expansion of cylindrical harmonics. We found that curl elements for Eq. (B2a) and divergence elements [93] for other equations of the system (B2) are the best choice for stable solutions.
For the wave equation Eq. (B2a), simulation domain radius R dom is defined via R dom = R + 500a 0 condition. R = r s N 1/3 e is the radius of the nanosphere, and, for dimers, it is the radius of one of the spheres. Perfectly matched layers (PML) are used in order to emulate an infinite domain and avoid unwanted reflections. PML thickness is set to t PML = 200a 0 for all the considered systems. Also, zero flux boundary condition is imposed on the electric field at the outer boundary of the PML. For Eqs. (B2b) -(B2d), simulations are done in a smaller domain, considering faster decay of variables P, F, and O compared to the electric field. The domain, as depicted in Fig. S2 (a) of SM [90], is a semicircle (consider the axial symmetry), for the nanospheres and, for the dimers, is the union of 2 circles centered at the centers of nanospheres. Moreover, to facilitate the calculations, only the "quarter" of the dimer is simulated with a corresponding perfect electric conductor condition at the intersection segment of 2 circles, as shown in Fig. S2 (b) of SM [90]. The radii for the circles is r dom = r + 25a 0 for QHT and QHT-PGS, but for QHT-PGSL, it is in the range r dom ≈ r + 12a 0 . The simulation domain is smaller for QHT-PGSL, because ν ≈ 1.12 decay slope is bigger in this case (see Sec. VI). Dirichlet boundary conditions P = 0, F = 0 and O = 0 are set on the simulation domain boundary. As stated in Section VI, mixed boundary condition with a fixed exponential decay can be used for QHT-PGSL, so that a very small simulation domain is enough for converged calculations.

S4. FITTING PROCEDURE TO EXTRACT THE OSCILLATOR STRENGTH
We fit the photoabsorption cross-section (of all QHT approaches and reference TD-DFT) with the function where γ is the broadening parameter and ω 0 the peak position. Note that f (ω) has a maximum at ω = ω 0 and f (ω 0 ) = 2/(πγ). This function integrates to 1 in the range 0 < ω < +∞ and originates from the Drude classical solution of a nanosphere, whose photoabsorption cross-section in CGS units is [2]: where we used that w 0 = w p / √ 3, ω p = 4πe 2 /m e (N/V ) and V = (4π/3)R 3 . For a jellium nanosphere of very large R, QHT and TD-DFT will reach the classical limit [3], thus the function f (ω) is the right function to be used in the fitting procedure, where we minimize: where σ m (ω) is the input absorption spectra (QHT or TDDFT). The error depends on three parameters f osc , ω 0 , and γ. For the range, we consider ω a = 1 eV whereas ω b is the energy position where other peaks/shoulders starts after the LSP resonance. In this way, only the LSP peak is included in the fit, whereas other high-energy peaks or shoulder are not considered. As an example, we show in Fig. S3 the QHT spectra for N e = 1074. The resulting fit is very accurate for all energies less than ω b . Other peaks above ω b are not included in the peaks, so that the oscillator strength of the first main peak can be obtained. For QHT, one can also do a fit with several f (ω), one for each peak, obtaining the same results. The situation complicates for the TD-DFT spectrum. In Fig. S4 a), we report the TD-DFT spectrum for N e = 6174 (i.e., the largest system considered) and the corresponding fit (green dashed-line) with Eq. (S32). The error in the fit is 2.1% in the fitting range and fitted f osc is 890. The fitting curve is not very accurate as the TD-DFT values (for low and high energies) are not well reproduced. In order to improve the fit, we considered an energy-dependent broadening [4], i.e., γ = g(ω) and to model the function g(ω), we inverted Eq. (S32) as a function on ω. For an arbitrary spectrum σ(ω) with f osc = ∞ 0 σ(ω)dω we define the function: It can be easily verified that if σ(ω) = f osc f (ω), then we obtain a constant function Γ(ω) = γ. In Fig. S4 b), we report the function Γ(ω) for the TD-DFT absorption spectra using Eq. (S35). There are some spikes at ω = ω 0 due to the presence vanishing first derivative at the main plasmon peak. For high-and low-energy, we see that Γ(ω) approaches exactly 2Γ 0 , where Γ 0 is the broadening used in TD-DFT calculation. Thus, we define a model g(ω) with the following limits g(ω ≤ ω low ) = 2Γ 0 , where ω low is close to ω a and g(ω > ω p ) = 2Γ 0 , where ω p is the plasma frequency. In between these limits, the function g(ω) grows up to the maximum values of g(ω 1 ) = γ, where ω 1 is fixed to be ω b i.e. the energy where the shoulder in the spectrum starts. If the TD-DFT absorption spectrum is fitted with such function, we obtain the red-curve in Fig. S4a) with an error reduced by a factor of 4 (0.6%) and also a reduced oscillator strength ( f osc = 872). The non-linear fit is done using four parameters: f osc , ω 0 , γ, and ω low . Note also that the function f (ω) has to be renormalized to integrate to 1 when a frequency dependent broadening is used. In Fig. S4 b) we report the function Γ(ω) obtained from the fitted curve. The agreement with the TD-DFT results is very good, for all energies before and after the plasmon shoulder.
In Table S1, we report all the parameters for the fits of the TD-DFT spectra. Note that for the fits, we compute the TD-DFT with a broadening of Γ 0 = 0.1 eV, to have a smooth curve, which is required especially for small spheres. The f osc of the LSP peak is thus very accurate (about 1%) for systems where the plasmon shoulder is clearly observed (i.e., for N e > 1000).  In order to improve the fit, we considered an energy-dependent broadening [4], i.e., γ = g(ω) and to model the function g(ω), we inverted Eq. (S32) as a function on ω. For an arbitrary spectrum σ(ω) with f osc = ∞ 0 σ(ω)dω we define the function: It can be easily verified that if σ(ω) = f osc f (ω), then we obtain a constant function Γ(ω) = γ. In Fig. S4 b), we report the function Γ(ω) for the TD-DFT absorption spectra using Eq. (S35). There are some spikes at ω = ω 0 due to the presence vanishing first derivative at the main plasmon peak. For high-and low-energy, we see that Γ(ω) approaches exactly 2Γ 0 , where Γ 0 is the broadening used in TD-DFT calculation. Thus, we define a model g(ω) with the following limits g(ω ≤ ω low ) = 2Γ 0 , where ω low is close to ω a and g(ω > ω p ) = 2Γ 0 , where ω p is the plasma frequency. In between these limits, the function g(ω) grows up to the maximum values of g(ω 1 ) = γ, where ω 1 is fixed to be ω b i.e. the energy where the shoulder in the spectrum starts. If the TD-DFT absorption spectrum is fitted with such function, we obtain the red-curve in Fig. S4a) with an error reduced by a factor of 4 (0.6%) and also a reduced oscillator strength ( f osc = 872). The non-linear fit is done using four parameters: f osc , ω 0 , γ, and ω low . Note also that the function f (ω) has to be renormalized to integrate to 1 when a frequency dependent broadening is used. In Fig. S4 b) we report the function Γ(ω) obtained from the fitted curve. The agreement with the TD-DFT results is very good, for all energies before and after the plasmon shoulder.
In Table S1, we report all the parameters for the fits of the TD-DFT spectra. Note that for the fits, we compute the TD-DFT with a broadening of Γ 0 = 0.1 eV, to have a smooth curve, which is required especially for small spheres. The f osc of the LSP peak is thus very accurate (about 1%) for systems where the plasmon shoulder is clearly observed (i.e., for N e > 1000).