A class of variational schemes for the hydrodynamic-electrodynamic model of lossless free electron gas in a quasi-neutral 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.In the past two decades, there have been impressive developments and significant advancement in applications of Surface Plasmon Polaritons (SPPs), bringing many new ideas into traditional electromagnetics and optics, such as the lithography beyond the diffraction limit, chip-scale photonic circuits, plasmonic metasurfaces, bio-photonics, etc. [1][2][3][4][5][6][7][8][9][10]. In the field of metal optics, plasmonics focuses on the collective motions of free electron gas in metal with self-consistent and external electromagnetic fields whose first-principle model is the classical particle-field theory [1][2][3]. Direct applications of the first-principle model in macroscopic simulations face many obstacles, such as the nonlinearity, the multi-scale, and the huge degrees-of-freedom. As a simplification, linearized phenomenological models, e.g., the Drude-Lorentz (DL) model, are widely used to describe macroscopic plasmonic phenomena [2,3]. In a mesoscopic context, kinetic and hydrodynamic descriptions are basic physical models of plasmonics, which involve both the dynamics of free electron gas and an electromagnetic field. Therefore, high-quality numerical schemes and simulations based on the hydrodynamic model are necessary in plasmonic research.The physics of SPPs can be described by the free electron gas model whose dynamical equations are hydrodynamic and Maxwell's equations [1]. For Maxwell's equations, many numerical schemes, such as the Finite-Difference Time-Domain (FDTD) method, the Finite Element (FE) method, and the Method of Moments (MoM) have been developed, that are widely used in modern electromagnetic engineering, Radio Frequency (RF) and microwave engineering, terahertz engineering, optical engineering, metamaterial design, accelerator design, fusion engineering, radio astrophysics, geophysics...

The relation between symmetries and local conservation laws, known as Noether's theorem, plays an important role in modern theoretical physics. As a discrete analog of the differentiable physical system, a good numerical scheme should admit the discrete local conservation laws and inherent mathematical structures. A class of variational schemes constructed for the hydrodynamic-electrodynamic model of lossless free-electron gas in a quasi-neutral background shows good properties in secular simulations of surface plasmon polaritons [Q. Chen et al., Phys. Rev. E 99, 023313 (2019)]. We show the discrete local conservation laws admitted by these schemes. Based on the gauge symmetry of the discrete action functional, a discrete charge conservation law is realized locally, which is consistent with the discrete Euler-Lagrange equations obtained from the variational schemes. Based on the discrete Euler-Lagrange equations, discrete local momentum and energy conservation laws are derived directly, which are rigorous in theory. The preservation of the discrete local conservation laws and Lagrangian symplectic structure ensure that the numerical scheme is correct in physics.