Accurate and Scalable Quantum Hydrodynamic Simulations of Plasmonic Nanostructures Within OFDFT

1. Introduction

Plasmonic nanostructure systems are at the forefront of nanoscience research and have found applications across a wide range of fields [1,2,3,4,5,6,7], including enhanced solar energy conversion [8,9], biomedical applications [10,11,12], surface-enhanced Raman scattering [13,14], plasmon rulers with ultrahigh sensitivity [15], plasmonic photocatalysis [16], plasmonic nanoantennas [17], sensors [18], plasmon lasers [19], and nano-optical tweezers [20], among others. As the size of nanodevices continues to shrink, quantum effects become increasingly significant and cannot be neglected [21,22,23,24]. The classical Drude model, based on the local response approximation, is no longer sufficient to accurately describe these systems, making the development of new simulation methods an urgent need.

Quantum hydrodynamic theory (QHT) offers a promising solution [25,26,27,28,29,30]. It provides results consistent with time-dependent density functional theory (TDDFT) for simple metals while requiring significantly fewer computational resources and supporting large-scale nanostructure modeling [25,26]. However, two major challenges must be addressed for its practical application. The first concerns the determination of the ground-state electron density, which serves as the critical input for optical response calculations. While density functional theory (DFT) or analytical models are typically used to compute ground-state density [25,26,27] , DFT becomes computationally prohibitive for large structures, and analytical models are limited to simple geometries such as spheres or infinite plates.To overcome these limitations, we choose using orbital-free (OF) DFT to calculate the ground-state density [28,29,30,31]. OFDFT can be viewed as a simplified, self-consistent Kohn–Sham-like equation considering only the lowest eigenvalue, making it suitable for implementation within the finite element method and enabling integration with QHT for efficient numerical calculations.

The second challenge lies in the accuracy and numerical stability of QHT, which depend on the form and parameters of the energy functional. A widely adopted approach is to use the Thomas Fermi–von Weizsäcker (TFvW) kinetic energy functional in combination with local density approximation (LDA) exchange–correlation (XC) terms [25,27,28,30]. However, TFvW often leads to non-convergent results, manifesting as spurious high-frequency peaks sensitive to the simulation domain size.To address this, a modified Laplacian-level kinetic energy functional, known as PGSLN, has been introduced and shown to significantly improve stability and accuracy when combined with model ground-state density [26]. To further extend its applicability, we propose the OF-PGSLN method, which integrates OFDFT-derived ground-state density with the PGSLN functional.

In this study, we first calibrate the energy functional form and parameters of OF-PGSLN by comparing the ground-state density against DFT and the optical response against TDDFT. We then demonstrate the capability of OF-PGSLN to model complex nanostructures by examining the optical response of nanodimers. In particular, for small interparticle gaps, OF-PGSLN captures significant differences in both the ground-state density and optical properties compared to the simple linear superposition model, which neglects interparticle interactions.In summary, we develop a QHT-based simulation framework that is both accurate and stable for arbitrary nanostructure modeling. The remaining part of this paper is organized as follows: Section 2 presents the theoretical framework; Section 3 details the DFT-based calibration of the ground-state energy functional; Section 4 calibrates the optical response functional using TDDFT and verifies generality using nanospheres of various sizes; Section 5 applies OF-PGSLN to nanodimers with different gap sizes and compares with the model ground-state approach; and Section 6 concludes the work.

2. Model and Method

The conventional linearized QHT response is governed by the following equations in the frequency domain [25,26,27,28,29,30]:

$$\nabla \times \nabla \times {\mathbf{E}}_{\mathrm{s}}-{\displaystyle \frac{{\omega }^2}{c^2}}{\mathbf{E}}_{\mathrm{s}}={\omega }^2{\mu }_0\mathbf{P},$$

(1a)

$${\displaystyle \frac{e{n}_0}{m_e}}\nabla {\left({\displaystyle \frac{\delta G\left[n\right]}{\delta n}}\right)}_1+({\omega }^2+i\gamma \omega )\mathbf{P}=-{\epsilon }_0{\omega }_p^2({\mathbf{E}}_{\mathrm{i}}+{\mathbf{E}}_{\mathrm{s}}).$$

(1b)

Here, ${\mathbf{E}}_{\mathrm{s}}$ is the scattered field for the free electron and ${\mathbf{E}}_{\mathrm{i}}$is the incident field. In this work, plane waves incident perpendicular to the z-axis are used. c is the speed of light in vacuum. ${\epsilon }_0$ and ${\mu }_0$ are the vacuum permittivity and permeability. $m_e$ and e are the electron mass and charge (in absolute value). $\gamma$ represents the phenomenological damping rate. In this paper, we choose $0.22\mathrm{eV}$, which can provide corresponding oscillator strength in line with TDDFT [26]. ${\omega }_p\left(\mathbf{r}\right)=\sqrt{e^2{n}_0/m_e{\epsilon }_0}$ is the plasma frequency with $n_0$ being the spatially dependent ground-state (equilibrium) electron density. ${\left(\delta G\left[n\right]/\delta n\right)}_1$ can be obtained by adopting a first-order perturbation approach where the perturbed density is taken as $n\left(r\right)=n_0\left(r\right)+n_1\left(r\right)$, with $n_1\left(r\right)={\displaystyle \frac{1}{e}}\nabla ·\mathbf{P}$ being a small (by assumption) first-order dynamic perturbation. Where the energy functional $G\left[n\right]$ is given by

$$G\left[n\right]=T_{\mathrm{s}}\left[n\right]+E_{\mathrm{XC}}\left[n\right].$$

(2)

Here, $E_{\mathrm{XC}}\left[n\right]$ is the XC functional, while $T_{\mathrm{s}}\left[n\right]$ is the noninteracting kinetic energy functional. Most recently [26], the noninteracting kinetic energy functional depending on the Laplacian of the electronic density has been proposed, which has the form

$$T_{\mathrm{s}}\left[n\right]=\int \tau (n,w,q)d^{3}r,$$

(3)

with $w={\left|\nabla n\right|}^2$ and $q={\nabla }^{2}n$. The kinetic energy density $\tau \left(n,w,q\right)$ is approximated as the sum of the von Weizsäcker (vW) ${\tau }^{\mathrm{vW}}\left(n,w\right)$, the Pauli-Gaussian (PG$\alpha$) ${\tau }^{\mathrm{PG}\alpha }\left(n,w\right)$, the Laplacian (L$\beta$) ${\tau }^{\mathrm{L}\beta }\left(n,q\right)$, and a modified term ${\tau }^{\mathrm{N}{q}_0}\left(n,q\right)$:

$$\tau (n,w,q)={\tau }^{\mathrm{vW}}(n,w)+{\tau }^{\mathrm{PG}\alpha }(n,w)+{\tau }^{\mathrm{L}\beta }(n,q)+{\tau }^{\mathrm{N}{q}_0}(n,q),$$

(4)

where

$${\tau }^{\mathrm{vW}}\left(n,w\right)=A{n}^{-1}w,$$

(5a)

$${\tau }^{\mathrm{PG}\alpha }\left(n,w\right)={\tau }^{\mathrm{TF}}\left(n\right)e^{-\alpha C{n}^{-8/3}w},$$

(5b)

$${\tau }^{\mathrm{L}\beta }\left(n,q\right)=\beta D{n}^{-5/3}{q}^2=\beta {\tau }^{\mathrm{TF}}\left(n\right)q_r^2,$$

(5c)

$${\tau }^{\mathrm{N}{q}_0}\left(n,q\right)=E{n}^{5/3}ln\left(1+q_r/q_0\right)=2{\tau }^{\mathrm{TF}}\beta q_0^{2}ln\left(1+q_r/q_0\right).$$

(5d)

In the above equations, ${\tau }^{\mathrm{TF}}\left(n\right)=3\left(E_h{a}_0^2\right){\left(3{\pi }^2\right)}^{2/3}{n}^{5/3}/10\equiv B{n}^{5/3}$ is the TF kinetic energy functional with the Hartree energ $E_h={\hslash }^2/\left(m_e{a}_0^2\right)$ and the coefficient $B=3\left(E_h{a}_0^2\right){\left(3{\pi }^2\right)}^{2/3}/10$. The other coefficients are $A={\lambda }_{\mathrm{W}}{E}_h{a}_0^2/8$, $C={\left(3{\pi }^2\right)}^{-2/3}/4$, $D=3{\left(3{\pi }^2\right)}^{-2/3}{E}_h{a}_0^2/160$, and $E=6\left(E_h{a}_0^2\right){\left(3{\pi }^2\right)}^{2/3}\beta q_0^2/10$. $q_r=3q/\left(40{\tau }^{\mathrm{TF}}\right)$ is the reduced Laplacian, which is very large for small density. It is emphasized that the Vw term near the metal surface plays a major role, and the parameter ${\lambda }_{\mathrm{W}}$ controls the magnitude of the spill out effect. The parameters $\alpha$ and $\beta$ are $\alpha =40/27$ and $\beta =0.25$. $q_0$ is a parameter, which can be used to tuning the energy position of the Bennett state. It is $q_0=700$ in order to have the same Bennett peak position as that obtained from the TDDFT calculations for a Na jellium nanosphere with 1074 electrons when using DFT dengsity [26].

The other undefined in Eq. (2) is the XC functional $E_{\mathrm{XC}}\left[n\right]$. As is stated in Ref. [28], an exact expression for $E_{\mathrm{XC}}\left[n\right]$ is not known in general, and it is not easy to calculate it numerically either. The most familiar one is the LDA, where no density gradients are considered. The Gunnarson and Lundqvist(GL) LDA XC functional is widely used [28]. The functional derivative with respect to the density n reads:

$${\displaystyle \frac{\delta E_{\mathrm{XC}}^{\mathrm{GL}}\left[n\right]}{\delta n}}=E_h\left[-a_0{\left({\displaystyle \frac{3}{\pi }}\right)}^{1/3}{n}^{1/3}-0.0333ln\left(1+{\displaystyle \frac{1}{x}}\right)\right],$$

(6)

Here, $x={(3/4\pi )}^{1/3}{n}^{-1/3}/11.4a_0$. Another LDA XC functional is the Perdew-Zunger(PZ) LDA parametrization, which is also used both for Na and Ag [25,26]. It reads:

$${\displaystyle \frac{\delta E_{\mathrm{XC}}^{\mathrm{PZ}}\left[n\right]}{\delta n}}=E_h\left[-a_0{\left({\displaystyle \frac{3}{\pi }}\right)}^{1/3}{n}^{1/3}+{\mu }_C\left[n\right]\right],$$

(7)

in which

$${\mu }_C\left[n\right]={\displaystyle \frac{{\alpha }_1+(7{\alpha }_1{\beta }_1/6)\sqrt{r_s}+(4{\alpha }_1{\beta }_2/3)r_s}{{(1+{\beta }_1\sqrt{r_s}+{\beta }_2{r}_s)}^2}}$$

(8)

with ${\alpha }_1=-0.1423,$${\beta }_1=1.0529$, ${\beta }_2=0.3334$ , and $r_s={(3/4\pi n)}^{1/3}/a_0$. In addition, Wigner (WG) LDA XC funcitomal adopted in the Ref. [27,30] on work function is also considered. As follows,

$${\displaystyle \frac{\delta E_{\mathrm{XC}}^{\mathrm{WG}}\left[n\right]}{\delta n}}=-{\displaystyle \frac{e^2}{{\epsilon }_0}}{\displaystyle \frac{0.0909a_0}{{\left(0.0625+7.8a_0{n}^{1/3}\right)}^2}}{n}^{{\displaystyle \frac{2}{3}}}+{\displaystyle \frac{e^2}{{\epsilon }_0}}\left({\displaystyle \frac{0.0467}{0.0625+7.8a_0{n}^{1/3}}}-0.0783\right)n^{1/3}$$

(9)

Since the XC potential does not have a direct microscopic representation, its explicit form remains ambiguous. In DFT, the choice of XC potential is crucial, as it significantly affects the predicted ground state properties. However, in the context of QHT, the influence of the XC potential—particularly on ground state characteristics—remains poorly understood, with limited systematic investigation to date. To address this gap, the present work aims to systematically examine how different forms of the XC potential impact ground state properties within the QHT framework.

To solve the scattering field formulation within the linearized hydrodynamic model (Equations (1)), one must first find the ground-state electronic density $n_0\left(r\right)$. It can be solved by the static equation of QHT as follows,

$$\nabla {\left({\displaystyle \frac{\delta G}{\delta n}}\right)}_0-e{E}_0=0$$

(10)

The static electric field $E_0=-\nabla {\varphi }_0$ is solved self-consistently with the Poissonequation:

$$\nabla ·{\epsilon }_0{\epsilon }_r\nabla {\varphi }_0=q_e\left(n_{+}-n_0\right)$$

(11)

Substituting Eq. (11) into the static equation yields,

$${\left({\displaystyle \frac{\delta G}{\delta n}}\right)}_0-q_e{\varphi }_0=\mu$$

(12)

where $\mu$ is the chemical potential. At this time, the vW kinetic energy functional can be regarded as an approximate kinetic energy operator, and the function of vW can be approximated with $\sqrt{n_0}$, and the approximate Schrödinger equation can be obtained,

$${\lambda }_W{\displaystyle \frac{{\hslash }^2}{2m}}{\nabla }^2\sqrt{n_0}+V_{eff}\sqrt{n_0}=\mu \sqrt{n_0}$$

(13)

Where $V_{eff}={\left({\displaystyle \frac{\delta T_{TF}}{\delta n}}\right)}_0+{\left({\displaystyle \frac{\delta E_{xc}}{\delta n}}\right)}_0+q_e{\varphi }_0$. The ground state density can be obtained by Equations (11) and (13), a method also known as OFDFT. And the QHT can also be reformulated within the Time-Dependent OFDFT [29]. For the ground state density, the use of the TFvW kinetic energy functional is sufficient to solve the problem.

3. Energy Functional Required for Accurate Ground State Density in OFDFT

First, it is necessary to determine the appropriate energy functional for calculating the ground state density. In this section, we compute the $n_0$ distribution of a jellium slab using QHT based on Equations (11) and (13), and compare it with DFT results. For the DFT calculations, the ground state density of an infinite slab in the $xy$ plane was obtained using the open-source software Octopus. The results are stable when the slab thickness in the z direction exceeds $L>5{\lambda }_F$, where ${\lambda }_F=3.27r_s$ [32,33].

Figure 2(a) illustrates how the coefficient ${\lambda }_{\mathrm{W}}$ affects the calculated work function when using three representative XC potentials: WG (black line), GL (red line), and PZ (blue line). The horizontal dashed line denotes the reference value obtained from the DFT calculation ($3.05\mathrm{eV}$). Taking the WG XC potential as an example, the work function increases from $2.63\mathrm{eV}$ to $3.57\mathrm{eV}$ as ${\lambda }_{\mathrm{W}}$ increases, intersecting the reference line at ${\lambda }_{\mathrm{W}}=0.43$, which is consistent with Ref. [27]. Similarly, the GL and PZ curves intersect the DFT reference line at ${\lambda }_{\mathrm{W}}=0.48$ and $0.60$, respectively. These results indicate that matching the DFT work function requires different ${\lambda }_{\mathrm{W}}$ values depending on the chosen XC potential. Figure 2(b) shows the effective potential $V_{\mathrm{eff}}$, normalized by the chemical potential $\mu$, for the respective ${\lambda }_{\mathrm{W}}$ values determined above. The effective potentials obtained from all three functionals converge around $50,\mathrm{a}.\mathrm{u}.$, and the corresponding work functions align well with the DFT result. Therefore, to reproduce the DFT-calculated work function, the optimal ${\lambda }_{\mathrm{W}}$ values for WG, GL, and PZ are 0.43, 0.48, and 0.60, respectively. The results demonstrate that the work function is significantly influenced by ${\lambda }_{\mathrm{W}}$: a larger ${\lambda }_{\mathrm{W}}$ leads to a stronger electron spill-out effect and an increase in the work function. Additionally, the choice of XC potential also contributes to variations in the work function.

Figure 1. Schematic diagram of the system. The slab is infinite in the $xy$ plane, modeled using a jellium background with $r_s=4$. The study focuses on the ground state density and work function in the z direction.

Figure 1. Schematic diagram of the system. The slab is infinite in the $xy$ plane, modeled using a jellium background with $r_s=4$. The study focuses on the ground state density and work function in the z direction.

Figure 2. (a) Work function as a function of the coefficient ${\lambda }_{\mathrm{W}}$ using WG (black), GL (red), and PZ (blue) XC potentials. The horizontal dashed line represents the DFT reference value. The intersection points with the DFT reference are at ${\lambda }_{\mathrm{W}}=0.43$, $0.48$, and $0.60$, respectively. (b) Effective potentials corresponding to these intersections, normalized by $\mu$.

Figure 2. (a) Work function as a function of the coefficient ${\lambda }_{\mathrm{W}}$ using WG (black), GL (red), and PZ (blue) XC potentials. The horizontal dashed line represents the DFT reference value. The intersection points with the DFT reference are at ${\lambda }_{\mathrm{W}}=0.43$, $0.48$, and $0.60$, respectively. (b) Effective potentials corresponding to these intersections, normalized by $\mu$.

The ground state electron density near the metal surface plays a critical role in excited-state calculations. Figure 3 presents the normalized ground state density at the surface of the infinite slab, using $n_b$ for normalization. Figure 3(a) shows the ground-state density computed using three XC functionals, each with a value of ${\lambda }_{\mathrm{W}}$ that reproduces the same work function as obtained by DFT. Specifically, the values used are ${\lambda }_{\mathrm{W}}=0.43$ for WG (black), ${\lambda }_{\mathrm{W}}=0.48$ for GL (red), and ${\lambda }_{\mathrm{W}}=0.60$ for PZ (blue). Although these functionals yield consistent work functions, their surface electron density differ significantly. The density derived from WG (${\lambda }_{\mathrm{W}}=0.43$) and GL (${\lambda }_{\mathrm{W}}=0.48$) decay more rapidly near the surface compared to the DFT result, while the PZ (${\lambda }_{\mathrm{W}}=0.60$) result shows close agreement with DFT. This suggests that the decay rate of the ground-state density is primarily governed by the coefficient ${\lambda }_{\mathrm{W}}$, which modulates the electron spill-out effect. A larger ${\lambda }_{\mathrm{W}}$ leads to greater electron spill-out and a slower decay of surface density.Figure 3(b) explores the impact of XC potentials on the ground state density using a fixed ${\lambda }_{\mathrm{W}}=0.60$ for all three functionals. The results show that, under the same ${\lambda }_{\mathrm{W}}$, the choice of XC potential has little influence on the surface electron density. Thus, the surface density is primarily determined by ${\lambda }_{\mathrm{W}}$ rather than the specific form of the XC potential.

These results collectively show that in QHT, the calculated work function is influenced by both ${\lambda }_{\mathrm{W}}$ and the XC potential, while the surface ground state density is predominantly affected by ${\lambda }_{\mathrm{W}}$ alone. Among the tested functionals, PZ with ${\lambda }_{\mathrm{W}}=0.60$ provides both an accurate work function and a ground state density profile that best agrees with DFT.

4. Energy Functional Required for Accurate Optical Response in QHT

In the previous section, we established that using the PZ with ${\lambda }_{\mathrm{W}}=0.60$ energy functional within QHT produces work functions and ground state density that are consistent with DFT results. In this section, TDDFT results are used as a reference to study the appropriate energy functionals for excited states.We calculate the absorption cross section (normalized by ${\sigma }_0=\pi R^2$) of a jellium sphere under plane wave illumination, as is illustrated in Figure 4. The absorption cross section is given by

$$\sigma \left(\omega \right)={\displaystyle \frac{\omega }{2I_0}}\int Im\left\{\mathbf{E}·{\mathbf{P}}^{*}\right\}dV,$$

(14)

where $\mathbf{E}$ and $\mathbf{P}$ are the electric field and polarization, respectively.

4.1. Optical Response Calibration by TDDFT

To extract the resonance features, we follow the spectral fitting procedure described in Ref. [26], where only the first localized surface plasmon (LSP) peak is fitted. As an example, Figure 5 shows the normalized absorption cross section for a nanosphere with 1074 electrons. The fitted curve accurately reproduces the first LSP response.

For excited-state calculations, we use the PGSLN kinetic energy functional. In practice, the LSP response inside the nanosphere is primarily governed by the TFVW component. Therefore, to accurately capture the LSP energy, it is sufficient to determine the appropriate ${\lambda }_{\mathrm{W}}$ parameter in the excited-state functional.

To compare with the TDDFT benchmark data and ensure consistency with previous work [26], we use 13 jellium nanospheres with electron numbers ranging from 338 to 6174. For calibration, we select three representative nanospheres: a small (1074 electrons), medium (2048 electrons), and large (6174 electrons) sphere. Figure 6(a) shows how the LSP energy varies with the coefficient ${\lambda }_{\mathrm{W}}$ for the three spheres. The results indicate that the LSP energy increases approximately linearly with ${\lambda }_{\mathrm{W}}$. For instance, in the 1074-electron sphere, increasing ${\lambda }_{\mathrm{W}}$ from 0.80 to 1.00 raises the LSP energy from $3.21\mathrm{eV}$ to $3.26\mathrm{eV}$, yielding a linear slope of $0.25\mathrm{eV}$. For the 2048 and 6174-electron spheres, the slopes are $0.21\mathrm{eV}$ and $0.29\mathrm{eV}$, respectively. These results demonstrate that the LSP energy is sensitive to ${\lambda }_{\mathrm{W}}$—a higher ${\lambda }_{\mathrm{W}}$ results in greater LSP energy due to enhanced electron spill-out. However, the slope of this increase is smaller for larger spheres, suggesting that the spill-out effect becomes less significant as the system size increases. Figure 6(b) presents the average relative error in LSP energy (compared to TDDFT results) for the three nanospheres as a function of ${\lambda }_{\mathrm{W}}$. The minimum error occurs between ${\lambda }_{\mathrm{W}}=0.88$ and $0.92$, with both yielding comparable accuracy. To maintain consistency with the ground state calibration, ${\lambda }_{\mathrm{W}}=0.90$ is selected for subsequent excited-state calculations.

4.2. LSP Energy and Corresponding Oscillator Strength of Different Sizes

In the previous section, the excited-state energy functional was determined using three representative nanospheres. To further validate the accuracy and generalizability of our method across a broader range of sizes, Figure 7 presents the LSP energies and corresponding oscillator strengths for nanospheres of various sizes, as is calculated using OF-PGSLN (red lines) and is compared with TDDFT results (black lines) from Ref. [26].

Figure 7(a) shows the LSP energies for nanospheres with electron numbers ranging from 338 to 6174. The OF-PGSLN results exhibit excellent agreement with TDDFT, with an average absolute error of $0.014\mathrm{eV}$. This error is slightly larger than that obtained using the DFT ground state density ($0.006\mathrm{eV}$), but consistent with the results based on the model ground state density reported in Ref. [26].Figure 7(b) displays the corresponding oscillator strengths, normalized to the classical local-response approximation value $f_{\mathrm{osc},\mathrm{Mie}}$. The average error across all 13 nanospheres is 0.021, which is slightly lower than the value obtained using the DFT ground state density (0.029), but higher than that of the model ground state density (0.008), as reported in Ref. [26].

These results demonstrate that the OF-PGSLN approach yields LSP energies and oscillator strengths in good agreement with TDDFT across a wide range of nanosphere sizes. In comparison with methods based on the DFT or model ground state density, OF-PGSLN makes a slightly larger error in oscillator strength than the model density, but performs better than the DFT-based method. For LSP energy, the accuracy of OF-PGSLN is comparable to that obtained using the model ground state density.In conclusion, both the ground and excited-state energy functionals in OF-PGSLN have been calibrated taking DFT and TDDFT results as references. Subsequent validation across nanospheres of varying sizes confirms that the LSP energy and oscillator strength predicted by OF-PGSLN exhibit smaller mean errors relative to TDDFT. A key advantage of the OF-PGSLN approach is its ability to bypass the computational complexity of DFT, while maintaining accuracy and enabling generalization to other physical systems.

5. Validity and limitations of Linear Superposition Approximation in Nanodimers

In the previous section, we demonstrated that the OF-PGSLN ground state density enables accurate predictions of both the LSP energy and oscillator strength in sodium nanospheres. We now extend this analysis to a dimer system composed of two sodium jellium spheres, each containing 1074 electrons. This system provides a platform for evaluating the validity of a linear superposition of ground state density across varying interparticle distances. Specifically, we assess the performance of the OF-derived ground state density in comparison with that obtained from a linear superposition of model ground state density for different gap sizes.

We begin with the small-gap regime, where strong quantum effects emerge due to significant wavefunction overlap. As shown in Figure 8(a), at $\mathrm{gap}=0.3\mathrm{nm}$, the absorption spectrum based on the model density (red line) exhibits a marked redshift in the first LSP peak relative to the spectrum obtained using the OF density (black line), with a difference of $56\mathrm{meV}$. To elucidate the source of this discrepancy, Figure 8(b) presents the ground state electron density at the center of the dimer. The OF result yields a minimum density of $1.5\times {10}^{-1}$, whereas the model result gives a significantly lower value of $9.3\times {10}^{-2}$. This underestimation arises from the model’s neglect of essential quantum mechanical effects—namely electron–electron repulsion, Fermi pressure, and quantum confinement—which become critical in the strongly coupled regime. These findings reveal the breakdown of the linear approximation at short distances and emphasize the necessity of using the OF ground state density when modeling closely spaced nanostructures.

In contrast, the large-gap regime corresponds to a situation in which the electron density of the two nanospheres remain largely non-overlapping, and mutual interactions are minimal. As illustrated in Figure 9(a), at $\mathrm{gap}=1.0\mathrm{nm}$, the absorption spectra calculated from both the model density (red line) and the OF density (black line) are nearly identical, differing by only $6\mathrm{meV}$ in LSP energy. This agreement is further supported by Figure 9(b), which shows the minimum electron densityat the dimer center: $6.38\times {10}^{-5}$ for the model density and $1.39\times {10}^{-4}$ for the OF density. The small discrepancy confirms the validity of the linear superposition approximation in the weakly coupled regime. Notably, this approximation offers a practical advantage: for systems comprising many nanostructures, individual ground state density can be precomputed and linearly combined, significantly reducing computational cost without sacrificing accuracy.

To quantitatively evaluate the performance of the linear approximation across different interparticle separations, Figure 10 plots the LSP energy error as a function of dimer gap. While the error remains below $10\mathrm{meV}$ for gaps larger than $0.8\mathrm{nm}$, it increases sharply below $0.6\mathrm{nm}$, saturating at approximately $54.2\mathrm{meV}$ in the strongly coupled regime. These results highlight the dual nature of the linear superposition approach: it fails under strong coupling but becomes a powerful and efficient modeling tool when the constituent nanostructures are sufficiently separated. This strategy is especially advantageous for simulating large nanoparticle assemblies, where reusing precomputed ground state density can drastically reduce computational overhead.

6. Conclusions

In this work, we developed OF-PGSLN, a QHT framework that integrates OFDFT for ground-state density with a Laplacian-level kinetic energy functional (PGSLN) for excited-state dynamics. The key conclusions are as follows:Using DFT as a benchmark, we determined that the PZ XC functional with ${\lambda }_{\mathrm{W}}=0.60$ yields the most accurate ground-state properties, including work functions and surface electron density. For excited states, calibrating against TDDFT, we found that ${\lambda }_{\mathrm{W}}=0.90$ enables precise reproduction of LSP resonances. Across sodium nanospheres with $338\le N_e\le 6174$, the average errors in resonance energy and oscillator strength are only $0.014\mathrm{eV}$ and $0.021$, respectively. We further applied the method to nanodimers, demonstrating that while the linear superposition model suffices at large gaps ($>0.8\mathrm{nm}$), it fails to capture critical density redistribution at small gaps ($<0.6\mathrm{nm}$), leading to significant errors—up to $56\mathrm{meV}$—which OF-PGSLN accurately resolves. By enabling stable, geometry-independent ground-state calculations and suppressing spurious spectral features through PGSLN, our approach provides a reliable and efficient platform for modeling quantum optical responses in complex nanostructures, with broad applicability to next-generation plasmonic and nanophotonic devices.