Self-Consistent Hydrodynamic Drude Model

• SC-HDM is a mesoscale model that integrates nonlocal electron dynamics, gradient corrections, and electron spill-out for accurate simulation of charge, current, and energy transport.
• It couples hydrodynamic conservation laws with Maxwell’s equations to capture plasmon dispersion, viscoelastic damping, and topological effects in metallic systems.
• The model underpins computational approaches such as finite-element and VIE methods, successfully validating experimental plasmonic and electronic transport phenomena.

The Self-Consistent Hydrodynamic Drude Model (SC-HDM) generalizes the conventional hydrodynamic Drude model to account for nonlocal electron dynamics, density-gradient corrections, electron spill-out, and viscoelastic exchange-correlation effects, yielding a self-consistent and physically grounded framework for charge, current, and energy transport in conducting systems. SC-HDM is formulated by coupling hydrodynamic conservation laws with approximations to the electronic energy functional, and is solved self-consistently with Maxwell’s equations or in kinetic response form, depending on context. The model reconciles the limitations of hard-wall boundary conditions and local Drude approximations, providing accurate predictions for transport, optical, and topological phenomena in metals, plasmonic nanostructures, and strongly correlated systems.

1. Fundamental Equations and Energy Functionals

SC-HDM is constructed from the continuum equations for charge density $n(\mathbf r,t)$, fluid velocity $\mathbf v(\mathbf r,t)$, and electromagnetic fields, supplemented by an internal energy functional $G[n]$. The kinetic energy is treated using a combination of the Thomas–Fermi (TF) term and the von Weizsäcker (vW) gradient correction: $$T_\text{TFW}[n] = T_\text{TF}[n] + T_\text{W}[n] = \frac{3}{10} \frac{\hbar^2}{m} (3\pi^2)^{2/3} \int n^{5/3} (\mathbf r) d^3 r + \frac{1}{72} \frac{\hbar^2}{m} \int \frac{|\nabla n(\mathbf r)|^2}{n(\mathbf r)} d^3 r$$ Exchange-correlation enters either via a local-density functional (e.g., Gunnarson–Lundqvist LDA) or via a viscoelastic stress tensor (Vignale–Kohn formalism) in more refined approaches (Ciracì, 2016). The governing equations for the hydrodynamic current $\mathbf J$ and continuity are (Toscano et al., 2014, Ciracì, 2016): $$\frac{\partial n}{\partial t} + \nabla \cdot \mathbf J = 0$$

$$\frac{\partial \mathbf J}{\partial t} = \frac{n e^2}{m}\mathbf E - \frac{e}{m} \mathbf J \times (\mathbf B + \nabla \times \mathbf A_\text{xc}) + \frac{n e}{m} \left( \nabla v_\text{xc}[n] - e \partial_t \mathbf A_\text{xc} \right) + \frac{e}{m} \nabla \cdot \Pi$$

where $\Pi$ contains TF and vW terms with appropriate functional derivatives.

2. Self-Consistent Linear Response and Nonlocal Effects

At the linear-response level, SC-HDM rigorously incorporates nonlocal response arising from gradient terms. The coupling to Maxwell's equations is established either via a direct finite-element implementation (Toscano et al., 2014, Ciracì, 2016) or, more recently, via a Volume Integral Equation (VIE) approach (Mystilidis et al., 28 Dec 2025), which exploits system symmetries for computational efficiency. The polarization density $P(\mathbf r)$ satisfies: $$\frac{e n_0(\mathbf r)}{m_e} \nabla \left( \delta G/\delta n \right)_1 + \omega(\omega+i\gamma) P(\mathbf r) = -\epsilon_0 \omega_p^2(\mathbf r) E(\mathbf r)$$ with local and nonlocal material response explicitly encoded through analytic $\gamma$-kernels that depend on $n_0(\mathbf r)$.

Electron spill-out is incorporated self-consistently through a smooth equilibrium density $n_0(\mathbf r)$, determined by the static Euler and Poisson equations, removing the need for hard-wall boundary conditions and allowing the natural emergence of induced charge decay (Toscano et al., 2014, Mystilidis et al., 28 Dec 2025).

3. Viscoelastic Exchange-Correlation and Plasmonic Damping

A distinguishing feature of SC-HDM is the inclusion of current-dependent exchange-correlation potentials that manifest as viscoelastic damping and additional $k$-dependent corrections to plasmon dispersion. In the Vignale–Kohn formalism (Ciracì, 2016), the XC vector potential is related to a generalized stress tensor: $$\sigma_{\mu\nu} = \tilde\eta(\omega, n)\left( \partial_\mu v_\nu + \partial_\nu v_\mu - \frac{2}{3} \delta_{\mu\nu} \nabla \cdot \mathbf v \right) + \tilde\zeta(\omega, n) \delta_{\mu\nu} \nabla \cdot \mathbf v$$ leading to broadened resonance linewidths and intrinsic nonlocal plasmon damping: $$\omega(k) \simeq \sqrt{\omega_p^2 + \beta^2 k^2} - i\left( \frac{\gamma}{2} + \frac{\zeta_\text{xc}}{2} k^2 \right)$$ This theoretical prediction is in quantitative agreement with ab initio TDDFT for nanospheres and tunneling gaps down to sub-nm (Ciracì, 2016).

4. Gradient Corrections, Drude Pole Renormalization, and Magnetic Response

SC-HDM extends beyond local Drude transport by retaining the first nontrivial gradient corrections to momentum conservation and energy flux (Goutéraux et al., 2023). In metals with weakly broken translation invariance, the constitutive relations for charge current and entropy flux at second order in gradients introduce new coefficients $\lambda_\rho, \lambda_s, \lambda_v$: $$j^i = \rho v^i - \sigma_0 \partial^i \mu - \alpha_0 \partial^i T + \ell^2 \lambda_\rho \partial^i \Phi^I \hat D\Phi_I + O(\ell^4)$$ These corrections contribute to the AC and DC thermoelectric conductivities as renormalized Drude weights and modify the cyclotron resonance in a magnetic field: $$\sigma(\omega) = \sigma_0 + \frac{(ρ + ℓ^2 λ_ρ)^2}{(χ_{ππ} + ℓ^2 λ_v)(Γ - iω)} + O(ℓ^2, ω)$$

$$ω_c = \frac{B (ρ + 2 ℓ^2 λ_ρ)}{χ_{ππ} + ℓ^2 λ_v} + O(ℓ^4)$$

Experimentally, these renormalizations resolve effective mass discrepancies between optical and thermodynamic measures in overdoped cuprates and explain Hall carrier density observations (Goutéraux et al., 2023).

5. Topological Properties and Boundary Conditions

SC-HDM possesses well-defined topological invariants in continuous systems, resolving the breakdowns inherent in local Drude models for high-wavenumber modes. Quantum-pressure-induced spatial dispersion furnishes a physical cutoff, restoring integer Chern numbers for bulk modes and robustly establishing the bulk-edge correspondence without ad hoc parameters (Pakniyat et al., 2021): $$\chi(k,\omega) = \frac{1}{1 - k^2/k_m^2}, \; k_m= \omega/\beta$$ The two TM bands admit $C_1=+1$, $C_2=-1$, with TE bands trivial, and the sum $C_1+C_2=0$. At interfaces, hydrodynamic boundary conditions enforce:

• Continuity of tangential $\mathbf E$ and normal displacement
• No normal penetration for the hydrodynamic current
• Continuity of pressure gradient at the surface

This formulation enables well-posed edge-mode spectrums and the physical realization of unidirectional topological SPPs.

6. Numerical and Computational Implementation

SC-HDM is typically implemented via finite-element solvers (e.g., COMSOL Multiphysics), solving nonlinear static and linear dynamic systems with full self-consistency in the electronic density and field responses (Toscano et al., 2014, Ciracì, 2016). The Volume Integral Equation method (Mystilidis et al., 28 Dec 2025) further accelerates these computations for spherical nanoparticles by transforming PDEs into efficiently solvable matrix equations with symmetry-adapted basis functions. Key benchmarks demonstrate that the VIE method for SC-HDM achieves simulation times comparable to local or simple hydrodynamic models, even as it captures nonlocal and spill-out physics.

Numerically, mesoscopic quantities such as Feibelman $d$-parameters, spatially resolved induced charge densities, and permittivity profiles are directly extractable from the SC-HDM framework, enabling direct comparisons to experiment and quantum chemistry benchmarks.

7. Model Variants: Halevi, GNOR, and Hybrid Response

The Halevi extension (Wegner et al., 2022) generalizes SC-HDM by introducing a frequency-dependent hydrodynamic velocity $\beta^2(\omega)$, accounting for the crossover from viscous to elastic electron-fluid response and generating a hybrid diffusion–wave current via auxiliary equations of motion: $$\partial_t \mathbf J_D + \gamma \mathbf J_D = -\frac{4 v_F^2}{15} \nabla\rho$$

$$\partial_t \mathbf J + \gamma \mathbf J = \varepsilon_0 \omega_p^2 \mathbf E - \beta_\text{HF}^2 \nabla\rho - \gamma (\mathbf J + \mathbf J_D)$$

This framework naturally predicts size-dependent damping and oscillatory contributions to the diffusive currents. Compared to GNOR, the diffusive length scale is set by kinetic theory and not by phenomenological constants, and the implementation in time-domain Maxwell solvers (FDTD, DGTD) is straightforward due to the limiting structure of the auxiliary variables.

8. Experimental Validation and Applications

SC-HDM predictions match experimental surface plasmon resonance (SPR) shifts for Na and Ag nanowires, reproduce Bennett resonances, and efficiently capture blueshifts/redshifts due to electron spill-out and screening effects (Toscano et al., 2014, Mystilidis et al., 28 Dec 2025). In strongly correlated electron systems, SC-HDM furnishes new diagnostics for transport measurements and resolves effective mass puzzles (Goutéraux et al., 2023). In quantum plasmonics, SC-HDM enables practical simulation of large NP arrays, gap structures, and topological band invariants with computational effort competitive with fully classical methods.

---

SC-HDM, by systematically encoding gradient, viscoelastic, spill-out, and topological effects within a unified kinetic–hydrodynamic framework, furnishes the state-of-the-art mesoscale model for electronic and photonic response in metals and plasmonic nanostructures, overcoming the limitations of previous local and phenomenological theories.