Variational schemes and geometric simulations for a hydrodynamic-electrodynamic model of surface plasmon polaritons

A class of variational schemes for the hydrodynamic-electrodynamic model of lossless free-electron gas in a quasineutral background is developed for high-quality simulations of surface plasmon polaritons. The Lagrangian density of lossless free-electron gas with a self-consistent electromagnetic field is established, and the dynamical equations with the associated constraints are obtained via a variational principle. Based on discrete exterior calculus, the action functional of this system is discretized and minimized to obtain the discrete dynamics. Newton-Raphson iteration and the biconjugate gradient stabilized method are equipped as a hybrid nonlinear-linear algebraic solver. Instead of discretizing the partial differential equations, the variational schemes have better numerical properties in secular simulations, as they preserve the discrete Lagrangian symplectic structure, gauge symmetry, and general energy-momentum density. Two numerical experiments were performed. The numerical results reproduce characteristic dispersion relations of bulk plasmons and surface plasmon polaritons, and the numerical errors of conserved quantities in all experiments are bounded by a small value after long-term simulations.

Figure 1

SPP configuration at an infinite metal-dielectric interface.

Figure 2

DEC-based discretization of the hydrodynamic-electrodynamic model in a rectangular lattice. The discrete spacelike submanifold is shown, and other types of discrete submanifolds can be given in the same way.

Figure 3

The complete iteration for the variational schemes of the hydrodynamic-electrodynamic model.

Figure 4

The evolution of the gauge field in the simulation of bulk plasmon oscillation.

Figure 5

Bulk plasmon dispersion relations: the numerical result (contour plot) vs the analytical one (solid line).

Figure 6

Numerical error of the Hamiltonian in the simulation: (a) the Hamiltonian of lossless free-electron gas; (b) the Hamiltonian of a self-consistent gauge field; and (c) the total Hamiltonian error. The Hamiltonian is in units of J.

Figure 7

The 2D SPP configuration used in the simulations.

Figure 8

SPP dispersion relations. Analytical vs numerical. Numerical results obtained at given wavelengths have good consistency with the analytical curve. The inset plots the wavelengths and decay lengths at the simulation region, which illustrates the subwavelength and energy localization of SPPs.

Figure 9

The slice of normalized SPP field $H_y$ obtained by simulations. The vacuum wavelength ${\lambda }_0=240$ nm. The insets plot the $X$ and $Z$ projections of $H_y$ at the region [35,65]$\times$[25,75]. It shows that the electromagnetic field of SPPs is strongly bounded by the interface, and the mode has a subwavelength spatial structure.